# | Line 1 | Line 1 | |
---|---|---|
1 | + | #include <cmath> |
2 | #include "Atom.hpp" | |
3 | #include "SRI.hpp" | |
4 | #include "AbstractClasses.hpp" | |
# | Line 32 | Line 33 | void NPTi::moveA() { | |
33 | ||
34 | void NPTi::moveA() { | |
35 | ||
36 | < | int i,j,k; |
36 | < | int atomIndex, aMatIndex; |
36 | > | int i, j; |
37 | DirectionalAtom* dAtom; | |
38 | < | double Tb[3]; |
39 | < | double ji[3]; |
38 | > | double Tb[3], ji[3]; |
39 | > | double A[3][3], I[3][3]; |
40 | > | double angle, mass; |
41 | > | double vel[3], pos[3], frc[3]; |
42 | > | |
43 | double rj[3]; | |
44 | double instaTemp, instaPress, instaVol; | |
45 | double tt2, tb2; | |
43 | – | double angle; |
46 | ||
47 | tt2 = tauThermostat * tauThermostat; | |
48 | tb2 = tauBarostat * tauBarostat; | |
# | Line 49 | Line 51 | void NPTi::moveA() { | |
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; |
62 | < | aMatIndex = i * 9; |
63 | < | |
61 | < | // velocity half step |
62 | < | for( j=atomIndex; j<(atomIndex+3); j++ ) |
63 | < | vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert |
64 | < | - vel[j]*(chi+eta)); |
61 | > | atoms[i]->getVel( vel ); |
62 | > | atoms[i]->getPos( pos ); |
63 | > | atoms[i]->getFrc( frc ); |
64 | ||
65 | < | // position whole step |
65 | > | mass = atoms[i]->getMass(); |
66 | ||
67 | < | rj[0] = pos[atomIndex]; |
68 | < | rj[1] = pos[atomIndex+1]; |
69 | < | rj[2] = pos[atomIndex+2]; |
70 | < | |
67 | > | for (j=0; j < 3; j++) { |
68 | > | vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta)); |
69 | > | rj[j] = pos[j]; |
70 | > | } |
71 | > | |
72 | > | atoms[i]->setVel( vel ); |
73 | > | |
74 | info->wrapVector(rj); | |
75 | ||
76 | < | pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]); |
77 | < | pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]); |
78 | < | pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]); |
79 | < | |
76 | > | for (j = 0; j < 3; j++) |
77 | > | pos[j] += dt * (vel[j] + eta*rj[j]); |
78 | > | |
79 | > | |
80 | > | atoms[i]->setPos( pos ); |
81 | > | |
82 | > | |
83 | if( atoms[i]->isDirectional() ){ | |
84 | ||
85 | dAtom = (DirectionalAtom *)atoms[i]; | |
86 | ||
87 | // get and convert the torque to body frame | |
88 | ||
89 | < | Tb[0] = dAtom->getTx(); |
85 | < | Tb[1] = dAtom->getTy(); |
86 | < | Tb[2] = dAtom->getTz(); |
87 | < | |
89 | > | dAtom->getTrq( Tb ); |
90 | dAtom->lab2Body( Tb ); | |
91 | ||
92 | // get the angular momentum, and propagate a half step | |
93 | ||
94 | < | ji[0] = dAtom->getJx(); |
95 | < | ji[1] = dAtom->getJy(); |
96 | < | ji[2] = dAtom->getJz(); |
94 | > | dAtom->getJ( ji ); |
95 | > | |
96 | > | for (j=0; j < 3; j++) |
97 | > | ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); |
98 | ||
96 | – | ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi); |
97 | – | ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi); |
98 | – | ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi); |
99 | – | |
99 | // use the angular velocities to propagate the rotation matrix a | |
100 | // full time step | |
101 | < | |
101 | > | |
102 | > | dAtom->getA(A); |
103 | > | dAtom->getI(I); |
104 | > | |
105 | // rotate about the x-axis | |
106 | < | angle = dt2 * ji[0] / dAtom->getIxx(); |
107 | < | this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); |
108 | < | |
106 | > | angle = dt2 * ji[0] / I[0][0]; |
107 | > | this->rotate( 1, 2, angle, ji, A ); |
108 | > | |
109 | // rotate about the y-axis | |
110 | < | angle = dt2 * ji[1] / dAtom->getIyy(); |
111 | < | this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); |
110 | > | angle = dt2 * ji[1] / I[1][1]; |
111 | > | this->rotate( 2, 0, angle, ji, A ); |
112 | ||
113 | // rotate about the z-axis | |
114 | < | angle = dt * ji[2] / dAtom->getIzz(); |
115 | < | this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] ); |
114 | > | angle = dt * ji[2] / I[2][2]; |
115 | > | this->rotate( 0, 1, angle, ji, A); |
116 | ||
117 | // rotate about the y-axis | |
118 | < | angle = dt2 * ji[1] / dAtom->getIyy(); |
119 | < | this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); |
118 | > | angle = dt2 * ji[1] / I[1][1]; |
119 | > | this->rotate( 2, 0, angle, ji, A ); |
120 | ||
121 | // rotate about the x-axis | |
122 | < | angle = dt2 * ji[0] / dAtom->getIxx(); |
123 | < | this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); |
122 | > | angle = dt2 * ji[0] / I[0][0]; |
123 | > | this->rotate( 1, 2, angle, ji, A ); |
124 | ||
125 | < | dAtom->setJx( ji[0] ); |
126 | < | dAtom->setJy( ji[1] ); |
127 | < | dAtom->setJz( ji[2] ); |
128 | < | } |
127 | < | |
125 | > | dAtom->setJ( ji ); |
126 | > | dAtom->setA( A ); |
127 | > | } |
128 | > | |
129 | } | |
130 | // Scale the box after all the positions have been moved: | |
131 | < | |
132 | < | info->scaleBox(exp(dt*eta)); |
133 | < | |
131 | > | |
132 | > | cerr << "eta = " << eta |
133 | > | << "; exp(dt*eta) = " << exp(eta*dt) << "\n"; |
134 | > | |
135 | > | info->scaleBox(exp(dt*eta)); |
136 | } | |
137 | ||
138 | void NPTi::moveB( void ){ | |
139 | < | int i,j,k; |
140 | < | int atomIndex; |
139 | > | |
140 | > | int i, j; |
141 | DirectionalAtom* dAtom; | |
142 | < | double Tb[3]; |
143 | < | double ji[3]; |
142 | > | double Tb[3], ji[3]; |
143 | > | double vel[3], frc[3]; |
144 | > | double mass; |
145 | > | |
146 | double instaTemp, instaPress, instaVol; | |
147 | double tt2, tb2; | |
148 | ||
# | Line 149 | Line 154 | void NPTi::moveB( void ){ | |
154 | instaVol = tStats->getVolume(); | |
155 | ||
156 | chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; | |
157 | < | eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2)); |
157 | > | eta += dt2 * ( instaVol * (instaPress - targetPressure) / |
158 | > | (p_convert*NkBT*tb2)); |
159 | ||
160 | for( i=0; i<nAtoms; i++ ){ | |
161 | < | atomIndex = i * 3; |
162 | < | |
161 | > | |
162 | > | atoms[i]->getVel( vel ); |
163 | > | atoms[i]->getFrc( frc ); |
164 | > | |
165 | > | mass = atoms[i]->getMass(); |
166 | > | |
167 | // velocity half step | |
168 | < | for( j=atomIndex; j<(atomIndex+3); j++ ) |
169 | < | for( j=atomIndex; j<(atomIndex+3); j++ ) |
160 | < | vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert |
161 | < | - vel[j]*(chi+eta)); |
168 | > | for (j=0; j < 3; j++) |
169 | > | vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta)); |
170 | ||
171 | + | atoms[i]->setVel( vel ); |
172 | + | |
173 | if( atoms[i]->isDirectional() ){ | |
174 | < | |
174 | > | |
175 | dAtom = (DirectionalAtom *)atoms[i]; | |
176 | < | |
177 | < | // get and convert the torque to body frame |
178 | < | |
179 | < | Tb[0] = dAtom->getTx(); |
170 | < | Tb[1] = dAtom->getTy(); |
171 | < | Tb[2] = dAtom->getTz(); |
172 | < | |
176 | > | |
177 | > | // get and convert the torque to body frame |
178 | > | |
179 | > | dAtom->getTrq( Tb ); |
180 | dAtom->lab2Body( Tb ); | |
181 | < | |
182 | < | // get the angular momentum, and complete the angular momentum |
183 | < | // half step |
184 | < | |
185 | < | ji[0] = dAtom->getJx(); |
186 | < | ji[1] = dAtom->getJy(); |
187 | < | ji[2] = dAtom->getJz(); |
188 | < | |
189 | < | ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi); |
183 | < | ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi); |
184 | < | ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi); |
185 | < | |
186 | < | dAtom->setJx( ji[0] ); |
187 | < | dAtom->setJy( ji[1] ); |
188 | < | dAtom->setJz( ji[2] ); |
181 | > | |
182 | > | // get the angular momentum, and propagate a half step |
183 | > | |
184 | > | dAtom->getJ( ji ); |
185 | > | |
186 | > | for (j=0; j < 3; j++) |
187 | > | ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); |
188 | > | |
189 | > | dAtom->setJ( ji ); |
190 | } | |
191 | } | |
192 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |