OOPSE/libmdtools/Verlet.cpp

#include <iostream>
#include <stdlib.h>

#include "Atom.hpp"
#include "SRI.hpp"
#include "Integrator.hpp"
#include "SimInfo.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "ExtendedSystem.hpp"

extern "C"{
  
  void v_constrain_a_( double &dt, int &n_atoms, double* mass, 
                       double* Rx, double* Ry, double* Rz, 
                       double* Vx, double* Vy, double* Vz, 
                       double* Fx, double* Fy, double* Fz, 
                       int &n_constrained, double *constr_sqr,
                       int* constr_i, int* constr_j, 
                       double &box_x, double &box_y, double &box_z );

  void v_constrain_b_( double &dt, int &n_atoms, double* mass, 
                       double* Rx, double* Ry, double* Rz, 
                       double* Vx, double* Vy, double* Vz, 
                       double* Fx, double* Fy, double* Fz, 
                       double &Kinetic, 
                       int &n_constrained, double *constr_sqr,
                       int* constr_i, int* constr_j, 
                       double &box_x, double &box_y, double &box_z );
}

  
Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){
  
  // get what information we need from the SimInfo object

  entry_plug = &info;
  myFF = the_ff;
  myES = the_es;
  
  c_natoms = info.n_atoms;
  c_atoms = info.atoms;
  nMols = info.n_mol;
  molecules = info.molecules;
  c_is_constrained = 0;
  c_box_x = info.box_x;
  c_box_y = info.box_y;
  c_box_z = info.box_z;
  
  // give a little love back to the SimInfo object

  if( info.the_integrator != NULL ) delete info.the_integrator;
  info.the_integrator = this;

  // the rest are initialization issues
  
  is_first = 1; // let the integrate method know when the first call is

  // mass array setup
  
  c_mass = new double[c_natoms];

  for(int i = 0; i < c_natoms; i++){
    c_mass[i] = c_atoms[i]->getMass();
  }
  
  // check for constraints

  Constraint *temp_con;
  Constraint *dummy_plug;
  temp_con = new Constraint[info.n_SRI];

  c_n_constrained = 0;
  int constrained = 0;
  SRI** theArray;
  for(int i = 0; i < nMols; i++){
    
    theArray = (SRI**) molecules[i].getMyBonds();
    for(int j=0; j<molecules[i].getNBonds(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        c_n_constrained++;
        constrained = 0;
      }
    }

    theArray = (SRI**) molecules[i].getMyBends();
    for(int j=0; j<molecules[i].getNBends(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        c_n_constrained++;
        constrained = 0;
      }
    }

    theArray = (SRI**) molecules[i].getMyTorsions();
    for(int j=0; j<molecules[i].getNTorsions(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        c_n_constrained++;
        constrained = 0;
      }
    }


  }

  if(c_n_constrained > 0){
    
    c_is_constrained = 1;
    c_constrained_i = new int[c_n_constrained];
    c_constrained_j = new int[c_n_constrained];
    c_constrained_dsqr = new double[c_n_constrained];
    
    for( int i = 0; i < c_n_constrained; i++){
      
      /* add 1 to the index for the fortran arrays. */
      
      c_constrained_i[i] = temp_con[i].get_a() + 1;
      c_constrained_j[i] = temp_con[i].get_b() + 1;
      c_constrained_dsqr[i] = temp_con[i].get_dsqr();
    }
  }
  
  delete[] temp_con;
}


Verlet::~Verlet(){
 
  if( c_is_constrained ){
    
    delete[] c_constrained_i;
    delete[] c_constrained_j;
    delete[] c_constrained_dsqr;
  }
  
  delete[] c_mass;
  c_mass = 0;
}


