src/restraints/Restraints.cpp

// Thermodynamic integration is not multiprocessor friendly right now

#include <iostream>
#include <stdlib.h>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <string>
#include <cstring>
#include <math.h>

using namespace std;

#include "restraints/Restraints.hpp"
#include "brains/SimInfo.hpp"
#include "utils/simError.h"
#include "io/basic_ifstrstream.hpp"

#define PI 3.14159265359
#define TWO_PI 6.28318530718

Restraints::Restraints(double lambdaVal, double lambdaExp){
  lambdaValue = lambdaVal;
  lambdaK = lambdaExp;
  vector<double> resConsts;
  const char *jolt = " \t\n;,";

#ifdef IS_MPI
  if(worldRank == 0 ){
#endif // is_mpi

    strcpy(springName, "HarmSpringConsts.txt");
    
    ifstream springs(springName);
    
    if (!springs) { 
      sprintf(painCave.errMsg,
              "Unable to open HarmSpringConsts.txt for reading, so the\n"
              "\tdefault spring constants will be loaded. If you want\n"
              "\tto specify spring constants, include a three line\n"
              "\tHarmSpringConsts.txt file in the execution directory.\n");
      painCave.severity = OOPSE_WARNING;
      painCave.isFatal = 0;
      simError();   
      
      // 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
    } else  {
      
      springs.getline(inLine,999,'\n');
      // the file is blank!
      if (springs.eof()){
      sprintf(painCave.errMsg,
              "HarmSpringConsts.txt file is not valid.\n"
              "\tThe file should contain four rows, the last three containing\n"
              "\ta label and the spring constant value. They should be listed\n"
              "\tin the following order: kDist (positional restrant), kTheta\n"
              "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
              "\trestraint: rotation about the z-axis).\n");
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();   
      }
      // read in spring constants and check to make sure it is a valid file
      springs.getline(inLine,999,'\n');
      while (!springs.eof()){
        if (NULL != inLine){
          token = strtok(inLine,jolt);
          token = strtok(NULL,jolt);
          if (NULL != token){
            strcpy(inValue,token);
            resConsts.push_back(atof(inValue));
          }
        }
        springs.getline(inLine,999,'\n');
      }
      if (resConsts.size() == 3){
        kDist = resConsts[0];
        kTheta = resConsts[1];
        kOmega = resConsts[2];
      }
      else {
        sprintf(painCave.errMsg,
                "HarmSpringConsts.txt file is not valid.\n"
                "\tThe file should contain four rows, the last three containing\n"
                "\ta label and the spring constant value. They should be listed\n"
                "\tin the following order: kDist (positional restrant), kTheta\n"
                "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
                "\trestraint: rotation about the z-axis).\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }
    }
#ifdef IS_MPI
  }
  
  MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
  
  sprintf( checkPointMsg,
           "Sucessfully opened and read spring file.\n");
  MPIcheckPoint();

#endif // is_mpi
  
  sprintf(painCave.errMsg,
          "The spring constants for thermodynamic integration are:\n"
          "\tkDist = %lf\n"
          "\tkTheta = %lf\n"
          "\tkOmega = %lf\n", kDist, kTheta, kOmega);
  painCave.severity = OOPSE_INFO;
  painCave.isFatal = 0;
  simError();   
}

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];

  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];

  normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z);
  ub0x = ub0x/normalize;
  ub0y = ub0y/normalize;
  ub0z = ub0z/normalize;

  // Theta is the dot product of the reference and new z-axes
  theta = acos(ub0z);

  return;
}

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;

  // translate the omega into a range between -PI and PI
  if (omega < -PI){
    tempOmega = omega / -TWO_PI;
    wholeTwoPis = floor(tempOmega);
    tempOmega = omega + TWO_PI*wholeTwoPis;
    if (tempOmega < -PI)
      omega = tempOmega + TWO_PI;
    else
      omega = tempOmega;
  }
  if (omega > PI){
    tempOmega = omega / TWO_PI;
    wholeTwoPis = floor(tempOmega);
    tempOmega = omega - TWO_PI*wholeTwoPis;
    if (tempOmega > PI)
      omega = tempOmega - TWO_PI;   
    else
      omega = tempOmega;
  }

  vb0x = sin(omega);
  vb0y = cos(omega);
  vb0z = 0.0;

  normalize = sqrt(vb0x*vb0x + vb0y*vb0y + vb0z*vb0z);
  vb0x = vb0x/normalize;
  vb0y = vb0y/normalize;
  vb0z = vb0z/normalize;

  return;
}

