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 423 by mmeineke, Thu Mar 27 20:12:15 2003 UTC vs.
Revision 469 by mmeineke, Mon Apr 7 20:06:31 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    molecules = entry_plug->molecules;
44    nMols = entry_plug->n_mol;
45  
# Line 48 | 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 <  }
52 <
53 <  
56 >  }  
57  
58    // check for constraints
59  
# Line 65 | Line 68 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
68    SRI** theArray;
69    for(int i = 0; i < nMols; i++){
70      
71 <    theArray = molecules[i].getMyBonds();
72 <    for(int j=0; j<molecules[i].getNbonds(); j++){
71 >    theArray = (SRI**) molecules[i].getMyBonds();
72 >    for(int j=0; j<molecules[i].getNBonds(); j++){
73        
74        constrained = theArray[j]->is_constrained();
75        
76        if(constrained){
77          
78          dummy_plug = theArray[j]->get_constraint();
79 <        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
80 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
81 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        c_n_constrained++;
83 >        n_constrained++;
84          constrained = 0;
85        }
86      }
87  
88 <    theArray = molecules[i].getMyBends();
89 <    for(int j=0; j<molecules[i].getNbends(); j++){
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[c_n_constrained].set_a( dummy_plug->get_a() );
97 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
98 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        c_n_constrained++;
100 >        n_constrained++;
101          constrained = 0;
102        }
103      }
104  
105 <    theArray = molecules[i].getMyTorsions();
106 <    for(int j=0; j<molecules[i].getNtorsions(); j++){
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[c_n_constrained].set_a( dummy_plug->get_a() );
114 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
115 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        c_n_constrained++;
117 >        n_constrained++;
118          constrained = 0;
119        }
120      }
# Line 183 | 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 );
191 <  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 324 | Line 335 | void Symplectic::integrate( void ){
335        
336        // calculate the forces
337        
338 <      myFF->doForces(calcPot, 0);
338 >      myFF->doForces(calcPot, calcStress);
339        
340        // move b
341  
# Line 395 | Line 406 | void Symplectic::integrate( void ){
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 >        calcStress = 0;
417 >      }
418      }
419    }
420    else{
# Line 494 | Line 512 | void Symplectic::integrate( void ){
512        
513        // calculate the forces
514        
515 <      myFF->doForces(calcPot,0);
515 >      myFF->doForces(calcPot,calcStress);
516        
517        for( i=0; i< nAtoms; i++ ){
518          
# Line 553 | Line 571 | void Symplectic::integrate( void ){
571          if( !(time % vel_n) ) tStats->velocitize();
572        }
573        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
574 <      if( !((time+1) % status_n) ) calcPot = 1;
575 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
574 >      if( !((time+1) % status_n) ) {
575 >        calcPot = 1;
576 >        calcStress = 1;
577 >      }
578 >      if( !(time % status_n) ){
579 >        e_out->writeStat( time * dt );
580 >        calcPot = 0;
581 >        calcStress = 0;
582 >      }      
583      }
584    }
585  
# Line 581 | Line 606 | void Symplectic::rotate( int axes1, int axes2, double
606  
607    for(i=0; i<3; i++){
608      for(j=0; j<3; j++){
609 <      tempA[i][j] = A[i][j];
609 >      tempA[j][i] = A[i][j];
610      }
611    }
612  
# Line 640 | Line 665 | void Symplectic::rotate( int axes1, int axes2, double
665      for(j=0; j<3; j++){
666        A[j][i] = 0.0;
667        for(k=0; k<3; k++){
668 <        A[j][i] += tempA[k][i] * rot[j][k];
668 >        A[j][i] += tempA[i][k] * rot[j][k];
669        }
670      }
671    }

Diff Legend

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