void Verlet::integrate( void ){
  
  int i, j; /* loop counters */
  int calcPot, calcStress;

  double kE;

  double *Rx = new double[c_natoms];
  double *Ry = new double[c_natoms];
  double *Rz = new double[c_natoms];

  double *Vx = new double[c_natoms];
  double *Vy = new double[c_natoms];
  double *Vz = new double[c_natoms];
  
  double *Fx = new double[c_natoms];
  double *Fy = new double[c_natoms];
  double *Fz = new double[c_natoms];
  
  int time;

  double press[9];

  double dt = entry_plug->dt;
  double runTime = entry_plug->run_time;
  double sampleTime = entry_plug->sampleTime;
  double statusTime = entry_plug->statusTime;
  double thermalTime = entry_plug->thermalTime;

  int n_loops  = (int)( runTime / dt );
  int sample_n = (int)( sampleTime / dt );
  int status_n = (int)( statusTime / dt );
  int vel_n    = (int)( thermalTime / dt );

  Thermo *tStats = new Thermo( entry_plug );

  StatWriter*  e_out    = new StatWriter( entry_plug );
  DumpWriter*  dump_out = new DumpWriter( entry_plug );

  // the first time integrate is called, the forces need to be initialized 
  
  myFF->doForces(1,1);
  
  if( entry_plug->setTemp ){
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );

  e_out->writeStat( 0.0 );

  calcPot = 0;

  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
    calcStress = 1;
  } else {
    calcStress = 0;
  }

  if( c_is_constrained ){
    for(i = 0; i < n_loops; i++){
      
      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
      
      // fill R, V, and F arrays and RATTLE in fortran
      
      for( j=0; j<c_natoms; j++ ){
        
        Rx[j] = c_atoms[j]->getX();
        Ry[j] = c_atoms[j]->getY();
        Rz[j] = c_atoms[j]->getZ();

        Vx[j] = c_atoms[j]->get_vx();
        Vy[j] = c_atoms[j]->get_vy();
        Vz[j] = c_atoms[j]->get_vz();

        Fx[j] = c_atoms[j]->getFx();
        Fy[j] = c_atoms[j]->getFy();
        Fz[j] = c_atoms[j]->getFz();

      }
      
      v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz, 
                      Fx, Fy, Fz,
                      c_n_constrained, c_constrained_dsqr, 
                      c_constrained_i, c_constrained_j,
                      c_box_x, c_box_y, c_box_z );

      for( j=0; j<c_natoms; j++ ){

        c_atoms[j]->setX(Rx[j]);
        c_atoms[j]->setY(Ry[j]);
        c_atoms[j]->setZ(Rz[j]);

        c_atoms[j]->set_vx(Vx[j]);
        c_atoms[j]->set_vy(Vy[j]);
        c_atoms[j]->set_vz(Vz[j]);
      }

      // calculate the forces
      
      myFF->doForces(calcPot,calcStress);
      
      // finish the constrain move ( same as above. )

      for( j=0; j<c_natoms; j++ ){

        Rx[j] = c_atoms[j]->getX();
        Ry[j] = c_atoms[j]->getY();
        Rz[j] = c_atoms[j]->getZ();

        Vx[j] = c_atoms[j]->get_vx();
        Vy[j] = c_atoms[j]->get_vy();
        Vz[j] = c_atoms[j]->get_vz();

        Fx[j] = c_atoms[j]->getFx();
        Fy[j] = c_atoms[j]->getFy();
        Fz[j] = c_atoms[j]->getFz();
      }
        

      v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz, 
                      Fx, Fy, Fz,
                      kE, c_n_constrained, c_constrained_dsqr, 
                      c_constrained_i, c_constrained_j,
                      c_box_x, c_box_y, c_box_z );
      
      for( j=0; j<c_natoms; j++ ){

        c_atoms[j]->setX(Rx[j]);
        c_atoms[j]->setY(Ry[j]);
        c_atoms[j]->setZ(Rz[j]);

        c_atoms[j]->set_vx(Vx[j]);
        c_atoms[j]->set_vy(Vy[j]);
        c_atoms[j]->set_vz(Vz[j]);
      }
      
      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
      
      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {        
        tStats->getPressureTensor(press);
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     press);
      }

      time = i + 1;
      
      if( entry_plug->setTemp ){
        if( !(time % vel_n) ) tStats->velocitize();
      }
      if( !(time % sample_n) ) dump_out->writeDump( time * dt );

      if( !((time+1) % status_n) ) {
        calcPot = 1;
        calcStress = 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
        else calcStress = 0;
      }
    }
  }
  else{
    for(i = 0; i < n_loops; i++){
 
      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
     
      move_a( dt );
      
      // calculate the forces
      
      myFF->doForces(calcPot,calcStress);
            
      // complete the verlet move
      
      move_b( dt );

      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
        tStats->getPressureTensor(press);
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     press);
      }

      time = i + 1;
      
      if( entry_plug->setTemp ){
        if( !(time % vel_n) ) tStats->velocitize();
      }
      if( !(time % sample_n) )  dump_out->writeDump( time * dt );
      if( !((time+1) % status_n) ) {
        calcPot = 1;
        calcStress = 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
        else calcStress = 0;
      }
    }
  }
  
  dump_out->writeFinal();

  delete dump_out;
  delete e_out;
  
}


