1 |
#include <math.h>
|
2 |
|
3 |
#include "primitives/Atom.hpp"
|
4 |
#include "primitives/SRI.hpp"
|
5 |
#include "primitives/AbstractClasses.hpp"
|
6 |
#include "brains/SimInfo.hpp"
|
7 |
#include "UseTheForce/ForceFields.hpp"
|
8 |
#include "brains/Thermo.hpp"
|
9 |
#include "io/ReadWrite.hpp"
|
10 |
#include "integrators/Integrator.hpp"
|
11 |
#include "utils/simError.h"
|
12 |
|
13 |
#ifdef IS_MPI
|
14 |
#include "brains/mpiSimulation.hpp"
|
15 |
#endif
|
16 |
|
17 |
|
18 |
// Basic isotropic thermostating and barostating via the Melchionna
|
19 |
// modification of the Hoover algorithm:
|
20 |
//
|
21 |
// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
|
22 |
// Molec. Phys., 78, 533.
|
23 |
//
|
24 |
// and
|
25 |
//
|
26 |
// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
|
27 |
|
28 |
template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
|
29 |
T( theInfo, the_ff )
|
30 |
{
|
31 |
GenericData* data;
|
32 |
DoubleGenericData * chiValue;
|
33 |
DoubleGenericData * integralOfChidtValue;
|
34 |
|
35 |
chiValue = NULL;
|
36 |
integralOfChidtValue = NULL;
|
37 |
|
38 |
chi = 0.0;
|
39 |
integralOfChidt = 0.0;
|
40 |
have_tau_thermostat = 0;
|
41 |
have_tau_barostat = 0;
|
42 |
have_target_temp = 0;
|
43 |
have_target_pressure = 0;
|
44 |
have_chi_tolerance = 0;
|
45 |
have_eta_tolerance = 0;
|
46 |
have_pos_iter_tolerance = 0;
|
47 |
|
48 |
// retrieve chi and integralOfChidt from simInfo
|
49 |
data = info->getProperty(CHIVALUE_ID);
|
50 |
if(data){
|
51 |
chiValue = dynamic_cast<DoubleGenericData*>(data);
|
52 |
}
|
53 |
|
54 |
data = info->getProperty(INTEGRALOFCHIDT_ID);
|
55 |
if(data){
|
56 |
integralOfChidtValue = dynamic_cast<DoubleGenericData*>(data);
|
57 |
}
|
58 |
|
59 |
// chi and integralOfChidt should appear by pair
|
60 |
if(chiValue && integralOfChidtValue){
|
61 |
chi = chiValue->getData();
|
62 |
integralOfChidt = integralOfChidtValue->getData();
|
63 |
}
|
64 |
|
65 |
oldPos = new double[3*integrableObjects.size()];
|
66 |
oldVel = new double[3*integrableObjects.size()];
|
67 |
oldJi = new double[3*integrableObjects.size()];
|
68 |
|
69 |
}
|
70 |
|
71 |
template<typename T> NPT<T>::~NPT() {
|
72 |
delete[] oldPos;
|
73 |
delete[] oldVel;
|
74 |
delete[] oldJi;
|
75 |
}
|
76 |
|
77 |
template<typename T> void NPT<T>::moveA() {
|
78 |
|
79 |
//new version of NPT
|
80 |
int i, j, k;
|
81 |
Vector3d Tb, ji;
|
82 |
double mass;
|
83 |
Vector3d vel;
|
84 |
Vector3d pos;
|
85 |
Vector3d frc;
|
86 |
Vector3d sc;
|
87 |
Vector3d COM;
|
88 |
|
89 |
instaTemp = tStats->getTemperature();
|
90 |
tStats->getPressureTensor( press );
|
91 |
instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
|
92 |
instaVol = tStats->getVolume();
|
93 |
|
94 |
tStats->getCOM(COM);
|
95 |
|
96 |
//evolve velocity half step
|
97 |
|
98 |
calcVelScale();
|
99 |
|
100 |
for( i=0; i<integrableObjects.size(); i++ ){
|
101 |
|
102 |
vel = integrableObjects[i]->getVel();
|
103 |
integrableObjects[i]->getFrc( frc );
|
104 |
|
105 |
mass = integrableObjects[i]->getMass();
|
106 |
|
107 |
getVelScaleA( sc, vel );
|
108 |
|
109 |
for (j=0; j < 3; j++) {
|
110 |
|
111 |
// velocity half step (use chi from previous step here):
|
112 |
vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
|
113 |
|
114 |
}
|
115 |
|
116 |
integrableObjects[i]->setVel( vel );
|
117 |
|
118 |
if( integrableObjects[i]->isDirectional() ){
|
119 |
|
120 |
// get and convert the torque to body frame
|
121 |
|
122 |
Tb = integrableObjects[i]->getTrq();
|
123 |
integrableObjects[i]->lab2Body( Tb );
|
124 |
|
125 |
// get the angular momentum, and propagate a half step
|
126 |
|
127 |
ji = integrableObjects[i]->getJ();
|
128 |
|
129 |
for (j=0; j < 3; j++)
|
130 |
ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
|
131 |
|
132 |
this->rotationPropagation( integrableObjects[i], ji );
|
133 |
|
134 |
integrableObjects[i]->setJ( ji );
|
135 |
}
|
136 |
}
|
137 |
|
138 |
// evolve chi and eta half step
|
139 |
|
140 |
evolveChiA();
|
141 |
evolveEtaA();
|
142 |
|
143 |
//calculate the integral of chidt
|
144 |
integralOfChidt += dt2*chi;
|
145 |
|
146 |
//save the old positions
|
147 |
for(i = 0; i < integrableObjects.size(); i++){
|
148 |
pos = integrableObjects[i]->getPos();
|
149 |
for(j = 0; j < 3; j++)
|
150 |
oldPos[i*3 + j] = pos[j];
|
151 |
}
|
152 |
|
153 |
//the first estimation of r(t+dt) is equal to r(t)
|
154 |
|
155 |
for(k = 0; k < 5; k ++){
|
156 |
|
157 |
for(i =0 ; i < integrableObjects.size(); i++){
|
158 |
|
159 |
vel = integrableObjects[i]->getVel();
|
160 |
pos = integrableObjects[i]->getPos();
|
161 |
|
162 |
this->getPosScale( pos, COM, i, sc );
|
163 |
|
164 |
for(j = 0; j < 3; j++)
|
165 |
pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
|
166 |
|
167 |
integrableObjects[i]->setPos( pos );
|
168 |
}
|
169 |
|
170 |
if(nConstrained)
|
171 |
constrainA();
|
172 |
}
|
173 |
|
174 |
|
175 |
// Scale the box after all the positions have been moved:
|
176 |
|
177 |
this->scaleSimBox();
|
178 |
}
|
179 |
|
180 |
template<typename T> void NPT<T>::moveB( void ){
|
181 |
|
182 |
//new version of NPT
|
183 |
int i, j, k;
|
184 |
Vector3d Tb;
|
185 |
Vector3d ji;
|
186 |
Vector3d sc;
|
187 |
Vector3d vel;
|
188 |
Vector3d frc;
|
189 |
double mass;
|
190 |
|
191 |
// Set things up for the iteration:
|
192 |
|
193 |
for( i=0; i<integrableObjects.size(); i++ ){
|
194 |
|
195 |
vel = integrableObjects[i]->getVel();
|
196 |
|
197 |
for (j=0; j < 3; j++)
|
198 |
oldVel[3*i + j] = vel[j];
|
199 |
|
200 |
if( integrableObjects[i]->isDirectional() ){
|
201 |
|
202 |
ji = integrableObjects[i]->getJ();
|
203 |
|
204 |
for (j=0; j < 3; j++)
|
205 |
oldJi[3*i + j] = ji[j];
|
206 |
|
207 |
}
|
208 |
}
|
209 |
|
210 |
// do the iteration:
|
211 |
|
212 |
instaVol = tStats->getVolume();
|
213 |
|
214 |
for (k=0; k < 4; k++) {
|
215 |
|
216 |
instaTemp = tStats->getTemperature();
|
217 |
instaPress = tStats->getPressure();
|
218 |
|
219 |
// evolve chi another half step using the temperature at t + dt/2
|
220 |
|
221 |
this->evolveChiB();
|
222 |
this->evolveEtaB();
|
223 |
this->calcVelScale();
|
224 |
|
225 |
for( i=0; i<integrableObjects.size(); i++ ){
|
226 |
|
227 |
integrableObjects[i]->getFrc( frc );
|
228 |
vel = integrableObjects[i]->getVel();
|
229 |
|
230 |
mass = integrableObjects[i]->getMass();
|
231 |
|
232 |
getVelScaleB( sc, i );
|
233 |
|
234 |
// velocity half step
|
235 |
for (j=0; j < 3; j++)
|
236 |
vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
|
237 |
|
238 |
integrableObjects[i]->setVel( vel );
|
239 |
|
240 |
if( integrableObjects[i]->isDirectional() ){
|
241 |
|
242 |
// get and convert the torque to body frame
|
243 |
|
244 |
Tb = integrableObjects[i]->getTrq();
|
245 |
integrableObjects[i]->lab2Body( Tb );
|
246 |
|
247 |
for (j=0; j < 3; j++)
|
248 |
ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
|
249 |
|
250 |
integrableObjects[i]->setJ( ji );
|
251 |
}
|
252 |
}
|
253 |
|
254 |
if(nConstrained)
|
255 |
constrainB();
|
256 |
|
257 |
if ( this->chiConverged() && this->etaConverged() ) break;
|
258 |
}
|
259 |
|
260 |
//calculate integral of chida
|
261 |
integralOfChidt += dt2*chi;
|
262 |
|
263 |
|
264 |
}
|
265 |
|
266 |
template<typename T> void NPT<T>::resetIntegrator() {
|
267 |
chi = 0.0;
|
268 |
T::resetIntegrator();
|
269 |
}
|
270 |
|
271 |
template<typename T> void NPT<T>::evolveChiA() {
|
272 |
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
|
273 |
oldChi = chi;
|
274 |
}
|
275 |
|
276 |
template<typename T> void NPT<T>::evolveChiB() {
|
277 |
|
278 |
prevChi = chi;
|
279 |
chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
|
280 |
}
|
281 |
|
282 |
template<typename T> bool NPT<T>::chiConverged() {
|
283 |
|
284 |
return ( fabs( prevChi - chi ) <= chiTolerance );
|
285 |
}
|
286 |
|
287 |
template<typename T> int NPT<T>::readyCheck() {
|
288 |
|
289 |
//check parent's readyCheck() first
|
290 |
if (T::readyCheck() == -1)
|
291 |
return -1;
|
292 |
|
293 |
// First check to see if we have a target temperature.
|
294 |
// Not having one is fatal.
|
295 |
|
296 |
if (!have_target_temp) {
|
297 |
sprintf( painCave.errMsg,
|
298 |
"NPT error: You can't use the NPT integrator\n"
|
299 |
" without a targetTemp!\n"
|
300 |
);
|
301 |
painCave.isFatal = 1;
|
302 |
simError();
|
303 |
return -1;
|
304 |
}
|
305 |
|
306 |
if (!have_target_pressure) {
|
307 |
sprintf( painCave.errMsg,
|
308 |
"NPT error: You can't use the NPT integrator\n"
|
309 |
" without a targetPressure!\n"
|
310 |
);
|
311 |
painCave.isFatal = 1;
|
312 |
simError();
|
313 |
return -1;
|
314 |
}
|
315 |
|
316 |
// We must set tauThermostat.
|
317 |
|
318 |
if (!have_tau_thermostat) {
|
319 |
sprintf( painCave.errMsg,
|
320 |
"NPT error: If you use the NPT\n"
|
321 |
" integrator, you must set tauThermostat.\n");
|
322 |
painCave.isFatal = 1;
|
323 |
simError();
|
324 |
return -1;
|
325 |
}
|
326 |
|
327 |
// We must set tauBarostat.
|
328 |
|
329 |
if (!have_tau_barostat) {
|
330 |
sprintf( painCave.errMsg,
|
331 |
"If you use the NPT integrator, you must set tauBarostat.\n");
|
332 |
painCave.severity = OOPSE_ERROR;
|
333 |
painCave.isFatal = 1;
|
334 |
simError();
|
335 |
return -1;
|
336 |
}
|
337 |
|
338 |
if (!have_chi_tolerance) {
|
339 |
sprintf( painCave.errMsg,
|
340 |
"Setting chi tolerance to 1e-6 in NPT integrator\n");
|
341 |
chiTolerance = 1e-6;
|
342 |
have_chi_tolerance = 1;
|
343 |
painCave.severity = OOPSE_INFO;
|
344 |
painCave.isFatal = 0;
|
345 |
simError();
|
346 |
}
|
347 |
|
348 |
if (!have_eta_tolerance) {
|
349 |
sprintf( painCave.errMsg,
|
350 |
"Setting eta tolerance to 1e-6 in NPT integrator");
|
351 |
etaTolerance = 1e-6;
|
352 |
have_eta_tolerance = 1;
|
353 |
painCave.severity = OOPSE_INFO;
|
354 |
painCave.isFatal = 0;
|
355 |
simError();
|
356 |
}
|
357 |
|
358 |
// We need NkBT a lot, so just set it here: This is the RAW number
|
359 |
// of integrableObjects, so no subtraction or addition of constraints or
|
360 |
// orientational degrees of freedom:
|
361 |
|
362 |
NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
|
363 |
|
364 |
// fkBT is used because the thermostat operates on more degrees of freedom
|
365 |
// than the barostat (when there are particles with orientational degrees
|
366 |
// of freedom).
|
367 |
|
368 |
fkBT = (double)(info->getNDF()) * kB * targetTemp;
|
369 |
|
370 |
tt2 = tauThermostat * tauThermostat;
|
371 |
tb2 = tauBarostat * tauBarostat;
|
372 |
|
373 |
return 1;
|
374 |
}
|