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;,";

  strcpy(springName, "HarmSpringConsts.txt");

  ifstream springs(springName);

  if (!springs) { 
    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 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
    return;
  }

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

  cout << "The Spring Constants are:\n\tkDist  = " << kDist << "\n\tkTheta = " << kTheta << "\n\tkOmega = " << kOmega << "\n";
}

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