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;
  double omega;

  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 you arestarting from the idealCrystal.in\n"
            "\treference configuration, so the omega values will all be set to\n"
            "\tzero. If this isn't the case, you should question 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:	1189
Committed:	Sat May 22 20:55:04 2004 UTC (20 years, 3 months ago) by chrisfen
File size:	10543 byte(s)
Log Message:	Fixed an error in Restraints.cpp... Too many arguements in a function call.
#	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	1187	double omega;
145	chrisfen	1180
146			tolerance = 5.72957795131e-7;
147
148			harmPotent = 0.0; // zero out the global harmonic potential variable
149
150			factor = 1 - pow(lambdaValue, lambdaK);
151
152	chrisfen	1187	for (i=0; i<nMol; i++){
153	chrisfen	1180	if (vecParticles[i]->isDirectional()){
154			vecParticles[i]->getPos(pos);
155			vecParticles[i]->getA(A);
156			Calc_rVal( pos, i );
157			Calc_body_thetaVal( A, i );
158	chrisfen	1189	omegaPass = vecParticles[i]->getZangle();
159			Calc_body_omegaVal( A, omegaPass );
160	chrisfen	1180
161			if (omega > PI \|\| omega < -PI)
162			cout << "oops... " << omega << "\n";
163
164			// first we calculate the derivatives
165			dVdrx = -kDist*delRx;
166			dVdry = -kDist*delRy;
167			dVdrz = -kDist*delRz;
168
169			// uTx... and vTx... are the body-fixed z and y unit vectors
170			uTx = 0.0;
171			uTy = 0.0;
172			uTz = 1.0;
173			vTx = 0.0;
174			vTy = 1.0;
175			vTz = 0.0;
176
177			dVdux = 0;
178			dVduy = 0;
179			dVduz = 0;
180			dVdvx = 0;
181			dVdvy = 0;
182			dVdvz = 0;
183
184			if (fabs(theta) > tolerance) {
185			dVdux = -(kThetatheta/sin(theta))ub0x;
186			dVduy = -(kThetatheta/sin(theta))ub0y;
187			dVduz = -(kThetatheta/sin(theta))ub0z;
188			}
189
190			if (fabs(omega) > tolerance) {
191			dVdvx = -(kOmegaomega/sin(omega))vb0x;
192			dVdvy = -(kOmegaomega/sin(omega))vb0y;
193			dVdvz = -(kOmegaomega/sin(omega))vb0z;
194			}
195
196			// next we calculate the restraint forces and torques
197			restraintFrc[0] = dVdrx;
198			restraintFrc[1] = dVdry;
199			restraintFrc[2] = dVdrz;
200			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
201
202			restraintTrq[0] = 0.0;
203			restraintTrq[1] = 0.0;
204			restraintTrq[2] = 0.0;
205
206			if (fabs(omega) > tolerance) {
207			restraintTrq[0] += 0.0;
208			restraintTrq[1] += 0.0;
209			restraintTrq[2] += vTy*dVdvx;
210			tempPotent += 0.5(kOmegaomega*omega);
211			}
212			if (fabs(theta) > tolerance) {
213			restraintTrq[0] += (uTz*dVduy);
214			restraintTrq[1] += -(uTz*dVdux);
215			restraintTrq[2] += 0.0;
216			tempPotent += 0.5(kThetatheta*theta);
217			}
218
219			for (j = 0; j < 3; j++) {
220			restraintFrc[j] *= factor;
221			restraintTrq[j] *= factor;
222			}
223
224			harmPotent += tempPotent;
225
226			// now we need to convert from body-fixed torques to space-fixed torques
227			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
228			+ A[2][0]*restraintTrq[2];
229			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
230			+ A[2][1]*restraintTrq[2];
231			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
232			+ A[2][2]*restraintTrq[2];
233
234			// now it's time to pass these temporary forces and torques
235			// to the total forces and torques
236			vecParticles[i]->addFrc(restraintFrc);
237			vecParticles[i]->addTrq(spaceTrq);
238			}
239			}
240
241			// and we can return the appropriately scaled potential energy
242			tempPotent = harmPotent * factor;
243			return tempPotent;
244			}
245
246	chrisfen	1187	void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
247	chrisfen	1180	double pos[3];
248			double A[3][3];
249			double RfromQ[3][3];
250			double quat0, quat1, quat2, quat3;
251			double dot;
252			// char *token;
253			// char fileName[200];
254			// char angleName[200];
255			// char inLine[1000];
256			// char inValue[200];
257			const char *delimit = " \t\n;,";
258
259			//open the idealCrystal.in file and zAngle.ang file
260			strcpy(fileName, "idealCrystal.in");
261			strcpy(angleName, "zAngle.ang");
262
263			ifstream crystalIn(fileName);
264			ifstream angleIn(angleName);
265
266			if (!crystalIn) {
267	chrisfen	1187	sprintf(painCave.errMsg,
268			"Restraints Error: Unable to open idealCrystal.in for reading.\n"
269			"\tMake sure a reference crystal file is in the current directory.\n");
270			painCave.isFatal = 1;
271			simError();
272
273	chrisfen	1180	return;
274			}
275
276			if (!angleIn) {
277	chrisfen	1187	sprintf(painCave.errMsg,
278			"Restraints Warning: The lack of a zAngle.ang file is mildly\n"
279			"\tunsettling... This means you arestarting from the idealCrystal.in\n"
280			"\treference configuration, so the omega values will all be set to\n"
281			"\tzero. If this isn't the case, you should question your results.\n");
282			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