ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPTfm.cpp
Revision: 767
Committed: Tue Sep 16 20:02:11 2003 UTC (21 years ago) by tim
File size: 20669 byte(s)
Log Message: 
fixed ecr grow in SimInfo

fixed conserved quantity in NPT (Still some small bug)

NPTi appears very stable.

File Contents

# User Rev Content
1 gezelter 617 #include <cmath>
2 gezelter 606 #include "Atom.hpp"
3     #include "Molecule.hpp"
4     #include "SRI.hpp"
5     #include "AbstractClasses.hpp"
6     #include "SimInfo.hpp"
7     #include "ForceFields.hpp"
8     #include "Thermo.hpp"
9     #include "ReadWrite.hpp"
10     #include "Integrator.hpp"
11     #include "simError.h"
12    
13 mmeineke 746
14 gezelter 606 // Basic non-isotropic thermostating and barostating via the Melchionna
15     // modification of the Hoover algorithm:
16     //
17     // Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
18     // Molec. Phys., 78, 533.
19     //
20     // and
21     //
22     // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
23    
24     // The NPTfm variant scales the molecular center-of-mass coordinates
25     // instead of the atomic coordinates
26    
27 tim 645 template<typename T> NPTfm<T>::NPTfm ( SimInfo *theInfo, ForceFields* the_ff):
28     T( theInfo, the_ff )
29 gezelter 606 {
30     int i, j;
31     chi = 0.0;
32 tim 763 integralOfChidt = 0.0;
33    
34 gezelter 606 for(i = 0; i < 3; i++)
35     for (j = 0; j < 3; j++)
36     eta[i][j] = 0.0;
37    
38     have_tau_thermostat = 0;
39     have_tau_barostat = 0;
40     have_target_temp = 0;
41     have_target_pressure = 0;
42     }
43    
44 tim 645 template<typename T> void NPTfm<T>::moveA() {
45 gezelter 606
46 tim 767 // int i, j, k;
47     // DirectionalAtom* dAtom;
48     // double Tb[3], ji[3];
49     // double A[3][3], I[3][3];
50     // double angle, mass;
51     // double vel[3], pos[3], frc[3];
52 gezelter 606
53 tim 767 // double rj[3];
54     // double instaTemp, instaPress, instaVol;
55     // double tt2, tb2;
56     // double sc[3];
57     // double eta2ij, smallScale, bigScale, offDiagMax;
58     // double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
59 gezelter 606
60 tim 767 // int nInMol;
61     // double rc[3];
62 gezelter 606
63 tim 767 // /*
64     // nMols = info->n_mol;
65     // myMolecules = info->molecules;
66 gezelter 606
67 tim 767 // tt2 = tauThermostat * tauThermostat;
68     // tb2 = tauBarostat * tauBarostat;
69 gezelter 606
70 tim 767 // instaTemp = tStats->getTemperature();
71     // tStats->getPressureTensor(press);
72     // instaVol = tStats->getVolume();
73 gezelter 606
74 tim 767 // // first evolve chi a half step
75 gezelter 606
76 tim 767 // chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
77 gezelter 606
78 tim 767 // for (i = 0; i < 3; i++ ) {
79     // for (j = 0; j < 3; j++ ) {
80     // if (i == j) {
81 gezelter 606
82 tim 767 // eta[i][j] += dt2 * instaVol *
83     // (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
84 gezelter 606
85 tim 767 // vScale[i][j] = eta[i][j] + chi;
86 gezelter 606
87 tim 767 // } else {
88 gezelter 606
89 tim 767 // eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
90 gezelter 606
91 tim 767 // vScale[i][j] = eta[i][j];
92 gezelter 606
93 tim 767 // }
94     // }
95     // }
96 gezelter 606
97    
98 tim 767 // for (i = 0; i < nMols; i++) {
99 gezelter 606
100 tim 767 // myMolecules[i].getCOM(rc);
101 gezelter 606
102 tim 767 // nInMol = myMolecules[i].getNAtoms();
103     // myAtoms = myMolecules[i].getMyAtoms();
104 gezelter 606
105 tim 767 // // find the minimum image coordinates of the molecular centers of mass:
106 gezelter 606
107 tim 767 // info->wrapVector(rc);
108 gezelter 606
109 tim 767 // for( j=0; j< nInMol; j++ ){
110 gezelter 606
111 tim 767 // if(myAtoms[j] != NULL) {
112 gezelter 606
113 tim 767 // myAtoms[j]->getVel( vel );
114     // myAtoms[j]->getPos( pos );
115     // myAtoms[j]->getFrc( frc );
116 gezelter 606
117 tim 767 // mass = myAtoms[j]->getMass();
118 gezelter 606
119 tim 767 // // velocity half step
120 gezelter 606
121 tim 767 // info->matVecMul3( vScale, vel, sc );
122 gezelter 606
123 tim 767 // for (k = 0; k < 3; k++)
124     // vel[k] += dt2 * ((frc[k] / mass) * eConvert - sc[k]);
125 gezelter 606
126 tim 767 // myAtoms[j]->setVel( vel );
127 gezelter 606
128 tim 767 // // position whole step
129 gezelter 606
130 tim 767 // info->matVecMul3( eta, rc, sc );
131 gezelter 606
132 tim 767 // for (k = 0; k < 3; k++ )
133     // pos[k] += dt * (vel[k] + sc[k]);
134 gezelter 617
135 tim 767 // myAtoms[j]->setPos( pos );
136 gezelter 606
137 tim 767 // if( myAtoms[j]->isDirectional() ){
138 gezelter 606
139 tim 767 // dAtom = (DirectionalAtom *)myAtoms[j];
140 gezelter 606
141 tim 767 // // get and convert the torque to body frame
142 gezelter 606
143 tim 767 // dAtom->getTrq( Tb );
144     // dAtom->lab2Body( Tb );
145 gezelter 606
146 tim 767 // // get the angular momentum, and propagate a half step
147 gezelter 606
148 tim 767 // dAtom->getJ( ji );
149 gezelter 606
150 tim 767 // for (k=0; k < 3; k++)
151     // ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
152 gezelter 606
153 tim 767 // // use the angular velocities to propagate the rotation matrix a
154     // // full time step
155 gezelter 606
156 tim 767 // dAtom->getA(A);
157     // dAtom->getI(I);
158 gezelter 606
159 tim 767 // // rotate about the x-axis
160     // angle = dt2 * ji[0] / I[0][0];
161     // this->rotate( 1, 2, angle, ji, A );
162 gezelter 606
163 tim 767 // // rotate about the y-axis
164     // angle = dt2 * ji[1] / I[1][1];
165     // this->rotate( 2, 0, angle, ji, A );
166 gezelter 606
167 tim 767 // // rotate about the z-axis
168     // angle = dt * ji[2] / I[2][2];
169     // this->rotate( 0, 1, angle, ji, A);
170 gezelter 606
171 tim 767 // // rotate about the y-axis
172     // angle = dt2 * ji[1] / I[1][1];
173     // this->rotate( 2, 0, angle, ji, A );
174 gezelter 606
175 tim 767 // // rotate about the x-axis
176     // angle = dt2 * ji[0] / I[0][0];
177     // this->rotate( 1, 2, angle, ji, A );
178 gezelter 606
179 tim 767 // dAtom->setJ( ji );
180     // dAtom->setA( A );
181     // }
182     // }
183     // }
184     // }
185 gezelter 606
186 tim 767 // // Scale the box after all the positions have been moved:
187 gezelter 606
188 tim 767 // // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
189     // // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
190 gezelter 606
191    
192 tim 767 // bigScale = 1.0;
193     // smallScale = 1.0;
194     // offDiagMax = 0.0;
195 gezelter 617
196 tim 767 // for(i=0; i<3; i++){
197     // for(j=0; j<3; j++){
198 gezelter 606
199 tim 767 // // Calculate the matrix Product of the eta array (we only need
200     // // the ij element right now):
201 gezelter 606
202 tim 767 // eta2ij = 0.0;
203     // for(k=0; k<3; k++){
204     // eta2ij += eta[i][k] * eta[k][j];
205     // }
206 gezelter 606
207 tim 767 // scaleMat[i][j] = 0.0;
208     // // identity matrix (see above):
209     // if (i == j) scaleMat[i][j] = 1.0;
210     // // Taylor expansion for the exponential truncated at second order:
211     // scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij;
212 gezelter 617
213 tim 767 // if (i != j)
214     // if (fabs(scaleMat[i][j]) > offDiagMax)
215     // offDiagMax = fabs(scaleMat[i][j]);
216     // }
217     // if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
218     // if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
219     // }
220 gezelter 606
221 tim 767 // if ((bigScale > 1.1) || (smallScale < 0.9)) {
222     // sprintf( painCave.errMsg,
223     // "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
224     // " Check your tauBarostat, as it is probably too small!\n\n"
225     // " scaleMat = [%lf\t%lf\t%lf]\n"
226     // " [%lf\t%lf\t%lf]\n"
227     // " [%lf\t%lf\t%lf]\n",
228     // scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
229     // scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
230     // scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
231     // painCave.isFatal = 1;
232     // simError();
233     // } else if (offDiagMax > 0.1) {
234     // sprintf( painCave.errMsg,
235     // "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
236     // " Check your tauBarostat, as it is probably too small!\n\n"
237     // " scaleMat = [%lf\t%lf\t%lf]\n"
238     // " [%lf\t%lf\t%lf]\n"
239     // " [%lf\t%lf\t%lf]\n",
240     // scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
241     // scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
242     // scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
243     // painCave.isFatal = 1;
244     // simError();
245     // } else {
246     // info->getBoxM(hm);
247     // info->matMul3(hm, scaleMat, hmnew);
248     // info->setBoxM(hmnew);
249     // }
250     // */
251    
252     // tt2 = tauThermostat * tauThermostat;
253     // tb2 = tauBarostat * tauBarostat;
254    
255     // instaTemp = tStats->getTemperature();
256     // tStats->getPressureTensor(press);
257     // instaVol = tStats->getVolume();
258    
259     // tStats->getCOM(COM);
260    
261     // //calculate scale factor of veloity
262     // for (i = 0; i < 3; i++ ) {
263     // for (j = 0; j < 3; j++ ) {
264     // vScale[i][j] = eta[i][j];
265    
266     // if (i == j) {
267     // vScale[i][j] += chi;
268     // }
269     // }
270     // }
271    
272    
273     // for (i = 0; i < nMols; i++) {
274    
275     // myMolecules[i].getCOM(rc);
276    
277     // nInMol = myMolecules[i].getNAtoms();
278     // myAtoms = myMolecules[i].getMyAtoms();
279    
280    
281     // for( j=0; j< nInMol; j++ ){
282    
283     // if(myAtoms[j] != NULL) {
284    
285     // myAtoms[j]->getVel( vel );
286     // myAtoms[j]->getFrc( frc );
287    
288     // mass = myAtoms[j]->getMass();
289    
290     // // velocity half step
291    
292     // info->matVecMul3( vScale, vel, sc );
293    
294     // for (k = 0; k < 3; k++)
295     // vel[k] += dt2 * ((frc[k] / mass) * eConvert - sc[k]);
296    
297     // myAtoms[j]->setVel( vel );
298    
299     // if( myAtoms[j]->isDirectional() ){
300    
301     // dAtom = (DirectionalAtom *)myAtoms[j];
302    
303     // // get and convert the torque to body frame
304    
305     // dAtom->getTrq( Tb );
306     // dAtom->lab2Body( Tb );
307    
308     // // get the angular momentum, and propagate a half step
309    
310     // dAtom->getJ( ji );
311    
312     // for (k=0; k < 3; k++)
313     // ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
314    
315     // // use the angular velocities to propagate the rotation matrix a
316     // // full time step
317    
318     // dAtom->getA(A);
319     // dAtom->getI(I);
320    
321     // // rotate about the x-axis
322     // angle = dt2 * ji[0] / I[0][0];
323     // this->rotate( 1, 2, angle, ji, A );
324    
325     // // rotate about the y-axis
326     // angle = dt2 * ji[1] / I[1][1];
327     // this->rotate( 2, 0, angle, ji, A );
328    
329     // // rotate about the z-axis
330     // angle = dt * ji[2] / I[2][2];
331     // this->rotate( 0, 1, angle, ji, A);
332    
333     // // rotate about the y-axis
334     // angle = dt2 * ji[1] / I[1][1];
335     // this->rotate( 2, 0, angle, ji, A );
336    
337     // // rotate about the x-axis
338     // angle = dt2 * ji[0] / I[0][0];
339     // this->rotate( 1, 2, angle, ji, A );
340    
341     // dAtom->setJ( ji );
342     // dAtom->setA( A );
343     // }
344     // }
345     // }
346     // }
347    
348    
349     // // advance chi half step
350     // chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
351    
352     // //calculate the integral of chidt
353     // integralOfChidt += dt2*chi;
354    
355     // //advance eta half step
356     // for(i = 0; i < 3; i ++)
357     // for(j = 0; j < 3; j++){
358     // if( i == j)
359     // eta[i][j] += dt2 * instaVol *
360     // (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
361     // else
362     // eta[i][j] += dt2 * instaVol * press[i][j] / ( NkBT*tb2);
363     // }
364    
365     // //save the old positions
366     // for(i = 0; i < nAtoms; i++){
367     // atoms[i]->getPos(pos);
368     // for(j = 0; j < 3; j++)
369     // oldPos[i*3 + j] = pos[j];
370     // }
371    
372     // //the first estimation of r(t+dt) is equal to r(t)
373    
374     // for(k = 0; k < 4; k ++){
375    
376     // for(i =0 ; i < nAtoms; i++){
377    
378     // atoms[i]->getVel(vel);
379     // atoms[i]->getPos(pos);
380    
381     // for(j = 0; j < 3; j++)
382     // rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
383    
384     // info->matVecMul3( eta, rj, sc );
385    
386     // for(j = 0; j < 3; j++)
387     // pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
388    
389     // atoms[i]->setPos( pos );
390    
391     // }
392    
393     // }
394    
395    
396     // // Scale the box after all the positions have been moved:
397    
398     // // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
399     // // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
400    
401     // bigScale = 1.0;
402     // smallScale = 1.0;
403     // offDiagMax = 0.0;
404    
405     // for(i=0; i<3; i++){
406     // for(j=0; j<3; j++){
407    
408     // // Calculate the matrix Product of the eta array (we only need
409     // // the ij element right now):
410    
411     // eta2ij = 0.0;
412     // for(k=0; k<3; k++){
413     // eta2ij += eta[i][k] * eta[k][j];
414     // }
415    
416     // scaleMat[i][j] = 0.0;
417     // // identity matrix (see above):
418     // if (i == j) scaleMat[i][j] = 1.0;
419     // // Taylor expansion for the exponential truncated at second order:
420     // scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij;
421    
422     // if (i != j)
423     // if (fabs(scaleMat[i][j]) > offDiagMax)
424     // offDiagMax = fabs(scaleMat[i][j]);
425     // }
426    
427     // if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
428     // if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
429     // }
430    
431     // if ((bigScale > 1.1) || (smallScale < 0.9)) {
432     // sprintf( painCave.errMsg,
433     // "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
434     // " Check your tauBarostat, as it is probably too small!\n\n"
435     // " scaleMat = [%lf\t%lf\t%lf]\n"
436     // " [%lf\t%lf\t%lf]\n"
437     // " [%lf\t%lf\t%lf]\n",
438     // scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
439     // scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
440     // scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
441     // painCave.isFatal = 1;
442     // simError();
443     // } else if (offDiagMax > 0.1) {
444     // sprintf( painCave.errMsg,
445     // "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
446     // " Check your tauBarostat, as it is probably too small!\n\n"
447     // " scaleMat = [%lf\t%lf\t%lf]\n"
448     // " [%lf\t%lf\t%lf]\n"
449     // " [%lf\t%lf\t%lf]\n",
450     // scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
451     // scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
452     // scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
453     // painCave.isFatal = 1;
454     // simError();
455     // } else {
456     // info->getBoxM(hm);
457     // info->matMul3(hm, scaleMat, hmnew);
458     // info->setBoxM(hmnew);
459     // }
460    
461    
462    
463 gezelter 606 }
464    
465 tim 645 template<typename T> void NPTfm<T>::moveB( void ){
466 gezelter 606
467 tim 767 // int i, j;
468     // DirectionalAtom* dAtom;
469     // double Tb[3], ji[3];
470     // double vel[3], frc[3];
471     // double mass;
472 gezelter 606
473 tim 767 // double instaTemp, instaPress, instaVol;
474     // double tt2, tb2;
475     // double sc[3];
476     // double press[3][3], vScale[3][3];
477     // double oldChi, prevChi;
478     // double oldEta[3][3], preEta[3][3], diffEta;
479 gezelter 606
480 tim 767 // /*
481     // tt2 = tauThermostat * tauThermostat;
482     // tb2 = tauBarostat * tauBarostat;
483    
484     // instaTemp = tStats->getTemperature();
485     // tStats->getPressureTensor(press);
486     // instaVol = tStats->getVolume();
487 gezelter 606
488 tim 767 // // first evolve chi a half step
489 gezelter 606
490 tim 767 // chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
491 gezelter 606
492 tim 767 // for (i = 0; i < 3; i++ ) {
493     // for (j = 0; j < 3; j++ ) {
494     // if (i == j) {
495 gezelter 606
496 tim 767 // eta[i][j] += dt2 * instaVol *
497     // (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
498 gezelter 606
499 tim 767 // vScale[i][j] = eta[i][j] + chi;
500 gezelter 606
501 tim 767 // } else {
502 gezelter 606
503 tim 767 // eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
504 gezelter 606
505 tim 767 // vScale[i][j] = eta[i][j];
506 gezelter 606
507 tim 767 // }
508     // }
509     // }
510 gezelter 606
511 tim 767 // for( i=0; i<nAtoms; i++ ){
512 gezelter 606
513 tim 767 // atoms[i]->getVel( vel );
514     // atoms[i]->getFrc( frc );
515 gezelter 606
516 tim 767 // mass = atoms[i]->getMass();
517 gezelter 606
518 tim 767 // // velocity half step
519 gezelter 606
520 tim 767 // info->matVecMul3( vScale, vel, sc );
521 gezelter 606
522 tim 767 // for (j = 0; j < 3; j++) {
523     // vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
524     // }
525 gezelter 606
526 tim 767 // atoms[i]->setVel( vel );
527 gezelter 606
528 tim 767 // if( atoms[i]->isDirectional() ){
529 gezelter 606
530 tim 767 // dAtom = (DirectionalAtom *)atoms[i];
531 gezelter 606
532 tim 767 // // get and convert the torque to body frame
533 gezelter 606
534 tim 767 // dAtom->getTrq( Tb );
535     // dAtom->lab2Body( Tb );
536 gezelter 606
537 tim 767 // // get the angular momentum, and propagate a half step
538 gezelter 606
539 tim 767 // dAtom->getJ( ji );
540 gezelter 606
541 tim 767 // for (j=0; j < 3; j++)
542     // ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
543 gezelter 606
544 tim 767 // dAtom->setJ( ji );
545 gezelter 606
546 tim 767 // }
547     // }
548     // */
549    
550     // tt2 = tauThermostat * tauThermostat;
551     // tb2 = tauBarostat * tauBarostat;
552    
553    
554     // // Set things up for the iteration:
555    
556     // oldChi = chi;
557    
558     // for(i = 0; i < 3; i++)
559     // for(j = 0; j < 3; j++)
560     // oldEta[i][j] = eta[i][j];
561    
562     // for( i=0; i<nAtoms; i++ ){
563    
564     // atoms[i]->getVel( vel );
565    
566     // for (j=0; j < 3; j++)
567     // oldVel[3*i + j] = vel[j];
568    
569     // if( atoms[i]->isDirectional() ){
570    
571     // dAtom = (DirectionalAtom *)atoms[i];
572    
573     // dAtom->getJ( ji );
574    
575     // for (j=0; j < 3; j++)
576     // oldJi[3*i + j] = ji[j];
577    
578     // }
579     // }
580    
581     // // do the iteration:
582    
583     // instaVol = tStats->getVolume();
584    
585     // for (k=0; k < 4; k++) {
586    
587     // instaTemp = tStats->getTemperature();
588     // tStats->getPressureTensor(press);
589    
590     // // evolve chi another half step using the temperature at t + dt/2
591    
592     // prevChi = chi;
593     // chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
594    
595     // for(i = 0; i < 3; i++)
596     // for(j = 0; j < 3; j++)
597     // preEta[i][j] = eta[i][j];
598    
599     // //advance eta half step and calculate scale factor for velocity
600     // for(i = 0; i < 3; i ++)
601     // for(j = 0; j < 3; j++){
602     // if( i == j){
603     // eta[i][j] = oldEta[i][j] + dt2 * instaVol *
604     // (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
605     // vScale[i][j] = eta[i][j] + chi;
606     // }
607     // else
608     // {
609     // eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
610     // vScale[i][j] = eta[i][j];
611     // }
612     // }
613    
614     // //advance velocity half step
615     // for( i=0; i<nAtoms; i++ ){
616    
617     // atoms[i]->getFrc( frc );
618     // atoms[i]->getVel(vel);
619    
620     // mass = atoms[i]->getMass();
621    
622     // info->matVecMul3( vScale, vel, sc );
623    
624     // for (j=0; j < 3; j++) {
625     // // velocity half step (use chi from previous step here):
626     // vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass) * eConvert - sc[j]);
627     // }
628    
629     // atoms[i]->setVel( vel );
630    
631     // if( atoms[i]->isDirectional() ){
632    
633     // dAtom = (DirectionalAtom *)atoms[i];
634    
635     // // get and convert the torque to body frame
636    
637     // dAtom->getTrq( Tb );
638     // dAtom->lab2Body( Tb );
639    
640     // for (j=0; j < 3; j++)
641     // ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
642    
643     // dAtom->setJ( ji );
644     // }
645     // }
646    
647    
648     // diffEta = 0;
649     // for(i = 0; i < 3; i++)
650     // diffEta += pow(preEta[i][i] - eta[i][i], 2);
651    
652     // if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
653     // break;
654     // }
655    
656     // //calculate integral of chida
657     // integralOfChidt += dt2*chi;
658 gezelter 606 }
659    
660 mmeineke 746 template<typename T> void NPTfm<T>::resetIntegrator() {
661     int i,j;
662    
663     chi = 0.0;
664    
665     for(i = 0; i < 3; i++)
666     for (j = 0; j < 3; j++)
667     eta[i][j] = 0.0;
668     }
669    
670 tim 645 template<typename T> int NPTfm<T>::readyCheck() {
671 tim 658
672     //check parent's readyCheck() first
673     if (T::readyCheck() == -1)
674     return -1;
675    
676 gezelter 606 // First check to see if we have a target temperature.
677     // Not having one is fatal.
678    
679     if (!have_target_temp) {
680     sprintf( painCave.errMsg,
681     "NPTfm error: You can't use the NPTfm integrator\n"
682     " without a targetTemp!\n"
683     );
684     painCave.isFatal = 1;
685     simError();
686     return -1;
687     }
688    
689     if (!have_target_pressure) {
690     sprintf( painCave.errMsg,
691     "NPTfm error: You can't use the NPTfm integrator\n"
692     " without a targetPressure!\n"
693     );
694     painCave.isFatal = 1;
695     simError();
696     return -1;
697     }
698    
699     // We must set tauThermostat.
700    
701     if (!have_tau_thermostat) {
702     sprintf( painCave.errMsg,
703     "NPTfm error: If you use the NPTfm\n"
704     " integrator, you must set tauThermostat.\n");
705     painCave.isFatal = 1;
706     simError();
707     return -1;
708     }
709    
710     // We must set tauBarostat.
711    
712     if (!have_tau_barostat) {
713     sprintf( painCave.errMsg,
714     "NPTfm error: If you use the NPTfm\n"
715     " integrator, you must set tauBarostat.\n");
716     painCave.isFatal = 1;
717     simError();
718     return -1;
719     }
720    
721     // We need NkBT a lot, so just set it here:
722    
723     NkBT = (double)info->ndf * kB * targetTemp;
724    
725     return 1;
726     }
727 tim 763
728     template<typename T> double NPTfm<T>::getConservedQuantity(void){
729    
730    
731 tim 767 // double conservedQuantity;
732     // double tb2;
733     // double trEta;
734     // double E_NPT;
735     // double U;
736     // double TS;
737     // double PV;
738     // double extra;
739 tim 763
740 tim 767 // U = tStats->getTotalE();
741 tim 763
742 tim 767 // TS = fkBT *
743     // (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
744 tim 763
745 tim 767 // PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
746    
747     // tb2 = tauBarostat * tauBarostat;
748    
749     // trEta = info->matTrace3(eta);
750    
751     // extra = (fkBT * tb2 * trEta * trEta / 2.0 ) / eConvert;
752    
753     // cout.width(8);
754     // cout.precision(8);
755 tim 763
756 tim 767 // cout << info->getTime() << "\t"
757     // << chi << "\t"
758     // << trEta << "\t"
759     // << U << "\t"
760     // << TS << "\t"
761     // << PV << "\t"
762     // << extra << "\t"
763     // << U+TS+PV+extra << endl;
764    
765     // conservedQuantity = U+TS+PV+extra;
766     // return conservedQuantity;
767 tim 763 }