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 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 , 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 , tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") )
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     tStats->getPressure());

      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;
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
      }
    }
  }
  else{
    for(i = 0; i < n_loops; i++){
 
      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt , 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 , tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") )
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     tStats->getPressure());

      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;
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 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:	475
Committed:	Tue Apr 8 12:44:18 2003 UTC (21 years, 3 months ago) by gezelter
File size:	11413 byte(s)
Log Message:	Changes to integrate the NVT and NPT ensembles
#	User	Rev	Content
1	mmeineke	377	#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	gezelter	466	#include "ExtendedSystem.hpp"
11	mmeineke	377
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	gezelter	466	Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){
34	mmeineke	377
35			// get what information we need from the SimInfo object
36
37			entry_plug = &info;
38			myFF = the_ff;
39	gezelter	466	myES = the_es;
40	mmeineke	423
41	mmeineke	377	c_natoms = info.n_atoms;
42			c_atoms = info.atoms;
43	mmeineke	423	nMols = info.n_mol;
44			molecules = info.molecules;
45	mmeineke	377	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	mmeineke	423	temp_con = new Constraint[info.n_SRI];
72	mmeineke	377
73			c_n_constrained = 0;
74			int constrained = 0;
75	mmeineke	423	SRI** theArray;
76			for(int i = 0; i < nMols; i++){
77	mmeineke	377
78	mmeineke	428	theArray = (SRI**) molecules[i].getMyBonds();
79			for(int j=0; j<molecules[i].getNBonds(); j++){
80	mmeineke	423
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	mmeineke	377
95	mmeineke	428	theArray = (SRI**) molecules[i].getMyBends();
96			for(int j=0; j<molecules[i].getNBends(); j++){
97	mmeineke	377
98	mmeineke	423	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	mmeineke	377
112	mmeineke	428	theArray = (SRI**) molecules[i].getMyTorsions();
113			for(int j=0; j<molecules[i].getNTorsions(); j++){
114	mmeineke	423
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	mmeineke	377	}
128	mmeineke	423
129
130	mmeineke	377	}
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	mmeineke	423
139	mmeineke	377	for( int i = 0; i < c_n_constrained; i++){
140
141			/* add 1 to the index for the fortran arrays. */
142	mmeineke	423
143	mmeineke	377	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	gezelter	468	int calcPot, calcStress;
171	mmeineke	377
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 dt = entry_plug->dt;
189			double runTime = entry_plug->run_time;
190			double sampleTime = entry_plug->sampleTime;
191			double statusTime = entry_plug->statusTime;
192			double thermalTime = entry_plug->thermalTime;
193
194			int n_loops = (int)( runTime / dt );
195			int sample_n = (int)( sampleTime / dt );
196			int status_n = (int)( statusTime / dt );
197			int vel_n = (int)( thermalTime / dt );
198
199			Thermo *tStats = new Thermo( entry_plug );
200
201			StatWriter* e_out = new StatWriter( entry_plug );
202			DumpWriter* dump_out = new DumpWriter( entry_plug );
203
204			// the first time integrate is called, the forces need to be initialized
205
206
207	gezelter	468	myFF->doForces(1,1);
208	mmeineke	377
209			if( entry_plug->setTemp ){
210			tStats->velocitize();
211			}
212
213			dump_out->writeDump( 0.0 );
214
215			e_out->writeStat( 0.0 );
216
217			calcPot = 0;
218
219	gezelter	475	if (!strcasecmp( entry_plug->ensemble, "NPT")) {
220			calcStress = 1;
221			} else {
222			calcStress = 0;
223			}
224
225	mmeineke	377	if( c_is_constrained ){
226			for(i = 0; i < n_loops; i++){
227
228	gezelter	471	if (!strcasecmp( entry_plug->ensemble, "NVT"))
229			myES->NoseHooverNVT( dt , tStats->getKinetic() );
230
231	mmeineke	377	// fill R, V, and F arrays and RATTLE in fortran
232
233			for( j=0; j<c_natoms; j++ ){
234
235			Rx[j] = c_atoms[j]->getX();
236			Ry[j] = c_atoms[j]->getY();
237			Rz[j] = c_atoms[j]->getZ();
238
239			Vx[j] = c_atoms[j]->get_vx();
240			Vy[j] = c_atoms[j]->get_vy();
241			Vz[j] = c_atoms[j]->get_vz();
242
243			Fx[j] = c_atoms[j]->getFx();
244			Fy[j] = c_atoms[j]->getFy();
245			Fz[j] = c_atoms[j]->getFz();
246
247			}
248
249			v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
250			Fx, Fy, Fz,
251			c_n_constrained, c_constrained_dsqr,
252			c_constrained_i, c_constrained_j,
253			c_box_x, c_box_y, c_box_z );
254
255			for( j=0; j<c_natoms; j++ ){
256
257			c_atoms[j]->setX(Rx[j]);
258			c_atoms[j]->setY(Ry[j]);
259			c_atoms[j]->setZ(Rz[j]);
260
261			c_atoms[j]->set_vx(Vx[j]);
262			c_atoms[j]->set_vy(Vy[j]);
263			c_atoms[j]->set_vz(Vz[j]);
264			}
265
266			// calculate the forces
267
268	gezelter	468	myFF->doForces(calcPot,calcStress);
269	mmeineke	377
270			// finish the constrain move ( same as above. )
271
272			for( j=0; j<c_natoms; j++ ){
273
274			Rx[j] = c_atoms[j]->getX();
275			Ry[j] = c_atoms[j]->getY();
276			Rz[j] = c_atoms[j]->getZ();
277
278			Vx[j] = c_atoms[j]->get_vx();
279			Vy[j] = c_atoms[j]->get_vy();
280			Vz[j] = c_atoms[j]->get_vz();
281
282			Fx[j] = c_atoms[j]->getFx();
283			Fy[j] = c_atoms[j]->getFy();
284			Fz[j] = c_atoms[j]->getFz();
285			}
286
287	gezelter	471
288	mmeineke	377	v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
289			Fx, Fy, Fz,
290			kE, c_n_constrained, c_constrained_dsqr,
291			c_constrained_i, c_constrained_j,
292			c_box_x, c_box_y, c_box_z );
293
294			for( j=0; j<c_natoms; j++ ){
295
296			c_atoms[j]->setX(Rx[j]);
297			c_atoms[j]->setY(Ry[j]);
298			c_atoms[j]->setZ(Rz[j]);
299
300			c_atoms[j]->set_vx(Vx[j]);
301			c_atoms[j]->set_vy(Vy[j]);
302			c_atoms[j]->set_vz(Vz[j]);
303			}
304
305	gezelter	471	if (!strcasecmp( entry_plug->ensemble, "NVT"))
306			myES->NoseHooverNVT( dt , tStats->getKinetic() );
307
308			if (!strcasecmp( entry_plug->ensemble, "NPT") )
309			myES->NoseHooverAndersonNPT( dt,
310			tStats->getKinetic(),
311			tStats->getPressure());
312
313	mmeineke	377	time = i + 1;
314
315			if( entry_plug->setTemp ){
316			if( !(time % vel_n) ) tStats->velocitize();
317			}
318			if( !(time % sample_n) ) dump_out->writeDump( time * dt );
319	gezelter	468	if( !((time+1) % status_n) ) {
320			calcPot = 1;
321	gezelter	475	// bitwise masking in case we need it for NPT
322			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
323	gezelter	468	}
324			if( !(time % status_n) ){
325			e_out->writeStat( time * dt );
326			calcPot = 0;
327	gezelter	475	// bitwise masking in case we need it for NPT
328			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
329	gezelter	468	}
330	mmeineke	377	}
331			}
332			else{
333			for(i = 0; i < n_loops; i++){
334	gezelter	471
335			if (!strcasecmp( entry_plug->ensemble, "NVT"))
336			myES->NoseHooverNVT( dt , tStats->getKinetic() );
337
338	mmeineke	377	move_a( dt );
339
340			// calculate the forces
341
342	gezelter	468	myFF->doForces(calcPot,calcStress);
343	mmeineke	377
344			// complete the verlet move
345
346			move_b( dt );
347
348	gezelter	471	if (!strcasecmp( entry_plug->ensemble, "NVT"))
349			myES->NoseHooverNVT( dt , tStats->getKinetic() );
350
351			if (!strcasecmp( entry_plug->ensemble, "NPT") )
352			myES->NoseHooverAndersonNPT( dt,
353			tStats->getKinetic(),
354			tStats->getPressure());
355
356	mmeineke	377	time = i + 1;
357
358			if( entry_plug->setTemp ){
359			if( !(time % vel_n) ) tStats->velocitize();
360			}
361			if( !(time % sample_n) ) dump_out->writeDump( time * dt );
362	gezelter	468	if( !((time+1) % status_n) ) {
363			calcPot = 1;
364	gezelter	475	// bitwise masking in case we need it for NPT
365			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
366	gezelter	468	}
367			if( !(time % status_n) ){
368			e_out->writeStat( time * dt );
369			calcPot = 0;
370	gezelter	475	// bitwise masking in case we need it for NPT
371			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
372	gezelter	468	}
373	mmeineke	377	}
374			}
375
376			dump_out->writeFinal();
377
378			delete dump_out;
379			delete e_out;
380
381			}
382
383
384			void Verlet::move_a(double dt){
385
386			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
387
388			double qx, qy, qz;
389			double vx, vy, vz;
390			int ma;
391			double h_dt = 0.5 * dt;
392			double h_dt2 = h_dt * dt;
393
394			for( ma = 0; ma < c_natoms; ma++){
395
396			qx = c_atoms[ma]->getX() + dt * c_atoms[ma]->get_vx() +
397			h_dt2 * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
398			qy = c_atoms[ma]->getY() + dt * c_atoms[ma]->get_vy() +
399			h_dt2 * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
400			qz = c_atoms[ma]->getZ() + dt * c_atoms[ma]->get_vz() +
401			h_dt2 * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
402
403			vx = c_atoms[ma]->get_vx() +
404			h_dt * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
405			vy = c_atoms[ma]->get_vy() +
406			h_dt * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
407			vz = c_atoms[ma]->get_vz() +
408			h_dt * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
409
410			c_atoms[ma]->setX(qx);
411			c_atoms[ma]->setY(qy);
412			c_atoms[ma]->setZ(qz);
413
414			c_atoms[ma]->set_vx(vx);
415			c_atoms[ma]->set_vy(vy);
416			c_atoms[ma]->set_vz(vz);
417			}
418			}
419
420			void Verlet::move_b( double dt ){
421
422			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
423
424			double vx, vy, vz;
425			int mb;
426			double h_dt = 0.5 * dt;
427
428
429			for( mb = 0; mb < c_natoms; mb++){
430
431			vx = c_atoms[mb]->get_vx() +
432			h_dt * c_atoms[mb]->getFx() * e_convert / c_atoms[mb]->getMass();
433			vy = c_atoms[mb]->get_vy() +
434			h_dt * c_atoms[mb]->getFy() * e_convert / c_atoms[mb]->getMass();
435			vz = c_atoms[mb]->get_vz() +
436			h_dt * c_atoms[mb]->getFz() * e_convert / c_atoms[mb]->getMass();
437
438			c_atoms[mb]->set_vx(vx);
439			c_atoms[mb]->set_vy(vy);
440			c_atoms[mb]->set_vz(vz);
441			}
442			}