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

  entry_plug = &info;

  c_natoms = info.n_atoms;
  c_atoms = info.atoms;
  c_sr_interactions = info.sr_interactions;
  longRange = info.longRange;
  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_b( double time_length, double dt, 
                            int n_bond_0, int n_bond_f,
                            int n_bend_0, int n_bend_f,
                            int n_torsion_0, int n_torsion_f,
                            bool do_bonds, bool do_bends, bool do_torsions,
                            bool do_LRI ){
  
//  double percent_tolerance = 0.001;
//  int max_iterations = 10000;

  int i, j; /* loop counters */
  double n_loops = time_length / dt;
  
  // the first time integrate is called, the forces need to be initialized 

  if(is_first){ 
    is_first = 0;
    
    for(i = 0; i < c_natoms; i++){
      c_atoms[i]->zeroForces();
    }
    
    if( do_bonds ){
      for(i = n_bond_0; i <= n_bond_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }

    if( do_bends ){
      for(i = n_bend_0; i <= n_bend_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }
    
    if( do_torsions ){
      for(i = n_torsion_0; i <= n_torsion_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }
 
    if( do_LRI ) longRange->calc_forces();
  }

  for(i = 0; i < n_loops; i++){
    
    move_a( dt );
    
    // calculate the forces
    
    for(j = 0; j < c_natoms; j++){
      c_atoms[j]->zeroForces();
    }
  

    if( do_bonds ){
      for(i = n_bond_0; i <= n_bond_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }

    if( do_bends ){
      for(i = n_bend_0; i <= n_bend_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }
    
    if( do_torsions ){
      for(i = n_torsion_0; i <= n_torsion_f; i++){
        c_sr_interactions[i]->calc_forces();
      }
    }

    if( do_LRI ) longRange->calc_forces();
    

    // complete the verlet move

    move_b( dt );
  }
}


void Verlet::integrate( void ){
  
  int i, j; /* loop counters */
  
  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];
  
  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 

  
  for(i = 0; i < c_natoms; i++){
    c_atoms[i]->zeroForces();
  }
  
  for(i = 0; i < c_n_SRI; i++){
    c_sr_interactions[i]->calc_forces();
  }
  
  longRange->calc_forces();
  
  if( entry_plug->setTemp ){
    tStats->velocitize();
  }
  
  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
      
      for(j = 0; j < c_natoms; j++){
        c_atoms[j]->zeroForces();
      }
      
      for(j = 0; j < c_n_SRI; j++){
        c_sr_interactions[j]->calc_forces();
      }
      
      longRange->calc_forces();
      
      // 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( entry_plug->setTemp ){
        if( !(i % vel_n) ) tStats->velocitize();
      }
      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
      if( !(i % status_n) ) e_out->writeStat( i * dt );

    }
  }
  else{
    for(i = 0; i < n_loops; i++){
      
      move_a( dt );
      
      // calculate the forces
      
      for(j = 0; j < c_natoms; j++){
        c_atoms[j]->zeroForces();
      }
      
      for(j = 0; j < c_n_SRI; j++){
        c_sr_interactions[j]->calc_forces();
      }
      
      longRange->calc_forces();
            
      // complete the verlet move
      
      move_b( dt );

      if( entry_plug->setTemp ){
        if( !(i % vel_n) ) tStats->velocitize();
      }
      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
      if( !(i % status_n) ) e_out->writeStat( i * dt );
    }
  }
  
  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;
  double v_sqr;
  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:	10
Committed:	Tue Jul 9 18:40:59 2002 UTC (21 years, 11 months ago) by mmeineke
Original Path:	*branches/mmeineke/mdtools/md_code/Verlet.cpp*
File size:	10560 byte(s)
Log Message:	everything you need to make libmdtools
#	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			Verlet::Verlet( SimInfo &info ){
34
35			// get what information we need from the SimInfo object
36
37			entry_plug = &info;
38
39			c_natoms = info.n_atoms;
40			c_atoms = info.atoms;
41			c_sr_interactions = info.sr_interactions;
42			longRange = info.longRange;
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
127			void Verlet::integrate_b( double time_length, double dt,
128			int n_bond_0, int n_bond_f,
129			int n_bend_0, int n_bend_f,
130			int n_torsion_0, int n_torsion_f,
131			bool do_bonds, bool do_bends, bool do_torsions,
132			bool do_LRI ){
133
134			// double percent_tolerance = 0.001;
135			// int max_iterations = 10000;
136
137			int i, j; /* loop counters */
138			double n_loops = time_length / dt;
139
140			// the first time integrate is called, the forces need to be initialized
141
142			if(is_first){
143			is_first = 0;
144
145			for(i = 0; i < c_natoms; i++){
146			c_atoms[i]->zeroForces();
147			}
148
149			if( do_bonds ){
150			for(i = n_bond_0; i <= n_bond_f; i++){
151			c_sr_interactions[i]->calc_forces();
152			}
153			}
154
155			if( do_bends ){
156			for(i = n_bend_0; i <= n_bend_f; i++){
157			c_sr_interactions[i]->calc_forces();
158			}
159			}
160
161			if( do_torsions ){
162			for(i = n_torsion_0; i <= n_torsion_f; i++){
163			c_sr_interactions[i]->calc_forces();
164			}
165			}
166
167			if( do_LRI ) longRange->calc_forces();
168			}
169
170			for(i = 0; i < n_loops; i++){
171
172			move_a( dt );
173
174			// calculate the forces
175
176			for(j = 0; j < c_natoms; j++){
177			c_atoms[j]->zeroForces();
178			}
179
180
181			if( do_bonds ){
182			for(i = n_bond_0; i <= n_bond_f; i++){
183			c_sr_interactions[i]->calc_forces();
184			}
185			}
186
187			if( do_bends ){
188			for(i = n_bend_0; i <= n_bend_f; i++){
189			c_sr_interactions[i]->calc_forces();
190			}
191			}
192
193			if( do_torsions ){
194			for(i = n_torsion_0; i <= n_torsion_f; i++){
195			c_sr_interactions[i]->calc_forces();
196			}
197			}
198
199			if( do_LRI ) longRange->calc_forces();
200
201
202			// complete the verlet move
203
204			move_b( dt );
205			}
206			}
207
208
209			void Verlet::integrate( void ){
210
211			int i, j; /* loop counters */
212
213			double kE;
214
215			double *Rx = new double[c_natoms];
216			double *Ry = new double[c_natoms];
217			double *Rz = new double[c_natoms];
218
219			double *Vx = new double[c_natoms];
220			double *Vy = new double[c_natoms];
221			double *Vz = new double[c_natoms];
222
223			double *Fx = new double[c_natoms];
224			double *Fy = new double[c_natoms];
225			double *Fz = new double[c_natoms];
226
227			double dt = entry_plug->dt;
228			double runTime = entry_plug->run_time;
229			double sampleTime = entry_plug->sampleTime;
230			double statusTime = entry_plug->statusTime;
231			double thermalTime = entry_plug->thermalTime;
232
233			int n_loops = (int)( runTime / dt );
234			int sample_n = (int)( sampleTime / dt );
235			int status_n = (int)( statusTime / dt );
236			int vel_n = (int)( thermalTime / dt );
237
238			Thermo *tStats = new Thermo( entry_plug );
239
240			StatWriter* e_out = new StatWriter( entry_plug );
241			DumpWriter* dump_out = new DumpWriter( entry_plug );
242
243			// the first time integrate is called, the forces need to be initialized
244
245
246			for(i = 0; i < c_natoms; i++){
247			c_atoms[i]->zeroForces();
248			}
249
250			for(i = 0; i < c_n_SRI; i++){
251			c_sr_interactions[i]->calc_forces();
252			}
253
254			longRange->calc_forces();
255
256			if( entry_plug->setTemp ){
257			tStats->velocitize();
258			}
259
260			if( c_is_constrained ){
261			for(i = 0; i < n_loops; i++){
262
263			// fill R, V, and F arrays and RATTLE in fortran
264
265			for( j=0; j<c_natoms; j++ ){
266
267			Rx[j] = c_atoms[j]->getX();
268			Ry[j] = c_atoms[j]->getY();
269			Rz[j] = c_atoms[j]->getZ();
270
271			Vx[j] = c_atoms[j]->get_vx();
272			Vy[j] = c_atoms[j]->get_vy();
273			Vz[j] = c_atoms[j]->get_vz();
274
275			Fx[j] = c_atoms[j]->getFx();
276			Fy[j] = c_atoms[j]->getFy();
277			Fz[j] = c_atoms[j]->getFz();
278
279			}
280
281			v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
282			Fx, Fy, Fz,
283			c_n_constrained, c_constrained_dsqr,
284			c_constrained_i, c_constrained_j,
285			c_box_x, c_box_y, c_box_z );
286
287			for( j=0; j<c_natoms; j++ ){
288
289			c_atoms[j]->setX(Rx[j]);
290			c_atoms[j]->setY(Ry[j]);
291			c_atoms[j]->setZ(Rz[j]);
292
293			c_atoms[j]->set_vx(Vx[j]);
294			c_atoms[j]->set_vy(Vy[j]);
295			c_atoms[j]->set_vz(Vz[j]);
296			}
297
298			// calculate the forces
299
300			for(j = 0; j < c_natoms; j++){
301			c_atoms[j]->zeroForces();
302			}
303
304			for(j = 0; j < c_n_SRI; j++){
305			c_sr_interactions[j]->calc_forces();
306			}
307
308			longRange->calc_forces();
309
310			// finish the constrain move ( same as above. )
311
312			for( j=0; j<c_natoms; j++ ){
313
314			Rx[j] = c_atoms[j]->getX();
315			Ry[j] = c_atoms[j]->getY();
316			Rz[j] = c_atoms[j]->getZ();
317
318			Vx[j] = c_atoms[j]->get_vx();
319			Vy[j] = c_atoms[j]->get_vy();
320			Vz[j] = c_atoms[j]->get_vz();
321
322			Fx[j] = c_atoms[j]->getFx();
323			Fy[j] = c_atoms[j]->getFy();
324			Fz[j] = c_atoms[j]->getFz();
325			}
326
327			v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
328			Fx, Fy, Fz,
329			kE, c_n_constrained, c_constrained_dsqr,
330			c_constrained_i, c_constrained_j,
331			c_box_x, c_box_y, c_box_z );
332
333			for( j=0; j<c_natoms; j++ ){
334
335			c_atoms[j]->setX(Rx[j]);
336			c_atoms[j]->setY(Ry[j]);
337			c_atoms[j]->setZ(Rz[j]);
338
339			c_atoms[j]->set_vx(Vx[j]);
340			c_atoms[j]->set_vy(Vy[j]);
341			c_atoms[j]->set_vz(Vz[j]);
342			}
343
344			if( entry_plug->setTemp ){
345			if( !(i % vel_n) ) tStats->velocitize();
346			}
347			if( !(i % sample_n) ) dump_out->writeDump( i * dt );
348			if( !(i % status_n) ) e_out->writeStat( i * dt );
349
350			}
351			}
352			else{
353			for(i = 0; i < n_loops; i++){
354
355			move_a( dt );
356
357			// calculate the forces
358
359			for(j = 0; j < c_natoms; j++){
360			c_atoms[j]->zeroForces();
361			}
362
363			for(j = 0; j < c_n_SRI; j++){
364			c_sr_interactions[j]->calc_forces();
365			}
366
367			longRange->calc_forces();
368
369			// complete the verlet move
370
371			move_b( dt );
372
373			if( entry_plug->setTemp ){
374			if( !(i % vel_n) ) tStats->velocitize();
375			}
376			if( !(i % sample_n) ) dump_out->writeDump( i * dt );
377			if( !(i % status_n) ) e_out->writeStat( i * dt );
378			}
379			}
380
381			dump_out->writeFinal();
382
383			delete dump_out;
384			delete e_out;
385
386			}
387
388
389			void Verlet::move_a(double dt){
390
391			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
392
393			double qx, qy, qz;
394			double vx, vy, vz;
395			int ma;
396			double h_dt = 0.5 * dt;
397			double h_dt2 = h_dt * dt;
398
399			for( ma = 0; ma < c_natoms; ma++){
400
401			qx = c_atoms[ma]->getX() + dt * c_atoms[ma]->get_vx() +
402			h_dt2 * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
403			qy = c_atoms[ma]->getY() + dt * c_atoms[ma]->get_vy() +
404			h_dt2 * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
405			qz = c_atoms[ma]->getZ() + dt * c_atoms[ma]->get_vz() +
406			h_dt2 * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
407
408			vx = c_atoms[ma]->get_vx() +
409			h_dt * c_atoms[ma]->getFx() * e_convert / c_atoms[ma]->getMass();
410			vy = c_atoms[ma]->get_vy() +
411			h_dt * c_atoms[ma]->getFy() * e_convert / c_atoms[ma]->getMass();
412			vz = c_atoms[ma]->get_vz() +
413			h_dt * c_atoms[ma]->getFz() * e_convert / c_atoms[ma]->getMass();
414
415			c_atoms[ma]->setX(qx);
416			c_atoms[ma]->setY(qy);
417			c_atoms[ma]->setZ(qz);
418
419			c_atoms[ma]->set_vx(vx);
420			c_atoms[ma]->set_vy(vy);
421			c_atoms[ma]->set_vz(vz);
422			}
423			}
424
425			void Verlet::move_b( double dt ){
426
427			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
428
429			double vx, vy, vz;
430			double v_sqr;
431			int mb;
432			double h_dt = 0.5 * dt;
433
434
435			for( mb = 0; mb < c_natoms; mb++){
436
437			vx = c_atoms[mb]->get_vx() +
438			h_dt * c_atoms[mb]->getFx() * e_convert / c_atoms[mb]->getMass();
439			vy = c_atoms[mb]->get_vy() +
440			h_dt * c_atoms[mb]->getFy() * e_convert / c_atoms[mb]->getMass();
441			vz = c_atoms[mb]->get_vz() +
442			h_dt * c_atoms[mb]->getFz() * e_convert / c_atoms[mb]->getMass();
443
444			c_atoms[mb]->set_vx(vx);
445			c_atoms[mb]->set_vy(vy);
446			c_atoms[mb]->set_vz(vz);
447			}
448			}