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 454 by gezelter, Fri Apr 4 01:57:11 2003 UTC vs.
Revision 475 by gezelter, Tue Apr 8 12:44:18 2003 UTC

# Line 5 | Line 5 | ExtendedSystem::ExtendedSystem( SimInfo &info ) {
5   #include "Thermo.hpp"
6   #include "ExtendedSystem.hpp"
7  
8 < ExtendedSystem::ExtendedSystem( SimInfo &info ) {
8 > ExtendedSystem::ExtendedSystem( SimInfo* the_entry_plug ) {
9  
10    // get what information we need from the SimInfo object
11    
12 <  entry_plug = &info;
13 <  nAtoms = info.n_atoms;
14 <  atoms = info.atoms;
15 <  nMols = info.n_mol;
16 <  molecules = info.molecules;
17 <  zeta = 0;
18 <
12 >  entry_plug = the_entry_plug;
13 >  zeta = 0.0;
14 >  epsilonDot = 0.0;
15   }
16  
17 < ExtendedSystem::~ExtendedSystem() {  
22 < }
17 > void ExtendedSystem::NoseHooverNVT( double dt, double ke ){
18  
24 –
25 – void ExtendedSystem::NoseHooverNVT( double dt ){
26 –
19    // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
20    
21    int i;
22 <  double kB, keconverter, NkBT, zetaScale, ke_temp;
22 >  double NkBT, zetaScale, ke_temp;
23    double vx, vy, vz, jx, jy, jz;
24 <  
25 <  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
26 <  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K
27 <  
28 <  ke_temp = getKinetic() * keconverter;
37 <  NkBT = (double)getNDF() * kB * targetTemp;
24 >  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
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  
30 <  // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin &
30 >  ke_temp = ke * e_convert;
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 <  zeta += dt*((ke_temp*2 - NkBT)/qmass);
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 <  
46 <  for(i = 0; i < n_atoms; 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 <    
49 <    atoms[i]->set_vx(vx - zetaScale * vx);
50 <    atoms[i]->set_vy(vy - zetaScale * vy);
51 <    atoms[i]->set_vz(vz - zetaScale * vz);
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( n_oriented ){
53 >  if( entry_plug->n_oriented ){
54      
55 <    for( i=0; i < n_atoms; i++ ){
55 >    for( i=0; i < entry_plug->n_atoms; i++ ){
56        
57        if( atoms[i]->isDirectional() ){
58          
# Line 65 | Line 62 | void ExtendedSystem::NoseHooverNVT( double dt ){
62          jy = dAtom->getJy();
63          jz = dAtom->getJz();
64          
65 <        dAtom->setJx( jx - zetaScale * jx);
66 <        dAtom->setJy( jy - zetaScale * jy);
67 <        dAtom->setJz( jz - zetaScale * jz);
65 >        dAtom->setJx(jx * (1.0 - zetaScale));
66 >        dAtom->setJy(jy * (1.0 - zetaScale));
67 >        dAtom->setJz(jz * (1.0 - zetaScale));
68        }
69      }  
70    }
71   }
72  
73  
74 < void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke,
75 <                                           double dt, double temp ) {
74 > void ExtendedSystem::NoseHooverAndersonNPT( double dt,
75 >                                            double ke,
76 >                                            double p_int ) {
77  
78    // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
79    // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500
80  
81 <  int i, j, degrees_freedom;
82 <  double pressure, dt, temp, pressure_units, epsilonScale;
83 <  double ke, kB, vxi, vyi, vzi, pressure_ext;
84 <  double boxx_old, boxy_old, boxz_old;
85 <  double keconverter, NkBT, zetaScale, ke_temp;
86 <  double jxi, jyi, jzi, scale;
81 >  double oldBox[3];
82 >  double newBox[3];
83 >  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
84 >  const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
85 >  const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to
86 >                                        // amu*Ang^2*fs^-2/K
87  
88 <  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
89 <  pressure_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
90 <  degrees_freedom = 6*nmol; // number of degrees of freedom for the system
91 <  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2/K
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 >  atoms = entry_plug->atoms;
94  
95 <  pressure_ext = pressure * pressure_units;
96 <  volume = boxx*boxy*boxz;
97 <  ke_temp = ke * keconverter;
98 <  NkBT = degrees_freedom*kB*temp;
95 >  p_ext = targetPressure * p_units;
96 >  p_mol = p_int;
97  
98 +  entry_plug->getBox(oldBox);
99 +  volume = oldBox[0]*oldBox[1]*oldBox[2];
100 +
101 +  ke_temp = ke * e_convert;
102 +  NkBT = (double)entry_plug->ndf * 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 +
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 * epsilon_dot), (1.0 / 3.0));
113 >  scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0));
114  
115 <  volume = volume * pow(scale, 3.0);
115 >  std::cerr << "p_mol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n";
116  
110 –  // perform affine transform to update positions with volume fluctuations
111 –  affine_transform( scale );
117  
118 <  // save old lengths and update box size
119 <  boxx_old = boxx;
120 <  boxy_old = boxy;
121 <  boxz_old = boxz;
118 >  newBox[0] = oldBox[0] * scale;
119 >  newBox[1] = oldBox[1] * scale;
120 >  newBox[2] = oldBox[2] * scale;
121 >  volume = newBox[0]*newBox[1]*newBox[2];
122  
123 <  boxx = boxx_old*scale;
119 <  boxy = boxy_old*scale;
120 <  boxz = boxz_old*scale;
123 >  entry_plug->setBox(newBox);
124  
125 +  // perform affine transform to update positions with volume fluctuations
126 +  this->AffineTransform( oldBox, newBox );
127 +
128    epsilonScale = epsilonDot * dt;
129  
130    // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
131    // qmass is set in the parameter file
132 <  zeta += dt * ( (ke_temp*2 - NkBT) / qmass );
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 < n_atoms; i++){
139 >  for(i = 0; i < entry_plug->n_atoms; i++){
140      
141      vx = atoms[i]->get_vx();
142      vy = atoms[i]->get_vy();
143      vz = atoms[i]->get_vz();
144      
145 <    atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale));
146 <    atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale));
147 <    atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale));
145 >    atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale));
146 >    atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale));
147 >    atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale));
148    }
149 <  if( n_oriented ){
149 >  if( entry_plug->n_oriented ){
150      
151 <    for( i=0; i < n_atoms; i++ ){
151 >    for( i=0; i < entry_plug->n_atoms; i++ ){
152        
153        if( atoms[i]->isDirectional() ){
154          
# Line 157 | Line 166 | void ExtendedSystem::AffineTransform( double scale ){
166    }
167   }
168  
169 < void ExtendedSystem::AffineTransform( double scale ){
169 > void ExtendedSystem::AffineTransform( double oldBox[3], double newBox[3] ){
170  
171    int i;
172 <  double boxx_old, boxy_old, boxz_old, percentScale;
173 <  double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi;
174 <  double[3] r;
172 >  double r[3];
173 >  double boxNum[3];
174 >  double percentScale[3];
175 >  double rxi, ryi, rzi;
176 >
177 >  molecules = entry_plug->molecules;
178      
179    // first determine the scaling factor from the box size change
180 <  percentScale = (boxx - boxx_old)/boxx_old;
180 >  percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0];
181 >  percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1];
182 >  percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2];
183    
184 +  std::cerr << "ps = " << percentScale[0] <<" " << percentScale[1] << " " << percentScale[2] << "\n";
185  
186 <  for (i=0; i < nMols; i++) {
186 >
187 >  for (i=0; i < entry_plug->n_mol; i++) {
188      
189 <    molecules[i]->getCOM(r);
189 >    molecules[i].getCOM(r);
190 >
191 >    // find the minimum image coordinates of the molecular centers of mass:    
192      
193 <    // find the minimum image coordinates of the molecular centers of mass:
194 <    
177 <    
178 <    boxx_num = boxx_old*copysign(1.0,r[0])*(double)(int)(fabs(r[0]/boxx_old)+0.5);
193 >    boxNum[0] = oldBox[0] * copysign(1.0,r[0]) *
194 >      (double)(int)(fabs(r[0]/oldBox[0]) + 0.5);
195  
196 <    boxx_num = boxx_old*dsign(1.0d0,rx(i))*int(abs(rx(i)/boxx_old)+0.5d0);
197 <    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);
196 >    boxNum[1] = oldBox[1] * copysign(1.0,r[1]) *
197 >      (double)(int)(fabs(r[1]/oldBox[1]) + 0.5);
198  
199 <    rxi = rx(i) - boxx_num;
200 <    ryi = ry(i) - boxy_num;
186 <    rzi = rz(i) - boxz_num;
199 >    boxNum[2] = oldBox[2] * copysign(1.0,r[2]) *
200 >      (double)(int)(fabs(r[2]/oldBox[2]) + 0.5);
201  
202 +
203 +    std::cerr << "boxNum = " << boxNum[0] << " " << boxNum[1] << " " << boxNum[2] << "\n";
204 +
205 +    rxi = r[0] - boxNum[0];
206 +    ryi = r[1] - boxNum[1];
207 +    rzi = r[2] - boxNum[2];
208 +
209      // update the minimum image coordinates using the scaling factor
210 <    rxi = rxi + rxi*percentScale;
211 <    ryi = ryi + ryi*percentScale;
212 <    rzi = rzi + rzi*percentScale;
210 >    rxi += rxi*percentScale[0];
211 >    ryi += ryi*percentScale[1];
212 >    rzi += rzi*percentScale[2];
213  
214 <    rx(i) = rxi + boxx_num;
215 <    ry(i) = ryi + boxy_num;
216 <    rz(i) = rzi + boxz_num;
214 >    r[0] = rxi + boxNum[0];
215 >    r[1] = ryi + boxNum[1];
216 >    r[2] = rzi + boxNum[2];
217 >
218 >    molecules[i].moveCOM(r);
219    }
220   }

Diff Legend

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