double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
  double pos[3];
  double A[3][3];
  double tolerance;
  double tempPotent;
  double factor;
  double spaceTrq[3];
  double omegaPass;

  tolerance = 5.72957795131e-7;

  harmPotent = 0.0;  // zero out the global harmonic potential variable

  factor = 1 - pow(lambdaValue, lambdaK);

  for (i=0; i<vecParticles.size(); i++){
    if (vecParticles[i]->isDirectional()){
      vecParticles[i]->getPos(pos);
      vecParticles[i]->getA(A);
      Calc_rVal( pos, i );
      Calc_body_thetaVal( A, i );
      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;
      uTz = 1.0;
      vTx = 0.0;
      vTy = 1.0;
      vTz = 0.0;

      dVdux = 0;
      dVduy = 0;
      dVduz = 0;
      dVdvx = 0;
      dVdvy = 0;
      dVdvz = 0;

      if (fabs(theta) > tolerance) {
        dVdux = -(kTheta*theta/sin(theta))*ub0x;
        dVduy = -(kTheta*theta/sin(theta))*ub0y;
        dVduz = -(kTheta*theta/sin(theta))*ub0z;
      }

      if (fabs(omega) > tolerance) {
        dVdvx = -(kOmega*omega/sin(omega))*vb0x;
        dVdvy = -(kOmega*omega/sin(omega))*vb0y;
        dVdvz = -(kOmega*omega/sin(omega))*vb0z;
      }

      // next we calculate the restraint forces and torques
      restraintFrc[0] = dVdrx;
      restraintFrc[1] = dVdry;
      restraintFrc[2] = dVdrz;
      tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);

      restraintTrq[0] = 0.0;
      restraintTrq[1] = 0.0;
      restraintTrq[2] = 0.0;

      if (fabs(omega) > tolerance) {
        restraintTrq[0] += 0.0;
        restraintTrq[1] += 0.0;
        restraintTrq[2] += vTy*dVdvx;
        tempPotent += 0.5*(kOmega*omega*omega);
      }
      if (fabs(theta) > tolerance) {
        restraintTrq[0] += (uTz*dVduy);
        restraintTrq[1] += -(uTz*dVdux);
        restraintTrq[2] += 0.0;
        tempPotent += 0.5*(kTheta*theta*theta);
      }

      for (j = 0; j < 3; j++) {
        restraintFrc[j] *= factor;
        restraintTrq[j] *= factor;
      }

      harmPotent += tempPotent;

      // now we need to convert from body-fixed torques to space-fixed torques
      spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1] 
        + A[2][0]*restraintTrq[2];
      spaceTrq[1] = A[0][1]*restraintTrq[0] + A[1][1]*restraintTrq[1] 
        + A[2][1]*restraintTrq[2];
      spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1] 
        + A[2][2]*restraintTrq[2];

      // now it's time to pass these temporary forces and torques
      // to the total forces and torques
      vecParticles[i]->addFrc(restraintFrc);
      vecParticles[i]->addTrq(spaceTrq);
    }
  }

  // and we can return the appropriately scaled potential energy
  tempPotent = harmPotent * factor;
  return tempPotent;
}