void Verlet::move_a(double dt){

  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
  
  double qx, qy, qz;
  double vx, vy, vz;
  int ma;
  double h_dt = 0.5 * dt;
  double h_dt2 = h_dt * dt;
  
  for( ma = 0; ma < c_natoms; ma++){

    qx = c_atoms[ma]->getX() + dt * c_atoms[ma]->get_vx() + 
      h_dt2 * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
    qy = c_atoms[ma]->getY() + dt * c_atoms[ma]->get_vy() + 
      h_dt2 * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
    qz = c_atoms[ma]->getZ() + dt * c_atoms[ma]->get_vz() + 
      h_dt2 * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();

    vx = c_atoms[ma]->get_vx() + 
      h_dt * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
    vy = c_atoms[ma]->get_vy() + 
      h_dt * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
    vz = c_atoms[ma]->get_vz() + 
      h_dt * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();

    c_atoms[ma]->setX(qx);
    c_atoms[ma]->setY(qy);
    c_atoms[ma]->setZ(qz);

    c_atoms[ma]->set_vx(vx);
    c_atoms[ma]->set_vy(vy);
    c_atoms[ma]->set_vz(vz);
  }
}
    
void Verlet::move_b( double dt ){

  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
  
  double vx, vy, vz;
  int mb;
  double h_dt = 0.5 * dt;
  
  
  for( mb = 0; mb < c_natoms; mb++){

    vx = c_atoms[mb]->get_vx() + 
      h_dt * c_atoms[mb]->getFx() * e_convert / c_atoms[mb]->getMass();
    vy = c_atoms[mb]->get_vy() + 
      h_dt * c_atoms[mb]->getFy() * e_convert / c_atoms[mb]->getMass();
    vz = c_atoms[mb]->get_vz() + 
      h_dt * c_atoms[mb]->getFz() * e_convert / c_atoms[mb]->getMass();

    c_atoms[mb]->set_vx(vx);
    c_atoms[mb]->set_vy(vy);
    c_atoms[mb]->set_vz(vz);
  }
}
Revision:	483
Committed:	Wed Apr 9 04:06:43 2003 UTC (21 years, 3 months ago) by gezelter
File size:	11273 byte(s)
Log Message:	fixes for NPT and NVT
#	Content
1	#include <iostream>
2	#include <stdlib.h>
3
4	#include "Atom.hpp"
5	#include "SRI.hpp"
6	#include "Integrator.hpp"
7	#include "SimInfo.hpp"
8	#include "Thermo.hpp"
9	#include "ReadWrite.hpp"
10	#include "ExtendedSystem.hpp"
11
12	extern "C"{
13
14	void v_constrain_a_( double &dt, int &n_atoms, double* mass,
15	double* Rx, double* Ry, double* Rz,
16	double* Vx, double* Vy, double* Vz,
17	double* Fx, double* Fy, double* Fz,
18	int &n_constrained, double *constr_sqr,
19	int* constr_i, int* constr_j,
20	double &box_x, double &box_y, double &box_z );
21
22	void v_constrain_b_( double &dt, int &n_atoms, double* mass,
23	double* Rx, double* Ry, double* Rz,
24	double* Vx, double* Vy, double* Vz,
25	double* Fx, double* Fy, double* Fz,
26	double &Kinetic,
27	int &n_constrained, double *constr_sqr,
28	int* constr_i, int* constr_j,
29	double &box_x, double &box_y, double &box_z );
30	}
31
32
33	Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){
34
35	// get what information we need from the SimInfo object
36
37	entry_plug = &info;
38	myFF = the_ff;
39	myES = the_es;
40
41	c_natoms = info.n_atoms;
42	c_atoms = info.atoms;
43	nMols = info.n_mol;
44	molecules = info.molecules;
45	c_is_constrained = 0;
46	c_box_x = info.box_x;
47	c_box_y = info.box_y;
48	c_box_z = info.box_z;
49
50	// give a little love back to the SimInfo object
51
52	if( info.the_integrator != NULL ) delete info.the_integrator;
53	info.the_integrator = this;
54
55	// the rest are initialization issues
56
57	is_first = 1; // let the integrate method know when the first call is
58
59	// mass array setup
60
61	c_mass = new double[c_natoms];
62
63	for(int i = 0; i < c_natoms; i++){
64	c_mass[i] = c_atoms[i]->getMass();
65	}
66
67	// check for constraints
68
69	Constraint *temp_con;
70	Constraint *dummy_plug;
71	temp_con = new Constraint[info.n_SRI];
72
73	c_n_constrained = 0;
74	int constrained = 0;
75	SRI** theArray;
76	for(int i = 0; i < nMols; i++){
77
78	theArray = (SRI**) molecules[i].getMyBonds();
79	for(int j=0; j<molecules[i].getNBonds(); j++){
80
81	constrained = theArray[j]->is_constrained();
82
83	if(constrained){
84
85	dummy_plug = theArray[j]->get_constraint();
86	temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
87	temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
88	temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
89
90	c_n_constrained++;
91	constrained = 0;
92	}
93	}
94
95	theArray = (SRI**) molecules[i].getMyBends();
96	for(int j=0; j<molecules[i].getNBends(); j++){
97
98	constrained = theArray[j]->is_constrained();
99
100	if(constrained){
101
102	dummy_plug = theArray[j]->get_constraint();
103	temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
104	temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
105	temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
106
107	c_n_constrained++;
108	constrained = 0;
109	}
110	}
111
112	theArray = (SRI**) molecules[i].getMyTorsions();
113	for(int j=0; j<molecules[i].getNTorsions(); j++){
114
115	constrained = theArray[j]->is_constrained();
116
117	if(constrained){
118
119	dummy_plug = theArray[j]->get_constraint();
120	temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
121	temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
122	temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
123
124	c_n_constrained++;
125	constrained = 0;
126	}
127	}
128
129
130	}
131
132	if(c_n_constrained > 0){
133
134	c_is_constrained = 1;
135	c_constrained_i = new int[c_n_constrained];
136	c_constrained_j = new int[c_n_constrained];
137	c_constrained_dsqr = new double[c_n_constrained];
138
139	for( int i = 0; i < c_n_constrained; i++){
140
141	/* add 1 to the index for the fortran arrays. */
142
143	c_constrained_i[i] = temp_con[i].get_a() + 1;
144	c_constrained_j[i] = temp_con[i].get_b() + 1;
145	c_constrained_dsqr[i] = temp_con[i].get_dsqr();
146	}
147	}
148
149	delete[] temp_con;
150	}
151
152
153	Verlet::~Verlet(){
154
155	if( c_is_constrained ){
156
157	delete[] c_constrained_i;
158	delete[] c_constrained_j;
159	delete[] c_constrained_dsqr;
160	}
161
162	delete[] c_mass;
163	c_mass = 0;
164	}
165
166
167	void Verlet::integrate( void ){
168
169	int i, j; /* loop counters */
170	int calcPot, calcStress;
171
172	double kE;
173
174	double *Rx = new double[c_natoms];
175	double *Ry = new double[c_natoms];
176	double *Rz = new double[c_natoms];
177
178	double *Vx = new double[c_natoms];
179	double *Vy = new double[c_natoms];
180	double *Vz = new double[c_natoms];
181
182	double *Fx = new double[c_natoms];
183	double *Fy = new double[c_natoms];
184	double *Fz = new double[c_natoms];
185
186	int time;
187
188	double press[9];
189
190	double dt = entry_plug->dt;
191	double runTime = entry_plug->run_time;
192	double sampleTime = entry_plug->sampleTime;
193	double statusTime = entry_plug->statusTime;
194	double thermalTime = entry_plug->thermalTime;
195
196	int n_loops = (int)( runTime / dt );
197	int sample_n = (int)( sampleTime / dt );
198	int status_n = (int)( statusTime / dt );
199	int vel_n = (int)( thermalTime / dt );
200
201	Thermo *tStats = new Thermo( entry_plug );
202
203	StatWriter* e_out = new StatWriter( entry_plug );
204	DumpWriter* dump_out = new DumpWriter( entry_plug );
205
206	// the first time integrate is called, the forces need to be initialized
207
208	myFF->doForces(1,1);
209
210	if( entry_plug->setTemp ){
211	tStats->velocitize();
212	}
213
214	dump_out->writeDump( 0.0 );
215
216	e_out->writeStat( 0.0 );
217
218	calcPot = 0;
219
220	if (!strcasecmp( entry_plug->ensemble, "NPT")) {
221	calcStress = 1;
222	} else {
223	calcStress = 0;
224	}
225
226	if( c_is_constrained ){
227	for(i = 0; i < n_loops; i++){
228
229	if (!strcasecmp( entry_plug->ensemble, "NVT"))
230	myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
231
232	// fill R, V, and F arrays and RATTLE in fortran
233
234	for( j=0; j<c_natoms; j++ ){
235
236	Rx[j] = c_atoms[j]->getX();
237	Ry[j] = c_atoms[j]->getY();
238	Rz[j] = c_atoms[j]->getZ();
239
240	Vx[j] = c_atoms[j]->get_vx();
241	Vy[j] = c_atoms[j]->get_vy();
242	Vz[j] = c_atoms[j]->get_vz();
243
244	Fx[j] = c_atoms[j]->getFx();
245	Fy[j] = c_atoms[j]->getFy();
246	Fz[j] = c_atoms[j]->getFz();
247
248	}
249
250	v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
251	Fx, Fy, Fz,
252	c_n_constrained, c_constrained_dsqr,
253	c_constrained_i, c_constrained_j,
254	c_box_x, c_box_y, c_box_z );
255
256	for( j=0; j<c_natoms; j++ ){
257
258	c_atoms[j]->setX(Rx[j]);
259	c_atoms[j]->setY(Ry[j]);
260	c_atoms[j]->setZ(Rz[j]);
261
262	c_atoms[j]->set_vx(Vx[j]);
263	c_atoms[j]->set_vy(Vy[j]);
264	c_atoms[j]->set_vz(Vz[j]);
265	}
266
267	// calculate the forces
268
269	myFF->doForces(calcPot,calcStress);
270
271	// finish the constrain move ( same as above. )
272
273	for( j=0; j<c_natoms; j++ ){
274
275	Rx[j] = c_atoms[j]->getX();
276	Ry[j] = c_atoms[j]->getY();
277	Rz[j] = c_atoms[j]->getZ();
278
279	Vx[j] = c_atoms[j]->get_vx();
280	Vy[j] = c_atoms[j]->get_vy();
281	Vz[j] = c_atoms[j]->get_vz();
282
283	Fx[j] = c_atoms[j]->getFx();
284	Fy[j] = c_atoms[j]->getFy();
285	Fz[j] = c_atoms[j]->getFz();
286	}
287
288
289	v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
290	Fx, Fy, Fz,
291	kE, c_n_constrained, c_constrained_dsqr,
292	c_constrained_i, c_constrained_j,
293	c_box_x, c_box_y, c_box_z );
294
295	for( j=0; j<c_natoms; j++ ){
296
297	c_atoms[j]->setX(Rx[j]);
298	c_atoms[j]->setY(Ry[j]);
299	c_atoms[j]->setZ(Rz[j]);
300
301	c_atoms[j]->set_vx(Vx[j]);
302	c_atoms[j]->set_vy(Vy[j]);
303	c_atoms[j]->set_vz(Vz[j]);
304	}
305
306	if (!strcasecmp( entry_plug->ensemble, "NVT"))
307	myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
308
309	if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
310	tStats->getPressureTensor(press);
311	myES->NoseHooverAndersonNPT( dt,
312	tStats->getKinetic(),
313	press);
314	}
315
316	time = i + 1;
317
318	if( entry_plug->setTemp ){
319	if( !(time % vel_n) ) tStats->velocitize();
320	}
321	if( !(time % sample_n) ) dump_out->writeDump( time * dt );
322
323	if( !((time+1) % status_n) ) {
324	calcPot = 1;
325	calcStress = 1;
326	}
327	if( !(time % status_n) ){
328	e_out->writeStat( time * dt );
329	calcPot = 0;
330	if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
331	else calcStress = 0;
332	}
333	}
334	}
335	else{
336	for(i = 0; i < n_loops; i++){
337
338	if (!strcasecmp( entry_plug->ensemble, "NVT"))
339	myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
340
341	move_a( dt );
342
343	// calculate the forces
344
345	myFF->doForces(calcPot,calcStress);
346
347	// complete the verlet move
348
349	move_b( dt );
350
351	if (!strcasecmp( entry_plug->ensemble, "NVT"))
352	myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
353
354	if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
355	tStats->getPressureTensor(press);
356	myES->NoseHooverAndersonNPT( dt,
357	tStats->getKinetic(),
358	press);
359	}
360
361	time = i + 1;
362
363	if( entry_plug->setTemp ){
364	if( !(time % vel_n) ) tStats->velocitize();
365	}
366	if( !(time % sample_n) ) dump_out->writeDump( time * dt );
367	if( !((time+1) % status_n) ) {
368	calcPot = 1;
369	calcStress = 1;
370	}
371	if( !(time % status_n) ){
372	e_out->writeStat( time * dt );
373	calcPot = 0;
374	if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
375	else calcStress = 0;
376	}
377	}
378	}
379
380	dump_out->writeFinal();
381
382	delete dump_out;
383	delete e_out;
384
385	}
386
387
388	void Verlet::move_a(double dt){
389
390	const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
391
392	double qx, qy, qz;
393	double vx, vy, vz;
394	int ma;
395	double h_dt = 0.5 * dt;
396	double h_dt2 = h_dt * dt;
397
398	for( ma = 0; ma < c_natoms; ma++){
399
400	qx = c_atoms[ma]->getX() + dt * c_atoms[ma]->get_vx() +
401	h_dt2 * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
402	qy = c_atoms[ma]->getY() + dt * c_atoms[ma]->get_vy() +
403	h_dt2 * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
404	qz = c_atoms[ma]->getZ() + dt * c_atoms[ma]->get_vz() +
405	h_dt2 * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
406
407	vx = c_atoms[ma]->get_vx() +
408	h_dt * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
409	vy = c_atoms[ma]->get_vy() +
410	h_dt * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
411	vz = c_atoms[ma]->get_vz() +
412	h_dt * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
413
414	c_atoms[ma]->setX(qx);
415	c_atoms[ma]->setY(qy);
416	c_atoms[ma]->setZ(qz);
417
418	c_atoms[ma]->set_vx(vx);
419	c_atoms[ma]->set_vy(vy);
420	c_atoms[ma]->set_vz(vz);
421	}
422	}
423
424	void Verlet::move_b( double dt ){
425
426	const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
427
428	double vx, vy, vz;
429	int mb;
430	double h_dt = 0.5 * dt;
431
432
433	for( mb = 0; mb < c_natoms; mb++){
434
435	vx = c_atoms[mb]->get_vx() +
436	h_dt * c_atoms[mb]->getFx() * e_convert / c_atoms[mb]->getMass();
437	vy = c_atoms[mb]->get_vy() +
438	h_dt * c_atoms[mb]->getFy() * e_convert / c_atoms[mb]->getMass();
439	vz = c_atoms[mb]->get_vz() +
440	h_dt * c_atoms[mb]->getFz() * e_convert / c_atoms[mb]->getMass();
441
442	c_atoms[mb]->set_vx(vx);
443	c_atoms[mb]->set_vy(vy);
444	c_atoms[mb]->set_vz(vz);
445	}
446	}