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 453 by gezelter, Thu Apr 3 23:49:20 2003 UTC vs.
Revision 471 by gezelter, Mon Apr 7 20:51:59 2003 UTC

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

Diff Legend

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