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 486 by mmeineke, Thu Apr 10 16:22:00 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 +  molecules = entry_plug->molecules;
42 +  nMols = entry_plug->n_mol;
43  
39 –  srInteractions = entry_plug->sr_interactions;
40 –  nSRI           =           entry_plug->n_SRI;
41 –
44    // give a little love back to the SimInfo object
45    
46    if( entry_plug->the_integrator != NULL ) delete entry_plug->the_integrator;
# Line 49 | Line 51 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
51    mass = new double[entry_plug->n_atoms];
52    for(int i = 0; i < entry_plug->n_atoms; i++){
53      mass[i] = entry_plug->atoms[i]->getMass();
54 <  }
53 <
54 >  }  
55  
56    // check for constraints
57  
# Line 58 | Line 59 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
59  
60    Constraint *temp_con;
61    Constraint *dummy_plug;
62 <  temp_con = new Constraint[nSRI];
62 >  temp_con = new Constraint[entry_plug->n_SRI];
63    n_constrained = 0;
64    int constrained = 0;
65    
66 <  for(int i = 0; i < nSRI; i++){
66 >  SRI** theArray;
67 >  for(int i = 0; i < nMols; i++){
68      
69 <    constrained = srInteractions[i]->is_constrained();
70 <    
69 <    if(constrained){
69 >    theArray = (SRI**) molecules[i].getMyBonds();
70 >    for(int j=0; j<molecules[i].getNBonds(); j++){
71        
72 <      dummy_plug = srInteractions[i]->get_constraint();
73 <      temp_con[n_constrained].set_a( dummy_plug->get_a() );
74 <      temp_con[n_constrained].set_b( dummy_plug->get_b() );
75 <      temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
72 >      constrained = theArray[j]->is_constrained();
73 >      
74 >      if(constrained){
75 >        
76 >        dummy_plug = theArray[j]->get_constraint();
77 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
78 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
79 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
80 >        
81 >        n_constrained++;
82 >        constrained = 0;
83 >      }
84 >    }
85  
86 <      n_constrained++;
87 <      constrained = 0;
86 >    theArray = (SRI**) molecules[i].getMyBends();
87 >    for(int j=0; j<molecules[i].getNBends(); j++){
88 >      
89 >      constrained = theArray[j]->is_constrained();
90 >      
91 >      if(constrained){
92 >        
93 >        dummy_plug = theArray[j]->get_constraint();
94 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
95 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
96 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
97 >        
98 >        n_constrained++;
99 >        constrained = 0;
100 >      }
101      }
102 +
103 +    theArray = (SRI**) molecules[i].getMyTorsions();
104 +    for(int j=0; j<molecules[i].getNTorsions(); j++){
105 +      
106 +      constrained = theArray[j]->is_constrained();
107 +      
108 +      if(constrained){
109 +        
110 +        dummy_plug = theArray[j]->get_constraint();
111 +        temp_con[n_constrained].set_a( dummy_plug->get_a() );
112 +        temp_con[n_constrained].set_b( dummy_plug->get_b() );
113 +        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
114 +        
115 +        n_constrained++;
116 +        constrained = 0;
117 +      }
118 +    }
119    }
120  
121    if(n_constrained > 0){
# Line 122 | Line 162 | void Symplectic::integrate( void ){
162    double dt2;                       // half the dt
163  
164    double vx, vy, vz;    // the velocities
165 < //  double vx2, vy2, vz2; // the square of the velocities
165 >  double vx2, vy2, vz2; // the square of the velocities
166    double rx, ry, rz;    // the postitions
167    
168    double ji[3];   // the body frame angular momentum
# Line 130 | Line 170 | void Symplectic::integrate( void ){
170    double Tb[3];   // torque in the body frame
171    double angle;   // the angle through which to rotate the rotation matrix
172    double A[3][3]; // the rotation matrix
173 +  double press[9];
174  
175    int time;
176  
# Line 144 | Line 185 | void Symplectic::integrate( void ){
185    int status_n = (int)( statusTime / dt );
186    int vel_n    = (int)( thermalTime / dt );
187  
188 <  int calcPot;
188 >  int calcPot, calcStress;
189  
190 <   Thermo *tStats = new Thermo( entry_plug );
190 >  Thermo *tStats;
191 >  StatWriter*  e_out;
192 >  DumpWriter*  dump_out;
193  
194 <  StatWriter*  e_out    = new StatWriter( entry_plug );
195 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
194 >  tStats   = new Thermo( entry_plug );
195 >  e_out    = new StatWriter( entry_plug );
196 >  dump_out = new DumpWriter( entry_plug );
197  
198    Atom** atoms = entry_plug->atoms;
199    DirectionalAtom* dAtom;
# Line 157 | Line 201 | void Symplectic::integrate( void ){
201  
202    // initialize the forces the before the first step
203  
204 <  myFF->doForces(1,0);
204 >  myFF->doForces(1,1);
205    
206    if( entry_plug->setTemp ){
207      
# Line 169 | Line 213 | void Symplectic::integrate( void ){
213    
214    calcPot = 0;
215  
216 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
217 +    calcStress = 1;
218 +  } else {
219 +    calcStress = 0;
220 +  }
221 +
222    if( n_constrained ){
223  
224      double *Rx = new double[nAtoms];
# Line 185 | Line 235 | void Symplectic::integrate( void ){
235      
236  
237      for( tl=0; tl < n_loops; tl++ ){
238 +
239 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
240 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
241        
242        for( j=0; j<nAtoms; j++ ){
243  
# Line 285 | Line 338 | void Symplectic::integrate( void ){
338        
339        // calculate the forces
340        
341 <      myFF->doForces(calcPot, 0);
341 >      myFF->doForces(calcPot, calcStress);
342        
343        // move b
344  
# Line 350 | Line 403 | void Symplectic::integrate( void ){
403          }
404        }
405      
406 +
407 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
408 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
409 +
410 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
411 +        tStats->getPressureTensor(press);
412 +        myES->NoseHooverAndersonNPT( dt,
413 +                                     tStats->getKinetic(),
414 +                                     press);
415 +      }
416 +
417        time = tl + 1;
418        
419        if( entry_plug->setTemp ){
420          if( !(time % vel_n) ) tStats->velocitize();
421        }
422        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
423 <      if( !((time+1) % status_n) ) calcPot = 1;
424 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
423 >      if( !((time+1) % status_n) ) {
424 >        calcPot = 1;
425 >        calcStress = 1;
426 >      }
427 >      if( !(time % status_n) ){
428 >        e_out->writeStat( time * dt );
429 >        calcPot = 0;
430 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
431 >        else calcStress = 0;
432 >      }
433      }
434    }
435    else{
# Line 367 | Line 439 | void Symplectic::integrate( void ){
439        kE = 0.0;
440        rot_kE= 0.0;
441        trans_kE = 0.0;
442 +
443 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
444 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
445        
446        for( i=0; i<nAtoms; i++ ){
447          
# Line 455 | Line 530 | void Symplectic::integrate( void ){
530        
531        // calculate the forces
532        
533 <      myFF->doForces(calcPot,0);
533 >      myFF->doForces(calcPot,calcStress);
534        
535        for( i=0; i< nAtoms; i++ ){
536          
# Line 472 | Line 547 | void Symplectic::integrate( void ){
547          atoms[i]->set_vy( vy );
548          atoms[i]->set_vz( vz );
549          
550 < //      vx2 = vx * vx;
551 < //      vy2 = vy * vy;
552 < //      vz2 = vz * vz;
550 >        vx2 = vx * vx;
551 >        vy2 = vy * vy;
552 >        vz2 = vz * vz;
553          
554          if( atoms[i]->isDirectional() ){
555  
# Line 506 | Line 581 | void Symplectic::integrate( void ){
581            dAtom->setJy( ji[1] );
582            dAtom->setJz( ji[2] );
583          }
584 +
585        }
586 <      
586 >    
587 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
588 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
589 >
590 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
591 >        tStats->getPressureTensor(press);
592 >        myES->NoseHooverAndersonNPT( dt,
593 >                                     tStats->getKinetic(),
594 >                                     press);
595 >      }
596 >  
597        time = tl + 1;
598        
599        if( entry_plug->setTemp ){
600          if( !(time % vel_n) ) tStats->velocitize();
601        }
602        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
603 <      if( !((time+1) % status_n) ) calcPot = 1;
604 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
603 >      if( !((time+1) % status_n) ) {
604 >        calcPot = 1;
605 >        calcStress = 1;
606 >      }
607 >      if( !(time % status_n) ){
608 >        e_out->writeStat( time * dt );
609 >        calcPot = 0;
610 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
611 >        else calcStress = 0;
612 >      }      
613      }
614    }
615  
# Line 542 | Line 636 | void Symplectic::rotate( int axes1, int axes2, double
636  
637    for(i=0; i<3; i++){
638      for(j=0; j<3; j++){
639 <      tempA[i][j] = A[i][j];
639 >      tempA[j][i] = A[i][j];
640      }
641    }
642  
# Line 601 | Line 695 | void Symplectic::rotate( int axes1, int axes2, double
695      for(j=0; j<3; j++){
696        A[j][i] = 0.0;
697        for(k=0; k<3; k++){
698 <        A[j][i] += tempA[k][i] * rot[j][k];
698 >        A[j][i] += tempA[i][k] * rot[j][k];
699        }
700      }
701    }

Diff Legend

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