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