src/restraints/Restraints.cpp

// Thermodynamic integration is not multiprocessor friendly right now

#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;
  vector<double> resConsts;
  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,
              "Unable to open HarmSpringConsts.txt for reading, so the\n"
              "\tdefault spring constants will be loaded. If you want\n"
              "\tto specify spring constants, include a three line\n"
              "\tHarmSpringConsts.txt file in the execution 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');
      // the file is blank!
      if (springs.eof()){
      sprintf(painCave.errMsg,
              "HarmSpringConsts.txt file is not valid.\n"
              "\tThe file should contain four rows, the last three containing\n"
              "\ta label and the spring constant value. They should be listed\n"
              "\tin the following order: kDist (positional restrant), kTheta\n"
              "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
              "\trestraint: rotation about the z-axis).\n");
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();   
      }
      // read in spring constants and check to make sure it is a valid file
      springs.getline(inLine,999,'\n');
      while (!springs.eof()){
        if (NULL != inLine){
          token = strtok(inLine,jolt);
          token = strtok(NULL,jolt);
          if (NULL != token){
            strcpy(inValue,token);
            resConsts.push_back(atof(inValue));
          }
        }
        springs.getline(inLine,999,'\n');
      }
      if (resConsts.size() == 3){
        kDist = resConsts[0];
        kTheta = resConsts[1];
        kOmega = resConsts[2];
      }
      else {
        sprintf(painCave.errMsg,
                "HarmSpringConsts.txt file is not valid.\n"
                "\tThe file should contain four rows, the last three containing\n"
                "\ta label and the spring constant value. They should be listed\n"
                "\tin the following order: kDist (positional restrant), kTheta\n"
                "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
                "\trestraint: rotation about the z-axis).\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();  
      }
    }
#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 );

      // 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){
  int idealSize;
  double pos[3];
  double A[3][3];
  double RfromQ[3][3];
  double quat0, quat1, quat2, quat3;
  double dot;
  vector<double> tempZangs;
  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);

  // check to see if these files are present in the execution directory
  if (!crystalIn) { 
    sprintf(painCave.errMsg,
            "Restraints Error: Unable to open idealCrystal.in for reading.\n"
            "\tMake sure a ref. crystal file is in the working directory.\n");
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();   
  }

  // it's not fatal to lack a zAngle.ang file, it just means you're starting
  // from the ideal crystal state
  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, the energy\n"
            "\tcalculations will be wrong.\n");
    painCave.severity = OOPSE_WARNING;
    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');
  // check to see if the crystal file is the same length as starting config.
  token = strtok(inLine,delimit);
  strcpy(inValue,token);
  idealSize = atoi(inValue);
  if (idealSize != vecParticles.size()) {
    sprintf(painCave.errMsg,
            "Restraints Error: Reference crystal file is not valid.\n"
            "\tMake sure the idealCrystal.in file is the same size as the\n"
            "\tstarting configuration. Using an incompatable crystal will\n"
            "\tlead to energy calculation failures.\n");
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();  
  }
  // else, the file is okay... let's continue
  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);
  }
  crystalIn.close();

  // now we read in the zAngle.ang file
  if (angleIn){
    angleIn.getline(inLine,999,'\n');
    angleIn.getline(inLine,999,'\n');
    while (!angleIn.eof()) {
      token = strtok(inLine,delimit);
      strcpy(inValue,token);
      tempZangs.push_back(atof(inValue));
      angleIn.getline(inLine,999,'\n');
    }

    // test to make sure the zAngle.ang file is the proper length
    if (tempZangs.size() == vecParticles.size())
      for (i=0; i<vecParticles.size(); i++)
        vecParticles[i]->setZangle(tempZangs[i]);
    else {
      sprintf(painCave.errMsg,
              "Restraints Error: the supplied zAngle file is not valid.\n"
              "\tMake sure the zAngle.ang file matches with the initial\n"
              "\tconfiguration (i.e. they're the same length). Using the wrong\n"
              "\tzAngle file will lead to errors in the energy calculations.\n");
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();  
    }
  }
  angleIn.close();

  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:	43
