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 43 by chrisfen, Wed Oct 6 20:01:07 2004 UTC vs.
Revision 221 by chrisfen, Tue Nov 23 22:48:31 2004 UTC

# Line 11 | Line 11 | using namespace std;
11  
12   using namespace std;
13  
14 < #include "Restraints.hpp"
15 < #include "SimInfo.hpp"
16 < #include "simError.h"
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  
# Line 119 | 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 188 | 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 201 | 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 <    if (vecParticles[i]->isDirectional()){
215 <      vecParticles[i]->getPos(pos);
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 >
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  
212 –      // first we calculate the derivatives
213 –      dVdrx = -kDist*delRx;
214 –      dVdry = -kDist*delRy;
215 –      dVdrz = -kDist*delRz;
216 –
304        // uTx... and vTx... are the body-fixed z and y unit vectors
305        uTx = 0.0;
306        uTy = 0.0;
# Line 222 | 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 241 | 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
245 <      restraintFrc[0] = dVdrx;
246 <      restraintFrc[1] = dVdry;
247 <      restraintFrc[2] = dVdrz;
248 <      tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
249 <
331 >      // next we calculate the restraint torques
332        restraintTrq[0] = 0.0;
333        restraintTrq[1] = 0.0;
334        restraintTrq[2] = 0.0;
# Line 264 | 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;
269 <        restraintTrq[j] *= factor;
270 <      }
349 >      // apply the lambda scaling factor to these torques
350 >      for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
351  
272 –      harmPotent += tempPotent;
273 –
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 279 | 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
283 <      // to the total forces and torques
284 <      vecParticles[i]->addFrc(restraintFrc);
360 >      // now pass this temporary torque vector to the total torque
361        vecParticles[i]->addTrq(spaceTrq);
362      }
287 –  }
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 <  int idealSize;
375 <  double pos[3];
297 <  double A[3][3];
298 <  double RfromQ[3][3];
299 <  double quat0, quat1, quat2, quat3;
300 <  double dot;
301 <  vector<double> tempZangs;
302 <  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
305 <  strcpy(fileName, "idealCrystal.in");
306 <  strcpy(angleName, "zAngle.ang");
307 <  
308 <  ifstream crystalIn(fileName);
309 <  ifstream angleIn(angleName);
377 >  char zOutName[200];
378  
379 <  // check to see if these files are present in the execution directory
312 <  if (!crystalIn) {
313 <    sprintf(painCave.errMsg,
314 <            "Restraints Error: Unable to open idealCrystal.in for reading.\n"
315 <            "\tMake sure a ref. crystal file is in the working directory.\n");
316 <    painCave.severity = OOPSE_ERROR;
317 <    painCave.isFatal = 1;
318 <    simError();  
319 <  }
379 >  std::cerr << nIntObj << " is the number of integrable objects\n";
380  
381 <  // it's not fatal to lack a zAngle.ang file, it just means you're starting
322 <  // from the ideal crystal state
323 <  if (!angleIn) {
324 <    sprintf(painCave.errMsg,
325 <            "Restraints Warning: The lack of a zAngle.ang file is mildly\n"
326 <            "\tunsettling... This means the simulation is starting from the\n"
327 <            "\tidealCrystal.in reference configuration, so the omega values\n"
328 <            "\twill all be set to zero. If this is not the case, the energy\n"
329 <            "\tcalculations will be wrong.\n");
330 <    painCave.severity = OOPSE_WARNING;
331 <    painCave.isFatal = 0;
332 <    simError();  
333 <  }
334 <
335 <  // A rather specific reader for OOPSE .eor files...
336 <  // Let's read in the perfect crystal file
337 <  crystalIn.getline(inLine,999,'\n');
338 <  // check to see if the crystal file is the same length as starting config.
339 <  token = strtok(inLine,delimit);
340 <  strcpy(inValue,token);
341 <  idealSize = atoi(inValue);
342 <  if (idealSize != vecParticles.size()) {
343 <    sprintf(painCave.errMsg,
344 <            "Restraints Error: Reference crystal file is not valid.\n"
345 <            "\tMake sure the idealCrystal.in file is the same size as the\n"
346 <            "\tstarting configuration. Using an incompatable crystal will\n"
347 <            "\tlead to energy calculation failures.\n");
348 <    painCave.severity = OOPSE_ERROR;
349 <    painCave.isFatal = 1;
350 <    simError();  
351 <  }
352 <  // else, the file is okay... let's continue
353 <  crystalIn.getline(inLine,999,'\n');
381 >  //#ifndef IS_MPI
382    
355 –  for (i=0; i<vecParticles.size(); i++) {
356 –    crystalIn.getline(inLine,999,'\n');
357 –    token = strtok(inLine,delimit);
358 –    token = strtok(NULL,delimit);
359 –    strcpy(inValue,token);
360 –    cofmPosX.push_back(atof(inValue));
361 –    token = strtok(NULL,delimit);
362 –    strcpy(inValue,token);
363 –    cofmPosY.push_back(atof(inValue));
364 –    token = strtok(NULL,delimit);
365 –    strcpy(inValue,token);
366 –    cofmPosZ.push_back(atof(inValue));
367 –    token = strtok(NULL,delimit);
368 –    token = strtok(NULL,delimit);
369 –    token = strtok(NULL,delimit);
370 –    token = strtok(NULL,delimit);
371 –    strcpy(inValue,token);
372 –    quat0 = atof(inValue);
373 –    token = strtok(NULL,delimit);
374 –    strcpy(inValue,token);
375 –    quat1 = atof(inValue);
376 –    token = strtok(NULL,delimit);
377 –    strcpy(inValue,token);
378 –    quat2 = atof(inValue);
379 –    token = strtok(NULL,delimit);
380 –    strcpy(inValue,token);
381 –    quat3 = atof(inValue);
382 –
383 –    // now build the rotation matrix and find the unit vectors
384 –    RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;
385 –    RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);
386 –    RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);
387 –    RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);
388 –    RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;
389 –    RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);
390 –    RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);
391 –    RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);
392 –    RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;
393 –    
394 –    normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
395 –                     + RfromQ[2][2]*RfromQ[2][2]);
396 –    uX0.push_back(RfromQ[2][0]/normalize);
397 –    uY0.push_back(RfromQ[2][1]/normalize);
398 –    uZ0.push_back(RfromQ[2][2]/normalize);
399 –
400 –    normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]
401 –                     + RfromQ[1][2]*RfromQ[1][2]);
402 –    vX0.push_back(RfromQ[1][0]/normalize);
403 –    vY0.push_back(RfromQ[1][1]/normalize);
404 –    vZ0.push_back(RfromQ[1][2]/normalize);
405 –  }
406 –  crystalIn.close();
407 –
408 –  // now we read in the zAngle.ang file
409 –  if (angleIn){
410 –    angleIn.getline(inLine,999,'\n');
411 –    angleIn.getline(inLine,999,'\n');
412 –    while (!angleIn.eof()) {
413 –      token = strtok(inLine,delimit);
414 –      strcpy(inValue,token);
415 –      tempZangs.push_back(atof(inValue));
416 –      angleIn.getline(inLine,999,'\n');
417 –    }
418 –
419 –    // test to make sure the zAngle.ang file is the proper length
420 –    if (tempZangs.size() == vecParticles.size())
421 –      for (i=0; i<vecParticles.size(); i++)
422 –        vecParticles[i]->setZangle(tempZangs[i]);
423 –    else {
424 –      sprintf(painCave.errMsg,
425 –              "Restraints Error: the supplied zAngle file is not valid.\n"
426 –              "\tMake sure the zAngle.ang file matches with the initial\n"
427 –              "\tconfiguration (i.e. they're the same length). Using the wrong\n"
428 –              "\tzAngle file will lead to errors in the energy calculations.\n");
429 –      painCave.severity = OOPSE_ERROR;
430 –      painCave.isFatal = 1;
431 –      simError();  
432 –    }
433 –  }
434 –  angleIn.close();
435 –
436 –  return;
437 – }
438 –
439 – void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
440 –
441 –  char zOutName[200];
442 –
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