src/restraints/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(double lambdaVal, double lambdaExp){
  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<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 );

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

  // now we can read in the zAngle.ang file
  if (angleIn){
    angleIn.getline(inLine,999,'\n');
    for (i=0; i<vecParticles.size(); 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<vecParticles.size(); i++) {
    angleOut << vecParticles[i]->getZangle() << "\n";
  }
  return;
}

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

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