void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
  int idealSize;
  double pos[3];
  double A[3][3];
  double RfromQ[3][3];
  double quat0, quat1, quat2, quat3;
  double dot;
  vector<double> tempZangs;
  const char *delimit = " \t\n;,";

  //open the idealCrystal.in file and zAngle.ang file
  strcpy(fileName, "idealCrystal.in");
  strcpy(angleName, "zAngle.ang");
  
  ifstrstream crystalIn(fileName);
  ifstrstream angleIn(angleName);

  // 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();   
  }

  // 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');
  
  for (i=0; i<vecParticles.size(); i++) {
    crystalIn.getline(inLine,999,'\n');
    token = strtok(inLine,delimit);
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    cofmPosX.push_back(atof(inValue));
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    cofmPosY.push_back(atof(inValue));
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    cofmPosZ.push_back(atof(inValue));
    token = strtok(NULL,delimit);
    token = strtok(NULL,delimit);
    token = strtok(NULL,delimit);
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    quat0 = atof(inValue);
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    quat1 = atof(inValue);
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    quat2 = atof(inValue);
    token = strtok(NULL,delimit);
    strcpy(inValue,token);
    quat3 = atof(inValue);

    // now build the rotation matrix and find the unit vectors
    RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;
    RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);
    RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);
    RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);
    RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;
    RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);
    RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);
    RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);
    RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;
    
    normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] 
                     + RfromQ[2][2]*RfromQ[2][2]);
    uX0.push_back(RfromQ[2][0]/normalize);
    uY0.push_back(RfromQ[2][1]/normalize);
    uZ0.push_back(RfromQ[2][2]/normalize);

    normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]
                     + RfromQ[1][2]*RfromQ[1][2]);
    vX0.push_back(RfromQ[1][0]/normalize);
    vY0.push_back(RfromQ[1][1]/normalize);
    vZ0.push_back(RfromQ[1][2]/normalize);
  }
  crystalIn.close();

  // now we read in the zAngle.ang file
  if (angleIn){
    angleIn.getline(inLine,999,'\n');
    angleIn.getline(inLine,999,'\n');
    while (!angleIn.eof()) {
      token = strtok(inLine,delimit);
      strcpy(inValue,token);
      tempZangs.push_back(atof(inValue));
      angleIn.getline(inLine,999,'\n');
    }

    // test to make sure the zAngle.ang file is the proper length
    if (tempZangs.size() == vecParticles.size())
      for (i=0; i<vecParticles.size(); i++)
        vecParticles[i]->setZangle(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<StuntDouble*> vecParticles){

  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; i<vecParticles.size(); i++) {
    angleOut << vecParticles[i]->getZangle() << "\n";
  }
  return;
}

double Restraints::getVharm(){
  return harmPotent;
}

