ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/integrators/NPTi.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-2.0/src/integrators/NPTi.cpp (file contents):
Revision 1837 by tim, Thu Dec 2 22:15:31 2004 UTC vs.
Revision 1897 by tim, Mon Dec 20 19:49:44 2004 UTC

# Line 9 | Line 9 | static IntegratorBuilder<NPTi>* NPTiCreator = new Inte
9  
10   namespace oopse {
11  
12 – static IntegratorBuilder<NPTi>* NPTiCreator = new IntegratorBuilder<NPTi>("NPTi");
13 –
12   // Basic isotropic thermostating and barostating via the Melchionna
13   // modification of the Hoover algorithm:
14   //
# Line 23 | Line 21 | NPTi::NPTi ( SimInfo *info) : NPT(info){
21  
22   NPTi::NPTi ( SimInfo *info) : NPT(info){
23  
26 –    Globals* globals = info_->getGlobals();
27 –
24   }
25  
26   void NPTi::evolveEtaA() {
# Line 56 | Line 52 | void NPTi::getPosScale(const Vector3d& pos, const Vect
52   void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
53                             int index, Vector3d& sc){
54      /**@todo*/
55 <    sc  = oldPos[index] + pos/2.0 -COM;
55 >    sc  = (oldPos[index] + pos)/2.0 -COM;
56      sc *= eta;
57   }
58  
# Line 89 | Line 85 | double NPTi::calcConservedQuantity(){
85  
86   double NPTi::calcConservedQuantity(){
87  
88 +    chi= currentSnapshot_->getChi();
89 +    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
90 +    loadEta();
91 +    // We need NkBT a lot, so just set it here: This is the RAW number
92 +    // of integrableObjects, so no subtraction or addition of constraints or
93 +    // orientational degrees of freedom:
94 +    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
95 +
96 +    // fkBT is used because the thermostat operates on more degrees of freedom
97 +    // than the barostat (when there are particles with orientational degrees
98 +    // of freedom).  
99 +    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
100 +    
101      double conservedQuantity;
102      double Energy;
103      double thermostat_kinetic;
# Line 111 | Line 120 | double NPTi::calcConservedQuantity(){
120      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
121          barostat_kinetic + barostat_potential;
122  
123 +    std::cout << "--------------------------------------------------------------" << std::endl;
124 +    std::cout << "time: " << currentSnapshot_->getTime() << std:: endl;
125 +    std::cout << "chi : " << chi << std::endl;
126 +    std::cout << "integralOfChidt : " << integralOfChidt << std::endl;
127 +    std::cout << "eta : " << eta << std::endl;    
128 +    std::cout << "NkBT: " << NkBT << std::endl;
129 +    std::cout << "fkBT: " << fkBT << std::endl;    
130 +    std::cout << "thermostat_kinetic : " << thermostat_kinetic<< std::endl;
131 +    std::cout << "thermostat_potential : " << thermostat_potential << std::endl;
132 +    std::cout << "barostat_kinetic : " << barostat_kinetic << std::endl;
133 +    std::cout << "barostat_potential : " << barostat_potential << std::endl;
134 +    std::cout << "Total Energy: " << Energy << std::endl;
135 +    std::cout << "Conserved Quantity: " << conservedQuantity <<std::endl;
136 +    std::cout << "--------------------------------------------------------------" << std::endl;
137      return conservedQuantity;
138   }
139  
140   void NPTi::loadEta() {
141      Mat3x3d etaMat = currentSnapshot_->getEta();
142      eta = etaMat(0,0);
143 <    if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || etaMat.isDiagonal()) {
144 <        sprintf( painCave.errMsg,
145 <                 "NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
146 <        painCave.isFatal = 1;
147 <        simError();
148 <    }
143 >    //if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || !etaMat.isDiagonal()) {
144 >    //    sprintf( painCave.errMsg,
145 >    //             "NPTi error: the diagonal elements of  eta matrix are not the same or etaMat is not a diagonal matrix");
146 >    //    painCave.isFatal = 1;
147 >    //    simError();
148 >    //}
149   }
150  
151   void NPTi::saveEta() {

Diff Legend

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