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 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 478 by gezelter, Tue Apr 8 14:39:40 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 | extern "C"{
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 <
40 <
39 >  myES = the_es;
40 >  
41    c_natoms = info.n_atoms;
42    c_atoms = info.atoms;
43 <  c_sr_interactions = info.sr_interactions;
44 <  c_n_SRI = info.n_SRI;
43 >  nMols = info.n_mol;
44 >  molecules = info.molecules;
45    c_is_constrained = 0;
46    c_box_x = info.box_x;
47    c_box_y = info.box_y;
# Line 67 | Line 68 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
68  
69    Constraint *temp_con;
70    Constraint *dummy_plug;
71 <  temp_con = new Constraint[c_n_SRI];
71 >  temp_con = new Constraint[info.n_SRI];
72  
73    c_n_constrained = 0;
74    int constrained = 0;
75 <  
76 <  for(int i = 0; i < c_n_SRI; i++){
75 >  SRI** theArray;
76 >  for(int i = 0; i < nMols; i++){
77      
78 <    constrained = c_sr_interactions[i]->is_constrained();
78 >    theArray = (SRI**) molecules[i].getMyBonds();
79 >    for(int j=0; j<molecules[i].getNBonds(); j++){
80 >      
81 >      constrained = theArray[j]->is_constrained();
82 >      
83 >      if(constrained){
84 >        
85 >        dummy_plug = theArray[j]->get_constraint();
86 >        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
87 >        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
88 >        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
89 >        
90 >        c_n_constrained++;
91 >        constrained = 0;
92 >      }
93 >    }
94  
95 <    if(constrained){
95 >    theArray = (SRI**) molecules[i].getMyBends();
96 >    for(int j=0; j<molecules[i].getNBends(); j++){
97        
98 <      dummy_plug = c_sr_interactions[i]->get_constraint();
99 <      temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
100 <      temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
101 <      temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
98 >      constrained = theArray[j]->is_constrained();
99 >      
100 >      if(constrained){
101 >        
102 >        dummy_plug = theArray[j]->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() );
106 >        
107 >        c_n_constrained++;
108 >        constrained = 0;
109 >      }
110 >    }
111  
112 <      c_n_constrained++;
113 <      constrained = 0;
112 >    theArray = (SRI**) molecules[i].getMyTorsions();
113 >    for(int j=0; j<molecules[i].getNTorsions(); j++){
114 >      
115 >      constrained = theArray[j]->is_constrained();
116 >      
117 >      if(constrained){
118 >        
119 >        dummy_plug = theArray[j]->get_constraint();
120 >        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
121 >        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
122 >        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
123 >        
124 >        c_n_constrained++;
125 >        constrained = 0;
126 >      }
127      }
128 +
129 +
130    }
131  
132    if(c_n_constrained > 0){
# Line 94 | Line 135 | Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){
135      c_constrained_i = new int[c_n_constrained];
136      c_constrained_j = new int[c_n_constrained];
137      c_constrained_dsqr = new double[c_n_constrained];
138 <
138 >    
139      for( int i = 0; i < c_n_constrained; i++){
140        
141        /* add 1 to the index for the fortran arrays. */
142 <
142 >      
143        c_constrained_i[i] = temp_con[i].get_a() + 1;
144        c_constrained_j[i] = temp_con[i].get_b() + 1;
145        c_constrained_dsqr[i] = temp_con[i].get_dsqr();
# Line 126 | 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 161 | Line 202 | void Verlet::integrate( void ){
202    DumpWriter*  dump_out = new DumpWriter( entry_plug );
203  
204    // the first time integrate is called, the forces need to be initialized
164 –
205    
206 <  myFF->doForces(1,0);
206 >  myFF->doForces(1,1);
207    
208    if( entry_plug->setTemp ){
209      tStats->velocitize();
# Line 175 | Line 215 | void Verlet::integrate( void ){
215  
216    calcPot = 0;
217  
218 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
219 +    calcStress = 1;
220 +  } else {
221 +    calcStress = 0;
222 +  }
223 +
224    if( c_is_constrained ){
225      for(i = 0; i < n_loops; i++){
226        
227 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
228 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
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 215 | Line 264 | void Verlet::integrate( void ){
264  
265        // calculate the forces
266        
267 <      myFF->doForces(calcPot,0);
267 >      myFF->doForces(calcPot,calcStress);
268        
269        // finish the constrain move ( same as above. )
270  
# Line 234 | Line 283 | void Verlet::integrate( void ){
283          Fz[j] = c_atoms[j]->getFz();
284        }
285          
286 +
287        v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
288                        Fx, Fy, Fz,
289                        kE, c_n_constrained, c_constrained_dsqr,
# Line 251 | Line 301 | void Verlet::integrate( void ){
301          c_atoms[j]->set_vz(Vz[j]);
302        }
303        
304 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
305 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
306 +      
307 +      if (!strcasecmp( entry_plug->ensemble, "NPT") )
308 +        myES->NoseHooverAndersonNPT( dt,
309 +                                     tStats->getKinetic(),
310 +                                     tStats->getPressure());
311 +
312        time = i + 1;
313        
314        if( entry_plug->setTemp ){
315          if( !(time % vel_n) ) tStats->velocitize();
316        }
317        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
318 <      if( !((time+1) % status_n) ) calcPot = 1;
319 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
318 >      if( !((time+1) % status_n) ) {
319 >        calcPot = 1;
320 >        // bitwise masking in case we need it for NPT
321 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
322 >      }
323 >      if( !(time % status_n) ){
324 >        e_out->writeStat( time * dt );
325 >        calcPot = 0;
326 >        // bitwise masking in case we need it for NPT
327 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
328 >      }
329      }
330    }
331    else{
332      for(i = 0; i < n_loops; i++){
333 <      
333 >
334 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
335 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
336 >    
337        move_a( dt );
338        
339        // calculate the forces
340        
341 <      myFF->doForces(calcPot,0);
341 >      myFF->doForces(calcPot,calcStress);
342              
343        // complete the verlet move
344        
345        move_b( dt );
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 +        myES->NoseHooverAndersonNPT( dt,
352 +                                     tStats->getKinetic(),
353 +                                     tStats->getPressure());
354 +
355        time = i + 1;
356        
357        if( entry_plug->setTemp ){
358          if( !(time % vel_n) ) tStats->velocitize();
359        }
360        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
361 <      if( !((time+1) % status_n) ) calcPot = 1;
362 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
361 >      if( !((time+1) % status_n) ) {
362 >        calcPot = 1;
363 >        // bitwise masking in case we need it for NPT
364 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
365 >      }
366 >      if( !(time % status_n) ){
367 >        e_out->writeStat( time * dt );
368 >        calcPot = 0;
369 >        // bitwise masking in case we need it for NPT
370 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
371 >      }
372      }
373    }
374    

Diff Legend

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