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"

#include "simError.h"

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,
                       int &isError );

  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,
                       int &isError);
}

  
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 );

  int isError;

  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++){
      
      
      // 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();

      }
      
      isError = 0;
      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,
                      isError);

      if( isError ){
        
        sprintf( painCave.errMsg,
                 "Constraint Failure in Fortran move A\n" );
        painCave.isFatal = 1;
        simError();
      }
      
#ifdef IS_MPI
      sprintf( checkPointMsg,
               "succesful return from move a.\n" );
      MPIcheckPoint();

#endif //is_mpi

      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() );
   
      // 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();
      }
        

      isError = 0;
      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, isError );

      if( isError ){
        
        sprintf( painCave.errMsg,
                 "Constraint Failure in Fortran move B\n" );
        painCave.isFatal = 1;
        simError();
      }

#ifdef IS_MPI
      sprintf( checkPointMsg,
               "succesful return from move B.\n" );
      MPIcheckPoint();

#endif //is_mpi      


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