OOPSE/libmdtools/Restraints.cpp

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

using namespace std;

#include "Restraints.hpp"
#include "SimInfo.hpp"
#include "simError.h"

#define PI 3.14159265359
#define TWO_PI 6.28318530718

Restraints::Restraints(int nMolInfo, double lambdaVal, double lambdaExp){
  nMol = nMolInfo;
  lambdaValue = lambdaVal;
  lambdaK = lambdaExp;

  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,
              "In Restraints: 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.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');
      springs.getline(inLine,999,'\n');
      token = strtok(inLine,jolt);
      token = strtok(NULL,jolt);
      strcpy(inValue,token);
      kDist = (atof(inValue));
      springs.getline(inLine,999,'\n');
      token = strtok(inLine,jolt);
      token = strtok(NULL,jolt);
      strcpy(inValue,token);
      kTheta = (atof(inValue));
      springs.getline(inLine,999,'\n');
      token = strtok(inLine,jolt);
      token = strtok(NULL,jolt);
      strcpy(inValue,token);
      kOmega = (atof(inValue));
      springs.close();
    }
#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<nMol; 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 );

      if (omega > PI || omega < -PI)
        cout << "oops... " << omega << "\n";

      // 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){
  double pos[3];
  double A[3][3];
  double RfromQ[3][3];
  double quat0, quat1, quat2, quat3;
  double dot;
//   char *token;
//   char fileName[200];
//   char angleName[200];
//   char inLine[1000];
//   char inValue[200];
  const char *delimit = " \t\n;,";

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

  if (!crystalIn) { 
    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) { 
    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...
  // Let's read in the perfect crystal file
  crystalIn.getline(inLine,999,'\n');
  crystalIn.getline(inLine,999,'\n');
  
  for (i=0; i<nMol; 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);
  }

  // now we can read in the zAngle.ang file
  if (angleIn){
    angleIn.getline(inLine,999,'\n');
    for (i=0; i<nMol; i++) {
      angleIn.getline(inLine,999,'\n');
      token = strtok(inLine,delimit);
      strcpy(inValue,token);
      vecParticles[i]->setZangle(atof(inValue));
    }
  }

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

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

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