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) { 
    cout << "Unable to open HarmSpringConsts.txt for reading.\n";

    // load place holder 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 << "Spring Constants: " << kDist << "\t" << kTheta << "\t" << 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_thetaVal(double matrix[3][3], int currentMol){
  uTx = matrix[2][0];
  uTy = matrix[2][1];
  uTz = matrix[2][2];

  normalize = sqrt(uTx*uTx + uTy*uTy + uTz*uTz);
  uTx = uTx/normalize;
  uTy = uTy/normalize;
  uTz = uTz/normalize;

  // Theta is the dot product of the reference and new z-axes
  theta = acos(uTx*uX0[currentMol]+uTy*uY0[currentMol]+uTz*uZ0[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_omegaVal(double matrix[3][3], int currentMol){
  double dot;

  uTx = matrix[2][0];
  uTy = matrix[2][1];
  uTz = matrix[2][2];
  vTx = matrix[1][0];
  vTy = matrix[1][1];
  vTz = matrix[1][2];

  normalize = sqrt(uTx*uTx + uTy*uTy + uTz*uTz);
  uTx = uTx/normalize;
  uTy = uTy/normalize;
  uTz = uTz/normalize;

  normalize = sqrt(vTx*vTx + vTy*vTy + vTz*vTz);
  vTx = vTx/normalize;
  vTy = vTy/normalize;
  vTz = vTz/normalize;

  dot = uTx * vX0[currentMol] + uTy * vY0[currentMol] + uTz * vZ0[currentMol];
  
  // Find the original y-axis vector projection on the current 
  // space-fixed xy-plane
  vProj0[0] = vX0[currentMol] - dot * uTx;
  vProj0[1] = vY0[currentMol] - dot * uTy;
  vProj0[2] = vZ0[currentMol] - dot * uTz;

  // Convert the projection to a unit vector
  vProjDist = sqrt(vProj0[0]*vProj0[0] + vProj0[1]*vProj0[1] 
                   + vProj0[2]*vProj0[2]);
  vProj0[0] = vProj0[0]/vProjDist;
  vProj0[1] = vProj0[1]/vProjDist;
  vProj0[2] = vProj0[2]/vProjDist;

  // Omega is the dot product of the new y-axis and the projection
  // of the reference y-axis on the current xy-plane
  omega = acos(vTx*vProj0[0] + vTy*vProj0[1] + vTz*vProj0[2]);

  return;
}

void Restraints::Calc_body_omegaVal(double matrix[3][3], int currentMol){
  double zRotator[3][3];
  double tempOmega;
  double wholeTwoPis;
  // Use the omega accumulated from the rotation propagation
  omega = zAngle[currentMol];

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

  //  atoms = atomPoint;

//   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

  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 );
      Calc_body_omegaVal( A, i );

      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(){
  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) { 
    cout << "Unable to open idealCrystal.in for reading.\n";
    return;
  }

  if (!angleIn) { 
    cout << "Unable to open zAngle.ang for reading.\n";
    cout << "The omega values are all assumed to be zero.\n";
  }

  // 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);
      zAngle[i] = (atof(inValue));
    }
  }

  return;
}

void Restraints::Determine_Lambda(){
//   double tempEps;

//   atoms = entry_plug->atoms;

//   if (!strcmp(atoms[0]->getType(),"SSD") || 
//       !strcmp(atoms[0]->getType(),"SSD_E") ||
//       !strcmp(atoms[0]->getType(),"SSD_RF") || 
//       !strcmp(atoms[0]->getType(),"SSD1")){
    
//     tempEps = atoms[0]->getEpsilon();
//     scaleLam = 1.0 - (tempEps/0.152); 
//   }
//   else if (!strcmp(atoms[0]->getType(),"O_TIP3P")){
//     tempEps = atoms[0]->getEpsilon();
//     scaleLam = 1.0 - (tempEps/0.1521);   
//   }
//   else if (!strcmp(atoms[0]->getType(),"O_TIP4P")){
//     tempEps = atoms[0]->getEpsilon();
//     scaleLam = 1.0 - (tempEps/0.1550);   
//   }
//   else if (!strcmp(atoms[0]->getType(),"O_TIP5P")){
//     tempEps = atoms[0]->getEpsilon();
//     scaleLam = 1.0 - (tempEps/0.16);   
//   }
//   else if (!strcmp(atoms[0]->getType(),"O_SPCE")){
//     tempEps = atoms[0]->getEpsilon();
//     scaleLam = 1.0 - (tempEps/0.15532);   
//   }
//   else
//     sprintf( painCave.errMsg,
//           "Error in setting the lambda scale: Restraints\n" );

//   if (fabs(scaleLam < 1e-9))
//       scaleLam = 0.0;
//   cout << "The scaleLam is " << scaleLam << "\n";
}

