# | Line 1 | Line 1 | |
---|---|---|
1 | #include <math.h> | |
2 | + | #include "Atom.hpp" |
3 | + | #include "Molecule.hpp" |
4 | + | #include "SimInfo.hpp" |
5 | + | #include "Thermo.hpp" |
6 | + | #include "ExtendedSystem.hpp" |
7 | + | #include "simError.h" |
8 | ||
9 | < | ExtendedSystem::ExtendedSystem( SimInfo &info ) { |
9 | > | ExtendedSystem::ExtendedSystem( SimInfo* the_entry_plug ) { |
10 | ||
11 | // get what information we need from the SimInfo object | |
12 | ||
13 | < | entry_plug = &info; |
14 | < | nAtoms = info.n_atoms; |
15 | < | atoms = info.atoms; |
16 | < | nMols = info.n_mol; |
17 | < | molecules = info.molecules; |
13 | > | entry_plug = the_entry_plug; |
14 | > | zeta = 0.0; |
15 | > | epsilonDot = 0.0; |
16 | > | epsilonDotX = 0.0; |
17 | > | epsilonDotY = 0.0; |
18 | > | epsilonDotZ = 0.0; |
19 | > | have_tau_thermostat = 0; |
20 | > | have_tau_barostat = 0; |
21 | > | have_target_temp = 0; |
22 | > | have_target_pressure = 0; |
23 | > | have_qmass = 0; |
24 | ||
25 | } | |
26 | ||
27 | < | ExtendedSystem::~ExtendedSystem() { |
16 | < | } |
27 | > | void ExtendedSystem::NoseHooverNVT( double dt, double ke ){ |
28 | ||
18 | – | |
19 | – | void ExtendedSystem::nose_hoover_nvt( double ke, double dt, double temp ){ |
20 | – | |
29 | // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 | |
30 | ||
31 | < | int i, j, degrees_freedom; |
32 | < | double ke, dt, temp, kB; |
33 | < | double keconverter, NkBT, zetaScale, ke_temp; |
34 | < | double vxi, vyi, vzi, jxi, jyi, jzi; |
35 | < | |
36 | < | degrees_freedom = 6*nmol; // number of degrees of freedom for the system |
37 | < | kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
30 | < | keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K |
31 | < | |
32 | < | ke_temp = ke * keconverter; |
33 | < | NkBT = degrees_freedom*kB*temp; |
31 | > | int i; |
32 | > | double NkBT, zetaScale, ke_temp; |
33 | > | double vx, vy, vz, jx, jy, jz; |
34 | > | const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
35 | > | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to |
36 | > | // amu*Ang^2*fs^-2/K |
37 | > | DirectionalAtom* dAtom; |
38 | ||
39 | < | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin & |
36 | < | // qmass is set in the parameter file |
37 | < | zeta = zeta + dt*((ke_temp*2 - NkBT)/qmass); |
38 | < | zetaScale = zeta * dt; |
39 | > | if (this->NVTready()) { |
40 | ||
41 | < | // perform thermostat scaling on linear velocities and angular momentum |
42 | < | for(i = 0, i < nmol; i++ ) { |
43 | < | vxi = vx(i)*zetaScale; |
44 | < | vyi = vy(i)*zetaScale; |
45 | < | vzi = vz(i)*zetaScale; |
46 | < | jxi = jx(i)*zetaScale; |
47 | < | jyi = jy(i)*zetaScale; |
48 | < | jzi = jz(i)*zetaScale; |
49 | < | |
50 | < | vx(i) = vx(i) - vxi; |
51 | < | vy(i) = vy(i) - vyi; |
52 | < | vz(i) = vz(i) - vzi; |
53 | < | jx(i) = jx(i) - jxi; |
54 | < | jy(i) = jy(i) - jyi; |
55 | < | jz(i) = jz(i) - jzi; |
41 | > | atoms = entry_plug->atoms; |
42 | > | |
43 | > | ke_temp = ke * e_convert; |
44 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
45 | > | |
46 | > | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
47 | > | // qmass is set in the parameter file |
48 | > | |
49 | > | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
50 | > | |
51 | > | zetaScale = zeta * dt; |
52 | > | |
53 | > | //std::cerr << "zetaScale = " << zetaScale << "\n"; |
54 | > | |
55 | > | // perform thermostat scaling on linear velocities and angular momentum |
56 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
57 | > | |
58 | > | vx = atoms[i]->get_vx(); |
59 | > | vy = atoms[i]->get_vy(); |
60 | > | vz = atoms[i]->get_vz(); |
61 | > | |
62 | > | atoms[i]->set_vx(vx * (1.0 - zetaScale)); |
63 | > | atoms[i]->set_vy(vy * (1.0 - zetaScale)); |
64 | > | atoms[i]->set_vz(vz * (1.0 - zetaScale)); |
65 | > | } |
66 | > | if( entry_plug->n_oriented ){ |
67 | > | |
68 | > | for( i=0; i < entry_plug->n_atoms; i++ ){ |
69 | > | |
70 | > | if( atoms[i]->isDirectional() ){ |
71 | > | |
72 | > | dAtom = (DirectionalAtom *)atoms[i]; |
73 | > | |
74 | > | jx = dAtom->getJx(); |
75 | > | jy = dAtom->getJy(); |
76 | > | jz = dAtom->getJz(); |
77 | > | |
78 | > | dAtom->setJx(jx * (1.0 - zetaScale)); |
79 | > | dAtom->setJy(jy * (1.0 - zetaScale)); |
80 | > | dAtom->setJz(jz * (1.0 - zetaScale)); |
81 | > | } |
82 | > | } |
83 | > | } |
84 | } | |
85 | } | |
86 | ||
87 | ||
88 | < | void ExtendedSystem::nose_hoover_anderson_npt(double pressure, double ke, double dt, |
89 | < | double temp ) { |
88 | > | void ExtendedSystem::NoseHooverAndersonNPT( double dt, |
89 | > | double ke, |
90 | > | double p_tensor[9] ) { |
91 | ||
92 | // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 | |
93 | // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500 | |
94 | ||
95 | < | int i, j, degrees_freedom; |
96 | < | double pressure, dt, temp, pressure_units, epsilonScale; |
97 | < | double ke, kB, vxi, vyi, vzi, pressure_ext; |
98 | < | double boxx_old, boxy_old, boxz_old; |
99 | < | double keconverter, NkBT, zetaScale, ke_temp; |
100 | < | double jxi, jyi, jzi, scale; |
95 | > | double oldBox[3]; |
96 | > | double newBox[3]; |
97 | > | const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
98 | > | const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1 |
99 | > | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to |
100 | > | // amu*Ang^2*fs^-2/K |
101 | ||
102 | < | kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
103 | < | pressure_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1 |
104 | < | degrees_freedom = 6*nmol; // number of degrees of freedom for the system |
105 | < | keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2/K |
102 | > | int i; |
103 | > | double p_ext, zetaScale, epsilonScale, scale, NkBT, ke_temp; |
104 | > | double volume, p_mol; |
105 | > | double vx, vy, vz, jx, jy, jz; |
106 | > | DirectionalAtom* dAtom; |
107 | ||
108 | < | pressure_ext = pressure * pressure_units; |
109 | < | volume = boxx*boxy*boxz; |
110 | < | ke_temp = ke * keconverter; |
111 | < | NkBT = degrees_freedom*kB*temp; |
108 | > | if (this->NPTready()) { |
109 | > | atoms = entry_plug->atoms; |
110 | > | |
111 | > | p_ext = targetPressure * p_units; |
112 | > | p_mol = (p_tensor[0] + p_tensor[4] + p_tensor[8])/3.0; |
113 | > | |
114 | > | entry_plug->getBox(oldBox); |
115 | > | volume = oldBox[0]*oldBox[1]*oldBox[2]; |
116 | > | |
117 | > | ke_temp = ke * e_convert; |
118 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
119 | > | |
120 | > | // propagate the strain rate |
121 | > | |
122 | > | epsilonDot += dt * ((p_mol - p_ext) * volume / |
123 | > | (tauBarostat*tauBarostat * kB * targetTemp) ); |
124 | > | |
125 | > | // determine the change in cell volume |
126 | > | scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0)); |
127 | > | //std::cerr << "pmol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n"; |
128 | > | |
129 | > | newBox[0] = oldBox[0] * scale; |
130 | > | newBox[1] = oldBox[1] * scale; |
131 | > | newBox[2] = oldBox[2] * scale; |
132 | > | volume = newBox[0]*newBox[1]*newBox[2]; |
133 | > | |
134 | > | entry_plug->setBox(newBox); |
135 | > | |
136 | > | // perform affine transform to update positions with volume fluctuations |
137 | > | this->AffineTransform( oldBox, newBox ); |
138 | > | |
139 | > | epsilonScale = epsilonDot * dt; |
140 | > | |
141 | > | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
142 | > | // qmass is set in the parameter file |
143 | > | |
144 | > | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
145 | > | zetaScale = zeta * dt; |
146 | > | |
147 | > | //std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale << "\n"; |
148 | > | |
149 | > | // apply barostating and thermostating to velocities and angular momenta |
150 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
151 | > | |
152 | > | vx = atoms[i]->get_vx(); |
153 | > | vy = atoms[i]->get_vy(); |
154 | > | vz = atoms[i]->get_vz(); |
155 | > | |
156 | > | atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale)); |
157 | > | atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale)); |
158 | > | atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale)); |
159 | > | } |
160 | > | if( entry_plug->n_oriented ){ |
161 | > | |
162 | > | for( i=0; i < entry_plug->n_atoms; i++ ){ |
163 | > | |
164 | > | if( atoms[i]->isDirectional() ){ |
165 | > | |
166 | > | dAtom = (DirectionalAtom *)atoms[i]; |
167 | > | |
168 | > | jx = dAtom->getJx(); |
169 | > | jy = dAtom->getJy(); |
170 | > | jz = dAtom->getJz(); |
171 | > | |
172 | > | dAtom->setJx( jx * (1.0 - zetaScale)); |
173 | > | dAtom->setJy( jy * (1.0 - zetaScale)); |
174 | > | dAtom->setJz( jz * (1.0 - zetaScale)); |
175 | > | } |
176 | > | } |
177 | > | } |
178 | > | } |
179 | > | } |
180 | ||
82 | – | // propogate the strain rate |
181 | ||
182 | < | epsilon_dot += dt * ( (p_mol - pressure_ext)*volume |
183 | < | / (tau_relax*tau_relax * kB * temp) ); |
182 | > | void ExtendedSystem::ConstantStress( double dt, |
183 | > | double ke, |
184 | > | double p_tensor[9] ) { |
185 | ||
186 | < | // determine the change in cell volume |
187 | < | scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0)); |
186 | > | double oldBox[3]; |
187 | > | double newBox[3]; |
188 | > | const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
189 | > | const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1 |
190 | > | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to |
191 | > | // amu*Ang^2*fs^-2/K |
192 | ||
193 | < | volume = volume * pow(scale, 3.0); |
193 | > | int i; |
194 | > | double p_ext, zetaScale, epsilonScale, scale, NkBT, ke_temp; |
195 | > | double pX_ext, pY_ext, pZ_ext; |
196 | > | double volume, p_mol; |
197 | > | double vx, vy, vz, jx, jy, jz; |
198 | > | DirectionalAtom* dAtom; |
199 | ||
200 | < | // perform affine transform to update positions with volume fluctuations |
201 | < | affine_transform( scale ); |
200 | > | if (this->NPTready()) { |
201 | > | atoms = entry_plug->atoms; |
202 | > | |
203 | > | p_ext = targetPressure * p_units; |
204 | ||
205 | < | // save old lengths and update box size |
206 | < | boxx_old = boxx; |
207 | < | boxy_old = boxy; |
208 | < | boxz_old = boxz; |
209 | < | |
210 | < | boxx = boxx_old*scale; |
211 | < | boxy = boxy_old*scale; |
212 | < | boxz = boxz_old*scale; |
213 | < | |
214 | < | epsilonScale = epsilon_dot * dt; |
205 | > | pX_ext = p_ext / 3.0; |
206 | > | pY_ext = p_ext / 3.0; |
207 | > | pZ_ext = p_ext / 3.0; |
208 | > | |
209 | > | entry_plug->getBox(oldBox); |
210 | > | volume = oldBox[0]*oldBox[1]*oldBox[2]; |
211 | > | |
212 | > | ke_temp = ke * e_convert; |
213 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
214 | > | |
215 | > | // propagate the strain rate |
216 | > | |
217 | > | epsilonDotX += dt * ((p_tensor[0] - pX_ext) * volume / |
218 | > | (tauBarostat*tauBarostat * kB * targetTemp) ); |
219 | > | epsilonDotY += dt * ((p_tensor[4] - pY_ext) * volume / |
220 | > | (tauBarostat*tauBarostat * kB * targetTemp) ); |
221 | > | epsilonDotZ += dt * ((p_tensor[8] - pZ_ext) * volume / |
222 | > | (tauBarostat*tauBarostat * kB * targetTemp) ); |
223 | > | |
224 | > | // determine the change in cell volume |
225 | ||
226 | < | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
227 | < | // qmass is set in the parameter file |
228 | < | zeta += dt * ( (ke_temp*2 - NkBT) / qmass ); |
229 | < | zetaScale = zeta * dt; |
230 | < | |
226 | > | //scale = pow( (1.0 + dt * 3.0 * (epsilonDot), (1.0 / 3.0)); |
227 | > | //std::cerr << "pmol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n"; |
228 | > | |
229 | > | newBox[0] = oldBox[0] * scale; |
230 | > | newBox[1] = oldBox[1] * scale; |
231 | > | newBox[2] = oldBox[2] * scale; |
232 | > | volume = newBox[0]*newBox[1]*newBox[2]; |
233 | > | |
234 | > | entry_plug->setBox(newBox); |
235 | > | |
236 | > | // perform affine transform to update positions with volume fluctuations |
237 | > | this->AffineTransform( oldBox, newBox ); |
238 | > | |
239 | > | epsilonScale = epsilonDot * dt; |
240 | > | |
241 | > | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
242 | > | // qmass is set in the parameter file |
243 | > | |
244 | > | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
245 | > | zetaScale = zeta * dt; |
246 | > | |
247 | > | //std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale << "\n"; |
248 | > | |
249 | // apply barostating and thermostating to velocities and angular momenta | |
250 | < | |
251 | < | for (i=0; i < nmol; i++) { |
252 | < | |
253 | < | vxi = vx(i)*epsilonScale; |
254 | < | vyi = vy(i)*epsilonScale; |
255 | < | vzi = vz(i)*epsilonScale; |
256 | < | vxi = vxi + vx(i)*zetaScale; |
257 | < | vyi = vyi + vy(i)*zetaScale; |
258 | < | vzi = vzi + vz(i)*zetaScale; |
259 | < | jxi = jx(i)*zetaScale; |
260 | < | jyi = jy(i)*zetaScale; |
261 | < | jzi = jz(i)*zetaScale; |
262 | < | |
263 | < | vx(i) = vx(i) - vxi; |
264 | < | vy(i) = vy(i) - vyi; |
265 | < | vz(i) = vz(i) - vzi; |
266 | < | jx(i) = jx(i) - jxi; |
267 | < | jy(i) = jy(i) - jyi; |
268 | < | jz(i) = jz(i) - jzi; |
250 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
251 | > | |
252 | > | vx = atoms[i]->get_vx(); |
253 | > | vy = atoms[i]->get_vy(); |
254 | > | vz = atoms[i]->get_vz(); |
255 | > | |
256 | > | atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale)); |
257 | > | atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale)); |
258 | > | atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale)); |
259 | > | } |
260 | > | if( entry_plug->n_oriented ){ |
261 | > | |
262 | > | for( i=0; i < entry_plug->n_atoms; i++ ){ |
263 | > | |
264 | > | if( atoms[i]->isDirectional() ){ |
265 | > | |
266 | > | dAtom = (DirectionalAtom *)atoms[i]; |
267 | > | |
268 | > | jx = dAtom->getJx(); |
269 | > | jy = dAtom->getJy(); |
270 | > | jz = dAtom->getJz(); |
271 | > | |
272 | > | dAtom->setJx( jx * (1.0 - zetaScale)); |
273 | > | dAtom->setJy( jy * (1.0 - zetaScale)); |
274 | > | dAtom->setJz( jz * (1.0 - zetaScale)); |
275 | > | } |
276 | > | } |
277 | > | } |
278 | } | |
132 | – | |
133 | – | |
134 | – | |
279 | } | |
280 | ||
281 | < | void ExtendedSystem::affine_transform( double scale ){ |
281 | > | void ExtendedSystem::AffineTransform( double oldBox[3], double newBox[3] ){ |
282 | ||
283 | int i; | |
284 | < | double boxx_old, boxy_old, boxz_old, percentScale; |
285 | < | double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi; |
286 | < | double[3] r; |
284 | > | double r[3]; |
285 | > | double boxNum[3]; |
286 | > | double percentScale[3]; |
287 | > | double delta[3]; |
288 | > | double rxi, ryi, rzi; |
289 | > | |
290 | > | molecules = entry_plug->molecules; |
291 | ||
292 | // first determine the scaling factor from the box size change | |
293 | < | percentScale = (boxx - boxx_old)/boxx_old; |
293 | > | percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0]; |
294 | > | percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1]; |
295 | > | percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2]; |
296 | ||
297 | + | for (i=0; i < entry_plug->n_mol; i++) { |
298 | + | |
299 | + | molecules[i].getCOM(r); |
300 | ||
301 | < | for (i=0; i < nMols; i++) { |
301 | > | // find the minimum image coordinates of the molecular centers of mass: |
302 | ||
303 | < | molecules[i]->getCOM(r); |
304 | < | |
152 | < | // find the minimum image coordinates of the molecular centers of mass: |
153 | < | |
154 | < | |
155 | < | boxx_num = boxx_old*copysign(1.0,r[0])*(double)(int)(fabs(r[0]/boxx_old)+0.5); |
303 | > | boxNum[0] = oldBox[0] * copysign(1.0,r[0]) * |
304 | > | (double)(int)(fabs(r[0]/oldBox[0]) + 0.5); |
305 | ||
306 | < | boxx_num = boxx_old*dsign(1.0d0,rx(i))*int(abs(rx(i)/boxx_old)+0.5d0); |
307 | < | boxy_num = boxy_old*dsign(1.0d0,ry(i))*int(abs(ry(i)/boxy_old)+0.5d0); |
159 | < | boxz_num = boxz_old*dsign(1.0d0,rz(i))*int(abs(rz(i)/boxz_old)+0.5d0); |
306 | > | boxNum[1] = oldBox[1] * copysign(1.0,r[1]) * |
307 | > | (double)(int)(fabs(r[1]/oldBox[1]) + 0.5); |
308 | ||
309 | < | rxi = rx(i) - boxx_num; |
310 | < | ryi = ry(i) - boxy_num; |
163 | < | rzi = rz(i) - boxz_num; |
309 | > | boxNum[2] = oldBox[2] * copysign(1.0,r[2]) * |
310 | > | (double)(int)(fabs(r[2]/oldBox[2]) + 0.5); |
311 | ||
312 | + | rxi = r[0] - boxNum[0]; |
313 | + | ryi = r[1] - boxNum[1]; |
314 | + | rzi = r[2] - boxNum[2]; |
315 | + | |
316 | // update the minimum image coordinates using the scaling factor | |
317 | < | rxi = rxi + rxi*percentScale; |
318 | < | ryi = ryi + ryi*percentScale; |
319 | < | rzi = rzi + rzi*percentScale; |
317 | > | rxi += rxi*percentScale[0]; |
318 | > | ryi += ryi*percentScale[1]; |
319 | > | rzi += rzi*percentScale[2]; |
320 | ||
321 | < | rx(i) = rxi + boxx_num; |
322 | < | ry(i) = ryi + boxy_num; |
323 | < | rz(i) = rzi + boxz_num; |
321 | > | delta[0] = r[0] - (rxi + boxNum[0]); |
322 | > | delta[1] = r[1] - (ryi + boxNum[1]); |
323 | > | delta[2] = r[2] - (rzi + boxNum[2]); |
324 | > | |
325 | > | molecules[i].moveCOM(delta); |
326 | } | |
327 | } | |
328 | + | |
329 | + | short int ExtendedSystem::NVTready() { |
330 | + | const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
331 | + | double NkBT; |
332 | + | |
333 | + | if (!have_target_temp) { |
334 | + | sprintf( painCave.errMsg, |
335 | + | "ExtendedSystem error: You can't use NVT without a targetTemp!\n" |
336 | + | ); |
337 | + | painCave.isFatal = 1; |
338 | + | simError(); |
339 | + | return -1; |
340 | + | } |
341 | + | |
342 | + | if (!have_qmass) { |
343 | + | if (have_tau_thermostat) { |
344 | + | |
345 | + | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
346 | + | std::cerr << "Setting qMass = " << tauThermostat * NkBT << "\n"; |
347 | + | this->setQmass(tauThermostat * NkBT); |
348 | + | |
349 | + | } else { |
350 | + | sprintf( painCave.errMsg, |
351 | + | "ExtendedSystem error: If you use the constant temperature\n" |
352 | + | " ensemble, you must set either tauThermostat or qMass.\n"); |
353 | + | painCave.isFatal = 1; |
354 | + | simError(); |
355 | + | } |
356 | + | } |
357 | + | |
358 | + | return 1; |
359 | + | } |
360 | + | |
361 | + | short int ExtendedSystem::NPTready() { |
362 | + | const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K |
363 | + | double NkBT; |
364 | + | |
365 | + | if (!have_target_temp) { |
366 | + | sprintf( painCave.errMsg, |
367 | + | "ExtendedSystem error: You can't use NPT without a targetTemp!\n" |
368 | + | ); |
369 | + | painCave.isFatal = 1; |
370 | + | simError(); |
371 | + | return -1; |
372 | + | } |
373 | + | |
374 | + | if (!have_target_pressure) { |
375 | + | sprintf( painCave.errMsg, |
376 | + | "ExtendedSystem error: You can't use NPT without a targetPressure!\n" |
377 | + | ); |
378 | + | painCave.isFatal = 1; |
379 | + | simError(); |
380 | + | return -1; |
381 | + | } |
382 | + | |
383 | + | if (!have_tau_barostat) { |
384 | + | sprintf( painCave.errMsg, |
385 | + | "ExtendedSystem error: If you use the NPT\n" |
386 | + | " ensemble, you must set tauBarostat.\n"); |
387 | + | painCave.isFatal = 1; |
388 | + | simError(); |
389 | + | } |
390 | + | |
391 | + | if (!have_qmass) { |
392 | + | if (have_tau_thermostat) { |
393 | + | |
394 | + | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
395 | + | std::cerr << "Setting qMass = " << tauThermostat * NkBT << "\n"; |
396 | + | this->setQmass(tauThermostat * NkBT); |
397 | + | |
398 | + | } else { |
399 | + | sprintf( painCave.errMsg, |
400 | + | "ExtendedSystem error: If you use the NPT\n" |
401 | + | " ensemble, you must set either tauThermostat or qMass.\n"); |
402 | + | painCave.isFatal = 1; |
403 | + | simError(); |
404 | + | } |
405 | + | } |
406 | + | return 1; |
407 | + | } |
408 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |