| 1 |
+ |
// Thermodynamic integration is not multiprocessor friendly right now |
| 2 |
+ |
|
| 3 |
|
#include <iostream> |
| 4 |
|
#include <stdlib.h> |
| 5 |
|
#include <cstdio> |
| 18 |
|
#define PI 3.14159265359 |
| 19 |
|
#define TWO_PI 6.28318530718 |
| 20 |
|
|
| 21 |
< |
Restraints::Restraints(int nMolInfo, double lambdaVal, double lambdaExp){ |
| 20 |
< |
nMol = nMolInfo; |
| 21 |
> |
Restraints::Restraints(double lambdaVal, double lambdaExp){ |
| 22 |
|
lambdaValue = lambdaVal; |
| 23 |
|
lambdaK = lambdaExp; |
| 24 |
< |
|
| 24 |
> |
vector<double> resConsts; |
| 25 |
|
const char *jolt = " \t\n;,"; |
| 26 |
|
|
| 27 |
|
#ifdef IS_MPI |
| 34 |
|
|
| 35 |
|
if (!springs) { |
| 36 |
|
sprintf(painCave.errMsg, |
| 37 |
< |
"In Restraints: Unable to open HarmSpringConsts.txt for reading.\n" |
| 38 |
< |
"\tDefault spring constants will be loaded. If you want to specify\n" |
| 39 |
< |
"\tspring constants, include a three line HarmSpringConsts.txt file\n" |
| 40 |
< |
"\tin the current directory.\n"); |
| 37 |
> |
"Unable to open HarmSpringConsts.txt for reading, so the\n" |
| 38 |
> |
"\tdefault spring constants will be loaded. If you want\n" |
| 39 |
> |
"\tto specify spring constants, include a three line\n" |
| 40 |
> |
"\tHarmSpringConsts.txt file in the execution directory.\n"); |
| 41 |
|
painCave.severity = OOPSE_WARNING; |
| 42 |
|
painCave.isFatal = 0; |
| 43 |
|
simError(); |
| 49 |
|
} else { |
| 50 |
|
|
| 51 |
|
springs.getline(inLine,999,'\n'); |
| 52 |
+ |
// the file is blank! |
| 53 |
+ |
if (springs.eof()){ |
| 54 |
+ |
sprintf(painCave.errMsg, |
| 55 |
+ |
"HarmSpringConsts.txt file is not valid.\n" |
| 56 |
+ |
"\tThe file should contain four rows, the last three containing\n" |
| 57 |
+ |
"\ta label and the spring constant value. They should be listed\n" |
| 58 |
+ |
"\tin the following order: kDist (positional restrant), kTheta\n" |
| 59 |
+ |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
| 60 |
+ |
"\trestraint: rotation about the z-axis).\n"); |
| 61 |
+ |
painCave.severity = OOPSE_ERROR; |
| 62 |
+ |
painCave.isFatal = 1; |
| 63 |
+ |
simError(); |
| 64 |
+ |
} |
| 65 |
+ |
// read in spring constants and check to make sure it is a valid file |
| 66 |
|
springs.getline(inLine,999,'\n'); |
| 67 |
< |
token = strtok(inLine,jolt); |
| 68 |
< |
token = strtok(NULL,jolt); |
| 69 |
< |
strcpy(inValue,token); |
| 70 |
< |
kDist = (atof(inValue)); |
| 71 |
< |
springs.getline(inLine,999,'\n'); |
| 72 |
< |
token = strtok(inLine,jolt); |
| 73 |
< |
token = strtok(NULL,jolt); |
| 74 |
< |
strcpy(inValue,token); |
| 75 |
< |
kTheta = (atof(inValue)); |
| 76 |
< |
springs.getline(inLine,999,'\n'); |
| 77 |
< |
token = strtok(inLine,jolt); |
| 78 |
< |
token = strtok(NULL,jolt); |
| 79 |
< |
strcpy(inValue,token); |
| 80 |
< |
kOmega = (atof(inValue)); |
| 81 |
< |
springs.close(); |
| 67 |
> |
while (!springs.eof()){ |
| 68 |
> |
if (NULL != inLine){ |
| 69 |
> |
token = strtok(inLine,jolt); |
| 70 |
> |
token = strtok(NULL,jolt); |
| 71 |
> |
if (NULL != token){ |
| 72 |
> |
strcpy(inValue,token); |
| 73 |
> |
resConsts.push_back(atof(inValue)); |
| 74 |
> |
} |
| 75 |
> |
} |
| 76 |
> |
springs.getline(inLine,999,'\n'); |
| 77 |
> |
} |
| 78 |
> |
if (resConsts.size() == 3){ |
| 79 |
> |
kDist = resConsts[0]; |
| 80 |
> |
kTheta = resConsts[1]; |
| 81 |
> |
kOmega = resConsts[2]; |
| 82 |
> |
} |
| 83 |
> |
else { |
| 84 |
> |
sprintf(painCave.errMsg, |
| 85 |
> |
"HarmSpringConsts.txt file is not valid.\n" |
| 86 |
> |
"\tThe file should contain four rows, the last three containing\n" |
| 87 |
> |
"\ta label and the spring constant value. They should be listed\n" |
| 88 |
> |
"\tin the following order: kDist (positional restrant), kTheta\n" |
| 89 |
> |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
| 90 |
> |
"\trestraint: rotation about the z-axis).\n"); |
| 91 |
> |
painCave.severity = OOPSE_ERROR; |
| 92 |
> |
painCave.isFatal = 1; |
| 93 |
> |
simError(); |
| 94 |
> |
} |
| 95 |
|
} |
| 96 |
|
#ifdef IS_MPI |
| 97 |
|
} |
| 200 |
|
|
| 201 |
|
factor = 1 - pow(lambdaValue, lambdaK); |
| 202 |
|
|
| 203 |
< |
for (i=0; i<nMol; i++){ |
| 203 |
> |
for (i=0; i<vecParticles.size(); i++){ |
| 204 |
|
if (vecParticles[i]->isDirectional()){ |
| 205 |
|
vecParticles[i]->getPos(pos); |
| 206 |
|
vecParticles[i]->getA(A); |
| 209 |
|
omegaPass = vecParticles[i]->getZangle(); |
| 210 |
|
Calc_body_omegaVal( A, omegaPass ); |
| 211 |
|
|
| 184 |
– |
if (omega > PI || omega < -PI) |
| 185 |
– |
cout << "oops... " << omega << "\n"; |
| 186 |
– |
|
| 212 |
|
// first we calculate the derivatives |
| 213 |
|
dVdrx = -kDist*delRx; |
| 214 |
|
dVdry = -kDist*delRy; |
| 292 |
|
} |
| 293 |
|
|
| 294 |
|
void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){ |
| 295 |
+ |
int idealSize; |
| 296 |
|
double pos[3]; |
| 297 |
|
double A[3][3]; |
| 298 |
|
double RfromQ[3][3]; |
| 299 |
|
double quat0, quat1, quat2, quat3; |
| 300 |
|
double dot; |
| 301 |
< |
// char *token; |
| 276 |
< |
// char fileName[200]; |
| 277 |
< |
// char angleName[200]; |
| 278 |
< |
// char inLine[1000]; |
| 279 |
< |
// char inValue[200]; |
| 301 |
> |
vector<double> tempZangs; |
| 302 |
|
const char *delimit = " \t\n;,"; |
| 303 |
|
|
| 304 |
|
//open the idealCrystal.in file and zAngle.ang file |
| 308 |
|
ifstream crystalIn(fileName); |
| 309 |
|
ifstream angleIn(angleName); |
| 310 |
|
|
| 311 |
+ |
// 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 reference crystal file is in the current directory.\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(); |
| 295 |
– |
|
| 296 |
– |
return; |
| 319 |
|
} |
| 320 |
|
|
| 321 |
+ |
// 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, you should\n" |
| 329 |
< |
"\tquestion your results.\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 |
|
} |
| 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'); |
| 354 |
|
|
| 355 |
< |
for (i=0; i<nMol; i++) { |
| 355 |
> |
for (i=0; i<vecParticles.size(); i++) { |
| 356 |
|
crystalIn.getline(inLine,999,'\n'); |
| 357 |
|
token = strtok(inLine,delimit); |
| 358 |
|
token = strtok(NULL,delimit); |
| 403 |
|
vY0.push_back(RfromQ[1][1]/normalize); |
| 404 |
|
vZ0.push_back(RfromQ[1][2]/normalize); |
| 405 |
|
} |
| 406 |
+ |
crystalIn.close(); |
| 407 |
|
|
| 408 |
< |
// now we can read in the zAngle.ang file |
| 408 |
> |
// now we read in the zAngle.ang file |
| 409 |
|
if (angleIn){ |
| 410 |
|
angleIn.getline(inLine,999,'\n'); |
| 411 |
< |
for (i=0; i<nMol; i++) { |
| 412 |
< |
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 |
< |
vecParticles[i]->setZangle(atof(inValue)); |
| 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 |
|
} |
| 444 |
|
|
| 445 |
|
ofstream angleOut(zOutName); |
| 446 |
|
angleOut << "This file contains the omega values for the .eor file\n"; |
| 447 |
< |
for (i=0; i<nMol; i++) { |
| 447 |
> |
for (i=0; i<vecParticles.size(); i++) { |
| 448 |
|
angleOut << vecParticles[i]->getZangle() << "\n"; |
| 449 |
|
} |
| 450 |
|
return; |
