# | Line 4 | Line 4 | |
---|---|---|
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; |
13 | < | nAtoms = info.n_atoms; |
14 | < | atoms = info.atoms; |
15 | < | nMols = info.n_mol; |
16 | < | 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 | ||
22 | – | ExtendedSystem::~ExtendedSystem() { |
23 | – | } |
24 | – | |
25 | – | |
27 | void ExtendedSystem::NoseHooverNVT( double dt, double ke ){ | |
28 | ||
29 | // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 | |
# | Line 33 | Line 34 | void ExtendedSystem::NoseHooverNVT( double dt, double | |
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 | < | |
37 | < | ke_temp = ke * e_convert; |
38 | < | NkBT = (double)getNDF() * kB * targetTemp; |
37 | > | DirectionalAtom* dAtom; |
38 | ||
39 | < | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
41 | < | // qmass is set in the parameter file |
39 | > | if (this->NVTready()) { |
40 | ||
41 | < | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
44 | < | zetaScale = zeta * dt; |
45 | < | |
46 | < | // perform thermostat scaling on linear velocities and angular momentum |
47 | < | for(i = 0; i < n_atoms; i++){ |
41 | > | atoms = entry_plug->atoms; |
42 | ||
43 | < | vx = atoms[i]->get_vx(); |
44 | < | vy = atoms[i]->get_vy(); |
51 | < | vz = atoms[i]->get_vz(); |
43 | > | ke_temp = ke * e_convert; |
44 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
45 | ||
46 | < | atoms[i]->set_vx(vx * (1.0 - zetaScale)); |
47 | < | atoms[i]->set_vy(vy * (1.0 - zetaScale)); |
55 | < | atoms[i]->set_vz(vz * (1.0 - zetaScale)); |
56 | < | } |
57 | < | if( n_oriented ){ |
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 | < | for( i=0; i < n_atoms; i++ ){ |
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 | < | if( atoms[i]->isDirectional() ){ |
59 | < | |
60 | < | dAtom = (DirectionalAtom *)atoms[i]; |
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 | < | jx = dAtom->getJx(); |
71 | < | jy = dAtom->getJy(); |
72 | < | jz = dAtom->getJz(); |
73 | < | |
74 | < | dAtom->setJx(jx * (1.0 - zetaScale)); |
75 | < | dAtom->setJy(jy * (1.0 - zetaScale)); |
76 | < | dAtom->setJz(jz * (1.0 - zetaScale)); |
77 | < | } |
78 | < | } |
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::NoseHooverAndersonNPT( double dt, | |
89 | double ke, | |
90 | < | double p_int ) { |
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 | |
# | Line 89 | Line 99 | void ExtendedSystem::NoseHooverAndersonNPT( double dt, | |
99 | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to | |
100 | // amu*Ang^2*fs^-2/K | |
101 | ||
102 | < | double p_ext; |
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 | < | p_ext = targetPressure * p_units; |
109 | < | p_mol = p_int * p_units; |
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 | ||
97 | – | getBox(oldBox); |
181 | ||
182 | < | volume = oldBox[0]*oldBox[1]*oldBox[2]; |
182 | > | void ExtendedSystem::ConstantStress( double dt, |
183 | > | double ke, |
184 | > | double p_tensor[9] ) { |
185 | ||
186 | < | ke_temp = ke * e_convert; |
187 | < | NkBT = (double)getNDF() * kB * targetTemp; |
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 | < | // propogate the strain rate |
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 | < | epsilonDot += dt * ((p_mol - p_ext) * volume / |
201 | < | (tauRelax*tauRelax * kB * targetTemp) ); |
200 | > | if (this->NPTready()) { |
201 | > | atoms = entry_plug->atoms; |
202 | > | |
203 | > | p_ext = targetPressure * p_units; |
204 | ||
205 | < | // determine the change in cell volume |
206 | < | scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0)); |
207 | < | |
208 | < | newBox[0] = oldBox[0] * scale; |
209 | < | newBox[1] = oldBox[1] * scale; |
210 | < | newBox[2] = oldBox[2] * scale; |
211 | < | volume = newBox[0]*newBox[1]*newBox[2]; |
212 | < | |
213 | < | // perform affine transform to update positions with volume fluctuations |
214 | < | this->AffineTransform( oldBox, newBox ); |
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 | < | epsilonScale = epsilonDot * dt; |
227 | < | |
122 | < | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin |
123 | < | // qmass is set in the parameter file |
124 | < | |
125 | < | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
126 | < | zetaScale = zeta * dt; |
127 | < | |
128 | < | // apply barostating and thermostating to velocities and angular momenta |
129 | < | for(i = 0; i < n_atoms; i++){ |
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 | < | vx = atoms[i]->get_vx(); |
230 | < | vy = atoms[i]->get_vy(); |
231 | < | vz = atoms[i]->get_vz(); |
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 | < | atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale)); |
136 | < | atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale)); |
137 | < | atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale)); |
138 | < | } |
139 | < | if( n_oriented ){ |
234 | > | entry_plug->setBox(newBox); |
235 | ||
236 | < | for( i=0; i < n_atoms; i++ ){ |
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 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
251 | ||
252 | < | if( atoms[i]->isDirectional() ){ |
253 | < | |
254 | < | dAtom = (DirectionalAtom *)atoms[i]; |
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 | < | jx = dAtom->getJx(); |
265 | < | jy = dAtom->getJy(); |
266 | < | jz = dAtom->getJz(); |
267 | < | |
268 | < | dAtom->setJx( jx * (1.0 - zetaScale)); |
269 | < | dAtom->setJy( jy * (1.0 - zetaScale)); |
270 | < | dAtom->setJz( jz * (1.0 - zetaScale)); |
271 | < | } |
272 | < | } |
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 | } | |
279 | } | |
280 | ||
# | Line 162 | Line 284 | void ExtendedSystem::AffineTransform( double oldBox[3] | |
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[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 < nMols; i++) { |
297 | > | for (i=0; i < entry_plug->n_mol; i++) { |
298 | ||
299 | < | molecules[i]->getCOM(r); |
300 | < | |
299 | > | molecules[i].getCOM(r); |
300 | > | |
301 | // find the minimum image coordinates of the molecular centers of mass: | |
302 | ||
303 | boxNum[0] = oldBox[0] * copysign(1.0,r[0]) * | |
# | Line 193 | Line 318 | void ExtendedSystem::AffineTransform( double oldBox[3] | |
318 | ryi += ryi*percentScale[1]; | |
319 | rzi += rzi*percentScale[2]; | |
320 | ||
321 | < | r[0] = rxi + boxNum[0]; |
322 | < | r[1] = ryi + boxNum[1]; |
323 | < | r[2] = rzi + boxNum[2]; |
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(r); |
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 |