mdtools/md_code/Verlet.cpp

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

#include "Atom.hpp"
#include "SRI.hpp"
#include "LRI.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;
  c_sr_interactions = info.sr_interactions;
  c_n_SRI = info.n_SRI;
  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[c_n_SRI];

  c_n_constrained = 0;
  int constrained = 0;
  
  for(int i = 0; i < c_n_SRI; i++){
    
    constrained = c_sr_interactions[i]->is_constrained();

    if(constrained){
      
      dummy_plug = c_sr_interactions[i]->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);
  
  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);
      
      // 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);
            
      // 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:	261
Committed:	Mon Feb 3 21:15:59 2003 UTC (21 years, 5 months ago) by chuckv
File size:	8722 byte(s)
Log Message:	We have working code today... MPI bug fixes to DumpWriter.
#	User	Rev	Content
1	mmeineke	10	#include <iostream>
2			#include <stdlib.h>
3
4			#include "Atom.hpp"
5			#include "SRI.hpp"
6			#include "LRI.hpp"
7			#include "Integrator.hpp"
8			#include "SimInfo.hpp"
9			#include "Thermo.hpp"
10			#include "ReadWrite.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	chuckv	249	Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
34	mmeineke	10
35			// get what information we need from the SimInfo object
36
37			entry_plug = &info;
38	chuckv	248	myFF = the_ff;
39	mmeineke	10
40	chuckv	248
41	mmeineke	10	c_natoms = info.n_atoms;
42			c_atoms = info.atoms;
43			c_sr_interactions = info.sr_interactions;
44			c_n_SRI = info.n_SRI;
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[c_n_SRI];
72
73			c_n_constrained = 0;
74			int constrained = 0;
75
76			for(int i = 0; i < c_n_SRI; i++){
77
78			constrained = c_sr_interactions[i]->is_constrained();
79
80			if(constrained){
81
82			dummy_plug = c_sr_interactions[i]->get_constraint();
83			temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
84			temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
85			temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
86
87			c_n_constrained++;
88			constrained = 0;
89			}
90			}
91
92			if(c_n_constrained > 0){
93
94			c_is_constrained = 1;
95			c_constrained_i = new int[c_n_constrained];
96			c_constrained_j = new int[c_n_constrained];
97			c_constrained_dsqr = new double[c_n_constrained];
98
99			for( int i = 0; i < c_n_constrained; i++){
100
101			/* add 1 to the index for the fortran arrays. */
102
103			c_constrained_i[i] = temp_con[i].get_a() + 1;
104			c_constrained_j[i] = temp_con[i].get_b() + 1;
105			c_constrained_dsqr[i] = temp_con[i].get_dsqr();
106			}
107			}
108
109			delete[] temp_con;
110			}
111
112
113			Verlet::~Verlet(){
114
115			if( c_is_constrained ){
116
117			delete[] c_constrained_i;
118			delete[] c_constrained_j;
119			delete[] c_constrained_dsqr;
120			}
121
122			delete[] c_mass;
123			c_mass = 0;
124			}
125
126
127			void Verlet::integrate( void ){
128
129			int i, j; /* loop counters */
130	chuckv	253	int calcPot;
131
132	mmeineke	10	double kE;
133
134			double *Rx = new double[c_natoms];
135			double *Ry = new double[c_natoms];
136			double *Rz = new double[c_natoms];
137
138			double *Vx = new double[c_natoms];
139			double *Vy = new double[c_natoms];
140			double *Vz = new double[c_natoms];
141
142			double *Fx = new double[c_natoms];
143			double *Fy = new double[c_natoms];
144			double *Fz = new double[c_natoms];
145
146	mmeineke	25	int time;
147
148	mmeineke	10	double dt = entry_plug->dt;
149			double runTime = entry_plug->run_time;
150			double sampleTime = entry_plug->sampleTime;
151			double statusTime = entry_plug->statusTime;
152			double thermalTime = entry_plug->thermalTime;
153
154			int n_loops = (int)( runTime / dt );
155			int sample_n = (int)( sampleTime / dt );
156			int status_n = (int)( statusTime / dt );
157			int vel_n = (int)( thermalTime / dt );
158
159			Thermo *tStats = new Thermo( entry_plug );
160
161			StatWriter* e_out = new StatWriter( entry_plug );
162			DumpWriter* dump_out = new DumpWriter( entry_plug );
163
164			// the first time integrate is called, the forces need to be initialized
165
166
167	chuckv	253	myFF->doForces(1);
168	mmeineke	10
169			if( entry_plug->setTemp ){
170			tStats->velocitize();
171			}
172
173	mmeineke	25	dump_out->writeDump( 0.0 );
174	chuckv	253
175	mmeineke	25	e_out->writeStat( 0.0 );
176	chuckv	261
177	chuckv	253	calcPot = 0;
178	mmeineke	25
179	mmeineke	10	if( c_is_constrained ){
180			for(i = 0; i < n_loops; i++){
181
182			// fill R, V, and F arrays and RATTLE in fortran
183
184			for( j=0; j<c_natoms; j++ ){
185
186			Rx[j] = c_atoms[j]->getX();
187			Ry[j] = c_atoms[j]->getY();
188			Rz[j] = c_atoms[j]->getZ();
189
190			Vx[j] = c_atoms[j]->get_vx();
191			Vy[j] = c_atoms[j]->get_vy();
192			Vz[j] = c_atoms[j]->get_vz();
193
194			Fx[j] = c_atoms[j]->getFx();
195			Fy[j] = c_atoms[j]->getFy();
196			Fz[j] = c_atoms[j]->getFz();
197
198			}
199
200			v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
201			Fx, Fy, Fz,
202			c_n_constrained, c_constrained_dsqr,
203			c_constrained_i, c_constrained_j,
204			c_box_x, c_box_y, c_box_z );
205
206			for( j=0; j<c_natoms; j++ ){
207
208			c_atoms[j]->setX(Rx[j]);
209			c_atoms[j]->setY(Ry[j]);
210			c_atoms[j]->setZ(Rz[j]);
211
212			c_atoms[j]->set_vx(Vx[j]);
213			c_atoms[j]->set_vy(Vy[j]);
214			c_atoms[j]->set_vz(Vz[j]);
215			}
216
217			// calculate the forces
218
219	chuckv	253	myFF->doForces(calcPot);
220	mmeineke	10
221			// finish the constrain move ( same as above. )
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			v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
239			Fx, Fy, Fz,
240			kE, c_n_constrained, c_constrained_dsqr,
241			c_constrained_i, c_constrained_j,
242			c_box_x, c_box_y, c_box_z );
243
244			for( j=0; j<c_natoms; j++ ){
245
246			c_atoms[j]->setX(Rx[j]);
247			c_atoms[j]->setY(Ry[j]);
248			c_atoms[j]->setZ(Rz[j]);
249
250			c_atoms[j]->set_vx(Vx[j]);
251			c_atoms[j]->set_vy(Vy[j]);
252			c_atoms[j]->set_vz(Vz[j]);
253			}
254
255	mmeineke	25	time = i + 1;
256
257	mmeineke	10	if( entry_plug->setTemp ){
258	mmeineke	25	if( !(time % vel_n) ) tStats->velocitize();
259	mmeineke	10	}
260	mmeineke	25	if( !(time % sample_n) ) dump_out->writeDump( time * dt );
261	chuckv	253	if( !((time+1) % status_n) ) calcPot = 1;
262			if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
263	mmeineke	10	}
264			}
265			else{
266			for(i = 0; i < n_loops; i++){
267
268			move_a( dt );
269
270			// calculate the forces
271
272	chuckv	253	myFF->doForces(calcPot);
273	mmeineke	10
274			// complete the verlet move
275
276			move_b( dt );
277
278	mmeineke	25	time = i + 1;
279
280	mmeineke	10	if( entry_plug->setTemp ){
281	mmeineke	25	if( !(time % vel_n) ) tStats->velocitize();
282	mmeineke	10	}
283	chuckv	253	if( !(time % sample_n) ) dump_out->writeDump( time * dt );
284			if( !((time+1) % status_n) ) calcPot = 1;
285			if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
286	mmeineke	10	}
287			}
288
289			dump_out->writeFinal();
290
291			delete dump_out;
292			delete e_out;
293
294			}
295
296
297			void Verlet::move_a(double dt){
298
299			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
300
301			double qx, qy, qz;
302			double vx, vy, vz;
303			int ma;
304			double h_dt = 0.5 * dt;
305			double h_dt2 = h_dt * dt;
306
307			for( ma = 0; ma < c_natoms; ma++){
308
309			qx = c_atoms[ma]->getX() + dt * c_atoms[ma]->get_vx() +
310			h_dt2 * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
311			qy = c_atoms[ma]->getY() + dt * c_atoms[ma]->get_vy() +
312			h_dt2 * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
313			qz = c_atoms[ma]->getZ() + dt * c_atoms[ma]->get_vz() +
314			h_dt2 * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
315
316			vx = c_atoms[ma]->get_vx() +
317			h_dt * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
318			vy = c_atoms[ma]->get_vy() +
319			h_dt * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
320			vz = c_atoms[ma]->get_vz() +
321			h_dt * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
322
323			c_atoms[ma]->setX(qx);
324			c_atoms[ma]->setY(qy);
325			c_atoms[ma]->setZ(qz);
326
327			c_atoms[ma]->set_vx(vx);
328			c_atoms[ma]->set_vy(vy);
329			c_atoms[ma]->set_vz(vz);
330			}
331			}
332
333			void Verlet::move_b( double dt ){
334
335			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
336
337			double vx, vy, vz;
338			int mb;
339			double h_dt = 0.5 * dt;
340
341
342			for( mb = 0; mb < c_natoms; mb++){
343
344			vx = c_atoms[mb]->get_vx() +
345			h_dt * c_atoms[mb]->getFx() * e_convert / c_atoms[mb]->getMass();
346			vy = c_atoms[mb]->get_vy() +
347			h_dt * c_atoms[mb]->getFy() * e_convert / c_atoms[mb]->getMass();
348			vz = c_atoms[mb]->get_vz() +
349			h_dt * c_atoms[mb]->getFz() * e_convert / c_atoms[mb]->getMass();
350
351			c_atoms[mb]->set_vx(vx);
352			c_atoms[mb]->set_vy(vy);
353			c_atoms[mb]->set_vz(vz);
354			}
355			}