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

Comparing:
branches/mmeineke/OOPSE/libmdtools/Symplectic.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Symplectic.cpp (file contents), Revision 468 by gezelter, Mon Apr 7 16:56:38 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  
# Line 144 | Line 184 | void Symplectic::integrate( void ){
184    int status_n = (int)( statusTime / dt );
185    int vel_n    = (int)( thermalTime / dt );
186  
187 <  int calcPot;
187 >  int calcPot, calcStress;
188  
189     Thermo *tStats = new Thermo( entry_plug );
190  
# Line 157 | Line 197 | void Symplectic::integrate( void ){
197  
198    // initialize the forces the before the first step
199  
200 <  myFF->doForces(1,0);
200 >  myFF->doForces(1,1);
201    
202    if( entry_plug->setTemp ){
203      
# Line 285 | Line 325 | void Symplectic::integrate( void ){
325        
326        // calculate the forces
327        
328 <      myFF->doForces(calcPot, 0);
328 >      myFF->doForces(calcPot, calcStress);
329        
330        // move b
331  
# Line 356 | Line 396 | void Symplectic::integrate( void ){
396          if( !(time % vel_n) ) tStats->velocitize();
397        }
398        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
399 <      if( !((time+1) % status_n) ) calcPot = 1;
400 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
399 >      if( !((time+1) % status_n) ) {
400 >        calcPot = 1;
401 >        calcStress = 1;
402 >      }
403 >      if( !(time % status_n) ){
404 >        e_out->writeStat( time * dt );
405 >        calcPot = 0;
406 >        calcStress = 0;
407 >      }
408      }
409    }
410    else{
# Line 455 | Line 502 | void Symplectic::integrate( void ){
502        
503        // calculate the forces
504        
505 <      myFF->doForces(calcPot,0);
505 >      myFF->doForces(calcPot,calcStress);
506        
507        for( i=0; i< nAtoms; i++ ){
508          
# Line 514 | Line 561 | void Symplectic::integrate( void ){
561          if( !(time % vel_n) ) tStats->velocitize();
562        }
563        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
564 <      if( !((time+1) % status_n) ) calcPot = 1;
565 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
564 >      if( !((time+1) % status_n) ) {
565 >        calcPot = 1;
566 >        calcStress = 1;
567 >      }
568 >      if( !(time % status_n) ){
569 >        e_out->writeStat( time * dt );
570 >        calcPot = 0;
571 >        calcStress = 0;
572 >      }      
573      }
574    }
575  
# Line 542 | Line 596 | void Symplectic::rotate( int axes1, int axes2, double
596  
597    for(i=0; i<3; i++){
598      for(j=0; j<3; j++){
599 <      tempA[i][j] = A[i][j];
599 >      tempA[j][i] = A[i][j];
600      }
601    }
602  
# Line 601 | Line 655 | void Symplectic::rotate( int axes1, int axes2, double
655      for(j=0; j<3; j++){
656        A[j][i] = 0.0;
657        for(k=0; k<3; k++){
658 <        A[j][i] += tempA[k][i] * rot[j][k];
658 >        A[j][i] += tempA[i][k] * rot[j][k];
659        }
660      }
661    }

Diff Legend

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