# | Line 10 | Line 10 | ExtendedSystem::ExtendedSystem( SimInfo* the_entry_plu | |
---|---|---|
10 | // get what information we need from the SimInfo object | |
11 | ||
12 | entry_plug = the_entry_plug; | |
13 | – | nAtoms = entry_plug->n_atoms; |
14 | – | atoms = entry_plug->atoms; |
15 | – | nMols = entry_plug->n_mol; |
16 | – | molecules = entry_plug->molecules; |
17 | – | nOriented = entry_plug->n_oriented; |
18 | – | ndf = entry_plug->ndf; |
13 | zeta = 0.0; | |
14 | epsilonDot = 0.0; | |
21 | – | |
15 | } | |
16 | ||
17 | void ExtendedSystem::NoseHooverNVT( double dt, double ke ){ | |
# | Line 32 | Line 25 | void ExtendedSystem::NoseHooverNVT( double dt, double | |
25 | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to | |
26 | // amu*Ang^2*fs^-2/K | |
27 | DirectionalAtom* dAtom; | |
28 | + | atoms = entry_plug->atoms; |
29 | ||
36 | – | |
30 | ke_temp = ke * e_convert; | |
31 | < | NkBT = (double)ndf * kB * targetTemp; |
31 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
32 | ||
33 | // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin | |
34 | // qmass is set in the parameter file | |
35 | ||
36 | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); | |
37 | + | |
38 | zetaScale = zeta * dt; | |
39 | + | |
40 | + | std::cerr << "zetaScale = " << zetaScale << "\n"; |
41 | ||
42 | // perform thermostat scaling on linear velocities and angular momentum | |
43 | < | for(i = 0; i < nAtoms; i++){ |
43 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
44 | ||
45 | vx = atoms[i]->get_vx(); | |
46 | vy = atoms[i]->get_vy(); | |
47 | vz = atoms[i]->get_vz(); | |
48 | < | |
48 | > | |
49 | atoms[i]->set_vx(vx * (1.0 - zetaScale)); | |
50 | atoms[i]->set_vy(vy * (1.0 - zetaScale)); | |
51 | atoms[i]->set_vz(vz * (1.0 - zetaScale)); | |
52 | } | |
53 | < | if( nOriented ){ |
53 | > | if( entry_plug->n_oriented ){ |
54 | ||
55 | < | for( i=0; i < nAtoms; i++ ){ |
55 | > | for( i=0; i < entry_plug->n_atoms; i++ ){ |
56 | ||
57 | if( atoms[i]->isDirectional() ){ | |
58 | ||
# | Line 89 | Line 85 | void ExtendedSystem::NoseHooverAndersonNPT( double dt, | |
85 | const double e_convert = 4.184e-4; // to convert ke from kcal/mol to | |
86 | // amu*Ang^2*fs^-2/K | |
87 | ||
88 | + | int i; |
89 | double p_ext, zetaScale, epsilonScale, scale, NkBT, ke_temp; | |
90 | double volume, p_mol; | |
91 | double vx, vy, vz, jx, jy, jz; | |
92 | DirectionalAtom* dAtom; | |
93 | < | int i; |
93 | > | atoms = entry_plug->atoms; |
94 | ||
95 | p_ext = targetPressure * p_units; | |
96 | < | p_mol = p_int * p_units; |
96 | > | p_mol = p_int; |
97 | ||
98 | entry_plug->getBox(oldBox); | |
102 | – | |
99 | volume = oldBox[0]*oldBox[1]*oldBox[2]; | |
100 | ||
101 | ke_temp = ke * e_convert; | |
102 | < | NkBT = (double)ndf * kB * targetTemp; |
102 | > | NkBT = (double)entry_plug->ndf * kB * targetTemp; |
103 | ||
104 | // propogate the strain rate | |
105 | ||
106 | epsilonDot += dt * ((p_mol - p_ext) * volume / | |
107 | (tauRelax*tauRelax * kB * targetTemp) ); | |
108 | ||
109 | + | |
110 | + | std::cerr << "dt = " << dt << " tauRelax = " << tauRelax << " kB = " << kB << "targetTemp = " << targetTemp << "\n"; |
111 | + | |
112 | // determine the change in cell volume | |
113 | scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0)); | |
114 | ||
115 | + | std::cerr << "p_mol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n"; |
116 | + | |
117 | + | |
118 | newBox[0] = oldBox[0] * scale; | |
119 | newBox[1] = oldBox[1] * scale; | |
120 | newBox[2] = oldBox[2] * scale; | |
# | Line 130 | Line 132 | void ExtendedSystem::NoseHooverAndersonNPT( double dt, | |
132 | ||
133 | zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass ); | |
134 | zetaScale = zeta * dt; | |
135 | + | |
136 | + | std::cerr << "zetaScale = " << zetaScale << "epsilonScale = " << epsilonScale << "\n"; |
137 | ||
138 | // apply barostating and thermostating to velocities and angular momenta | |
139 | < | for(i = 0; i < nAtoms; i++){ |
139 | > | for(i = 0; i < entry_plug->n_atoms; i++){ |
140 | ||
141 | vx = atoms[i]->get_vx(); | |
142 | vy = atoms[i]->get_vy(); | |
# | Line 142 | Line 146 | void ExtendedSystem::NoseHooverAndersonNPT( double dt, | |
146 | atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale)); | |
147 | atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale)); | |
148 | } | |
149 | < | if( nOriented ){ |
149 | > | if( entry_plug->n_oriented ){ |
150 | ||
151 | < | for( i=0; i < nAtoms; i++ ){ |
151 | > | for( i=0; i < entry_plug->n_atoms; i++ ){ |
152 | ||
153 | if( atoms[i]->isDirectional() ){ | |
154 | ||
# | Line 168 | Line 172 | void ExtendedSystem::AffineTransform( double oldBox[3] | |
172 | double r[3]; | |
173 | double boxNum[3]; | |
174 | double percentScale[3]; | |
175 | + | double delta[3]; |
176 | double rxi, ryi, rzi; | |
177 | + | |
178 | + | molecules = entry_plug->molecules; |
179 | ||
180 | // first determine the scaling factor from the box size change | |
181 | percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0]; | |
182 | percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1]; | |
183 | percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2]; | |
184 | ||
185 | < | for (i=0; i < nMols; i++) { |
185 | > | for (i=0; i < entry_plug->n_mol; i++) { |
186 | ||
187 | molecules[i].getCOM(r); | |
188 | < | |
188 | > | |
189 | // find the minimum image coordinates of the molecular centers of mass: | |
190 | ||
191 | boxNum[0] = oldBox[0] * copysign(1.0,r[0]) * | |
# | Line 199 | Line 206 | void ExtendedSystem::AffineTransform( double oldBox[3] | |
206 | ryi += ryi*percentScale[1]; | |
207 | rzi += rzi*percentScale[2]; | |
208 | ||
209 | < | r[0] = rxi + boxNum[0]; |
210 | < | r[1] = ryi + boxNum[1]; |
211 | < | r[2] = rzi + boxNum[2]; |
209 | > | delta[0] = r[0] - (rxi + boxNum[0]); |
210 | > | delta[1] = r[1] - (ryi + boxNum[1]); |
211 | > | delta[2] = r[2] - (rzi + boxNum[2]); |
212 | ||
213 | < | molecules[i].moveCOM(r); |
213 | > | molecules[i].moveCOM(delta); |
214 | } | |
215 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |