ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/restraints/Restraints.cpp
(Generate patch)

Comparing trunk/src/restraints/Restraints.cpp (file contents):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 221 by chrisfen, Tue Nov 23 22:48:31 2004 UTC

# Line 1 | Line 1
1 + // Thermodynamic integration is not multiprocessor friendly right now
2 +
3   #include <iostream>
4   #include <stdlib.h>
5   #include <cstdio>
# Line 12 | Line 14 | using namespace std;
14   #include "restraints/Restraints.hpp"
15   #include "brains/SimInfo.hpp"
16   #include "utils/simError.h"
17 + #include "io/basic_ifstrstream.hpp"
18  
19 + #ifdef IS_MPI
20 + #include<mpi.h>
21 + #include "brains/mpiSimulation.hpp"
22 + #endif // is_mpi
23 +
24   #define PI 3.14159265359
25   #define TWO_PI 6.28318530718
26  
27   Restraints::Restraints(double lambdaVal, double lambdaExp){
28    lambdaValue = lambdaVal;
29    lambdaK = lambdaExp;
30 <
30 >  vector<double> resConsts;
31    const char *jolt = " \t\n;,";
32  
33   #ifdef IS_MPI
# Line 32 | Line 40 | Restraints::Restraints(double lambdaVal, double lambda
40      
41      if (!springs) {
42        sprintf(painCave.errMsg,
43 <              "In Restraints: Unable to open HarmSpringConsts.txt for reading.\n"
44 <              "\tDefault spring constants will be loaded. If you want to specify\n"
45 <              "\tspring constants, include a three line HarmSpringConsts.txt file\n"
46 <              "\tin the current directory.\n");
43 >              "Unable to open HarmSpringConsts.txt for reading, so the\n"
44 >              "\tdefault spring constants will be loaded. If you want\n"
45 >              "\tto specify spring constants, include a three line\n"
46 >              "\tHarmSpringConsts.txt file in the execution directory.\n");
47        painCave.severity = OOPSE_WARNING;
48        painCave.isFatal = 0;
49        simError();  
# Line 47 | Line 55 | Restraints::Restraints(double lambdaVal, double lambda
55      } else  {
56        
57        springs.getline(inLine,999,'\n');
58 +      // the file is blank!
59 +      if (springs.eof()){
60 +      sprintf(painCave.errMsg,
61 +              "HarmSpringConsts.txt file is not valid.\n"
62 +              "\tThe file should contain four rows, the last three containing\n"
63 +              "\ta label and the spring constant value. They should be listed\n"
64 +              "\tin the following order: kDist (positional restrant), kTheta\n"
65 +              "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
66 +              "\trestraint: rotation about the z-axis).\n");
67 +      painCave.severity = OOPSE_ERROR;
68 +      painCave.isFatal = 1;
69 +      simError();  
70 +      }
71 +      // read in spring constants and check to make sure it is a valid file
72        springs.getline(inLine,999,'\n');
73 <      token = strtok(inLine,jolt);
74 <      token = strtok(NULL,jolt);
75 <      strcpy(inValue,token);
76 <      kDist = (atof(inValue));
77 <      springs.getline(inLine,999,'\n');
78 <      token = strtok(inLine,jolt);
79 <      token = strtok(NULL,jolt);
80 <      strcpy(inValue,token);
81 <      kTheta = (atof(inValue));
82 <      springs.getline(inLine,999,'\n');
83 <      token = strtok(inLine,jolt);
84 <      token = strtok(NULL,jolt);
85 <      strcpy(inValue,token);
86 <      kOmega = (atof(inValue));
87 <      springs.close();
73 >      while (!springs.eof()){
74 >        if (NULL != inLine){
75 >          token = strtok(inLine,jolt);
76 >          token = strtok(NULL,jolt);
77 >          if (NULL != token){
78 >            strcpy(inValue,token);
79 >            resConsts.push_back(atof(inValue));
80 >          }
81 >        }
82 >        springs.getline(inLine,999,'\n');
83 >      }
84 >      if (resConsts.size() == 3){
85 >        kDist = resConsts[0];
86 >        kTheta = resConsts[1];
87 >        kOmega = resConsts[2];
88 >      }
89 >      else {
90 >        sprintf(painCave.errMsg,
91 >                "HarmSpringConsts.txt file is not valid.\n"
92 >                "\tThe file should contain four rows, the last three containing\n"
93 >                "\ta label and the spring constant value. They should be listed\n"
94 >                "\tin the following order: kDist (positional restrant), kTheta\n"
95 >                "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
96 >                "\trestraint: rotation about the z-axis).\n");
97 >        painCave.severity = OOPSE_ERROR;
98 >        painCave.isFatal = 1;
99 >        simError();  
100 >      }
101      }
102   #ifdef IS_MPI
103    }
# Line 90 | Line 125 | Restraints::~Restraints(){
125   Restraints::~Restraints(){
126   }
127  
128 < void Restraints::Calc_rVal(double position[3], int currentMol){
129 <  delRx = position[0] - cofmPosX[currentMol];
130 <  delRy = position[1] - cofmPosY[currentMol];
131 <  delRz = position[2] - cofmPosZ[currentMol];
128 > void Restraints::Calc_rVal(double position[3], double refPosition[3]){
129 >  delRx = position[0] - refPosition[0];
130 >  delRy = position[1] - refPosition[1];
131 >  delRz = position[2] - refPosition[2];
132  
133    return;
134   }
135  
136 < void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
137 <  ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol]
138 <    + matrix[0][2]*uZ0[currentMol];
139 <  ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol]
140 <    + matrix[1][2]*uZ0[currentMol];
141 <  ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol]
142 <    + matrix[2][2]*uZ0[currentMol];
136 > void Restraints::Calc_body_thetaVal(double matrix[3][3], double refUnit[3]){
137 >  ub0x = matrix[0][0]*refUnit[0] + matrix[0][1]*refUnit[1]
138 >    + matrix[0][2]*refUnit[2];
139 >  ub0y = matrix[1][0]*refUnit[0] + matrix[1][1]*refUnit[1]
140 >    + matrix[1][2]*refUnit[2];
141 >  ub0z = matrix[2][0]*refUnit[0] + matrix[2][1]*refUnit[1]
142 >    + matrix[2][2]*refUnit[2];
143  
144    normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z);
145    ub0x = ub0x/normalize;
# Line 159 | Line 194 | double Restraints::Calc_Restraint_Forces(vector<StuntD
194   double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
195    double pos[3];
196    double A[3][3];
197 +  double refPos[3];
198 +  double refVec[3];
199    double tolerance;
200    double tempPotent;
201    double factor;
202    double spaceTrq[3];
203    double omegaPass;
204 +  GenericData* data;
205 +  DoubleGenericData* doubleData;
206  
207    tolerance = 5.72957795131e-7;
208  
# Line 172 | Line 211 | double Restraints::Calc_Restraint_Forces(vector<StuntD
211    factor = 1 - pow(lambdaValue, lambdaK);
212  
213    for (i=0; i<vecParticles.size(); i++){
214 +    // obtain the current and reference positions
215 +    vecParticles[i]->getPos(pos);
216 +
217 +    data = vecParticles[i]->getProperty("refPosX");
218 +    if (data){
219 +      doubleData = dynamic_cast<DoubleGenericData*>(data);
220 +      if (!doubleData){
221 +        cerr << "Can't obtain refPosX from StuntDouble\n";
222 +        return 0.0;
223 +      }
224 +      else refPos[0] = doubleData->getData();
225 +    }
226 +    data = vecParticles[i]->getProperty("refPosY");
227 +    if (data){
228 +      doubleData = dynamic_cast<DoubleGenericData*>(data);
229 +      if (!doubleData){
230 +        cerr << "Can't obtain refPosY from StuntDouble\n";
231 +        return 0.0;
232 +      }
233 +      else refPos[1] = doubleData->getData();
234 +    }
235 +    data = vecParticles[i]->getProperty("refPosZ");
236 +    if (data){
237 +      doubleData = dynamic_cast<DoubleGenericData*>(data);
238 +      if (!doubleData){
239 +        cerr << "Can't obtain refPosZ from StuntDouble\n";
240 +        return 0.0;
241 +      }
242 +      else refPos[2] = doubleData->getData();
243 +    }
244 +
245 +    // calculate the displacement
246 +    Calc_rVal( pos, refPos );
247 +
248 +    // calculate the derivatives
249 +    dVdrx = -kDist*delRx;
250 +    dVdry = -kDist*delRy;
251 +    dVdrz = -kDist*delRz;
252 +    
253 +    // next we calculate the restraint forces
254 +    restraintFrc[0] = dVdrx;
255 +    restraintFrc[1] = dVdry;
256 +    restraintFrc[2] = dVdrz;
257 +    tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
258 +
259 +    // apply the lambda scaling factor to the forces
260 +    for (j = 0; j < 3; j++) restraintFrc[j] *= factor;
261 +
262 +    // and add the temporary force to the total force
263 +    vecParticles[i]->addFrc(restraintFrc);
264 +
265 +    // if the particle is directional, we accumulate the rot. restraints
266      if (vecParticles[i]->isDirectional()){
267 <      vecParticles[i]->getPos(pos);
267 >
268 >      // get the current rotation matrix and reference vector
269        vecParticles[i]->getA(A);
270 <      Calc_rVal( pos, i );
271 <      Calc_body_thetaVal( A, i );
270 >      
271 >      data = vecParticles[i]->getProperty("refVectorX");
272 >      if (data){
273 >        doubleData = dynamic_cast<DoubleGenericData*>(data);
274 >        if (!doubleData){
275 >          cerr << "Can't obtain refVectorX from StuntDouble\n";
276 >          return 0.0;
277 >        }
278 >        else refVec[0] = doubleData->getData();
279 >      }
280 >      data = vecParticles[i]->getProperty("refVectorY");
281 >      if (data){
282 >        doubleData = dynamic_cast<DoubleGenericData*>(data);
283 >        if (!doubleData){
284 >          cerr << "Can't obtain refVectorY from StuntDouble\n";
285 >          return 0.0;
286 >        }
287 >        else refVec[1] = doubleData->getData();
288 >      }
289 >      data = vecParticles[i]->getProperty("refVectorZ");
290 >      if (data){
291 >        doubleData = dynamic_cast<DoubleGenericData*>(data);
292 >        if (!doubleData){
293 >          cerr << "Can't obtain refVectorZ from StuntDouble\n";
294 >          return 0.0;
295 >        }
296 >        else refVec[2] = doubleData->getData();
297 >      }
298 >      
299 >      // calculate the theta and omega displacements
300 >      Calc_body_thetaVal( A, refVec );
301        omegaPass = vecParticles[i]->getZangle();
302        Calc_body_omegaVal( A, omegaPass );
303  
183 –      if (omega > PI || omega < -PI)
184 –        cout << "oops... " << omega << "\n";
185 –
186 –      // first we calculate the derivatives
187 –      dVdrx = -kDist*delRx;
188 –      dVdry = -kDist*delRy;
189 –      dVdrz = -kDist*delRz;
190 –
304        // uTx... and vTx... are the body-fixed z and y unit vectors
305        uTx = 0.0;
306        uTy = 0.0;
# Line 196 | Line 309 | double Restraints::Calc_Restraint_Forces(vector<StuntD
309        vTy = 1.0;
310        vTz = 0.0;
311  
312 <      dVdux = 0;
313 <      dVduy = 0;
314 <      dVduz = 0;
315 <      dVdvx = 0;
316 <      dVdvy = 0;
317 <      dVdvz = 0;
312 >      dVdux = 0.0;
313 >      dVduy = 0.0;
314 >      dVduz = 0.0;
315 >      dVdvx = 0.0;
316 >      dVdvy = 0.0;
317 >      dVdvz = 0.0;
318  
319        if (fabs(theta) > tolerance) {
320          dVdux = -(kTheta*theta/sin(theta))*ub0x;
# Line 215 | Line 328 | double Restraints::Calc_Restraint_Forces(vector<StuntD
328          dVdvz = -(kOmega*omega/sin(omega))*vb0z;
329        }
330  
331 <      // next we calculate the restraint forces and torques
219 <      restraintFrc[0] = dVdrx;
220 <      restraintFrc[1] = dVdry;
221 <      restraintFrc[2] = dVdrz;
222 <      tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
223 <
331 >      // next we calculate the restraint torques
332        restraintTrq[0] = 0.0;
333        restraintTrq[1] = 0.0;
334        restraintTrq[2] = 0.0;
# Line 238 | Line 346 | double Restraints::Calc_Restraint_Forces(vector<StuntD
346          tempPotent += 0.5*(kTheta*theta*theta);
347        }
348  
349 <      for (j = 0; j < 3; j++) {
350 <        restraintFrc[j] *= factor;
243 <        restraintTrq[j] *= factor;
244 <      }
349 >      // apply the lambda scaling factor to these torques
350 >      for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
351  
246 –      harmPotent += tempPotent;
247 –
352        // now we need to convert from body-fixed torques to space-fixed torques
353        spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1]
354          + A[2][0]*restraintTrq[2];
# Line 253 | Line 357 | double Restraints::Calc_Restraint_Forces(vector<StuntD
357        spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1]
358          + A[2][2]*restraintTrq[2];
359  
360 <      // now it's time to pass these temporary forces and torques
257 <      // to the total forces and torques
258 <      vecParticles[i]->addFrc(restraintFrc);
360 >      // now pass this temporary torque vector to the total torque
361        vecParticles[i]->addTrq(spaceTrq);
362      }
261 –  }
363  
364 <  // and we can return the appropriately scaled potential energy
364 >    // update the total harmonic potential with this object's contribution
365 >    harmPotent += tempPotent;
366 >  }
367 >  
368 >  // we can finish by returning the appropriately scaled potential energy
369    tempPotent = harmPotent * factor;
370    return tempPotent;
371   }
372  
373 < void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
374 <  double pos[3];
375 <  double A[3][3];
271 <  double RfromQ[3][3];
272 <  double quat0, quat1, quat2, quat3;
273 <  double dot;
274 < //   char *token;
275 < //   char fileName[200];
276 < //   char angleName[200];
277 < //   char inLine[1000];
278 < //   char inValue[200];
279 <  const char *delimit = " \t\n;,";
373 > void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles,
374 >                                   int currTime,
375 >                                   int nIntObj){
376  
377 <  //open the idealCrystal.in file and zAngle.ang file
282 <  strcpy(fileName, "idealCrystal.in");
283 <  strcpy(angleName, "zAngle.ang");
284 <  
285 <  ifstream crystalIn(fileName);
286 <  ifstream angleIn(angleName);
377 >  char zOutName[200];
378  
379 <  if (!crystalIn) {
289 <    sprintf(painCave.errMsg,
290 <            "Restraints Error: Unable to open idealCrystal.in for reading.\n"
291 <            "\tMake sure a reference crystal file is in the current directory.\n");
292 <    painCave.isFatal = 1;
293 <    simError();  
294 <    
295 <    return;
296 <  }
379 >  std::cerr << nIntObj << " is the number of integrable objects\n";
380  
381 <  if (!angleIn) {
299 <    sprintf(painCave.errMsg,
300 <            "Restraints Warning: The lack of a zAngle.ang file is mildly\n"
301 <            "\tunsettling... This means the simulation is starting from the\n"
302 <            "\tidealCrystal.in reference configuration, so the omega values\n"
303 <            "\twill all be set to zero. If this is not the case, you should\n"
304 <            "\tquestion your results.\n");
305 <    painCave.isFatal = 0;
306 <    simError();  
307 <  }
308 <
309 <  // A rather specific reader for OOPSE .eor files...
310 <  // Let's read in the perfect crystal file
311 <  crystalIn.getline(inLine,999,'\n');
312 <  crystalIn.getline(inLine,999,'\n');
381 >  //#ifndef IS_MPI
382    
314 –  for (i=0; i<vecParticles.size(); i++) {
315 –    crystalIn.getline(inLine,999,'\n');
316 –    token = strtok(inLine,delimit);
317 –    token = strtok(NULL,delimit);
318 –    strcpy(inValue,token);
319 –    cofmPosX.push_back(atof(inValue));
320 –    token = strtok(NULL,delimit);
321 –    strcpy(inValue,token);
322 –    cofmPosY.push_back(atof(inValue));
323 –    token = strtok(NULL,delimit);
324 –    strcpy(inValue,token);
325 –    cofmPosZ.push_back(atof(inValue));
326 –    token = strtok(NULL,delimit);
327 –    token = strtok(NULL,delimit);
328 –    token = strtok(NULL,delimit);
329 –    token = strtok(NULL,delimit);
330 –    strcpy(inValue,token);
331 –    quat0 = atof(inValue);
332 –    token = strtok(NULL,delimit);
333 –    strcpy(inValue,token);
334 –    quat1 = atof(inValue);
335 –    token = strtok(NULL,delimit);
336 –    strcpy(inValue,token);
337 –    quat2 = atof(inValue);
338 –    token = strtok(NULL,delimit);
339 –    strcpy(inValue,token);
340 –    quat3 = atof(inValue);
341 –
342 –    // now build the rotation matrix and find the unit vectors
343 –    RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;
344 –    RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);
345 –    RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);
346 –    RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);
347 –    RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;
348 –    RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);
349 –    RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);
350 –    RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);
351 –    RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;
352 –    
353 –    normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
354 –                     + RfromQ[2][2]*RfromQ[2][2]);
355 –    uX0.push_back(RfromQ[2][0]/normalize);
356 –    uY0.push_back(RfromQ[2][1]/normalize);
357 –    uZ0.push_back(RfromQ[2][2]/normalize);
358 –
359 –    normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]
360 –                     + RfromQ[1][2]*RfromQ[1][2]);
361 –    vX0.push_back(RfromQ[1][0]/normalize);
362 –    vY0.push_back(RfromQ[1][1]/normalize);
363 –    vZ0.push_back(RfromQ[1][2]/normalize);
364 –  }
365 –
366 –  // now we can read in the zAngle.ang file
367 –  if (angleIn){
368 –    angleIn.getline(inLine,999,'\n');
369 –    for (i=0; i<vecParticles.size(); i++) {
370 –      angleIn.getline(inLine,999,'\n');
371 –      token = strtok(inLine,delimit);
372 –      strcpy(inValue,token);
373 –      vecParticles[i]->setZangle(atof(inValue));
374 –    }
375 –  }
376 –
377 –  return;
378 – }
379 –
380 – void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
381 –
382 –  char zOutName[200];
383 –
383    strcpy(zOutName,"zAngle.ang");
384 <
384 >  
385    ofstream angleOut(zOutName);
386 <  angleOut << "This file contains the omega values for the .eor file\n";
386 >  angleOut << currTime << ": omega values at this time\n";
387    for (i=0; i<vecParticles.size(); i++) {
388      angleOut << vecParticles[i]->getZangle() << "\n";
389    }
390 +
391    return;
392   }
393  

Diff Legend

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