1 |
+ |
#include <cmath> |
2 |
|
#include "Atom.hpp" |
3 |
|
#include "SRI.hpp" |
4 |
|
#include "AbstractClasses.hpp" |
43 |
|
double tt2, tb2; |
44 |
|
double angle; |
45 |
|
|
46 |
+ |
|
47 |
|
tt2 = tauThermostat * tauThermostat; |
48 |
|
tb2 = tauBarostat * tauBarostat; |
49 |
|
|
51 |
|
instaPress = tStats->getPressure(); |
52 |
|
instaVol = tStats->getVolume(); |
53 |
|
|
54 |
< |
// first evolve chi a half step |
54 |
> |
// first evolve chi a half step |
55 |
|
|
56 |
|
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
57 |
< |
eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2)); |
57 |
> |
eta += dt2 * ( instaVol * (instaPress - targetPressure) / |
58 |
> |
(p_convert*NkBT*tb2)); |
59 |
|
|
60 |
|
for( i=0; i<nAtoms; i++ ){ |
61 |
|
atomIndex = i * 3; |
68 |
|
|
69 |
|
// position whole step |
70 |
|
|
71 |
< |
for( j=atomIndex; j<(atomIndex+3); j=j+3 ) { |
72 |
< |
rj[0] = pos[j]; |
73 |
< |
rj[1] = pos[j+1]; |
74 |
< |
rj[2] = pos[j+2]; |
71 |
> |
rj[0] = pos[atomIndex]; |
72 |
> |
rj[1] = pos[atomIndex+1]; |
73 |
> |
rj[2] = pos[atomIndex+2]; |
74 |
> |
|
75 |
> |
info->wrapVector(rj); |
76 |
|
|
77 |
< |
info->wrapVector(rj); |
78 |
< |
|
79 |
< |
pos[j] += dt * (vel[j] + eta*rj[0]); |
76 |
< |
pos[j+1] += dt * (vel[j+1] + eta*rj[1]); |
77 |
< |
pos[j+2] += dt * (vel[j+2] + eta*rj[2]); |
78 |
< |
} |
79 |
< |
|
80 |
< |
// Scale the box after all the positions have been moved: |
81 |
< |
|
82 |
< |
info->scaleBox(exp(dt*eta)); |
77 |
> |
pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]); |
78 |
> |
pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]); |
79 |
> |
pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]); |
80 |
|
|
81 |
|
if( atoms[i]->isDirectional() ){ |
82 |
|
|
129 |
|
} |
130 |
|
|
131 |
|
} |
132 |
+ |
// Scale the box after all the positions have been moved: |
133 |
+ |
|
134 |
+ |
cerr << "eta = " << eta |
135 |
+ |
<< "; exp(dt*eta) = " << exp(eta*dt) << "\n"; |
136 |
+ |
|
137 |
+ |
info->scaleBox(exp(dt*eta)); |
138 |
+ |
|
139 |
|
} |
140 |
|
|
141 |
|
void NPTi::moveB( void ){ |
155 |
|
instaVol = tStats->getVolume(); |
156 |
|
|
157 |
|
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
158 |
< |
eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2)); |
158 |
> |
eta += dt2 * ( instaVol * (instaPress - targetPressure) / |
159 |
> |
(p_convert*NkBT*tb2)); |
160 |
|
|
161 |
|
for( i=0; i<nAtoms; i++ ){ |
162 |
|
atomIndex = i * 3; |