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

Comparing:
branches/mmeineke/OOPSE/libmdtools/Verlet.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Verlet.cpp (file contents), Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

# Line 7 | Line 7
7   #include "SimInfo.hpp"
8   #include "Thermo.hpp"
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 16 | 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 25 | 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    
37 < Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
37 > Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){
38    
39    // get what information we need from the SimInfo object
40  
41    entry_plug = &info;
42    myFF = the_ff;
43 <
44 <
43 >  myES = the_es;
44 >  
45    c_natoms = info.n_atoms;
46    c_atoms = info.atoms;
47 <  c_sr_interactions = info.sr_interactions;
48 <  c_n_SRI = info.n_SRI;
47 >  nMols = info.n_mol;
48 >  molecules = info.molecules;
49    c_is_constrained = 0;
50    c_box_x = info.box_x;
51    c_box_y = info.box_y;
# Line 67 | Line 72 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
72  
73    Constraint *temp_con;
74    Constraint *dummy_plug;
75 <  temp_con = new Constraint[c_n_SRI];
75 >  temp_con = new Constraint[info.n_SRI];
76  
77    c_n_constrained = 0;
78    int constrained = 0;
79 <  
80 <  for(int i = 0; i < c_n_SRI; i++){
79 >  SRI** theArray;
80 >  for(int i = 0; i < nMols; i++){
81      
82 <    constrained = c_sr_interactions[i]->is_constrained();
82 >    theArray = (SRI**) molecules[i].getMyBonds();
83 >    for(int j=0; j<molecules[i].getNBonds(); j++){
84 >      
85 >      constrained = theArray[j]->is_constrained();
86 >      
87 >      if(constrained){
88 >        
89 >        dummy_plug = theArray[j]->get_constraint();
90 >        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
91 >        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
92 >        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
93 >        
94 >        c_n_constrained++;
95 >        constrained = 0;
96 >      }
97 >    }
98  
99 <    if(constrained){
99 >    theArray = (SRI**) molecules[i].getMyBends();
100 >    for(int j=0; j<molecules[i].getNBends(); j++){
101        
102 <      dummy_plug = c_sr_interactions[i]->get_constraint();
103 <      temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
104 <      temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
105 <      temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
102 >      constrained = theArray[j]->is_constrained();
103 >      
104 >      if(constrained){
105 >        
106 >        dummy_plug = theArray[j]->get_constraint();
107 >        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
108 >        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
109 >        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
110 >        
111 >        c_n_constrained++;
112 >        constrained = 0;
113 >      }
114 >    }
115  
116 <      c_n_constrained++;
117 <      constrained = 0;
116 >    theArray = (SRI**) molecules[i].getMyTorsions();
117 >    for(int j=0; j<molecules[i].getNTorsions(); j++){
118 >      
119 >      constrained = theArray[j]->is_constrained();
120 >      
121 >      if(constrained){
122 >        
123 >        dummy_plug = theArray[j]->get_constraint();
124 >        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
125 >        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
126 >        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
127 >        
128 >        c_n_constrained++;
129 >        constrained = 0;
130 >      }
131      }
132 +
133 +
134    }
135  
136    if(c_n_constrained > 0){
# Line 94 | Line 139 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
139      c_constrained_i = new int[c_n_constrained];
140      c_constrained_j = new int[c_n_constrained];
141      c_constrained_dsqr = new double[c_n_constrained];
142 <
142 >    
143      for( int i = 0; i < c_n_constrained; i++){
144        
145        /* add 1 to the index for the fortran arrays. */
146 <
146 >      
147        c_constrained_i[i] = temp_con[i].get_a() + 1;
148        c_constrained_j[i] = temp_con[i].get_b() + 1;
149        c_constrained_dsqr[i] = temp_con[i].get_dsqr();
# Line 126 | Line 171 | void Verlet::integrate( void ){
171   void Verlet::integrate( void ){
172    
173    int i, j; /* loop counters */
174 <  int calcPot;
174 >  int calcPot, calcStress;
175  
176    double kE;
177  
# Line 144 | 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 155 | Line 202 | void Verlet::integrate( void ){
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  
209    StatWriter*  e_out    = new StatWriter( entry_plug );
210    DumpWriter*  dump_out = new DumpWriter( entry_plug );
211  
212    // the first time integrate is called, the forces need to be initialized
164 –
213    
214 <  myFF->doForces(1,0);
214 >  myFF->doForces(1,1);
215    
216    if( entry_plug->setTemp ){
217      tStats->velocitize();
# Line 175 | Line 223 | void Verlet::integrate( void ){
223  
224    calcPot = 0;
225  
226 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
227 +    calcStress = 1;
228 +  } else {
229 +    calcStress = 0;
230 +  }
231 +
232    if( c_is_constrained ){
233      for(i = 0; i < n_loops; i++){
234        
235 +      
236        // fill R, V, and F arrays and RATTLE in fortran
237 <
237 >      
238        for( j=0; j<c_natoms; j++ ){
239 <
239 >        
240          Rx[j] = c_atoms[j]->getX();
241          Ry[j] = c_atoms[j]->getY();
242          Rz[j] = c_atoms[j]->getZ();
# Line 196 | 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 213 | 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,0);
293 >      myFF->doForces(calcPot,calcStress);
294        
295        // finish the constrain move ( same as above. )
296  
# Line 234 | Line 309 | void Verlet::integrate( void ){
309          Fz[j] = c_atoms[j]->getFz();
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 251 | Line 344 | void Verlet::integrate( void ){
344          c_atoms[j]->set_vz(Vz[j]);
345        }
346        
347 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
348 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
349 +      
350 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {        
351 +        tStats->getPressureTensor(press);
352 +        myES->NoseHooverAndersonNPT( dt,
353 +                                     tStats->getKinetic(),
354 +                                     press);
355 +      }
356 +
357        time = i + 1;
358        
359        if( entry_plug->setTemp ){
360          if( !(time % vel_n) ) tStats->velocitize();
361        }
362        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
363 <      if( !((time+1) % status_n) ) calcPot = 1;
364 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
363 >
364 >      if( !((time+1) % status_n) ) {
365 >        calcPot = 1;
366 >        calcStress = 1;
367 >      }
368 >      if( !(time % status_n) ){
369 >        e_out->writeStat( time * dt );
370 >        calcPot = 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 <      
378 >
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        
387 <      myFF->doForces(calcPot,0);
387 >      myFF->doForces(calcPot,calcStress);
388              
389        // complete the verlet move
390        
391        move_b( dt );
392  
393 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
394 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
395 +
396 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
397 +        tStats->getPressureTensor(press);
398 +        myES->NoseHooverAndersonNPT( dt,
399 +                                     tStats->getKinetic(),
400 +                                     press);
401 +      }
402 +
403        time = i + 1;
404        
405        if( entry_plug->setTemp ){
406          if( !(time % vel_n) ) tStats->velocitize();
407        }
408        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
409 <      if( !((time+1) % status_n) ) calcPot = 1;
410 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
409 >      if( !((time+1) % status_n) ) {
410 >        calcPot = 1;
411 >        calcStress = 1;
412 >      }
413 >      if( !(time % status_n) ){
414 >        e_out->writeStat( time * dt );
415 >        calcPot = 0;
416 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
417 >        else calcStress = 0;
418 >      }
419      }
420    }
421    

Diff Legend

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