ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/Verlet.cpp
(Generate patch)

Comparing:
branches/mmeineke/mdtools/md_code/Verlet.cpp (file contents), Revision 10 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
trunk/mdtools/md_code/Verlet.cpp (file contents), Revision 249 by chuckv, Mon Jan 27 21:28:19 2003 UTC

# Line 30 | Line 30 | Verlet::Verlet( SimInfo &info ){
30   }
31  
32    
33 < Verlet::Verlet( SimInfo &info ){
33 > Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
34    
35    // get what information we need from the SimInfo object
36  
37    entry_plug = &info;
38 +  myFF = the_ff;
39  
40 +
41    c_natoms = info.n_atoms;
42    c_atoms = info.atoms;
43    c_sr_interactions = info.sr_interactions;
42  longRange = info.longRange;
44    c_n_SRI = info.n_SRI;
45    c_is_constrained = 0;
46    c_box_x = info.box_x;
# Line 123 | Line 124 | Verlet::~Verlet(){
124   }
125  
126  
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
127   void Verlet::integrate( void ){
128    
129    int i, j; /* loop counters */
# Line 224 | Line 142 | void Verlet::integrate( void ){
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;
# Line 243 | Line 163 | void Verlet::integrate( void ){
163    // the first time integrate is called, the forces need to be initialized
164  
165    
166 <  for(i = 0; i < c_natoms; i++){
247 <    c_atoms[i]->zeroForces();
248 <  }
166 >  myFF->doForces();
167    
250  for(i = 0; i < c_n_SRI; i++){
251    c_sr_interactions[i]->calc_forces();
252  }
253  
254  longRange->calc_forces();
255  
168    if( entry_plug->setTemp ){
169      tStats->velocitize();
170    }
171    
172 +  dump_out->writeDump( 0.0 );
173 +  e_out->writeStat( 0.0 );
174 +
175    if( c_is_constrained ){
176      for(i = 0; i < n_loops; i++){
177        
# Line 297 | Line 212 | void Verlet::integrate( void ){
212  
213        // calculate the forces
214        
215 <      for(j = 0; j < c_natoms; j++){
301 <        c_atoms[j]->zeroForces();
302 <      }
215 >      myFF->doForces();
216        
304      for(j = 0; j < c_n_SRI; j++){
305        c_sr_interactions[j]->calc_forces();
306      }
307      
308      longRange->calc_forces();
309      
217        // finish the constrain move ( same as above. )
218  
219        for( j=0; j<c_natoms; j++ ){
# Line 341 | Line 248 | void Verlet::integrate( void ){
248          c_atoms[j]->set_vz(Vz[j]);
249        }
250        
251 +      time = i + 1;
252 +      
253        if( entry_plug->setTemp ){
254 <        if( !(i % vel_n) ) tStats->velocitize();
254 >        if( !(time % vel_n) ) tStats->velocitize();
255        }
256 <      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
257 <      if( !(i % status_n) ) e_out->writeStat( i * dt );
349 <
256 >      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
257 >      if( !(time % status_n) ) e_out->writeStat( time * dt );
258      }
259    }
260    else{
# Line 356 | Line 264 | void Verlet::integrate( void ){
264        
265        // calculate the forces
266        
267 <      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();
267 >      myFF->doForces();
268              
269        // complete the verlet move
270        
271        move_b( dt );
272  
273 +      time = i + 1;
274 +      
275        if( entry_plug->setTemp ){
276 <        if( !(i % vel_n) ) tStats->velocitize();
276 >        if( !(time % vel_n) ) tStats->velocitize();
277        }
278 <      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
279 <      if( !(i % status_n) ) e_out->writeStat( i * dt );
278 >      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
279 >      if( !(time % status_n) ) e_out->writeStat( time * dt );
280      }
281    }
282    
# Line 427 | Line 329 | void Verlet::move_b( double dt ){
329    const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
330    
331    double vx, vy, vz;
430  double v_sqr;
332    int mb;
333    double h_dt = 0.5 * dt;
334    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines