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

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 482 by chuckv, Tue Apr 8 22:38:43 2003 UTC

# Line 5 | Line 5
5   #include "Thermo.hpp"
6   #include "ReadWrite.hpp"
7   #include "ForceFields.hpp"
8 + #include "ExtendedSystem.hpp"
9   #include "simError.h"
10  
11   extern "C"{
# Line 30 | Line 31 | extern "C"{
31  
32  
33  
34 < Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff ){
34 > Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff,
35 >                        ExtendedSystem* the_es ){
36    entry_plug = the_entry_plug;
37    myFF = the_ff;
38 +  myES = the_es;
39    isFirst = 1;
40  
41 +  std::cerr<< "calling symplectic constructor\n";
42  
43 <  srInteractions = entry_plug->sr_interactions;
44 <  nSRI           =           entry_plug->n_SRI;
43 >  molecules = entry_plug->molecules;
44 >  nMols = entry_plug->n_mol;
45  
46    // give a little love back to the SimInfo object
47    
# Line 49 | Line 53 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
53    mass = new double[entry_plug->n_atoms];
54    for(int i = 0; i < entry_plug->n_atoms; i++){
55      mass[i] = entry_plug->atoms[i]->getMass();
56 <  }
53 <
56 >  }  
57  
58    // check for constraints
59  
# Line 58 | Line 61 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
61  
62    Constraint *temp_con;
63    Constraint *dummy_plug;
64 <  temp_con = new Constraint[nSRI];
64 >  temp_con = new Constraint[entry_plug->n_SRI];
65    n_constrained = 0;
66    int constrained = 0;
67    
68 <  for(int i = 0; i < nSRI; i++){
68 >  SRI** theArray;
69 >  for(int i = 0; i < nMols; i++){
70      
71 <    constrained = srInteractions[i]->is_constrained();
72 <    
69 <    if(constrained){
71 >    theArray = (SRI**) molecules[i].getMyBonds();
72 >    for(int j=0; j<molecules[i].getNBonds(); j++){
73        
74 <      dummy_plug = srInteractions[i]->get_constraint();
75 <      temp_con[n_constrained].set_a( dummy_plug->get_a() );
76 <      temp_con[n_constrained].set_b( dummy_plug->get_b() );
77 <      temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
74 >      constrained = theArray[j]->is_constrained();
75 >      
76 >      if(constrained){
77 >        
78 >        dummy_plug = theArray[j]->get_constraint();
79 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
80 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
81 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
82 >        
83 >        n_constrained++;
84 >        constrained = 0;
85 >      }
86 >    }
87  
88 <      n_constrained++;
89 <      constrained = 0;
88 >    theArray = (SRI**) molecules[i].getMyBends();
89 >    for(int j=0; j<molecules[i].getNBends(); j++){
90 >      
91 >      constrained = theArray[j]->is_constrained();
92 >      
93 >      if(constrained){
94 >        
95 >        dummy_plug = theArray[j]->get_constraint();
96 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
97 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
98 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
99 >        
100 >        n_constrained++;
101 >        constrained = 0;
102 >      }
103      }
104 +
105 +    theArray = (SRI**) molecules[i].getMyTorsions();
106 +    for(int j=0; j<molecules[i].getNTorsions(); j++){
107 +      
108 +      constrained = theArray[j]->is_constrained();
109 +      
110 +      if(constrained){
111 +        
112 +        dummy_plug = theArray[j]->get_constraint();
113 +        temp_con[n_constrained].set_a( dummy_plug->get_a() );
114 +        temp_con[n_constrained].set_b( dummy_plug->get_b() );
115 +        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
116 +        
117 +        n_constrained++;
118 +        constrained = 0;
119 +      }
120 +    }
121    }
122  
123    if(n_constrained > 0){
# Line 122 | Line 164 | void Symplectic::integrate( void ){
164    double dt2;                       // half the dt
165  
166    double vx, vy, vz;    // the velocities
167 < //  double vx2, vy2, vz2; // the square of the velocities
167 >  double vx2, vy2, vz2; // the square of the velocities
168    double rx, ry, rz;    // the postitions
169    
170    double ji[3];   // the body frame angular momentum
# Line 144 | Line 186 | void Symplectic::integrate( void ){
186    int status_n = (int)( statusTime / dt );
187    int vel_n    = (int)( thermalTime / dt );
188  
189 <  int calcPot;
189 >  int calcPot, calcStress;
190  
191 <   Thermo *tStats = new Thermo( entry_plug );
191 >  Thermo *tStats;
192 >  StatWriter*  e_out;
193 >  DumpWriter*  dump_out;
194  
195 <  StatWriter*  e_out    = new StatWriter( entry_plug );
152 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
195 >  std::cerr << "about to call new thermo\n";
196  
197 +  tStats   = new Thermo( entry_plug );
198 +  e_out    = new StatWriter( entry_plug );
199 +
200 +  std::cerr << "calling dumpWriter \n";
201 +  dump_out = new DumpWriter( entry_plug );
202 +  std::cerr << "called dumpWriter \n";
203 +
204    Atom** atoms = entry_plug->atoms;
205    DirectionalAtom* dAtom;
206    dt2 = 0.5 * dt;
207  
208    // initialize the forces the before the first step
209  
210 <  myFF->doForces(1,0);
210 >  myFF->doForces(1,1);
211    
212    if( entry_plug->setTemp ){
213      
# Line 169 | Line 219 | void Symplectic::integrate( void ){
219    
220    calcPot = 0;
221  
222 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
223 +    calcStress = 1;
224 +  } else {
225 +    calcStress = 0;
226 +  }
227 +
228    if( n_constrained ){
229  
230      double *Rx = new double[nAtoms];
# Line 185 | Line 241 | void Symplectic::integrate( void ){
241      
242  
243      for( tl=0; tl < n_loops; tl++ ){
244 +
245 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
246 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
247        
248        for( j=0; j<nAtoms; j++ ){
249  
# Line 285 | Line 344 | void Symplectic::integrate( void ){
344        
345        // calculate the forces
346        
347 <      myFF->doForces(calcPot, 0);
347 >      myFF->doForces(calcPot, calcStress);
348        
349        // move b
350  
# Line 350 | Line 409 | void Symplectic::integrate( void ){
409          }
410        }
411      
412 +
413 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
414 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
415 +
416 +      if (!strcasecmp( entry_plug->ensemble, "NPT") )
417 +        myES->NoseHooverAndersonNPT( dt,
418 +                                     tStats->getKinetic(),
419 +                                     tStats->getPressure());
420 +
421        time = tl + 1;
422        
423        if( entry_plug->setTemp ){
424          if( !(time % vel_n) ) tStats->velocitize();
425        }
426        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
427 <      if( !((time+1) % status_n) ) calcPot = 1;
428 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
427 >      if( !((time+1) % status_n) ) {
428 >        calcPot = 1;
429 >        // bitwise masking in case we need it for NPT
430 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
431 >      }
432 >      if( !(time % status_n) ){
433 >        e_out->writeStat( time * dt );
434 >        calcPot = 0;
435 >        // bitwise masking in case we need it for NPT
436 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
437 >      }
438      }
439    }
440    else{
# Line 367 | Line 444 | void Symplectic::integrate( void ){
444        kE = 0.0;
445        rot_kE= 0.0;
446        trans_kE = 0.0;
447 +
448 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
449 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
450        
451        for( i=0; i<nAtoms; i++ ){
452          
# Line 455 | Line 535 | void Symplectic::integrate( void ){
535        
536        // calculate the forces
537        
538 <      myFF->doForces(calcPot,0);
538 >      myFF->doForces(calcPot,calcStress);
539        
540        for( i=0; i< nAtoms; i++ ){
541          
# Line 472 | Line 552 | void Symplectic::integrate( void ){
552          atoms[i]->set_vy( vy );
553          atoms[i]->set_vz( vz );
554          
555 < //      vx2 = vx * vx;
556 < //      vy2 = vy * vy;
557 < //      vz2 = vz * vz;
555 >        vx2 = vx * vx;
556 >        vy2 = vy * vy;
557 >        vz2 = vz * vz;
558          
559          if( atoms[i]->isDirectional() ){
560  
# Line 505 | Line 585 | void Symplectic::integrate( void ){
585            dAtom->setJx( ji[0] );
586            dAtom->setJy( ji[1] );
587            dAtom->setJz( ji[2] );
588 <        }
588 >        }
589        }
590 <      
590 >    
591 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
592 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
593 >
594 >      if (!strcasecmp( entry_plug->ensemble, "NPT") )
595 >        myES->NoseHooverAndersonNPT( dt,
596 >                                     tStats->getKinetic(),
597 >                                     tStats->getPressure());
598 >  
599        time = tl + 1;
600        
601        if( entry_plug->setTemp ){
602          if( !(time % vel_n) ) tStats->velocitize();
603        }
604        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605 <      if( !((time+1) % status_n) ) calcPot = 1;
606 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
605 >      if( !((time+1) % status_n) ) {
606 >        calcPot = 1;
607 >        // bitwise masking in case we need it for NPT
608 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
609 >      }
610 >      if( !(time % status_n) ){
611 >        e_out->writeStat( time * dt );
612 >        calcPot = 0;
613 >        // bitwise masking in case we need it for NPT
614 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
615 >      }      
616      }
617    }
618  
# Line 542 | Line 639 | void Symplectic::rotate( int axes1, int axes2, double
639  
640    for(i=0; i<3; i++){
641      for(j=0; j<3; j++){
642 <      tempA[i][j] = A[i][j];
642 >      tempA[j][i] = A[i][j];
643      }
644    }
645  
# Line 601 | Line 698 | void Symplectic::rotate( int axes1, int axes2, double
698      for(j=0; j<3; j++){
699        A[j][i] = 0.0;
700        for(k=0; k<3; k++){
701 <        A[j][i] += tempA[k][i] * rot[j][k];
701 >        A[j][i] += tempA[i][k] * rot[j][k];
702        }
703      }
704    }

Diff Legend

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