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

Comparing trunk/OOPSE/libmdtools/Verlet.cpp (file contents):
Revision 423 by mmeineke, Thu Mar 27 20:12:15 2003 UTC vs.
Revision 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC

# Line 7 | Line 7
7   #include "SimInfo.hpp"
8   #include "Thermo.hpp"
9   #include "ReadWrite.hpp"
10 + #include "ExtendedSystem.hpp"
11  
12   extern "C"{
13    
# Line 29 | Line 30 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
30   }
31  
32    
33 < Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
33 > Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){
34    
35    // get what information we need from the SimInfo object
36  
37    entry_plug = &info;
38    myFF = the_ff;
39 <
39 >  myES = the_es;
40    
41    c_natoms = info.n_atoms;
42    c_atoms = info.atoms;
# Line 74 | Line 75 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
75    SRI** theArray;
76    for(int i = 0; i < nMols; i++){
77      
78 <    theArray = molecules[i].getMyBonds();
79 <    for(int j=0; j<molecules[i].getNbonds(); j++){
78 >    theArray = (SRI**) molecules[i].getMyBonds();
79 >    for(int j=0; j<molecules[i].getNBonds(); j++){
80        
81        constrained = theArray[j]->is_constrained();
82        
# Line 91 | Line 92 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
92        }
93      }
94  
95 <    theArray = molecules[i].getMyBends();
96 <    for(int j=0; j<molecules[i].getNbends(); j++){
95 >    theArray = (SRI**) molecules[i].getMyBends();
96 >    for(int j=0; j<molecules[i].getNBends(); j++){
97        
98        constrained = theArray[j]->is_constrained();
99        
# Line 108 | Line 109 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
109        }
110      }
111  
112 <    theArray = molecules[i].getMyTorsions();
113 <    for(int j=0; j<molecules[i].getNtorsions(); j++){
112 >    theArray = (SRI**) molecules[i].getMyTorsions();
113 >    for(int j=0; j<molecules[i].getNTorsions(); j++){
114        
115        constrained = theArray[j]->is_constrained();
116        
# Line 166 | Line 167 | void Verlet::integrate( void ){
167   void Verlet::integrate( void ){
168    
169    int i, j; /* loop counters */
170 <  int calcPot;
170 >  int calcPot, calcStress;
171  
172    double kE;
173  
# Line 184 | Line 185 | void Verlet::integrate( void ){
185    
186    int time;
187  
188 +  double press[9];
189 +
190    double dt = entry_plug->dt;
191    double runTime = entry_plug->run_time;
192    double sampleTime = entry_plug->sampleTime;
# Line 201 | Line 204 | void Verlet::integrate( void ){
204    DumpWriter*  dump_out = new DumpWriter( entry_plug );
205  
206    // the first time integrate is called, the forces need to be initialized
204 –
207    
208 <  myFF->doForces(1,0);
208 >  myFF->doForces(1,1);
209    
210    if( entry_plug->setTemp ){
211      tStats->velocitize();
# Line 214 | Line 216 | void Verlet::integrate( void ){
216    e_out->writeStat( 0.0 );
217  
218    calcPot = 0;
219 +
220 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
221 +    calcStress = 1;
222 +  } else {
223 +    calcStress = 0;
224 +  }
225  
226    if( c_is_constrained ){
227      for(i = 0; i < n_loops; i++){
228        
229 +      
230        // fill R, V, and F arrays and RATTLE in fortran
231 <
231 >      
232        for( j=0; j<c_natoms; j++ ){
233 <
233 >        
234          Rx[j] = c_atoms[j]->getX();
235          Ry[j] = c_atoms[j]->getY();
236          Rz[j] = c_atoms[j]->getZ();
# Line 253 | Line 262 | void Verlet::integrate( void ){
262          c_atoms[j]->set_vz(Vz[j]);
263        }
264  
265 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
266 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
267 +  
268        // calculate the forces
269        
270 <      myFF->doForces(calcPot,0);
270 >      myFF->doForces(calcPot,calcStress);
271        
272        // finish the constrain move ( same as above. )
273  
# Line 274 | Line 286 | void Verlet::integrate( void ){
286          Fz[j] = c_atoms[j]->getFz();
287        }
288          
289 +
290        v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
291                        Fx, Fy, Fz,
292                        kE, c_n_constrained, c_constrained_dsqr,
# Line 291 | Line 304 | void Verlet::integrate( void ){
304          c_atoms[j]->set_vz(Vz[j]);
305        }
306        
307 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
308 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
309 +      
310 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {        
311 +        tStats->getPressureTensor(press);
312 +        myES->NoseHooverAndersonNPT( dt,
313 +                                     tStats->getKinetic(),
314 +                                     press);
315 +      }
316 +
317        time = i + 1;
318        
319        if( entry_plug->setTemp ){
320          if( !(time % vel_n) ) tStats->velocitize();
321        }
322        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
323 <      if( !((time+1) % status_n) ) calcPot = 1;
324 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
323 >
324 >      if( !((time+1) % status_n) ) {
325 >        calcPot = 1;
326 >        calcStress = 1;
327 >      }
328 >      if( !(time % status_n) ){
329 >        e_out->writeStat( time * dt );
330 >        calcPot = 0;
331 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
332 >        else calcStress = 0;
333 >      }
334      }
335    }
336    else{
337      for(i = 0; i < n_loops; i++){
338 <      
338 >
339 >    
340        move_a( dt );
341 +
342 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
343 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
344        
345        // calculate the forces
346        
347 <      myFF->doForces(calcPot,0);
347 >      myFF->doForces(calcPot,calcStress);
348              
349        // complete the verlet move
350        
351        move_b( dt );
352  
353 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
354 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
355 +
356 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
357 +        tStats->getPressureTensor(press);
358 +        myES->NoseHooverAndersonNPT( dt,
359 +                                     tStats->getKinetic(),
360 +                                     press);
361 +      }
362 +
363        time = i + 1;
364        
365        if( entry_plug->setTemp ){
366          if( !(time % vel_n) ) tStats->velocitize();
367        }
368        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
369 <      if( !((time+1) % status_n) ) calcPot = 1;
370 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
369 >      if( !((time+1) % status_n) ) {
370 >        calcPot = 1;
371 >        calcStress = 1;
372 >      }
373 >      if( !(time % status_n) ){
374 >        e_out->writeStat( time * dt );
375 >        calcPot = 0;
376 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
377 >        else calcStress = 0;
378 >      }
379      }
380    }
381    

Diff Legend

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