OOPSE/libmdtools/ExtendedSystem.cpp

#include <math.h>
#include "Atom.hpp"
#include "Molecule.hpp"
#include "SimInfo.hpp"
#include "Thermo.hpp"
#include "ExtendedSystem.hpp"

ExtendedSystem::ExtendedSystem( SimInfo &info ) {

  // get what information we need from the SimInfo object
  
  entry_plug = &info;
  nAtoms = info.n_atoms;
  atoms = info.atoms;
  nMols = info.n_mol;
  molecules = info.molecules;
  zeta = 0;

}

ExtendedSystem::~ExtendedSystem() {  
}


void ExtendedSystem::NoseHooverNVT( double dt ){

  // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
  
  int i;
  double kB, keconverter, NkBT, zetaScale, ke_temp;
  double vx, vy, vz, jx, jy, jz;
  
  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K
  
  ke_temp = getKinetic() * keconverter;
  NkBT = (double)getNDF() * kB * targetTemp;

  // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin &
  // qmass is set in the parameter file
  zeta += dt*((ke_temp*2 - NkBT)/qmass);
  zetaScale = zeta * dt;

  // perform thermostat scaling on linear velocities and angular momentum
  
  for(i = 0; i < n_atoms; i++){
    
    vx = atoms[i]->get_vx();
    vy = atoms[i]->get_vy();
    vz = atoms[i]->get_vz();
    
    atoms[i]->set_vx(vx - zetaScale * vx);
    atoms[i]->set_vy(vy - zetaScale * vy);
    atoms[i]->set_vz(vz - zetaScale * vz);
  }
  if( n_oriented ){
    
    for( i=0; i < n_atoms; i++ ){
      
      if( atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *)atoms[i];
        
        jx = dAtom->getJx();
        jy = dAtom->getJy();
        jz = dAtom->getJz();
        
        dAtom->setJx( jx - zetaScale * jx);
        dAtom->setJy( jy - zetaScale * jy);
        dAtom->setJz( jz - zetaScale * jz);
      }
    }   
  }
}


void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke, 
                                           double dt, double temp ) {

  // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 
  // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500

  int i, j, degrees_freedom;
  double pressure, dt, temp, pressure_units, epsilonScale;
  double ke, kB, vxi, vyi, vzi, pressure_ext;
  double boxx_old, boxy_old, boxz_old;
  double keconverter, NkBT, zetaScale, ke_temp;
  double jxi, jyi, jzi, scale;

  kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
  pressure_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
  degrees_freedom = 6*nmol; // number of degrees of freedom for the system
  keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2/K

  pressure_ext = pressure * pressure_units;
  volume = boxx*boxy*boxz;
  ke_temp = ke * keconverter;
  NkBT = degrees_freedom*kB*temp;

  // propogate the strain rate

  epsilon_dot +=  dt * ( (p_mol - pressure_ext)*volume 
                         / (tau_relax*tau_relax * kB * targetTemp) );

  // determine the change in cell volume
  scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0));

  volume = volume * pow(scale, 3.0);

  // perform affine transform to update positions with volume fluctuations
  affine_transform( scale );

  // save old lengths and update box size
  boxx_old = boxx;
  boxy_old = boxy;
  boxz_old = boxz;

  boxx = boxx_old*scale;
  boxy = boxy_old*scale;
  boxz = boxz_old*scale;

  epsilonScale = epsilonDot * dt;

  // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin 
  // qmass is set in the parameter file
  zeta += dt * ( (ke_temp*2 - NkBT) / qmass );
  zetaScale = zeta * dt;
  
  // apply barostating and thermostating to velocities and angular momenta
  for(i = 0; i < n_atoms; i++){
    
    vx = atoms[i]->get_vx();
    vy = atoms[i]->get_vy();
    vz = atoms[i]->get_vz();
    
    atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale));
    atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale));
    atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale));
  }
  if( n_oriented ){
    
    for( i=0; i < n_atoms; i++ ){
      
      if( atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *)atoms[i];
        
        jx = dAtom->getJx();
        jy = dAtom->getJy();
        jz = dAtom->getJz();
        
        dAtom->setJx( jx * (1.0 - zetaScale));
        dAtom->setJy( jy * (1.0 - zetaScale));
        dAtom->setJz( jz * (1.0 - zetaScale));
      }
    }   
  }
}

void ExtendedSystem::AffineTransform( double scale ){

  int i;
  double boxx_old, boxy_old, boxz_old, percentScale;
  double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi;
  double[3] r;
    
  // first determine the scaling factor from the box size change
  percentScale = (boxx - boxx_old)/boxx_old;
  

  for (i=0; i < nMols; i++) {
    
    molecules[i]->getCOM(r);
    
    // find the minimum image coordinates of the molecular centers of mass:
    
    
    boxx_num = boxx_old*copysign(1.0,r[0])*(double)(int)(fabs(r[0]/boxx_old)+0.5);

    boxx_num = boxx_old*dsign(1.0d0,rx(i))*int(abs(rx(i)/boxx_old)+0.5d0);
    boxy_num = boxy_old*dsign(1.0d0,ry(i))*int(abs(ry(i)/boxy_old)+0.5d0);
    boxz_num = boxz_old*dsign(1.0d0,rz(i))*int(abs(rz(i)/boxz_old)+0.5d0);

    rxi = rx(i) - boxx_num;
    ryi = ry(i) - boxy_num;
    rzi = rz(i) - boxz_num;

    // update the minimum image coordinates using the scaling factor
    rxi = rxi + rxi*percentScale;
    ryi = ryi + ryi*percentScale;
    rzi = rzi + rzi*percentScale;

    rx(i) = rxi + boxx_num;
    ry(i) = ryi + boxy_num;
    rz(i) = rzi + boxz_num;
  }
}
Revision:	454
Committed:	Fri Apr 4 01:57:11 2003 UTC (21 years, 3 months ago) by gezelter
File size:	5438 byte(s)
Log Message:	changes for extended system code
#	User	Rev	Content
1	gezelter	453	#include <math.h>
2	gezelter	454	#include "Atom.hpp"
3			#include "Molecule.hpp"
4			#include "SimInfo.hpp"
5			#include "Thermo.hpp"
6			#include "ExtendedSystem.hpp"
7	gezelter	453
8			ExtendedSystem::ExtendedSystem( SimInfo &info ) {
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	gezelter	454	zeta = 0;
18	gezelter	453
19			}
20
21			ExtendedSystem::~ExtendedSystem() {
22			}
23
24
25	gezelter	454	void ExtendedSystem::NoseHooverNVT( double dt ){
26	gezelter	453
27			// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
28
29	gezelter	454	int i;
30			double kB, keconverter, NkBT, zetaScale, ke_temp;
31			double vx, vy, vz, jx, jy, jz;
32	gezelter	453
33			kB = 8.31451e-7; // boltzmann constant in amuAng^2fs^-2/K
34			keconverter = 4.184e-4; // to convert ke from kcal/mol to amuAng^2fs^-2 / K
35
36	gezelter	454	ke_temp = getKinetic() * keconverter;
37			NkBT = (double)getNDF() * kB * targetTemp;
38	gezelter	453
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	gezelter	454	zeta += dt((ke_temp2 - NkBT)/qmass);
42	gezelter	453	zetaScale = zeta * dt;
43
44			// perform thermostat scaling on linear velocities and angular momentum
45	gezelter	454
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	gezelter	453	}
56	gezelter	454	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	gezelter	453	}
75
76
77	gezelter	454	void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke,
78			double dt, double temp ) {
79	gezelter	453
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
82
83			int i, j, degrees_freedom;
84			double pressure, dt, temp, pressure_units, epsilonScale;
85			double ke, kB, vxi, vyi, vzi, pressure_ext;
86			double boxx_old, boxy_old, boxz_old;
87			double keconverter, NkBT, zetaScale, ke_temp;
88			double jxi, jyi, jzi, scale;
89
90			kB = 8.31451e-7; // boltzmann constant in amuAng^2fs^-2/K
91			pressure_units = 6.10192996e-9; // converts atm to amufs^-2Ang^-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 amuAng^2fs^-2/K
94
95			pressure_ext = pressure * pressure_units;
96			volume = boxxboxyboxz;
97			ke_temp = ke * keconverter;
98			NkBT = degrees_freedomkBtemp;
99
100			// propogate the strain rate
101
102			epsilon_dot += dt * ( (p_mol - pressure_ext)*volume
103	gezelter	454	/ (tau_relaxtau_relax kB * targetTemp) );
104	gezelter	453
105			// determine the change in cell volume
106			scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0));
107
108			volume = volume * pow(scale, 3.0);
109
110			// perform affine transform to update positions with volume fluctuations
111			affine_transform( scale );
112
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
122	gezelter	454	epsilonScale = epsilonDot * dt;
123	gezelter	453
124			// advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
125			// qmass is set in the parameter file
126			zeta += dt * ( (ke_temp*2 - NkBT) / qmass );
127			zetaScale = zeta * dt;
128
129			// apply barostating and thermostating to velocities and angular momenta
130	gezelter	454	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	gezelter	453	}
140	gezelter	454	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	gezelter	453	}
159
160	gezelter	454	void ExtendedSystem::AffineTransform( double scale ){
161	gezelter	453
162			int i;
163			double boxx_old, boxy_old, boxz_old, percentScale;
164			double boxx_num, boxy_num, boxz_num, rxi, ryi, rzi;
165			double[3] r;
166
167			// first determine the scaling factor from the box size change
168			percentScale = (boxx - boxx_old)/boxx_old;
169
170
171			for (i=0; i < nMols; i++) {
172
173			molecules[i]->getCOM(r);
174
175			// find the minimum image coordinates of the molecular centers of mass:
176
177
178			boxx_num = boxx_oldcopysign(1.0,r[0])(double)(int)(fabs(r[0]/boxx_old)+0.5);
179
180			boxx_num = boxx_olddsign(1.0d0,rx(i))int(abs(rx(i)/boxx_old)+0.5d0);
181			boxy_num = boxy_olddsign(1.0d0,ry(i))int(abs(ry(i)/boxy_old)+0.5d0);
182			boxz_num = boxz_olddsign(1.0d0,rz(i))int(abs(rz(i)/boxz_old)+0.5d0);
183
184			rxi = rx(i) - boxx_num;
185			ryi = ry(i) - boxy_num;
186			rzi = rz(i) - boxz_num;
187
188			// update the minimum image coordinates using the scaling factor
189			rxi = rxi + rxi*percentScale;
190			ryi = ryi + ryi*percentScale;
191			rzi = rzi + rzi*percentScale;
192
193			rx(i) = rxi + boxx_num;
194			ry(i) = ryi + boxy_num;
195			rz(i) = rzi + boxz_num;
196			}
197			}