ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/OOPSE/libmdtools/NPTf.cpp
(Generate patch)

Comparing branches/new-templateless/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 848, Fri Oct 31 18:28:53 2003 UTC vs.
Revision 852 by mmeineke, Thu Nov 6 18:20:47 2003 UTC

# Line 1 | Line 1
1 + #include <stdlib.h>
2   #include <math.h>
3 + #include <string.h>
4 +
5   #include "Atom.hpp"
6   #include "SRI.hpp"
7   #include "AbstractClasses.hpp"
# Line 23 | Line 26
26   //
27   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28  
29 < template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
30 <  T( theInfo, the_ff )
29 > NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
30 >  NPT( theInfo, the_ff )
31   {
32    GenericData* data;
33 <  DoubleArrayData * etaValue;
31 <  vector<double> etaArray;
33 >  double *etaArray;
34    int i,j;
35  
36    for(i = 0; i < 3; i++){
# Line 39 | Line 41 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
41      }
42    }
43  
44 <    // retrieve eta array from simInfo if it exists
45 <    data = info->getProperty(ETAVALUE_ID);
46 <    if(data){
47 <      etaValue = dynamic_cast<DoubleArrayData*>(data);
48 <
49 <      if(etaValue){
50 <        etaArray = etaValue->getData();
51 <
52 <        for(i = 0; i < 3; i++){
53 <          for (j = 0; j < 3; j++){
54 <            eta[i][j] = etaArray[3*i+j];
55 <            oldEta[i][j] = eta[i][j];
56 <          }
57 <        }
58 <
57 <      }
44 >  // retrieve eta array from simInfo if it exists
45 >  data = info->getProperty(ETAVALUE_ID);
46 >  if(data != NULL){
47 >    
48 >    int test = data->getDarray(etaArray);
49 >    
50 >    if( test == 9 ){
51 >      
52 >      for(i = 0; i < 3; i++){
53 >        for (j = 0; j < 3; j++){
54 >          eta[i][j] = etaArray[3*i+j];
55 >          oldEta[i][j] = eta[i][j];
56 >        }
57 >      }    
58 >      delete[] etaArray;
59      }
60 <
60 >    else
61 >      std::cerr << "NPTf error: etaArray is not length 9 (actual = " << test
62 >                << ").\n"
63 >                << "            Simulation wil proceed with eta = 0;\n";
64 >  }
65   }
66  
67 < template<typename T> NPTf<T>::~NPTf() {
67 > NPTf::~NPTf() {
68  
69    // empty for now
70   }
71  
72 < template<typename T> void NPTf<T>::resetIntegrator() {
72 > void NPTf::resetIntegrator() {
73  
74    int i, j;
75  
# Line 72 | Line 77 | template<typename T> void NPTf<T>::resetIntegrator() {
77      for (j = 0; j < 3; j++)
78        eta[i][j] = 0.0;
79  
80 <  T::resetIntegrator();
80 >  NPT::resetIntegrator();
81   }
82  
83 < template<typename T> void NPTf<T>::evolveEtaA() {
83 > void NPTf::evolveEtaA() {
84  
85    int i, j;
86  
# Line 94 | Line 99 | template<typename T> void NPTf<T>::evolveEtaA() {
99        oldEta[i][j] = eta[i][j];
100   }
101  
102 < template<typename T> void NPTf<T>::evolveEtaB() {
102 > void NPTf::evolveEtaB() {
103  
104    int i,j;
105  
# Line 114 | Line 119 | template<typename T> void NPTf<T>::evolveEtaB() {
119    }
120   }
121  
122 < template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
122 > void NPTf::getVelScaleA(double sc[3], double vel[3]) {
123    int i,j;
124    double vScale[3][3];
125  
# Line 131 | Line 136 | template<typename T> void NPTf<T>::getVelScaleA(double
136    info->matVecMul3( vScale, vel, sc );
137   }
138  
139 < template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
139 > void NPTf::getVelScaleB(double sc[3], int index ){
140    int i,j;
141    double myVel[3];
142    double vScale[3][3];
143  
144 + //   std::cerr << "velScaleB chi = " << chi << "\n";
145 +
146    for (i = 0; i < 3; i++ ) {
147      for (j = 0; j < 3; j++ ) {
148        vScale[i][j] = eta[i][j];
# Line 152 | Line 159 | template<typename T> void NPTf<T>::getVelScaleB(double
159    info->matVecMul3( vScale, myVel, sc );
160   }
161  
162 < template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
163 <                                               int index, double sc[3]){
162 > void NPTf::getPosScale(double pos[3], double COM[3],
163 >                       int index, double sc[3]){
164    int j;
165    double rj[3];
166  
# Line 163 | Line 170 | template<typename T> void NPTf<T>::getPosScale(double
170    info->matVecMul3( eta, rj, sc );
171   }
172  
173 < template<typename T> void NPTf<T>::scaleSimBox( void ){
173 > void NPTf::scaleSimBox( void ){
174  
175    int i,j,k;
176    double scaleMat[3][3];
# Line 208 | Line 215 | template<typename T> void NPTf<T>::scaleSimBox( void )
215      if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
216    }
217  
218 <  if ((bigScale > 1.1) || (smallScale < 0.9)) {
218 >  if ((bigScale > 1.01) || (smallScale < 0.99)) {
219      sprintf( painCave.errMsg,
220 <             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
220 >             "NPTf error: Attempting a Box scaling of more than 1 percent.\n"
221               " Check your tauBarostat, as it is probably too small!\n\n"
222               " scaleMat = [%lf\t%lf\t%lf]\n"
223               "            [%lf\t%lf\t%lf]\n"
# Line 220 | Line 227 | template<typename T> void NPTf<T>::scaleSimBox( void )
227               scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
228      painCave.isFatal = 1;
229      simError();
230 <  } else if (offDiagMax > 0.1) {
230 >  } else if (offDiagMax > 0.01) {
231      sprintf( painCave.errMsg,
232 <             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
232 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
233               " Check your tauBarostat, as it is probably too small!\n\n"
234               " scaleMat = [%lf\t%lf\t%lf]\n"
235               "            [%lf\t%lf\t%lf]\n"
# Line 239 | Line 246 | template<typename T> void NPTf<T>::scaleSimBox( void )
246    }
247   }
248  
249 < template<typename T> bool NPTf<T>::etaConverged() {
249 > bool NPTf::etaConverged() {
250    int i;
251    double diffEta, sumEta;
252  
# Line 252 | Line 259 | template<typename T> bool NPTf<T>::etaConverged() {
259    return ( diffEta <= etaTolerance );
260   }
261  
262 < template<typename T> double NPTf<T>::getConservedQuantity(void){
262 > double NPTf::getConservedQuantity(void){
263  
264    double conservedQuantity;
265    double totalEnergy;
# Line 283 | Line 290 | template<typename T> double NPTf<T>::getConservedQuant
290    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
291      barostat_kinetic + barostat_potential;
292  
286 – //   cout.width(8);
287 – //   cout.precision(8);
288 –
289 – //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
290 – //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
291 – //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
292 –
293    return conservedQuantity;
294  
295   }
296  
297 < template<typename T> string NPTf<T>::getAdditionalParameters(void){
298 <  string parameters;
299 <  const int BUFFERSIZE = 2000; // size of the read buffer
300 <  char buffer[BUFFERSIZE];
297 > char* NPTf::getAdditionalParameters(void){
298  
299 <  sprintf(buffer,"\t%lf\t%lf;", chi, integralOfChidt);
300 <  parameters += buffer;
299 >  sprintf(addParamBuffer,
300 >          "\t%G\t%G;"
301 >          "\t%G\t%G\t%G;"
302 >          "\t%G\t%G\t%G;"
303 >          "\t%G\t%G\t%G;",
304 >          chi, integralOfChidt,
305 >          eta[0][0], eta[0][1], eta[0][2],
306 >          eta[1][0], eta[1][1], eta[1][2],
307 >          eta[2][0], eta[2][1], eta[2][2]
308 >          );
309  
310 <  for(int i = 0; i < 3; i++){
306 <    sprintf(buffer,"\t%lf\t%lf\t%lf;", eta[3*i], eta[3*i+1], eta[3*i+2]);
307 <    parameters += buffer;
308 <  }
309 <
310 <  return parameters;
311 <
310 >  return addParamBuffer;
311   }

Diff Legend

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