--- trunk/OOPSE-2.0/src/restraints/Restraints.cpp 2004/11/11 21:47:25 1731 +++ trunk/OOPSE-2.0/src/restraints/Restraints.cpp 2004/11/23 22:48:31 1772 @@ -16,6 +16,11 @@ using namespace std; #include "utils/simError.h" #include "io/basic_ifstrstream.hpp" +#ifdef IS_MPI +#include +#include "brains/mpiSimulation.hpp" +#endif // is_mpi + #define PI 3.14159265359 #define TWO_PI 6.28318530718 @@ -120,21 +125,21 @@ void Restraints::Calc_rVal(double position[3], int cur Restraints::~Restraints(){ } -void Restraints::Calc_rVal(double position[3], int currentMol){ - delRx = position[0] - cofmPosX[currentMol]; - delRy = position[1] - cofmPosY[currentMol]; - delRz = position[2] - cofmPosZ[currentMol]; +void Restraints::Calc_rVal(double position[3], double refPosition[3]){ + delRx = position[0] - refPosition[0]; + delRy = position[1] - refPosition[1]; + delRz = position[2] - refPosition[2]; return; } -void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){ - ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol] - + matrix[0][2]*uZ0[currentMol]; - ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol] - + matrix[1][2]*uZ0[currentMol]; - ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol] - + matrix[2][2]*uZ0[currentMol]; +void Restraints::Calc_body_thetaVal(double matrix[3][3], double refUnit[3]){ + ub0x = matrix[0][0]*refUnit[0] + matrix[0][1]*refUnit[1] + + matrix[0][2]*refUnit[2]; + ub0y = matrix[1][0]*refUnit[0] + matrix[1][1]*refUnit[1] + + matrix[1][2]*refUnit[2]; + ub0z = matrix[2][0]*refUnit[0] + matrix[2][1]*refUnit[1] + + matrix[2][2]*refUnit[2]; normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z); ub0x = ub0x/normalize; @@ -189,11 +194,15 @@ double Restraints::Calc_Restraint_Forces(vector vecParticles){ double pos[3]; double A[3][3]; + double refPos[3]; + double refVec[3]; double tolerance; double tempPotent; double factor; double spaceTrq[3]; double omegaPass; + GenericData* data; + DoubleGenericData* doubleData; tolerance = 5.72957795131e-7; @@ -202,19 +211,96 @@ double Restraints::Calc_Restraint_Forces(vectorisDirectional()){ - vecParticles[i]->getPos(pos); + // obtain the current and reference positions + vecParticles[i]->getPos(pos); + + data = vecParticles[i]->getProperty("refPosX"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refPosX from StuntDouble\n"; + return 0.0; + } + else refPos[0] = doubleData->getData(); + } + data = vecParticles[i]->getProperty("refPosY"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refPosY from StuntDouble\n"; + return 0.0; + } + else refPos[1] = doubleData->getData(); + } + data = vecParticles[i]->getProperty("refPosZ"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refPosZ from StuntDouble\n"; + return 0.0; + } + else refPos[2] = doubleData->getData(); + } + + // calculate the displacement + Calc_rVal( pos, refPos ); + + // calculate the derivatives + dVdrx = -kDist*delRx; + dVdry = -kDist*delRy; + dVdrz = -kDist*delRz; + + // next we calculate the restraint forces + restraintFrc[0] = dVdrx; + restraintFrc[1] = dVdry; + restraintFrc[2] = dVdrz; + tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz); + + // apply the lambda scaling factor to the forces + for (j = 0; j < 3; j++) restraintFrc[j] *= factor; + + // and add the temporary force to the total force + vecParticles[i]->addFrc(restraintFrc); + + // if the particle is directional, we accumulate the rot. restraints + if (vecParticles[i]->isDirectional()){ + + // get the current rotation matrix and reference vector vecParticles[i]->getA(A); - Calc_rVal( pos, i ); - Calc_body_thetaVal( A, i ); + + data = vecParticles[i]->getProperty("refVectorX"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refVectorX from StuntDouble\n"; + return 0.0; + } + else refVec[0] = doubleData->getData(); + } + data = vecParticles[i]->getProperty("refVectorY"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refVectorY from StuntDouble\n"; + return 0.0; + } + else refVec[1] = doubleData->getData(); + } + data = vecParticles[i]->getProperty("refVectorZ"); + if (data){ + doubleData = dynamic_cast(data); + if (!doubleData){ + cerr << "Can't obtain refVectorZ from StuntDouble\n"; + return 0.0; + } + else refVec[2] = doubleData->getData(); + } + + // calculate the theta and omega displacements + Calc_body_thetaVal( A, refVec ); omegaPass = vecParticles[i]->getZangle(); Calc_body_omegaVal( A, omegaPass ); - // first we calculate the derivatives - dVdrx = -kDist*delRx; - dVdry = -kDist*delRy; - dVdrz = -kDist*delRz; - // uTx... and vTx... are the body-fixed z and y unit vectors uTx = 0.0; uTy = 0.0; @@ -223,12 +309,12 @@ double Restraints::Calc_Restraint_Forces(vector tolerance) { dVdux = -(kTheta*theta/sin(theta))*ub0x; @@ -242,12 +328,7 @@ double Restraints::Calc_Restraint_Forces(vectoraddFrc(restraintFrc); + // now pass this temporary torque vector to the total torque vecParticles[i]->addTrq(spaceTrq); } - } - // and we can return the appropriately scaled potential energy + // update the total harmonic potential with this object's contribution + harmPotent += tempPotent; + } + + // we can finish by returning the appropriately scaled potential energy tempPotent = harmPotent * factor; return tempPotent; } -void Restraints::Store_Init_Info(vector vecParticles){ - int idealSize; - double pos[3]; - double A[3][3]; - double RfromQ[3][3]; - double quat0, quat1, quat2, quat3; - double dot; - vector tempZangs; - const char *delimit = " \t\n;,"; +void Restraints::Write_zAngle_File(vector vecParticles, + int currTime, + int nIntObj){ - //open the idealCrystal.in file and zAngle.ang file - strcpy(fileName, "idealCrystal.in"); - strcpy(angleName, "zAngle.ang"); - - ifstrstream crystalIn(fileName); - ifstrstream angleIn(angleName); + char zOutName[200]; - // check to see if these files are present in the execution directory - if (!crystalIn) { - sprintf(painCave.errMsg, - "Restraints Error: Unable to open idealCrystal.in for reading.\n" - "\tMake sure a ref. crystal file is in the working directory.\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } + std::cerr << nIntObj << " is the number of integrable objects\n"; - // it's not fatal to lack a zAngle.ang file, it just means you're starting - // from the ideal crystal state - if (!angleIn) { - sprintf(painCave.errMsg, - "Restraints Warning: The lack of a zAngle.ang file is mildly\n" - "\tunsettling... This means the simulation is starting from the\n" - "\tidealCrystal.in reference configuration, so the omega values\n" - "\twill all be set to zero. If this is not the case, the energy\n" - "\tcalculations will be wrong.\n"); - painCave.severity = OOPSE_WARNING; - painCave.isFatal = 0; - simError(); - } - - // A rather specific reader for OOPSE .eor files... - // Let's read in the perfect crystal file - crystalIn.getline(inLine,999,'\n'); - // check to see if the crystal file is the same length as starting config. - token = strtok(inLine,delimit); - strcpy(inValue,token); - idealSize = atoi(inValue); - if (idealSize != vecParticles.size()) { - sprintf(painCave.errMsg, - "Restraints Error: Reference crystal file is not valid.\n" - "\tMake sure the idealCrystal.in file is the same size as the\n" - "\tstarting configuration. Using an incompatable crystal will\n" - "\tlead to energy calculation failures.\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } - // else, the file is okay... let's continue - crystalIn.getline(inLine,999,'\n'); + //#ifndef IS_MPI - for (i=0; isetZangle(tempZangs[i]); - else { - sprintf(painCave.errMsg, - "Restraints Error: the supplied zAngle file is not valid.\n" - "\tMake sure the zAngle.ang file matches with the initial\n" - "\tconfiguration (i.e. they're the same length). Using the wrong\n" - "\tzAngle file will lead to errors in the energy calculations.\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } - } - angleIn.close(); - - return; -} - -void Restraints::Write_zAngle_File(vector vecParticles){ - - char zOutName[200]; - strcpy(zOutName,"zAngle.ang"); - + ofstream angleOut(zOutName); - angleOut << "This file contains the omega values for the .eor file\n"; + angleOut << currTime << ": omega values at this time\n"; for (i=0; igetZangle() << "\n"; } + return; }