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"

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 ){
  
  // get what information we need from the SimInfo object

  entry_plug = &info;
  myFF = the_ff;

  
  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;

  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,0);
  
  if( entry_plug->setTemp ){
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );

  e_out->writeStat( 0.0 );

  calcPot = 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();

      }
      
      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,0);
      
      // 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]);
      }
      
      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;
      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
    }
  }
  else{
    for(i = 0; i < n_loops; i++){
      
      move_a( dt );
      
      // calculate the forces
      
      myFF->doForces(calcPot,0);
            
      // complete the verlet move
      
      move_b( dt );

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