void Restraints::Write_zAngle_File(){

  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 << zAngle[i] << "\n";

  return;
}

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

Revision:	1180
Committed:	Thu May 20 20:24:07 2004 UTC (20 years, 1 month ago) by chrisfen
File size:	12621 byte(s)
Log Message:	Several additions... Restraints.cpp and .hpp were included for restraining particles in thermodynamic integration. By including these, changes were made in Integrator, SimInfo, ForceFeilds, SimSetup, StatWriter, and possibly some other files. Two bass keywords were also added for performing thermodynamic integration: a lambda value one and a k power one.
#	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			cout << "Unable to open HarmSpringConsts.txt for reading.\n";
32
33			// load place holder spring constants
34			kDist = 6; // spring constant in units of kcal/(mol*ang^2)
35			kTheta = 7.5; // in units of kcal/mol
36			kOmega = 13.5; // in units of kcal/mol
37			return;
38			}
39
40			springs.getline(inLine,999,'\n');
41			springs.getline(inLine,999,'\n');
42			token = strtok(inLine,jolt);
43			token = strtok(NULL,jolt);
44			strcpy(inValue,token);
45			kDist = (atof(inValue));
46			springs.getline(inLine,999,'\n');
47			token = strtok(inLine,jolt);
48			token = strtok(NULL,jolt);
49			strcpy(inValue,token);
50			kTheta = (atof(inValue));
51			springs.getline(inLine,999,'\n');
52			token = strtok(inLine,jolt);
53			token = strtok(NULL,jolt);
54			strcpy(inValue,token);
55			kOmega = (atof(inValue));
56			springs.close();
57
58			cout << "Spring Constants: " << kDist << "\t" << kTheta << "\t" << kOmega << "\n";
59			}
60
61			Restraints::~Restraints(){
62			}
63
64			void Restraints::Calc_rVal(double position[3], int currentMol){
65			delRx = position[0] - cofmPosX[currentMol];
66			delRy = position[1] - cofmPosY[currentMol];
67			delRz = position[2] - cofmPosZ[currentMol];
68
69			return;
70			}
71
72			void Restraints::Calc_thetaVal(double matrix[3][3], int currentMol){
73			uTx = matrix[2][0];
74			uTy = matrix[2][1];
75			uTz = matrix[2][2];
76
77			normalize = sqrt(uTxuTx + uTyuTy + uTz*uTz);
78			uTx = uTx/normalize;
79			uTy = uTy/normalize;
80			uTz = uTz/normalize;
81
82			// Theta is the dot product of the reference and new z-axes
83			theta = acos(uTxuX0[currentMol]+uTyuY0[currentMol]+uTz*uZ0[currentMol]);
84
85			return;
86			}
87
88			void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
89			ub0x = matrix[0][0]uX0[currentMol] + matrix[0][1]uY0[currentMol]
90			+ matrix[0][2]*uZ0[currentMol];
91			ub0y = matrix[1][0]uX0[currentMol] + matrix[1][1]uY0[currentMol]
92			+ matrix[1][2]*uZ0[currentMol];
93			ub0z = matrix[2][0]uX0[currentMol] + matrix[2][1]uY0[currentMol]
94			+ matrix[2][2]*uZ0[currentMol];
95
96			normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
97			ub0x = ub0x/normalize;
98			ub0y = ub0y/normalize;
99			ub0z = ub0z/normalize;
100
101			// Theta is the dot product of the reference and new z-axes
102			theta = acos(ub0z);
103
104			return;
105			}
106
107			void Restraints::Calc_omegaVal(double matrix[3][3], int currentMol){
108			double dot;
109
110			uTx = matrix[2][0];
111			uTy = matrix[2][1];
112			uTz = matrix[2][2];
113			vTx = matrix[1][0];
114			vTy = matrix[1][1];
115			vTz = matrix[1][2];
116
117			normalize = sqrt(uTxuTx + uTyuTy + uTz*uTz);
118			uTx = uTx/normalize;
119			uTy = uTy/normalize;
120			uTz = uTz/normalize;
121
122			normalize = sqrt(vTxvTx + vTyvTy + vTz*vTz);
123			vTx = vTx/normalize;
124			vTy = vTy/normalize;
125			vTz = vTz/normalize;
126
127			dot = uTx * vX0[currentMol] + uTy * vY0[currentMol] + uTz * vZ0[currentMol];
128
129			// Find the original y-axis vector projection on the current
130			// space-fixed xy-plane
131			vProj0[0] = vX0[currentMol] - dot * uTx;
132			vProj0[1] = vY0[currentMol] - dot * uTy;
133			vProj0[2] = vZ0[currentMol] - dot * uTz;
134
135			// Convert the projection to a unit vector
136			vProjDist = sqrt(vProj0[0]vProj0[0] + vProj0[1]vProj0[1]
137			+ vProj0[2]*vProj0[2]);
138			vProj0[0] = vProj0[0]/vProjDist;
139			vProj0[1] = vProj0[1]/vProjDist;
140			vProj0[2] = vProj0[2]/vProjDist;
141
142			// Omega is the dot product of the new y-axis and the projection
143			// of the reference y-axis on the current xy-plane
144			omega = acos(vTxvProj0[0] + vTyvProj0[1] + vTz*vProj0[2]);
145
146			return;
147			}
148
149			void Restraints::Calc_body_omegaVal(double matrix[3][3], int currentMol){
150			double zRotator[3][3];
151			double tempOmega;
152			double wholeTwoPis;
153			// Use the omega accumulated from the rotation propagation
154			omega = zAngle[currentMol];
155
156			// translate the omega into a range between -PI and PI
157			if (omega < -PI){
158			tempOmega = omega / -TWO_PI;
159			wholeTwoPis = floor(tempOmega);
160			tempOmega = omega + TWO_PI*wholeTwoPis;
161			if (tempOmega < -PI)
162			omega = tempOmega + TWO_PI;
163			else
164			omega = tempOmega;
165			}
166			if (omega > PI){
167			tempOmega = omega / TWO_PI;
168			wholeTwoPis = floor(tempOmega);
169			tempOmega = omega - TWO_PI*wholeTwoPis;
170			if (tempOmega > PI)
171			omega = tempOmega - TWO_PI;
172			else
173			omega = tempOmega;
174			}
175
176			vb0x = sin(omega);
177			vb0y = cos(omega);
178			vb0z = 0.0;
179
180			normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
181			vb0x = vb0x/normalize;
182			vb0y = vb0y/normalize;
183			vb0z = vb0z/normalize;
184
185			return;
186			}
187
188			double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
189			double pos[3];
190			double A[3][3];
191			double tolerance;
192			double tempPotent;
193			double factor;
194			double spaceTrq[3];
195
196			// atoms = atomPoint;
197
198			// kDist = 6; // spring constant in units of kcal/(mol*ang^2)
199			// kTheta = 7.5; // in units of kcal/mol
200			// kOmega = 13.5; // in units of kcal/mol
201
202			tolerance = 5.72957795131e-7;
203
204			harmPotent = 0.0; // zero out the global harmonic potential variable
205
206			factor = 1 - pow(lambdaValue, lambdaK);
207
208			for (i=0; i<vecParticles.size(); i++){
209			if (vecParticles[i]->isDirectional()){
210			vecParticles[i]->getPos(pos);
211			vecParticles[i]->getA(A);
212			Calc_rVal( pos, i );
213			Calc_body_thetaVal( A, i );
214			Calc_body_omegaVal( A, i );
215
216			if (omega > PI \|\| omega < -PI)
217			cout << "oops... " << omega << "\n";
218
219			// first we calculate the derivatives
220			dVdrx = -kDist*delRx;
221			dVdry = -kDist*delRy;
222			dVdrz = -kDist*delRz;
223
224			// uTx... and vTx... are the body-fixed z and y unit vectors
225			uTx = 0.0;
226			uTy = 0.0;
227			uTz = 1.0;
228			vTx = 0.0;
229			vTy = 1.0;
230			vTz = 0.0;
231
232			dVdux = 0;
233			dVduy = 0;
234			dVduz = 0;
235			dVdvx = 0;
236			dVdvy = 0;
237			dVdvz = 0;
238
239			if (fabs(theta) > tolerance) {
240			dVdux = -(kThetatheta/sin(theta))ub0x;
241			dVduy = -(kThetatheta/sin(theta))ub0y;
242			dVduz = -(kThetatheta/sin(theta))ub0z;
243			}
244
245			if (fabs(omega) > tolerance) {
246			dVdvx = -(kOmegaomega/sin(omega))vb0x;
247			dVdvy = -(kOmegaomega/sin(omega))vb0y;
248			dVdvz = -(kOmegaomega/sin(omega))vb0z;
249			}
250
251			// next we calculate the restraint forces and torques
252			restraintFrc[0] = dVdrx;
253			restraintFrc[1] = dVdry;
254			restraintFrc[2] = dVdrz;
255			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
256
257			restraintTrq[0] = 0.0;
258			restraintTrq[1] = 0.0;
259			restraintTrq[2] = 0.0;
260
261			if (fabs(omega) > tolerance) {
262			restraintTrq[0] += 0.0;
263			restraintTrq[1] += 0.0;
264			restraintTrq[2] += vTy*dVdvx;
265			tempPotent += 0.5(kOmegaomega*omega);
266			}
267			if (fabs(theta) > tolerance) {
268			restraintTrq[0] += (uTz*dVduy);
269			restraintTrq[1] += -(uTz*dVdux);
270			restraintTrq[2] += 0.0;
271			tempPotent += 0.5(kThetatheta*theta);
272			}
273
274			for (j = 0; j < 3; j++) {
275			restraintFrc[j] *= factor;
276			restraintTrq[j] *= factor;
277			}
278
279			harmPotent += tempPotent;
280
281			// now we need to convert from body-fixed torques to space-fixed torques
282			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
283			+ A[2][0]*restraintTrq[2];
284			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
285			+ A[2][1]*restraintTrq[2];
286			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
287			+ A[2][2]*restraintTrq[2];
288
289			// now it's time to pass these temporary forces and torques
290			// to the total forces and torques
291			vecParticles[i]->addFrc(restraintFrc);
292			vecParticles[i]->addTrq(spaceTrq);
293			}
294			}
295
296			// and we can return the appropriately scaled potential energy
297			tempPotent = harmPotent * factor;
298			return tempPotent;
299			}
300
301			void Restraints::Store_Init_Info(){
302			double pos[3];
303			double A[3][3];
304			double RfromQ[3][3];
305			double quat0, quat1, quat2, quat3;
306			double dot;
307			// char *token;
308			// char fileName[200];
309			// char angleName[200];
310			// char inLine[1000];
311			// char inValue[200];
312			const char *delimit = " \t\n;,";
313
314			//open the idealCrystal.in file and zAngle.ang file
315			strcpy(fileName, "idealCrystal.in");
316			strcpy(angleName, "zAngle.ang");
317
318			ifstream crystalIn(fileName);
319			ifstream angleIn(angleName);
320
321			if (!crystalIn) {
322			cout << "Unable to open idealCrystal.in for reading.\n";
323			return;
324			}
325
326			if (!angleIn) {
327			cout << "Unable to open zAngle.ang for reading.\n";
328			cout << "The omega values are all assumed to be zero.\n";
329			}
330
331			// A rather specific reader for OOPSE .eor files...
332			// Let's read in the perfect crystal file
333			crystalIn.getline(inLine,999,'\n');
334			crystalIn.getline(inLine,999,'\n');
335
336			for (i=0; i<nMol; i++) {
337			crystalIn.getline(inLine,999,'\n');
338			token = strtok(inLine,delimit);
339			token = strtok(NULL,delimit);
340			strcpy(inValue,token);
341			cofmPosX.push_back(atof(inValue));
342			token = strtok(NULL,delimit);
343			strcpy(inValue,token);
344			cofmPosY.push_back(atof(inValue));
345			token = strtok(NULL,delimit);
346			strcpy(inValue,token);
347			cofmPosZ.push_back(atof(inValue));
348			token = strtok(NULL,delimit);
349			token = strtok(NULL,delimit);
350			token = strtok(NULL,delimit);
351			token = strtok(NULL,delimit);
352			strcpy(inValue,token);
353			quat0 = atof(inValue);
354			token = strtok(NULL,delimit);
355			strcpy(inValue,token);
356			quat1 = atof(inValue);
357			token = strtok(NULL,delimit);
358			strcpy(inValue,token);
359			quat2 = atof(inValue);
360			token = strtok(NULL,delimit);
361			strcpy(inValue,token);
362			quat3 = atof(inValue);
363
364			// now build the rotation matrix and find the unit vectors
365			RfromQ[0][0] = quat0quat0 + quat1quat1 - quat2quat2 - quat3quat3;
366			RfromQ[0][1] = 2(quat1quat2 + quat0*quat3);
367			RfromQ[0][2] = 2(quat1quat3 - quat0*quat2);
368			RfromQ[1][0] = 2(quat1quat2 - quat0*quat3);
369			RfromQ[1][1] = quat0quat0 - quat1quat1 + quat2quat2 - quat3quat3;
370			RfromQ[1][2] = 2(quat2quat3 + quat0*quat1);
371			RfromQ[2][0] = 2(quat1quat3 + quat0*quat2);
372			RfromQ[2][1] = 2(quat2quat3 - quat0*quat1);
373			RfromQ[2][2] = quat0quat0 - quat1quat1 - quat2quat2 + quat3quat3;
374
375			normalize = sqrt(RfromQ[2][0]RfromQ[2][0] + RfromQ[2][1]RfromQ[2][1]
376			+ RfromQ[2][2]*RfromQ[2][2]);
377			uX0.push_back(RfromQ[2][0]/normalize);
378			uY0.push_back(RfromQ[2][1]/normalize);
379			uZ0.push_back(RfromQ[2][2]/normalize);
380
381			normalize = sqrt(RfromQ[1][0]RfromQ[1][0] + RfromQ[1][1]RfromQ[1][1]
382			+ RfromQ[1][2]*RfromQ[1][2]);
383			vX0.push_back(RfromQ[1][0]/normalize);
384			vY0.push_back(RfromQ[1][1]/normalize);
385			vZ0.push_back(RfromQ[1][2]/normalize);
386			}
387
388			// now we can read in the zAngle.ang file
389			if (angleIn){
390			angleIn.getline(inLine,999,'\n');
391			for (i=0; i<nMol; i++) {
392			angleIn.getline(inLine,999,'\n');
393			token = strtok(inLine,delimit);
394			strcpy(inValue,token);
395			zAngle[i] = (atof(inValue));
396			}
397			}
398
399			return;
400			}
401
402			void Restraints::Determine_Lambda(){
403			// double tempEps;
404
405			// atoms = entry_plug->atoms;
406
407			// if (!strcmp(atoms[0]->getType(),"SSD") \|\|
408			// !strcmp(atoms[0]->getType(),"SSD_E") \|\|
409			// !strcmp(atoms[0]->getType(),"SSD_RF") \|\|
410			// !strcmp(atoms[0]->getType(),"SSD1")){
411
412			// tempEps = atoms[0]->getEpsilon();
413			// scaleLam = 1.0 - (tempEps/0.152);
414			// }
415			// else if (!strcmp(atoms[0]->getType(),"O_TIP3P")){
416			// tempEps = atoms[0]->getEpsilon();
417			// scaleLam = 1.0 - (tempEps/0.1521);
418			// }
419			// else if (!strcmp(atoms[0]->getType(),"O_TIP4P")){
420			// tempEps = atoms[0]->getEpsilon();
421			// scaleLam = 1.0 - (tempEps/0.1550);
422			// }
423			// else if (!strcmp(atoms[0]->getType(),"O_TIP5P")){
424			// tempEps = atoms[0]->getEpsilon();
425			// scaleLam = 1.0 - (tempEps/0.16);
426			// }
427			// else if (!strcmp(atoms[0]->getType(),"O_SPCE")){
428			// tempEps = atoms[0]->getEpsilon();
429			// scaleLam = 1.0 - (tempEps/0.15532);
430			// }
431			// else
432			// sprintf( painCave.errMsg,
433			// "Error in setting the lambda scale: Restraints\n" );
434
435			// if (fabs(scaleLam < 1e-9))
436			// scaleLam = 0.0;
437			// cout << "The scaleLam is " << scaleLam << "\n";
438			}
439
440			void Restraints::Write_zAngle_File(){
441
442			char zOutName[200];
443
444			strcpy(zOutName,"zAngle.ang");
445
446			ofstream angleOut(zOutName);
447			angleOut << "This file contains the omega values for the .eor file\n";
448			for (i=0; i<nMol; i++)
449			angleOut << zAngle[i] << "\n";
450
451			return;
452			}
453
454			double Restraints::getVharm(){
455			return harmPotent;
456			}
457