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 478 by gezelter, Tue Apr 8 14:39:40 2003 UTC vs.
Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

# Line 9 | Line 9
9   #include "ReadWrite.hpp"
10   #include "ExtendedSystem.hpp"
11  
12 + #include "simError.h"
13 +
14   extern "C"{
15    
16    void v_constrain_a_( double &dt, int &n_atoms, double* mass,
# Line 17 | Line 19 | extern "C"{
19                         double* Fx, double* Fy, double* Fz,
20                         int &n_constrained, double *constr_sqr,
21                         int* constr_i, int* constr_j,
22 <                       double &box_x, double &box_y, double &box_z );
22 >                       double &box_x, double &box_y, double &box_z,
23 >                       int &isError );
24  
25    void v_constrain_b_( double &dt, int &n_atoms, double* mass,
26                         double* Rx, double* Ry, double* Rz,
# Line 26 | Line 29 | extern "C"{
29                         double &Kinetic,
30                         int &n_constrained, double *constr_sqr,
31                         int* constr_i, int* constr_j,
32 <                       double &box_x, double &box_y, double &box_z );
32 >                       double &box_x, double &box_y, double &box_z,
33 >                       int &isError);
34   }
35  
36    
# Line 185 | Line 189 | void Verlet::integrate( void ){
189    
190    int time;
191  
192 +  double press[9];
193 +
194    double dt = entry_plug->dt;
195    double runTime = entry_plug->run_time;
196    double sampleTime = entry_plug->sampleTime;
# Line 195 | Line 201 | void Verlet::integrate( void ){
201    int sample_n = (int)( sampleTime / dt );
202    int status_n = (int)( statusTime / dt );
203    int vel_n    = (int)( thermalTime / dt );
204 +
205 +  int isError;
206  
207    Thermo *tStats = new Thermo( entry_plug );
208  
# Line 224 | Line 232 | void Verlet::integrate( void ){
232    if( c_is_constrained ){
233      for(i = 0; i < n_loops; i++){
234        
227 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
228 –        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
235        
236        // fill R, V, and F arrays and RATTLE in fortran
237        
# Line 245 | Line 251 | void Verlet::integrate( void ){
251  
252        }
253        
254 +      isError = 0;
255        v_constrain_a_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
256                        Fx, Fy, Fz,
257                        c_n_constrained, c_constrained_dsqr,
258                        c_constrained_i, c_constrained_j,
259 <                      c_box_x, c_box_y, c_box_z );
259 >                      c_box_x, c_box_y, c_box_z,
260 >                      isError);
261 >
262 >      if( isError ){
263 >        
264 >        sprintf( painCave.errMsg,
265 >                 "Constraint Failure in Fortran move A\n" );
266 >        painCave.isFatal = 1;
267 >        simError();
268 >      }
269 >      
270 > #ifdef IS_MPI
271 >      sprintf( checkPointMsg,
272 >               "succesful return from move a.\n" );
273 >      MPIcheckPoint();
274  
275 + #endif //is_mpi
276 +
277        for( j=0; j<c_natoms; j++ ){
278  
279          c_atoms[j]->setX(Rx[j]);
# Line 262 | Line 285 | void Verlet::integrate( void ){
285          c_atoms[j]->set_vz(Vz[j]);
286        }
287  
288 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
289 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
290 +  
291        // calculate the forces
292        
293        myFF->doForces(calcPot,calcStress);
# Line 284 | Line 310 | void Verlet::integrate( void ){
310        }
311          
312  
313 +      isError = 0;
314        v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
315                        Fx, Fy, Fz,
316                        kE, c_n_constrained, c_constrained_dsqr,
317                        c_constrained_i, c_constrained_j,
318 <                      c_box_x, c_box_y, c_box_z );
319 <      
318 >                      c_box_x, c_box_y, c_box_z, isError );
319 >
320 >      if( isError ){
321 >        
322 >        sprintf( painCave.errMsg,
323 >                 "Constraint Failure in Fortran move B\n" );
324 >        painCave.isFatal = 1;
325 >        simError();
326 >      }
327 >
328 > #ifdef IS_MPI
329 >      sprintf( checkPointMsg,
330 >               "succesful return from move B.\n" );
331 >      MPIcheckPoint();
332 >
333 > #endif //is_mpi      
334 >
335 >
336        for( j=0; j<c_natoms; j++ ){
337  
338          c_atoms[j]->setX(Rx[j]);
# Line 304 | Line 347 | void Verlet::integrate( void ){
347        if (!strcasecmp( entry_plug->ensemble, "NVT"))
348          myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
349        
350 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
350 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {        
351 >        tStats->getPressureTensor(press);
352          myES->NoseHooverAndersonNPT( dt,
353                                       tStats->getKinetic(),
354 <                                     tStats->getPressure());
354 >                                     press);
355 >      }
356  
357        time = i + 1;
358        
# Line 315 | Line 360 | void Verlet::integrate( void ){
360          if( !(time % vel_n) ) tStats->velocitize();
361        }
362        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
363 +
364        if( !((time+1) % status_n) ) {
365          calcPot = 1;
366 <        // bitwise masking in case we need it for NPT
321 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
366 >        calcStress = 1;
367        }
368        if( !(time % status_n) ){
369          e_out->writeStat( time * dt );
370          calcPot = 0;
371 <        // bitwise masking in case we need it for NPT
372 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
371 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
372 >        else calcStress = 0;
373        }
374      }
375    }
376    else{
377      for(i = 0; i < n_loops; i++){
378  
334 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
335 –        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
379      
380        move_a( dt );
381 +
382 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
383 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
384        
385        // calculate the forces
386        
# Line 347 | Line 393 | void Verlet::integrate( void ){
393        if (!strcasecmp( entry_plug->ensemble, "NVT"))
394          myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
395  
396 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
396 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
397 >        tStats->getPressureTensor(press);
398          myES->NoseHooverAndersonNPT( dt,
399                                       tStats->getKinetic(),
400 <                                     tStats->getPressure());
400 >                                     press);
401 >      }
402  
403        time = i + 1;
404        
# Line 360 | Line 408 | void Verlet::integrate( void ){
408        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
409        if( !((time+1) % status_n) ) {
410          calcPot = 1;
411 <        // bitwise masking in case we need it for NPT
364 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
411 >        calcStress = 1;
412        }
413        if( !(time % status_n) ){
414          e_out->writeStat( time * dt );
415          calcPot = 0;
416 <        // bitwise masking in case we need it for NPT
417 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
416 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
417 >        else calcStress = 0;
418        }
419      }
420    }

Diff Legend

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