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; |
50 |
|
|
51 |
|
zetaScale = zeta * dt; |
52 |
|
|
53 |
< |
std::cerr << "zetaScale = " << zetaScale << "\n"; |
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++){ |
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 |
109 |
|
atoms = entry_plug->atoms; |
110 |
|
|
111 |
|
p_ext = targetPressure * p_units; |
112 |
< |
p_mol = p_int * p_units; |
113 |
< |
|
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 |
< |
// propogate the strain rate |
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; |
144 |
|
zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); |
145 |
|
zetaScale = zeta * dt; |
146 |
|
|
147 |
< |
std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale << "\n"; |
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 |
> |
|
181 |
> |
|
182 |
> |
void ExtendedSystem::ConstantStress( double dt, |
183 |
> |
double ke, |
184 |
> |
double p_tensor[9] ) { |
185 |
> |
|
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 |
> |
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 |
> |
if (this->NPTready()) { |
201 |
> |
atoms = entry_plug->atoms; |
202 |
|
|
203 |
+ |
p_ext = targetPressure * p_units; |
204 |
+ |
|
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 |
+ |
//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 |
|
for(i = 0; i < entry_plug->n_atoms; i++){ |
251 |
|
|
355 |
|
} |
356 |
|
} |
357 |
|
|
358 |
< |
return 0; |
358 |
> |
return 1; |
359 |
|
} |
360 |
|
|
361 |
|
short int ExtendedSystem::NPTready() { |
403 |
|
simError(); |
404 |
|
} |
405 |
|
} |
406 |
< |
return 0; |
406 |
> |
return 1; |
407 |
|
} |
408 |
|
|