OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
NPAT.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "integrators/NPAT.hpp"
46
47#include "brains/SimInfo.hpp"
48#include "brains/Thermo.hpp"
49#include "integrators/IntegratorCreator.hpp"
51#include "utils/Constants.hpp"
52#include "utils/simError.h"
53
54namespace OpenMD {
55
56 NPAT::NPAT(SimInfo* info) : NPT(info) {
57 Globals* simParams = info_->getSimParams();
58
59 // Default value of privilegedAxis is "z"
60 if (simParams->getPrivilegedAxis() == "x")
61 axis_ = 0;
62 else if (simParams->getPrivilegedAxis() == "y")
63 axis_ = 1;
64 else if (simParams->getPrivilegedAxis() == "z")
65 axis_ = 2;
66 }
67
68 void NPAT::evolveEtaA() {
69 eta(axis_, axis_) +=
70 dt2 * instaVol *
71 (press(axis_, axis_) - targetPressure / Constants::pressureConvert) /
72 (NkBT * tb2);
73 oldEta_ = eta;
74 }
75
76 void NPAT::evolveEtaB() {
77 prevEta_ = eta;
78 eta(axis_, axis_) = oldEta_(axis_, axis_) +
79 dt2 * instaVol *
80 (press(axis_, axis_) -
81 targetPressure / Constants::pressureConvert) /
82 (NkBT * tb2);
83 }
84
85 void NPAT::calcVelScale() {
86 for (int i = 0; i < 3; i++) {
87 for (int j = 0; j < 3; j++) {
88 vScale_(i, j) = eta(i, j);
89
90 if (i == j) { vScale_(i, j) += thermostat.first; }
91 }
92 }
93 }
94
95 void NPAT::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
96 sc = vScale_ * vel;
97 }
98
99 void NPAT::getVelScaleB(Vector3d& sc, int index) {
100 sc = vScale_ * oldVel[index];
101 }
102
103 void NPAT::getPosScale(const Vector3d& pos, const Vector3d& COM, int index,
104 Vector3d& sc) {
105 /**@todo */
106 Vector3d rj = (oldPos[index] + pos) / (RealType)2.0 - COM;
107 sc = eta * rj;
108 }
109
110 void NPAT::scaleSimBox() {
111 Mat3x3d scaleMat;
112
113 for (int i = 0; i < 3; i++) {
114 for (int j = 0; j < 3; j++) {
115 scaleMat(i, j) = 0.0;
116 if (i == j) { scaleMat(i, j) = 1.0; }
117 }
118 }
119
120 scaleMat(axis_, axis_) = exp(dt * eta(axis_, axis_));
121 Mat3x3d hmat = snap->getHmat();
122 hmat = hmat * scaleMat;
123 snap->setHmat(hmat);
124 }
125
126 bool NPAT::etaConverged() {
127 int i;
128 RealType diffEta, sumEta;
129
130 sumEta = 0;
131 for (i = 0; i < 3; i++) {
132 sumEta += pow(prevEta_(i, i) - eta(i, i), 2);
133 }
134
135 diffEta = sqrt(sumEta / 3.0);
136
137 return (diffEta <= etaTolerance);
138 }
139
140 RealType NPAT::calcConservedQuantity() {
141 thermostat = snap->getThermostat();
142 loadEta();
143
144 // We need NkBT a lot, so just set it here: This is the RAW number
145 // of integrableObjects, so no subtraction or addition of constraints or
146 // orientational degrees of freedom:
147 NkBT = info_->getNGlobalIntegrableObjects() * Constants::kB * targetTemp;
148
149 // fkBT is used because the thermostat operates on more degrees of freedom
150 // than the barostat (when there are particles with orientational degrees
151 // of freedom).
152 fkBT = info_->getNdf() * Constants::kB * targetTemp;
153
154 RealType conservedQuantity;
155 RealType totalEnergy;
156 RealType thermostat_kinetic;
157 RealType thermostat_potential;
158 RealType barostat_kinetic;
159 RealType barostat_potential;
160 RealType trEta;
161
162 totalEnergy = thermo.getTotalEnergy();
163
164 thermostat_kinetic = fkBT * tt2 * thermostat.first * thermostat.first /
165 (2.0 * Constants::energyConvert);
166
167 thermostat_potential = fkBT * thermostat.second / Constants::energyConvert;
168
169 SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
170 trEta = tmp.trace();
171
172 barostat_kinetic = NkBT * tb2 * trEta / (2.0 * Constants::energyConvert);
173
174 barostat_potential =
175 (targetPressure * thermo.getVolume() / Constants::pressureConvert) /
176 Constants::energyConvert;
177
178 conservedQuantity = totalEnergy + thermostat_kinetic +
179 thermostat_potential + barostat_kinetic +
180 barostat_potential;
181
182 return conservedQuantity;
183 }
184
185 void NPAT::loadEta() {
186 eta = snap->getBarostat();
187
188 // if (!eta.isDiagonal()) {
189 // snprintf( painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
190 // "NPAT error: the diagonal elements of eta matrix are not the
191 // same or etaMat is not a diagonal matrix");
192 // painCave.isFatal = 1;
193 // simError();
194 //}
195 }
196
197 void NPAT::saveEta() { snap->setBarostat(eta); }
198
199} // namespace OpenMD
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
int getNdf()
Returns the number of degrees of freedom.
Definition SimInfo.hpp:220
int getNGlobalIntegrableObjects()
Returns the total number of integrable objects (total number of rigid bodies plus the total number of...
Definition SimInfo.hpp:139
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
void setHmat(const Mat3x3d &m)
Sets the H-Matrix.
Definition Snapshot.cpp:217
Real trace() const
Returns the trace of this matrix.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.