Committed:	Wed Oct 6 20:01:07 2004 UTC (20 years, 6 months ago) by chrisfen
File size:	13734 byte(s)
Log Message:	made the file reading in restraints more intelligent - no longer gives uninformative seg faults
#	User	Rev	Content
1	chrisfen	43	// Thermodynamic integration is not multiprocessor friendly right now
2
3	gezelter	2	#include <iostream>
4			#include <stdlib.h>
5			#include <cstdio>
6			#include <fstream>
7			#include <iomanip>
8			#include <string>
9			#include <cstring>
10			#include <math.h>
11
12			using namespace std;
13
14	chrisfen	43	#include "Restraints.hpp"
15			#include "SimInfo.hpp"
16			#include "simError.h"
17	gezelter	2
18			#define PI 3.14159265359
19			#define TWO_PI 6.28318530718
20
21			Restraints::Restraints(double lambdaVal, double lambdaExp){
22			lambdaValue = lambdaVal;
23			lambdaK = lambdaExp;
24	chrisfen	43	vector<double> resConsts;
25	gezelter	2	const char *jolt = " \t\n;,";
26
27			#ifdef IS_MPI
28			if(worldRank == 0 ){
29			#endif // is_mpi
30
31			strcpy(springName, "HarmSpringConsts.txt");
32
33			ifstream springs(springName);
34
35			if (!springs) {
36			sprintf(painCave.errMsg,
37	chrisfen	43	"Unable to open HarmSpringConsts.txt for reading, so the\n"
38			"\tdefault spring constants will be loaded. If you want\n"
39			"\tto specify spring constants, include a three line\n"
40			"\tHarmSpringConsts.txt file in the execution directory.\n");
41	gezelter	2	painCave.severity = OOPSE_WARNING;
42			painCave.isFatal = 0;
43			simError();
44
45			// load default spring constants
46			kDist = 6; // spring constant in units of kcal/(mol*ang^2)
47			kTheta = 7.5; // in units of kcal/mol
48			kOmega = 13.5; // in units of kcal/mol
49			} else {
50
51			springs.getline(inLine,999,'\n');
52	chrisfen	43	// the file is blank!
53			if (springs.eof()){
54			sprintf(painCave.errMsg,
55			"HarmSpringConsts.txt file is not valid.\n"
56			"\tThe file should contain four rows, the last three containing\n"
57			"\ta label and the spring constant value. They should be listed\n"
58			"\tin the following order: kDist (positional restrant), kTheta\n"
59			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
60			"\trestraint: rotation about the z-axis).\n");
61			painCave.severity = OOPSE_ERROR;
62			painCave.isFatal = 1;
63			simError();
64			}
65			// read in spring constants and check to make sure it is a valid file
66	gezelter	2	springs.getline(inLine,999,'\n');
67	chrisfen	43	while (!springs.eof()){
68			if (NULL != inLine){
69			token = strtok(inLine,jolt);
70			token = strtok(NULL,jolt);
71			if (NULL != token){
72			strcpy(inValue,token);
73			resConsts.push_back(atof(inValue));
74			}
75			}
76			springs.getline(inLine,999,'\n');
77			}
78			if (resConsts.size() == 3){
79			kDist = resConsts[0];
80			kTheta = resConsts[1];
81			kOmega = resConsts[2];
82			}
83			else {
84			sprintf(painCave.errMsg,
85			"HarmSpringConsts.txt file is not valid.\n"
86			"\tThe file should contain four rows, the last three containing\n"
87			"\ta label and the spring constant value. They should be listed\n"
88			"\tin the following order: kDist (positional restrant), kTheta\n"
89			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
90			"\trestraint: rotation about the z-axis).\n");
91			painCave.severity = OOPSE_ERROR;
92			painCave.isFatal = 1;
93			simError();
94			}
95	gezelter	2	}
96			#ifdef IS_MPI
97			}
98
99			MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
100			MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
101			MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
102
103			sprintf( checkPointMsg,
104			"Sucessfully opened and read spring file.\n");
105			MPIcheckPoint();
106
107			#endif // is_mpi
108
109			sprintf(painCave.errMsg,
110			"The spring constants for thermodynamic integration are:\n"
111			"\tkDist = %lf\n"
112			"\tkTheta = %lf\n"
113			"\tkOmega = %lf\n", kDist, kTheta, kOmega);
114			painCave.severity = OOPSE_INFO;
115			painCave.isFatal = 0;
116			simError();
117			}
118
119			Restraints::~Restraints(){
120			}
121
122			void Restraints::Calc_rVal(double position[3], int currentMol){
123			delRx = position[0] - cofmPosX[currentMol];
124			delRy = position[1] - cofmPosY[currentMol];
125			delRz = position[2] - cofmPosZ[currentMol];
126
127			return;
128			}
129
130			void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
131			ub0x = matrix[0][0]uX0[currentMol] + matrix[0][1]uY0[currentMol]
132			+ matrix[0][2]*uZ0[currentMol];
133			ub0y = matrix[1][0]uX0[currentMol] + matrix[1][1]uY0[currentMol]
134			+ matrix[1][2]*uZ0[currentMol];
135			ub0z = matrix[2][0]uX0[currentMol] + matrix[2][1]uY0[currentMol]
136			+ matrix[2][2]*uZ0[currentMol];
137
138			normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
139			ub0x = ub0x/normalize;
140			ub0y = ub0y/normalize;
141			ub0z = ub0z/normalize;
142
143			// Theta is the dot product of the reference and new z-axes
144			theta = acos(ub0z);
145
146			return;
147			}
148
149			void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){
150			double zRotator[3][3];
151			double tempOmega;
152			double wholeTwoPis;
153			// Use the omega accumulated from the rotation propagation
154			omega = zAngle;
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			double omegaPass;
196
197			tolerance = 5.72957795131e-7;
198
199			harmPotent = 0.0; // zero out the global harmonic potential variable
200
201			factor = 1 - pow(lambdaValue, lambdaK);
202
203			for (i=0; i<vecParticles.size(); i++){
204			if (vecParticles[i]->isDirectional()){
205			vecParticles[i]->getPos(pos);
206			vecParticles[i]->getA(A);
207			Calc_rVal( pos, i );
208			Calc_body_thetaVal( A, i );
209			omegaPass = vecParticles[i]->getZangle();
210			Calc_body_omegaVal( A, omegaPass );
211
212			// first we calculate the derivatives
213			dVdrx = -kDist*delRx;
214			dVdry = -kDist*delRy;
215			dVdrz = -kDist*delRz;
216
217			// uTx... and vTx... are the body-fixed z and y unit vectors
218			uTx = 0.0;
219			uTy = 0.0;
220			uTz = 1.0;
221			vTx = 0.0;
222			vTy = 1.0;
223			vTz = 0.0;
224
225			dVdux = 0;
226			dVduy = 0;
227			dVduz = 0;
228			dVdvx = 0;
229			dVdvy = 0;
230			dVdvz = 0;
231
232			if (fabs(theta) > tolerance) {
233			dVdux = -(kThetatheta/sin(theta))ub0x;
234			dVduy = -(kThetatheta/sin(theta))ub0y;
235			dVduz = -(kThetatheta/sin(theta))ub0z;
236			}
237
238			if (fabs(omega) > tolerance) {
239			dVdvx = -(kOmegaomega/sin(omega))vb0x;
240			dVdvy = -(kOmegaomega/sin(omega))vb0y;
241			dVdvz = -(kOmegaomega/sin(omega))vb0z;
242			}
243
244			// next we calculate the restraint forces and torques
245			restraintFrc[0] = dVdrx;
246			restraintFrc[1] = dVdry;
247			restraintFrc[2] = dVdrz;
248			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
249
250			restraintTrq[0] = 0.0;
251			restraintTrq[1] = 0.0;
252			restraintTrq[2] = 0.0;
253
254			if (fabs(omega) > tolerance) {
255			restraintTrq[0] += 0.0;
256			restraintTrq[1] += 0.0;
257			restraintTrq[2] += vTy*dVdvx;
258			tempPotent += 0.5(kOmegaomega*omega);
259			}
260			if (fabs(theta) > tolerance) {
261			restraintTrq[0] += (uTz*dVduy);
262			restraintTrq[1] += -(uTz*dVdux);
263			restraintTrq[2] += 0.0;
264			tempPotent += 0.5(kThetatheta*theta);
265			}
266
267			for (j = 0; j < 3; j++) {
268			restraintFrc[j] *= factor;
269			restraintTrq[j] *= factor;
270			}
271
272			harmPotent += tempPotent;
273
274			// now we need to convert from body-fixed torques to space-fixed torques
275			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
276			+ A[2][0]*restraintTrq[2];
277			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
278			+ A[2][1]*restraintTrq[2];
279			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
280			+ A[2][2]*restraintTrq[2];
281
282			// now it's time to pass these temporary forces and torques
283			// to the total forces and torques
284			vecParticles[i]->addFrc(restraintFrc);
285			vecParticles[i]->addTrq(spaceTrq);
286			}
287			}
288
289			// and we can return the appropriately scaled potential energy
290			tempPotent = harmPotent * factor;
291			return tempPotent;
292			}
293
294			void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
295	chrisfen	43	int idealSize;
296	gezelter	2	double pos[3];
297			double A[3][3];
298			double RfromQ[3][3];
299			double quat0, quat1, quat2, quat3;
300			double dot;
301	chrisfen	43	vector<double> tempZangs;
302	gezelter	2	const char *delimit = " \t\n;,";
303
304			//open the idealCrystal.in file and zAngle.ang file
305			strcpy(fileName, "idealCrystal.in");
306			strcpy(angleName, "zAngle.ang");
307
308			ifstream crystalIn(fileName);
309			ifstream angleIn(angleName);
310
311	chrisfen	43	// check to see if these files are present in the execution directory
312	gezelter	2	if (!crystalIn) {
313			sprintf(painCave.errMsg,
314			"Restraints Error: Unable to open idealCrystal.in for reading.\n"
315	chrisfen	43	"\tMake sure a ref. crystal file is in the working directory.\n");
316			painCave.severity = OOPSE_ERROR;
317	gezelter	2	painCave.isFatal = 1;
318			simError();
319			}
320
321	chrisfen	43	// it's not fatal to lack a zAngle.ang file, it just means you're starting
322			// from the ideal crystal state
323	gezelter	2	if (!angleIn) {
324			sprintf(painCave.errMsg,
325			"Restraints Warning: The lack of a zAngle.ang file is mildly\n"
326			"\tunsettling... This means the simulation is starting from the\n"
327			"\tidealCrystal.in reference configuration, so the omega values\n"
328	chrisfen	43	"\twill all be set to zero. If this is not the case, the energy\n"
329			"\tcalculations will be wrong.\n");
330			painCave.severity = OOPSE_WARNING;
331	gezelter	2	painCave.isFatal = 0;
332			simError();
333			}
334
335			// A rather specific reader for OOPSE .eor files...
336			// Let's read in the perfect crystal file
337			crystalIn.getline(inLine,999,'\n');
338	chrisfen	43	// check to see if the crystal file is the same length as starting config.
339			token = strtok(inLine,delimit);
340			strcpy(inValue,token);
341			idealSize = atoi(inValue);
342			if (idealSize != vecParticles.size()) {
343			sprintf(painCave.errMsg,
344			"Restraints Error: Reference crystal file is not valid.\n"
345			"\tMake sure the idealCrystal.in file is the same size as the\n"
346			"\tstarting configuration. Using an incompatable crystal will\n"
347			"\tlead to energy calculation failures.\n");
348			painCave.severity = OOPSE_ERROR;
349			painCave.isFatal = 1;
350			simError();
351			}
352			// else, the file is okay... let's continue
353	gezelter	2	crystalIn.getline(inLine,999,'\n');
354
355			for (i=0; i<vecParticles.size(); i++) {
356			crystalIn.getline(inLine,999,'\n');
357			token = strtok(inLine,delimit);
358			token = strtok(NULL,delimit);
359			strcpy(inValue,token);
360			cofmPosX.push_back(atof(inValue));
361			token = strtok(NULL,delimit);
362			strcpy(inValue,token);
363			cofmPosY.push_back(atof(inValue));
364			token = strtok(NULL,delimit);
365			strcpy(inValue,token);
366			cofmPosZ.push_back(atof(inValue));
367			token = strtok(NULL,delimit);
368			token = strtok(NULL,delimit);
369			token = strtok(NULL,delimit);
370			token = strtok(NULL,delimit);
371			strcpy(inValue,token);
372			quat0 = atof(inValue);
373			token = strtok(NULL,delimit);
374			strcpy(inValue,token);
375			quat1 = atof(inValue);
376			token = strtok(NULL,delimit);
377			strcpy(inValue,token);
378			quat2 = atof(inValue);
379			token = strtok(NULL,delimit);
380			strcpy(inValue,token);
381			quat3 = atof(inValue);
382
383			// now build the rotation matrix and find the unit vectors
384			RfromQ[0][0] = quat0quat0 + quat1quat1 - quat2quat2 - quat3quat3;
385			RfromQ[0][1] = 2(quat1quat2 + quat0*quat3);
386			RfromQ[0][2] = 2(quat1quat3 - quat0*quat2);
387			RfromQ[1][0] = 2(quat1quat2 - quat0*quat3);
388			RfromQ[1][1] = quat0quat0 - quat1quat1 + quat2quat2 - quat3quat3;
389			RfromQ[1][2] = 2(quat2quat3 + quat0*quat1);
390			RfromQ[2][0] = 2(quat1quat3 + quat0*quat2);
391			RfromQ[2][1] = 2(quat2quat3 - quat0*quat1);
392			RfromQ[2][2] = quat0quat0 - quat1quat1 - quat2quat2 + quat3quat3;
393
394			normalize = sqrt(RfromQ[2][0]RfromQ[2][0] + RfromQ[2][1]RfromQ[2][1]
395			+ RfromQ[2][2]*RfromQ[2][2]);
396			uX0.push_back(RfromQ[2][0]/normalize);
397			uY0.push_back(RfromQ[2][1]/normalize);
398			uZ0.push_back(RfromQ[2][2]/normalize);
399
400			normalize = sqrt(RfromQ[1][0]RfromQ[1][0] + RfromQ[1][1]RfromQ[1][1]
401			+ RfromQ[1][2]*RfromQ[1][2]);
402			vX0.push_back(RfromQ[1][0]/normalize);
403			vY0.push_back(RfromQ[1][1]/normalize);
404			vZ0.push_back(RfromQ[1][2]/normalize);
405			}
406	chrisfen	43	crystalIn.close();
407	gezelter	2
408	chrisfen	43	// now we read in the zAngle.ang file
409	gezelter	2	if (angleIn){
410			angleIn.getline(inLine,999,'\n');
411	chrisfen	43	angleIn.getline(inLine,999,'\n');
412			while (!angleIn.eof()) {
413	gezelter	2	token = strtok(inLine,delimit);
414			strcpy(inValue,token);
415	chrisfen	43	tempZangs.push_back(atof(inValue));
416			angleIn.getline(inLine,999,'\n');
417	gezelter	2	}
418	chrisfen	43
419			// test to make sure the zAngle.ang file is the proper length
420			if (tempZangs.size() == vecParticles.size())
421			for (i=0; i<vecParticles.size(); i++)
422			vecParticles[i]->setZangle(tempZangs[i]);
423			else {
424			sprintf(painCave.errMsg,
425			"Restraints Error: the supplied zAngle file is not valid.\n"
426			"\tMake sure the zAngle.ang file matches with the initial\n"
427			"\tconfiguration (i.e. they're the same length). Using the wrong\n"
428			"\tzAngle file will lead to errors in the energy calculations.\n");
429			painCave.severity = OOPSE_ERROR;
430			painCave.isFatal = 1;
431			simError();
432			}
433	gezelter	2	}
434	chrisfen	43	angleIn.close();
435	gezelter	2
436			return;
437			}
438
439			void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
440
441			char zOutName[200];
442
443			strcpy(zOutName,"zAngle.ang");
444
445			ofstream angleOut(zOutName);
446			angleOut << "This file contains the omega values for the .eor file\n";
447			for (i=0; i<vecParticles.size(); i++) {
448			angleOut << vecParticles[i]->getZangle() << "\n";
449			}
450			return;
451			}
452
453			double Restraints::getVharm(){
454			return harmPotent;
455			}
456