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 454 by gezelter, Fri Apr 4 01:57:11 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 ) {
9  
# Line 9 | Line 14 | ExtendedSystem::ExtendedSystem( SimInfo &info ) {
14    atoms = info.atoms;
15    nMols = info.n_mol;
16    molecules = info.molecules;
17 +  zeta = 0;
18  
19   }
20  
# Line 16 | Line 22 | ExtendedSystem::~ExtendedSystem() {  
22   }
23  
24  
25 < void ExtendedSystem::nose_hoover_nvt( double ke, double dt, double temp ){
25 > void ExtendedSystem::NoseHooverNVT( double dt ){
26  
27    // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
28    
29 <  int i, j, degrees_freedom;
30 <  double ke, dt, temp, kB;
31 <  double keconverter, NkBT, zetaScale, ke_temp;
26 <  double vxi, vyi, vzi, jxi, jyi, jzi;
29 >  int i;
30 >  double kB, keconverter, NkBT, zetaScale, ke_temp;
31 >  double vx, vy, vz, jx, jy, jz;
32    
28 –  degrees_freedom = 6*nmol; // number of degrees of freedom for the system
33    kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
34    keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K
35    
36 <  ke_temp = ke * keconverter;
37 <  NkBT = degrees_freedom*kB*temp;
36 >  ke_temp = getKinetic() * keconverter;
37 >  NkBT = (double)getNDF() * 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 >  zeta += dt*((ke_temp*2 - NkBT)/qmass);
42    zetaScale = zeta * dt;
43  
44    // perform thermostat scaling on linear velocities and angular momentum
45 <  for(i = 0, i < nmol; i++ ) {
46 <    vxi = vx(i)*zetaScale;
47 <    vyi = vy(i)*zetaScale;
48 <    vzi = vz(i)*zetaScale;
49 <    jxi = jx(i)*zetaScale;
50 <    jyi = jy(i)*zetaScale;
51 <    jzi = jz(i)*zetaScale;
52 <
53 <    vx(i) = vx(i) - vxi;
54 <    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;
45 >  
46 >  for(i = 0; i < n_atoms; i++){
47 >    
48 >    vx = atoms[i]->get_vx();
49 >    vy = atoms[i]->get_vy();
50 >    vz = atoms[i]->get_vz();
51 >    
52 >    atoms[i]->set_vx(vx - zetaScale * vx);
53 >    atoms[i]->set_vy(vy - zetaScale * vy);
54 >    atoms[i]->set_vz(vz - zetaScale * vz);
55    }
56 +  if( n_oriented ){
57 +    
58 +    for( i=0; i < n_atoms; i++ ){
59 +      
60 +      if( atoms[i]->isDirectional() ){
61 +        
62 +        dAtom = (DirectionalAtom *)atoms[i];
63 +        
64 +        jx = dAtom->getJx();
65 +        jy = dAtom->getJy();
66 +        jz = dAtom->getJz();
67 +        
68 +        dAtom->setJx( jx - zetaScale * jx);
69 +        dAtom->setJy( jy - zetaScale * jy);
70 +        dAtom->setJz( jz - zetaScale * jz);
71 +      }
72 +    }  
73 +  }
74   }
75  
76  
77 < void ExtendedSystem::nose_hoover_anderson_npt(double pressure, double ke, double dt,
78 <                                   double temp ) {
77 > void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke,
78 >                                           double dt, double temp ) {
79  
80    // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
81    // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500
# Line 82 | Line 100 | void ExtendedSystem::nose_hoover_anderson_npt(double p
100    // propogate the strain rate
101  
102    epsilon_dot +=  dt * ( (p_mol - pressure_ext)*volume
103 <                         / (tau_relax*tau_relax * kB * temp) );
103 >                         / (tau_relax*tau_relax * kB * targetTemp) );
104  
105    // determine the change in cell volume
106    scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0));
# Line 101 | Line 119 | void ExtendedSystem::nose_hoover_anderson_npt(double p
119    boxy = boxy_old*scale;
120    boxz = boxz_old*scale;
121  
122 <  epsilonScale = epsilon_dot * dt;
122 >  epsilonScale = epsilonDot * dt;
123  
124    // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
125    // qmass is set in the parameter file
# Line 109 | Line 127 | void ExtendedSystem::nose_hoover_anderson_npt(double p
127    zetaScale = zeta * dt;
128    
129    // apply barostating and thermostating to velocities and angular momenta
130 <  
131 <  for (i=0; i < nmol; i++) {
132 <
133 <    vxi = vx(i)*epsilonScale;
134 <    vyi = vy(i)*epsilonScale;
135 <    vzi = vz(i)*epsilonScale;
136 <    vxi = vxi + vx(i)*zetaScale;
137 <    vyi = vyi + vy(i)*zetaScale;
138 <    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;
130 >  for(i = 0; i < n_atoms; i++){
131 >    
132 >    vx = atoms[i]->get_vx();
133 >    vy = atoms[i]->get_vy();
134 >    vz = atoms[i]->get_vz();
135 >    
136 >    atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale));
137 >    atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale));
138 >    atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale));
139    }
140 <  
141 <
142 <
140 >  if( n_oriented ){
141 >    
142 >    for( i=0; i < n_atoms; i++ ){
143 >      
144 >      if( atoms[i]->isDirectional() ){
145 >        
146 >        dAtom = (DirectionalAtom *)atoms[i];
147 >        
148 >        jx = dAtom->getJx();
149 >        jy = dAtom->getJy();
150 >        jz = dAtom->getJz();
151 >        
152 >        dAtom->setJx( jx * (1.0 - zetaScale));
153 >        dAtom->setJy( jy * (1.0 - zetaScale));
154 >        dAtom->setJz( jz * (1.0 - zetaScale));
155 >      }
156 >    }  
157 >  }
158   }
159  
160 < void ExtendedSystem::affine_transform( double scale ){
160 > void ExtendedSystem::AffineTransform( double scale ){
161  
162    int i;
163    double boxx_old, boxy_old, boxz_old, percentScale;

Diff Legend

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