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

Comparing trunk/mdtools/md_code/Verlet.cpp (file contents):
Revision 11 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
Revision 253 by chuckv, Thu Jan 30 15:20:21 2003 UTC

# Line 30 | Line 30 | extern "C"{
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 120 | Line 121 | Verlet::~Verlet(){
121    
122    delete[] c_mass;
123    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 –  }
124   }
125  
126  
127   void Verlet::integrate( void ){
128    
129    int i, j; /* loop counters */
130 <  
130 >  int calcPot;
131 >
132    double kE;
133  
134    double *Rx = new double[c_natoms];
# Line 224 | Line 143 | void Verlet::integrate( void ){
143    double *Fy = new double[c_natoms];
144    double *Fz = new double[c_natoms];
145    
146 +  int time;
147 +
148    double dt = entry_plug->dt;
149    double runTime = entry_plug->run_time;
150    double sampleTime = entry_plug->sampleTime;
# Line 243 | Line 164 | void Verlet::integrate( void ){
164    // the first time integrate is called, the forces need to be initialized
165  
166    
167 <  for(i = 0; i < c_natoms; i++){
247 <    c_atoms[i]->zeroForces();
248 <  }
167 >  myFF->doForces(1);
168    
250 –  for(i = 0; i < c_n_SRI; i++){
251 –    c_sr_interactions[i]->calc_forces();
252 –  }
253 –  
254 –  longRange->calc_forces();
255 –  
169    if( entry_plug->setTemp ){
170      tStats->velocitize();
171    }
172    
173 +  dump_out->writeDump( 0.0 );
174 +
175 +  e_out->writeStat( 0.0 );
176 +  calcPot = 0;
177 +
178    if( c_is_constrained ){
179      for(i = 0; i < n_loops; i++){
180        
# Line 297 | Line 215 | void Verlet::integrate( void ){
215  
216        // calculate the forces
217        
218 <      for(j = 0; j < c_natoms; j++){
301 <        c_atoms[j]->zeroForces();
302 <      }
218 >      myFF->doForces(calcPot);
219        
304 –      for(j = 0; j < c_n_SRI; j++){
305 –        c_sr_interactions[j]->calc_forces();
306 –      }
307 –      
308 –      longRange->calc_forces();
309 –      
220        // finish the constrain move ( same as above. )
221  
222        for( j=0; j<c_natoms; j++ ){
# Line 341 | Line 251 | void Verlet::integrate( void ){
251          c_atoms[j]->set_vz(Vz[j]);
252        }
253        
254 +      time = i + 1;
255 +      
256        if( entry_plug->setTemp ){
257 <        if( !(i % vel_n) ) tStats->velocitize();
257 >        if( !(time % vel_n) ) tStats->velocitize();
258        }
259 <      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
260 <      if( !(i % status_n) ) e_out->writeStat( i * dt );
261 <
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{
# Line 356 | Line 268 | void Verlet::integrate( void ){
268        
269        // calculate the forces
270        
271 <      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();
271 >      myFF->doForces(calcPot);
272              
273        // complete the verlet move
274        
275        move_b( dt );
276  
277 +      time = i + 1;
278 +      
279        if( entry_plug->setTemp ){
280 <        if( !(i % vel_n) ) tStats->velocitize();
280 >        if( !(time % vel_n) ) tStats->velocitize();
281        }
282 <      if( !(i % sample_n) ) dump_out->writeDump( i * dt );
283 <      if( !(i % status_n) ) e_out->writeStat( i * dt );
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    
# Line 427 | Line 334 | void Verlet::move_b( double dt ){
334    const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
335    
336    double vx, vy, vz;
430 –  double v_sqr;
337    int mb;
338    double h_dt = 0.5 * dt;
339    

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines