--- trunk/OOPSE/libmdtools/Restraints.cpp 2004/05/20 20:24:07 1180 +++ trunk/OOPSE/libmdtools/Restraints.cpp 2004/05/24 21:23:36 1194 @@ -28,9 +28,15 @@ Restraints::Restraints(int nMolInfo, double lambdaVal, ifstream springs(springName); if (!springs) { - cout << "Unable to open HarmSpringConsts.txt for reading.\n"; + sprintf(painCave.errMsg, + "Restraints Warning: Unable to open HarmSpringConsts.txt for reading.\n" + "\tDefault spring constants will be loaded. If you want to specify\n" + "\tspring constants, include a three line HarmSpringConsts.txt file\n" + "\tin the current directory.\n"); + painCave.isFatal = 0; + simError(); - // load place holder spring constants + // load default spring constants kDist = 6; // spring constant in units of kcal/(mol*ang^2) kTheta = 7.5; // in units of kcal/mol kOmega = 13.5; // in units of kcal/mol @@ -55,7 +61,7 @@ Restraints::Restraints(int nMolInfo, double lambdaVal, kOmega = (atof(inValue)); springs.close(); - cout << "Spring Constants: " << kDist << "\t" << kTheta << "\t" << kOmega << "\n"; + cout << "The Spring Constants are:\n\tkDist = " << kDist << "\n\tkTheta = " << kTheta << "\n\tkOmega = " << kOmega << "\n"; } Restraints::~Restraints(){ @@ -69,22 +75,6 @@ void Restraints::Calc_rVal(double position[3], int cur return; } -void Restraints::Calc_thetaVal(double matrix[3][3], int currentMol){ - uTx = matrix[2][0]; - uTy = matrix[2][1]; - uTz = matrix[2][2]; - - normalize = sqrt(uTx*uTx + uTy*uTy + uTz*uTz); - uTx = uTx/normalize; - uTy = uTy/normalize; - uTz = uTz/normalize; - - // Theta is the dot product of the reference and new z-axes - theta = acos(uTx*uX0[currentMol]+uTy*uY0[currentMol]+uTz*uZ0[currentMol]); - - 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]; @@ -104,54 +94,12 @@ void Restraints::Calc_body_thetaVal(double matrix[3][3 return; } -void Restraints::Calc_omegaVal(double matrix[3][3], int currentMol){ - double dot; - - uTx = matrix[2][0]; - uTy = matrix[2][1]; - uTz = matrix[2][2]; - vTx = matrix[1][0]; - vTy = matrix[1][1]; - vTz = matrix[1][2]; - - normalize = sqrt(uTx*uTx + uTy*uTy + uTz*uTz); - uTx = uTx/normalize; - uTy = uTy/normalize; - uTz = uTz/normalize; - - normalize = sqrt(vTx*vTx + vTy*vTy + vTz*vTz); - vTx = vTx/normalize; - vTy = vTy/normalize; - vTz = vTz/normalize; - - dot = uTx * vX0[currentMol] + uTy * vY0[currentMol] + uTz * vZ0[currentMol]; - - // Find the original y-axis vector projection on the current - // space-fixed xy-plane - vProj0[0] = vX0[currentMol] - dot * uTx; - vProj0[1] = vY0[currentMol] - dot * uTy; - vProj0[2] = vZ0[currentMol] - dot * uTz; - - // Convert the projection to a unit vector - vProjDist = sqrt(vProj0[0]*vProj0[0] + vProj0[1]*vProj0[1] - + vProj0[2]*vProj0[2]); - vProj0[0] = vProj0[0]/vProjDist; - vProj0[1] = vProj0[1]/vProjDist; - vProj0[2] = vProj0[2]/vProjDist; - - // Omega is the dot product of the new y-axis and the projection - // of the reference y-axis on the current xy-plane - omega = acos(vTx*vProj0[0] + vTy*vProj0[1] + vTz*vProj0[2]); - - return; -} - -void Restraints::Calc_body_omegaVal(double matrix[3][3], int currentMol){ +void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){ double zRotator[3][3]; double tempOmega; double wholeTwoPis; // Use the omega accumulated from the rotation propagation - omega = zAngle[currentMol]; + omega = zAngle; // translate the omega into a range between -PI and PI if (omega < -PI){ @@ -192,12 +140,7 @@ double Restraints::Calc_Restraint_Forces(vectorisDirectional()){ vecParticles[i]->getPos(pos); vecParticles[i]->getA(A); Calc_rVal( pos, i ); Calc_body_thetaVal( A, i ); - Calc_body_omegaVal( A, i ); + omegaPass = vecParticles[i]->getZangle(); + Calc_body_omegaVal( A, omegaPass ); if (omega > PI || omega < -PI) cout << "oops... " << omega << "\n"; @@ -298,7 +242,7 @@ double Restraints::Calc_Restraint_Forces(vector vecParticles){ double pos[3]; double A[3][3]; double RfromQ[3][3]; @@ -319,13 +263,24 @@ void Restraints::Store_Init_Info(){ ifstream angleIn(angleName); if (!crystalIn) { - cout << "Unable to open idealCrystal.in for reading.\n"; + sprintf(painCave.errMsg, + "Restraints Error: Unable to open idealCrystal.in for reading.\n" + "\tMake sure a reference crystal file is in the current directory.\n"); + painCave.isFatal = 1; + simError(); + return; } if (!angleIn) { - cout << "Unable to open zAngle.ang for reading.\n"; - cout << "The omega values are all assumed to be zero.\n"; + 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, you should\n" + "\tquestion your results.\n"); + painCave.isFatal = 0; + simError(); } // A rather specific reader for OOPSE .eor files... @@ -392,62 +347,24 @@ void Restraints::Store_Init_Info(){ angleIn.getline(inLine,999,'\n'); token = strtok(inLine,delimit); strcpy(inValue,token); - zAngle[i] = (atof(inValue)); + vecParticles[i]->setZangle(atof(inValue)); } } return; } -void Restraints::Determine_Lambda(){ -// double tempEps; +void Restraints::Write_zAngle_File(vector vecParticles){ -// atoms = entry_plug->atoms; - -// if (!strcmp(atoms[0]->getType(),"SSD") || -// !strcmp(atoms[0]->getType(),"SSD_E") || -// !strcmp(atoms[0]->getType(),"SSD_RF") || -// !strcmp(atoms[0]->getType(),"SSD1")){ - -// tempEps = atoms[0]->getEpsilon(); -// scaleLam = 1.0 - (tempEps/0.152); -// } -// else if (!strcmp(atoms[0]->getType(),"O_TIP3P")){ -// tempEps = atoms[0]->getEpsilon(); -// scaleLam = 1.0 - (tempEps/0.1521); -// } -// else if (!strcmp(atoms[0]->getType(),"O_TIP4P")){ -// tempEps = atoms[0]->getEpsilon(); -// scaleLam = 1.0 - (tempEps/0.1550); -// } -// else if (!strcmp(atoms[0]->getType(),"O_TIP5P")){ -// tempEps = atoms[0]->getEpsilon(); -// scaleLam = 1.0 - (tempEps/0.16); -// } -// else if (!strcmp(atoms[0]->getType(),"O_SPCE")){ -// tempEps = atoms[0]->getEpsilon(); -// scaleLam = 1.0 - (tempEps/0.15532); -// } -// else -// sprintf( painCave.errMsg, -// "Error in setting the lambda scale: Restraints\n" ); - -// if (fabs(scaleLam < 1e-9)) -// scaleLam = 0.0; -// cout << "The scaleLam is " << scaleLam << "\n"; -} - -void Restraints::Write_zAngle_File(){ - char zOutName[200]; strcpy(zOutName,"zAngle.ang"); ofstream angleOut(zOutName); angleOut << "This file contains the omega values for the .eor file\n"; - for (i=0; igetZangle() << "\n"; + } return; }