ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ExtendedSystem.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/ExtendedSystem.cpp (file contents):
Revision 456 by gezelter, Fri Apr 4 01:57:11 2003 UTC vs.
Revision 457 by gezelter, Fri Apr 4 19:16:11 2003 UTC

# Line 14 | Line 14 | ExtendedSystem::ExtendedSystem( SimInfo &info ) {
14    atoms = info.atoms;
15    nMols = info.n_mol;
16    molecules = info.molecules;
17 <  zeta = 0;
17 >  zeta = 0.0;
18 >  epsilonDot = 0.0;
19  
20   }
21  
# Line 22 | Line 23 | ExtendedSystem::~ExtendedSystem() {  
23   }
24  
25  
26 < void ExtendedSystem::NoseHooverNVT( double dt ){
26 > void ExtendedSystem::NoseHooverNVT( double dt, double ke ){
27  
28    // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
29    
30    int i;
31 <  double kB, keconverter, NkBT, zetaScale, ke_temp;
31 >  double NkBT, zetaScale, ke_temp;
32    double vx, vy, vz, jx, jy, jz;
33 <  
34 <  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
35 <  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K
36 <  
37 <  ke_temp = getKinetic() * keconverter;
33 >  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
34 >  const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to
35 >                                        // amu*Ang^2*fs^-2/K
36 >    
37 >  ke_temp = ke * e_convert;
38    NkBT = (double)getNDF() * kB * targetTemp;
39  
40 <  // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin &
40 >  // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
41    // qmass is set in the parameter file
42 <  zeta += dt*((ke_temp*2 - NkBT)/qmass);
42 >
43 >  zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
44    zetaScale = zeta * dt;
45  
46    // perform thermostat scaling on linear velocities and angular momentum
45 –  
47    for(i = 0; i < n_atoms; i++){
48      
49      vx = atoms[i]->get_vx();
50      vy = atoms[i]->get_vy();
51      vz = atoms[i]->get_vz();
52      
53 <    atoms[i]->set_vx(vx - zetaScale * vx);
54 <    atoms[i]->set_vy(vy - zetaScale * vy);
55 <    atoms[i]->set_vz(vz - zetaScale * vz);
53 >    atoms[i]->set_vx(vx * (1.0 - zetaScale));
54 >    atoms[i]->set_vy(vy * (1.0 - zetaScale));
55 >    atoms[i]->set_vz(vz * (1.0 - zetaScale));
56    }
57    if( n_oriented ){
58      
# Line 65 | Line 66 | void ExtendedSystem::NoseHooverNVT( double dt ){
66          jy = dAtom->getJy();
67          jz = dAtom->getJz();
68          
69 <        dAtom->setJx( jx - zetaScale * jx);
70 <        dAtom->setJy( jy - zetaScale * jy);
71 <        dAtom->setJz( jz - zetaScale * jz);
69 >        dAtom->setJx(jx * (1.0 - zetaScale));
70 >        dAtom->setJy(jy * (1.0 - zetaScale));
71 >        dAtom->setJz(jz * (1.0 - zetaScale));
72        }
73      }  
74    }
75   }
76  
77  
78 < void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke,
79 <                                           double dt, double temp ) {
78 > void ExtendedSystem::NoseHooverAndersonNPT( double dt,
79 >                                            double ke,
80 >                                            double p_int ) {
81  
82    // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
83    // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500
84  
85 <  int i, j, degrees_freedom;
86 <  double pressure, dt, temp, pressure_units, epsilonScale;
87 <  double ke, kB, vxi, vyi, vzi, pressure_ext;
88 <  double boxx_old, boxy_old, boxz_old;
89 <  double keconverter, NkBT, zetaScale, ke_temp;
90 <  double jxi, jyi, jzi, scale;
85 >  double oldBox[3];
86 >  double newBox[3];
87 >  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
88 >  const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
89 >  const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to
90 >                                        // amu*Ang^2*fs^-2/K
91  
92 <  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
91 <  pressure_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
92 <  degrees_freedom = 6*nmol; // number of degrees of freedom for the system
93 <  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2/K
92 >  double p_ext;
93  
94 <  pressure_ext = pressure * pressure_units;
95 <  volume = boxx*boxy*boxz;
97 <  ke_temp = ke * keconverter;
98 <  NkBT = degrees_freedom*kB*temp;
94 >  p_ext = targetPressure * p_units;
95 >  p_mol = p_int * p_units;
96  
97 +  getBox(oldBox);
98 +
99 +  volume = oldBox[0]*oldBox[1]*oldBox[2];
100 +
101 +  ke_temp = ke * e_convert;
102 +  NkBT = (double)getNDF() * kB * targetTemp;
103 +
104    // propogate the strain rate
105  
106 <  epsilon_dot +=  dt * ( (p_mol - pressure_ext)*volume
107 <                         / (tau_relax*tau_relax * kB * targetTemp) );
106 >  epsilonDot +=  dt * ((p_mol - p_ext) * volume /
107 >                       (tauRelax*tauRelax * kB * targetTemp) );
108  
109    // determine the change in cell volume
110 <  scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0));
110 >  scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0));
111  
112 <  volume = volume * pow(scale, 3.0);
112 >  newBox[0] = oldBox[0] * scale;
113 >  newBox[1] = oldBox[1] * scale;
114 >  newBox[2] = oldBox[2] * scale;
115 >  volume = newBox[0]*newBox[1]*newBox[2];
116  
117    // perform affine transform to update positions with volume fluctuations
118 <  affine_transform( scale );
118 >  this->AffineTransform( oldBox, newBox );
119  
113 –  // save old lengths and update box size
114 –  boxx_old = boxx;
115 –  boxy_old = boxy;
116 –  boxz_old = boxz;
117 –
118 –  boxx = boxx_old*scale;
119 –  boxy = boxy_old*scale;
120 –  boxz = boxz_old*scale;
121 –
120    epsilonScale = epsilonDot * dt;
121  
122    // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
123    // qmass is set in the parameter file
124 <  zeta += dt * ( (ke_temp*2 - NkBT) / qmass );
124 >
125 >  zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
126    zetaScale = zeta * dt;
127    
128    // apply barostating and thermostating to velocities and angular momenta
# Line 133 | Line 132 | void ExtendedSystem::NoseHooverAndersonNPT(double pres
132      vy = atoms[i]->get_vy();
133      vz = atoms[i]->get_vz();
134      
135 <    atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale));
136 <    atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale));
137 <    atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale));
135 >    atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale));
136 >    atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale));
137 >    atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale));
138    }
139    if( n_oriented ){
140      
# Line 157 | Line 156 | void ExtendedSystem::NoseHooverAndersonNPT(double pres
156    }
157   }
158  
159 < void ExtendedSystem::AffineTransform( double scale ){
159 > void ExtendedSystem::AffineTransform( double oldBox[3], double newBox[3] ){
160  
161    int i;
162 <  double boxx_old, boxy_old, boxz_old, percentScale;
163 <  double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi;
164 <  double[3] r;
162 >  double r[3];
163 >  double boxNum[3];
164 >  double percentScale[3];
165 >  double rxi, ryi, rzi;
166      
167    // first determine the scaling factor from the box size change
168 <  percentScale = (boxx - boxx_old)/boxx_old;
168 >  percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0];
169 >  percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1];
170 >  percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2];
171    
170 –
172    for (i=0; i < nMols; i++) {
173      
174      molecules[i]->getCOM(r);
175      
176 <    // find the minimum image coordinates of the molecular centers of mass:
176 >    // find the minimum image coordinates of the molecular centers of mass:    
177      
178 <    
179 <    boxx_num = boxx_old*copysign(1.0,r[0])*(double)(int)(fabs(r[0]/boxx_old)+0.5);
178 >    boxNum[0] = oldBox[0] * copysign(1.0,r[0]) *
179 >      (double)(int)(fabs(r[0]/oldBox[0]) + 0.5);
180  
181 <    boxx_num = boxx_old*dsign(1.0d0,rx(i))*int(abs(rx(i)/boxx_old)+0.5d0);
182 <    boxy_num = boxy_old*dsign(1.0d0,ry(i))*int(abs(ry(i)/boxy_old)+0.5d0);
182 <    boxz_num = boxz_old*dsign(1.0d0,rz(i))*int(abs(rz(i)/boxz_old)+0.5d0);
181 >    boxNum[1] = oldBox[1] * copysign(1.0,r[1]) *
182 >      (double)(int)(fabs(r[1]/oldBox[1]) + 0.5);
183  
184 <    rxi = rx(i) - boxx_num;
185 <    ryi = ry(i) - boxy_num;
186 <    rzi = rz(i) - boxz_num;
184 >    boxNum[2] = oldBox[2] * copysign(1.0,r[2]) *
185 >      (double)(int)(fabs(r[2]/oldBox[2]) + 0.5);
186  
187 +    rxi = r[0] - boxNum[0];
188 +    ryi = r[1] - boxNum[1];
189 +    rzi = r[2] - boxNum[2];
190 +
191      // update the minimum image coordinates using the scaling factor
192 <    rxi = rxi + rxi*percentScale;
193 <    ryi = ryi + ryi*percentScale;
194 <    rzi = rzi + rzi*percentScale;
192 >    rxi += rxi*percentScale[0];
193 >    ryi += ryi*percentScale[1];
194 >    rzi += rzi*percentScale[2];
195  
196 <    rx(i) = rxi + boxx_num;
197 <    ry(i) = ryi + boxy_num;
198 <    rz(i) = rzi + boxz_num;
196 >    r[0] = rxi + boxNum[0];
197 >    r[1] = ryi + boxNum[1];
198 >    r[2] = rzi + boxNum[2];
199 >
200 >    molecules[i]->moveCOM(r);
201    }
202   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines