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

Comparing branches/new_design/OOPSE-3.0/src/integrators/NPTxyz.cpp (file contents):
Revision 1866 by tim, Sun Dec 5 17:09:27 2004 UTC vs.
Revision 1867 by tim, Tue Dec 7 23:08:14 2004 UTC

# Line 22 | Line 22 | double NPTxyz::calcConservedQuantity(){
22      
23   double NPTxyz::calcConservedQuantity(){
24  
25 <  double conservedQuantity;
26 <  double totalEnergy;
27 <  double thermostat_kinetic;
28 <  double thermostat_potential;
29 <  double barostat_kinetic;
30 <  double barostat_potential;
31 <  double trEta;
25 >    // We need NkBT a lot, so just set it here: This is the RAW number
26 >    // of integrableObjects, so no subtraction or addition of constraints or
27 >    // orientational degrees of freedom:
28 >    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
29  
30 <  totalEnergy = thermo.getTotalE();
30 >    // fkBT is used because the thermostat operates on more degrees of freedom
31 >    // than the barostat (when there are particles with orientational degrees
32 >    // of freedom).  
33 >    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;        
34  
35 <  thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
35 >    double conservedQuantity;
36 >    double totalEnergy;
37 >    double thermostat_kinetic;
38 >    double thermostat_potential;
39 >    double barostat_kinetic;
40 >    double barostat_potential;
41 >    double trEta;
42  
43 <  thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
43 >    totalEnergy = thermo.getTotalE();
44  
45 +    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
46 +
47 +    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
48 +
49      SquareMatrix<double, 3> tmp = eta.transpose() * eta;
50      trEta = tmp.trace();
51  
52 <  barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
52 >    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
53  
54 <  barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
54 >    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
55  
56 <  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
57 <    barostat_kinetic + barostat_potential;
56 >    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
57 >        barostat_kinetic + barostat_potential;
58  
59  
60 <  return conservedQuantity;
60 >    return conservedQuantity;
61  
62   }
63  
64      
65   void NPTxyz::scaleSimBox(){
66  
67 <  int i,j,k;
68 <  Mat3x3d scaleMat;
69 <  double eta2ij, scaleFactor;
70 <  double bigScale, smallScale, offDiagMax;
71 <  Mat3x3d hm;
72 <  Mat3x3d hmnew;
67 >    int i,j,k;
68 >    Mat3x3d scaleMat;
69 >    double eta2ij, scaleFactor;
70 >    double bigScale, smallScale, offDiagMax;
71 >    Mat3x3d hm;
72 >    Mat3x3d hmnew;
73  
74  
75  
# Line 68 | Line 78 | void NPTxyz::scaleSimBox(){
78    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
79    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
80  
81 <  bigScale = 1.0;
82 <  smallScale = 1.0;
83 <  offDiagMax = 0.0;
81 >    bigScale = 1.0;
82 >    smallScale = 1.0;
83 >    offDiagMax = 0.0;
84  
85 <  for(i=0; i<3; i++){
86 <    for(j=0; j<3; j++){
87 <      scaleMat(i, j) = 0.0;
88 <      if(i==j) scaleMat(i, j) = 1.0;
85 >    for(i=0; i<3; i++){
86 >        for(j=0; j<3; j++){
87 >            scaleMat(i, j) = 0.0;
88 >            if(i==j) {
89 >                scaleMat(i, j) = 1.0;
90 >            }
91 >        }
92      }
80 –  }
93  
94 <  for(i=0;i<3;i++){
94 >    for(i=0;i<3;i++){
95  
96      // calculate the scaleFactors
97  
98 <    scaleFactor = exp(dt*eta(i, i));
98 >        scaleFactor = exp(dt*eta(i, i));
99  
100 <    scaleMat(i, i) = scaleFactor;
100 >        scaleMat(i, i) = scaleFactor;
101  
102 <    if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
103 <    if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
104 <  }
102 >        if (scaleMat(i, i) > bigScale) {
103 >            bigScale = scaleMat(i, i);
104 >        }
105 >        
106 >        if (scaleMat(i, i) < smallScale) {
107 >            smallScale = scaleMat(i, i);
108 >        }
109 >    }
110  
111 <  if ((bigScale > 1.1) || (smallScale < 0.9)) {
112 <    sprintf( painCave.errMsg,
113 <             "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
114 <             " Check your tauBarostat, as it is probably too small!\n\n"
115 <             " scaleMat = [%lf\t%lf\t%lf]\n"
116 <             "            [%lf\t%lf\t%lf]\n"
117 <             "            [%lf\t%lf\t%lf]\n",
118 <             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
119 <             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
120 <             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
121 <    painCave.isFatal = 1;
122 <    simError();
123 <  } else {
111 >    if ((bigScale > 1.1) || (smallScale < 0.9)) {
112 >        sprintf( painCave.errMsg,
113 >            "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
114 >            " Check your tauBarostat, as it is probably too small!\n\n"
115 >            " scaleMat = [%lf\t%lf\t%lf]\n"
116 >            "            [%lf\t%lf\t%lf]\n"
117 >            "            [%lf\t%lf\t%lf]\n",
118 >        scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
119 >        scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
120 >        scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
121 >        painCave.isFatal = 1;
122 >        simError();
123 >    } else {
124  
125          Mat3x3d hmat = currentSnapshot_->getHmat();
126          hmat = hmat *scaleMat;
127          currentSnapshot_->setHmat(hmat);
128 <  }
128 >    }
129   }
130  
131   void NPTxyz::loadEta() {

Diff Legend

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