Revision:	1731
Committed:	Thu Nov 11 21:47:25 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	13800 byte(s)
Log Message:	working towards parallel restraints
#	User	Rev	Content
1	chrisfen	1534	// Thermodynamic integration is not multiprocessor friendly right now
2
3	gezelter	1490	#include <iostream>
4			#include <stdlib.h>
5			#include <cstdio>
6			#include <fstream>
7			#include <iomanip>
8			#include <string>
9			#include <cstring>
10			#include <math.h>
11
12			using namespace std;
13
14	gezelter	1537	#include "restraints/Restraints.hpp"
15			#include "brains/SimInfo.hpp"
16			#include "utils/simError.h"
17	chrisfen	1731	#include "io/basic_ifstrstream.hpp"
18	gezelter	1490
19			#define PI 3.14159265359
20			#define TWO_PI 6.28318530718
21
22			Restraints::Restraints(double lambdaVal, double lambdaExp){
23			lambdaValue = lambdaVal;
24			lambdaK = lambdaExp;
25	chrisfen	1534	vector<double> resConsts;
26	gezelter	1490	const char *jolt = " \t\n;,";
27
28			#ifdef IS_MPI
29			if(worldRank == 0 ){
30			#endif // is_mpi
31
32			strcpy(springName, "HarmSpringConsts.txt");
33
34			ifstream springs(springName);
35
36			if (!springs) {
37			sprintf(painCave.errMsg,
38	chrisfen	1534	"Unable to open HarmSpringConsts.txt for reading, so the\n"
39			"\tdefault spring constants will be loaded. If you want\n"
40			"\tto specify spring constants, include a three line\n"
41			"\tHarmSpringConsts.txt file in the execution directory.\n");
42	gezelter	1490	painCave.severity = OOPSE_WARNING;
43			painCave.isFatal = 0;
44			simError();
45
46			// load default spring constants
47			kDist = 6; // spring constant in units of kcal/(mol*ang^2)
48			kTheta = 7.5; // in units of kcal/mol
49			kOmega = 13.5; // in units of kcal/mol
50			} else {
51
52			springs.getline(inLine,999,'\n');
53	chrisfen	1534	// the file is blank!
54			if (springs.eof()){
55			sprintf(painCave.errMsg,
56			"HarmSpringConsts.txt file is not valid.\n"
57			"\tThe file should contain four rows, the last three containing\n"
58			"\ta label and the spring constant value. They should be listed\n"
59			"\tin the following order: kDist (positional restrant), kTheta\n"
60			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
61			"\trestraint: rotation about the z-axis).\n");
62			painCave.severity = OOPSE_ERROR;
63			painCave.isFatal = 1;
64			simError();
65			}
66			// read in spring constants and check to make sure it is a valid file
67	gezelter	1490	springs.getline(inLine,999,'\n');
68	chrisfen	1534	while (!springs.eof()){
69			if (NULL != inLine){
70			token = strtok(inLine,jolt);
71			token = strtok(NULL,jolt);
72			if (NULL != token){
73			strcpy(inValue,token);
74			resConsts.push_back(atof(inValue));
75			}
76			}
77			springs.getline(inLine,999,'\n');
78			}
79			if (resConsts.size() == 3){
80			kDist = resConsts[0];
81			kTheta = resConsts[1];
82			kOmega = resConsts[2];
83			}
84			else {
85			sprintf(painCave.errMsg,
86			"HarmSpringConsts.txt file is not valid.\n"
87			"\tThe file should contain four rows, the last three containing\n"
88			"\ta label and the spring constant value. They should be listed\n"
89			"\tin the following order: kDist (positional restrant), kTheta\n"
90			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
91			"\trestraint: rotation about the z-axis).\n");
92			painCave.severity = OOPSE_ERROR;
93			painCave.isFatal = 1;
94			simError();
95			}
96	gezelter	1490	}
97			#ifdef IS_MPI
98			}
99
100			MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
101			MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
102			MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
103
104			sprintf( checkPointMsg,
105			"Sucessfully opened and read spring file.\n");
106			MPIcheckPoint();
107
108			#endif // is_mpi
109
110			sprintf(painCave.errMsg,
111			"The spring constants for thermodynamic integration are:\n"
112			"\tkDist = %lf\n"
113			"\tkTheta = %lf\n"
114			"\tkOmega = %lf\n", kDist, kTheta, kOmega);
115			painCave.severity = OOPSE_INFO;
116			painCave.isFatal = 0;
117			simError();
118			}
119
120			Restraints::~Restraints(){
121			}
122
123			void Restraints::Calc_rVal(double position[3], int currentMol){
124			delRx = position[0] - cofmPosX[currentMol];
125			delRy = position[1] - cofmPosY[currentMol];
126			delRz = position[2] - cofmPosZ[currentMol];
127
128			return;
129			}
130
131			void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
132			ub0x = matrix[0][0]uX0[currentMol] + matrix[0][1]uY0[currentMol]
133			+ matrix[0][2]*uZ0[currentMol];
134			ub0y = matrix[1][0]uX0[currentMol] + matrix[1][1]uY0[currentMol]
135			+ matrix[1][2]*uZ0[currentMol];
136			ub0z = matrix[2][0]uX0[currentMol] + matrix[2][1]uY0[currentMol]
137			+ matrix[2][2]*uZ0[currentMol];
138
139			normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
140			ub0x = ub0x/normalize;
141			ub0y = ub0y/normalize;
142			ub0z = ub0z/normalize;
143
144			// Theta is the dot product of the reference and new z-axes
145			theta = acos(ub0z);
146
147			return;
148			}
149
150			void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){
151			double zRotator[3][3];
152			double tempOmega;
153			double wholeTwoPis;
154			// Use the omega accumulated from the rotation propagation
155			omega = zAngle;
156
157			// translate the omega into a range between -PI and PI
158			if (omega < -PI){
159			tempOmega = omega / -TWO_PI;
160			wholeTwoPis = floor(tempOmega);
161			tempOmega = omega + TWO_PI*wholeTwoPis;
162			if (tempOmega < -PI)
163			omega = tempOmega + TWO_PI;
164			else
165			omega = tempOmega;
166			}
167			if (omega > PI){
168			tempOmega = omega / TWO_PI;
169			wholeTwoPis = floor(tempOmega);
170			tempOmega = omega - TWO_PI*wholeTwoPis;
171			if (tempOmega > PI)
172			omega = tempOmega - TWO_PI;
173			else
174			omega = tempOmega;
175			}
176
177			vb0x = sin(omega);
178			vb0y = cos(omega);
179			vb0z = 0.0;
180
181			normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
182			vb0x = vb0x/normalize;
183			vb0y = vb0y/normalize;
184			vb0z = vb0z/normalize;
185
186			return;
187			}
188
189			double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
190			double pos[3];
191			double A[3][3];
192			double tolerance;
193			double tempPotent;
194			double factor;
195			double spaceTrq[3];
196			double omegaPass;
197
198			tolerance = 5.72957795131e-7;
199
200			harmPotent = 0.0; // zero out the global harmonic potential variable
201
202			factor = 1 - pow(lambdaValue, lambdaK);
203
204			for (i=0; i<vecParticles.size(); i++){
205			if (vecParticles[i]->isDirectional()){
206			vecParticles[i]->getPos(pos);
207			vecParticles[i]->getA(A);
208			Calc_rVal( pos, i );
209			Calc_body_thetaVal( A, i );
210			omegaPass = vecParticles[i]->getZangle();
211			Calc_body_omegaVal( A, omegaPass );
212
213			// first we calculate the derivatives
214			dVdrx = -kDist*delRx;
215			dVdry = -kDist*delRy;
216			dVdrz = -kDist*delRz;
217
218			// uTx... and vTx... are the body-fixed z and y unit vectors
219			uTx = 0.0;
220			uTy = 0.0;
221			uTz = 1.0;
222			vTx = 0.0;
223			vTy = 1.0;
224			vTz = 0.0;
225
226			dVdux = 0;
227			dVduy = 0;
228			dVduz = 0;
229			dVdvx = 0;
230			dVdvy = 0;
231			dVdvz = 0;
232
233			if (fabs(theta) > tolerance) {
234			dVdux = -(kThetatheta/sin(theta))ub0x;
235			dVduy = -(kThetatheta/sin(theta))ub0y;
236			dVduz = -(kThetatheta/sin(theta))ub0z;
237			}
238
239			if (fabs(omega) > tolerance) {
240			dVdvx = -(kOmegaomega/sin(omega))vb0x;
241			dVdvy = -(kOmegaomega/sin(omega))vb0y;
242			dVdvz = -(kOmegaomega/sin(omega))vb0z;
243			}
244
245			// next we calculate the restraint forces and torques
246			restraintFrc[0] = dVdrx;
247			restraintFrc[1] = dVdry;
248			restraintFrc[2] = dVdrz;
249			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
250
251			restraintTrq[0] = 0.0;
252			restraintTrq[1] = 0.0;
253			restraintTrq[2] = 0.0;
254
255			if (fabs(omega) > tolerance) {
256			restraintTrq[0] += 0.0;
257			restraintTrq[1] += 0.0;
258			restraintTrq[2] += vTy*dVdvx;
259			tempPotent += 0.5(kOmegaomega*omega);
260			}
261			if (fabs(theta) > tolerance) {
262			restraintTrq[0] += (uTz*dVduy);
263			restraintTrq[1] += -(uTz*dVdux);
264			restraintTrq[2] += 0.0;
265			tempPotent += 0.5(kThetatheta*theta);
266			}
267
268			for (j = 0; j < 3; j++) {
269			restraintFrc[j] *= factor;
270			restraintTrq[j] *= factor;
271			}
272
273			harmPotent += tempPotent;
274
275			// now we need to convert from body-fixed torques to space-fixed torques
276			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
277			+ A[2][0]*restraintTrq[2];
278			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
279			+ A[2][1]*restraintTrq[2];
280			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
281			+ A[2][2]*restraintTrq[2];
282
283			// now it's time to pass these temporary forces and torques
284			// to the total forces and torques
285			vecParticles[i]->addFrc(restraintFrc);
286			vecParticles[i]->addTrq(spaceTrq);
287			}
288			}
289
290			// and we can return the appropriately scaled potential energy
291			tempPotent = harmPotent * factor;
292			return tempPotent;
293			}
294
295			void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
296	chrisfen	1534	int idealSize;
297	gezelter	1490	double pos[3];
298			double A[3][3];
299			double RfromQ[3][3];
300			double quat0, quat1, quat2, quat3;
301			double dot;
302	chrisfen	1534	vector<double> tempZangs;
303	gezelter	1490	const char *delimit = " \t\n;,";
304
305			//open the idealCrystal.in file and zAngle.ang file
306			strcpy(fileName, "idealCrystal.in");
307			strcpy(angleName, "zAngle.ang");
308
309	chrisfen	1731	ifstrstream crystalIn(fileName);
310			ifstrstream angleIn(angleName);
311	gezelter	1490
312	chrisfen	1534	// check to see if these files are present in the execution directory
313	gezelter	1490	if (!crystalIn) {
314			sprintf(painCave.errMsg,
315			"Restraints Error: Unable to open idealCrystal.in for reading.\n"
316	chrisfen	1534	"\tMake sure a ref. crystal file is in the working directory.\n");
317			painCave.severity = OOPSE_ERROR;
318	gezelter	1490	painCave.isFatal = 1;
319			simError();
320			}
321
322	chrisfen	1534	// it's not fatal to lack a zAngle.ang file, it just means you're starting
323			// from the ideal crystal state
324	gezelter	1490	if (!angleIn) {
325			sprintf(painCave.errMsg,
326			"Restraints Warning: The lack of a zAngle.ang file is mildly\n"
327			"\tunsettling... This means the simulation is starting from the\n"
328			"\tidealCrystal.in reference configuration, so the omega values\n"
329	chrisfen	1534	"\twill all be set to zero. If this is not the case, the energy\n"
330			"\tcalculations will be wrong.\n");
331			painCave.severity = OOPSE_WARNING;
332	gezelter	1490	painCave.isFatal = 0;
333			simError();
334			}
335
336			// A rather specific reader for OOPSE .eor files...
337			// Let's read in the perfect crystal file
338			crystalIn.getline(inLine,999,'\n');
339	chrisfen	1534	// check to see if the crystal file is the same length as starting config.
340			token = strtok(inLine,delimit);
341			strcpy(inValue,token);
342			idealSize = atoi(inValue);
343			if (idealSize != vecParticles.size()) {
344			sprintf(painCave.errMsg,
345			"Restraints Error: Reference crystal file is not valid.\n"
346			"\tMake sure the idealCrystal.in file is the same size as the\n"
347			"\tstarting configuration. Using an incompatable crystal will\n"
348			"\tlead to energy calculation failures.\n");
349			painCave.severity = OOPSE_ERROR;
350			painCave.isFatal = 1;
351			simError();
352			}
353			// else, the file is okay... let's continue
354	gezelter	1490	crystalIn.getline(inLine,999,'\n');
355
356			for (i=0; i<vecParticles.size(); i++) {
357			crystalIn.getline(inLine,999,'\n');
358			token = strtok(inLine,delimit);
359			token = strtok(NULL,delimit);
360			strcpy(inValue,token);
361			cofmPosX.push_back(atof(inValue));
362			token = strtok(NULL,delimit);
363			strcpy(inValue,token);
364			cofmPosY.push_back(atof(inValue));
365			token = strtok(NULL,delimit);
366			strcpy(inValue,token);
367			cofmPosZ.push_back(atof(inValue));
368			token = strtok(NULL,delimit);
369			token = strtok(NULL,delimit);
370			token = strtok(NULL,delimit);
371			token = strtok(NULL,delimit);
372			strcpy(inValue,token);
373			quat0 = atof(inValue);
374			token = strtok(NULL,delimit);
375			strcpy(inValue,token);
376			quat1 = atof(inValue);
377			token = strtok(NULL,delimit);
378			strcpy(inValue,token);
379			quat2 = atof(inValue);
380			token = strtok(NULL,delimit);
381			strcpy(inValue,token);
382			quat3 = atof(inValue);
383
384			// now build the rotation matrix and find the unit vectors
385			RfromQ[0][0] = quat0quat0 + quat1quat1 - quat2quat2 - quat3quat3;
386			RfromQ[0][1] = 2(quat1quat2 + quat0*quat3);
387			RfromQ[0][2] = 2(quat1quat3 - quat0*quat2);
388			RfromQ[1][0] = 2(quat1quat2 - quat0*quat3);
389			RfromQ[1][1] = quat0quat0 - quat1quat1 + quat2quat2 - quat3quat3;
390			RfromQ[1][2] = 2(quat2quat3 + quat0*quat1);
391			RfromQ[2][0] = 2(quat1quat3 + quat0*quat2);
392			RfromQ[2][1] = 2(quat2quat3 - quat0*quat1);
393			RfromQ[2][2] = quat0quat0 - quat1quat1 - quat2quat2 + quat3quat3;
394
395			normalize = sqrt(RfromQ[2][0]RfromQ[2][0] + RfromQ[2][1]RfromQ[2][1]
396			+ RfromQ[2][2]*RfromQ[2][2]);
397			uX0.push_back(RfromQ[2][0]/normalize);
398			uY0.push_back(RfromQ[2][1]/normalize);
399			uZ0.push_back(RfromQ[2][2]/normalize);
400
401			normalize = sqrt(RfromQ[1][0]RfromQ[1][0] + RfromQ[1][1]RfromQ[1][1]
402			+ RfromQ[1][2]*RfromQ[1][2]);
403			vX0.push_back(RfromQ[1][0]/normalize);
404			vY0.push_back(RfromQ[1][1]/normalize);
405			vZ0.push_back(RfromQ[1][2]/normalize);
406			}
407	chrisfen	1534	crystalIn.close();
408	gezelter	1490
409	chrisfen	1534	// now we read in the zAngle.ang file
410	gezelter	1490	if (angleIn){
411			angleIn.getline(inLine,999,'\n');
412	chrisfen	1534	angleIn.getline(inLine,999,'\n');
413			while (!angleIn.eof()) {
414	gezelter	1490	token = strtok(inLine,delimit);
415			strcpy(inValue,token);
416	chrisfen	1534	tempZangs.push_back(atof(inValue));
417			angleIn.getline(inLine,999,'\n');
418	gezelter	1490	}
419	chrisfen	1534
420			// test to make sure the zAngle.ang file is the proper length
421			if (tempZangs.size() == vecParticles.size())
422			for (i=0; i<vecParticles.size(); i++)
423			vecParticles[i]->setZangle(tempZangs[i]);
424			else {
425			sprintf(painCave.errMsg,
426			"Restraints Error: the supplied zAngle file is not valid.\n"
427			"\tMake sure the zAngle.ang file matches with the initial\n"
428			"\tconfiguration (i.e. they're the same length). Using the wrong\n"
429			"\tzAngle file will lead to errors in the energy calculations.\n");
430			painCave.severity = OOPSE_ERROR;
431			painCave.isFatal = 1;
432			simError();
433			}
434	gezelter	1490	}
435	chrisfen	1534	angleIn.close();
436	gezelter	1490
437			return;
438			}
439
440			void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
441
442			char zOutName[200];
443
444			strcpy(zOutName,"zAngle.ang");
445
446			ofstream angleOut(zOutName);
447			angleOut << "This file contains the omega values for the .eor file\n";
448			for (i=0; i<vecParticles.size(); i++) {
449			angleOut << vecParticles[i]->getZangle() << "\n";
450			}
451			return;
452			}
453
454			double Restraints::getVharm(){
455			return harmPotent;
456			}
457