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/Restraints.hpp"

#include "brains/SimInfo.hpp"

#include "utils/simError.h"


#define PI 3.14159265359

#define TWO_PI 6.28318530718


Restraints::Restraints(double lambdaVal, double lambdaExp){

  lambdaValue = lambdaVal;

  lambdaK = lambdaExp;

   std::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()){

      pos = vecParticles[i]->getPos();

      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;

   std::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:	1826
Committed:	Thu Dec 2 05:04:20 2004 UTC (19 years, 8 months ago) by tim
File size:	14229 byte(s)
Log Message:	adding qulified name prefix std
#	User	Rev	Content
1	tim	1826	// Thermodynamic integration is not multiprocessor friendly right now
2
3
4
5			#include <iostream>
6
7			#include <stdlib.h>
8
9			#include <cstdio>
10
11			#include <fstream>
12
13			#include <iomanip>
14
15			#include <string>
16
17			#include <cstring>
18
19			#include <math.h>
20
21
22
23			using namespace std;
24
25
26
27			#include "restraints/Restraints.hpp"
28
29			#include "brains/SimInfo.hpp"
30
31			#include "utils/simError.h"
32
33
34
35			#define PI 3.14159265359
36
37			#define TWO_PI 6.28318530718
38
39
40
41			Restraints::Restraints(double lambdaVal, double lambdaExp){
42
43			lambdaValue = lambdaVal;
44
45			lambdaK = lambdaExp;
46
47			std::vector<double> resConsts;
48
49			const char *jolt = " \t\n;,";
50
51
52
53			#ifdef IS_MPI
54
55			if(worldRank == 0 ){
56
57			#endif // is_mpi
58
59
60
61			strcpy(springName, "HarmSpringConsts.txt");
62
63
64
65			ifstream springs(springName);
66
67
68
69			if (!springs) {
70
71			sprintf(painCave.errMsg,
72
73			"Unable to open HarmSpringConsts.txt for reading, so the\n"
74
75			"\tdefault spring constants will be loaded. If you want\n"
76
77			"\tto specify spring constants, include a three line\n"
78
79			"\tHarmSpringConsts.txt file in the execution directory.\n");
80
81			painCave.severity = OOPSE_WARNING;
82
83			painCave.isFatal = 0;
84
85			simError();
86
87
88
89			// load default spring constants
90
91			kDist = 6; // spring constant in units of kcal/(mol*ang^2)
92
93			kTheta = 7.5; // in units of kcal/mol
94
95			kOmega = 13.5; // in units of kcal/mol
96
97			} else {
98
99
100
101			springs.getline(inLine,999,'\n');
102
103			// the file is blank!
104
105			if (springs.eof()){
106
107			sprintf(painCave.errMsg,
108
109			"HarmSpringConsts.txt file is not valid.\n"
110
111			"\tThe file should contain four rows, the last three containing\n"
112
113			"\ta label and the spring constant value. They should be listed\n"
114
115			"\tin the following order: kDist (positional restrant), kTheta\n"
116
117			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
118
119			"\trestraint: rotation about the z-axis).\n");
120
121			painCave.severity = OOPSE_ERROR;
122
123			painCave.isFatal = 1;
124
125			simError();
126
127			}
128
129			// read in spring constants and check to make sure it is a valid file
130
131			springs.getline(inLine,999,'\n');
132
133			while (!springs.eof()){
134
135			if (NULL != inLine){
136
137			token = strtok(inLine,jolt);
138
139			token = strtok(NULL,jolt);
140
141			if (NULL != token){
142
143			strcpy(inValue,token);
144
145			resConsts.push_back(atof(inValue));
146
147			}
148
149			}
150
151			springs.getline(inLine,999,'\n');
152
153			}
154
155			if (resConsts.size() == 3){
156
157			kDist = resConsts[0];
158
159			kTheta = resConsts[1];
160
161			kOmega = resConsts[2];
162
163			}
164
165			else {
166
167			sprintf(painCave.errMsg,
168
169			"HarmSpringConsts.txt file is not valid.\n"
170
171			"\tThe file should contain four rows, the last three containing\n"
172
173			"\ta label and the spring constant value. They should be listed\n"
174
175			"\tin the following order: kDist (positional restrant), kTheta\n"
176
177			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
178
179			"\trestraint: rotation about the z-axis).\n");
180
181			painCave.severity = OOPSE_ERROR;
182
183			painCave.isFatal = 1;
184
185			simError();
186
187			}
188
189			}
190
191			#ifdef IS_MPI
192
193			}
194
195
196
197			MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
198
199			MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
200
201			MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
202
203
204
205			sprintf( checkPointMsg,
206
207			"Sucessfully opened and read spring file.\n");
208
209			MPIcheckPoint();
210
211
212
213			#endif // is_mpi
214
215
216
217			sprintf(painCave.errMsg,
218
219			"The spring constants for thermodynamic integration are:\n"
220
221			"\tkDist = %lf\n"
222
223			"\tkTheta = %lf\n"
224
225			"\tkOmega = %lf\n", kDist, kTheta, kOmega);
226
227			painCave.severity = OOPSE_INFO;
228
229			painCave.isFatal = 0;
230
231			simError();
232
233			}
234
235
236
237			Restraints::~Restraints(){
238
239			}
240
241
242
243			void Restraints::Calc_rVal(double position[3], int currentMol){
244
245			delRx = position[0] - cofmPosX[currentMol];
246
247			delRy = position[1] - cofmPosY[currentMol];
248
249			delRz = position[2] - cofmPosZ[currentMol];
250
251
252
253			return;
254
255			}
256
257
258
259			void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
260
261			ub0x = matrix[0][0]uX0[currentMol] + matrix[0][1]uY0[currentMol]
262
263			+ matrix[0][2]*uZ0[currentMol];
264
265			ub0y = matrix[1][0]uX0[currentMol] + matrix[1][1]uY0[currentMol]
266
267			+ matrix[1][2]*uZ0[currentMol];
268
269			ub0z = matrix[2][0]uX0[currentMol] + matrix[2][1]uY0[currentMol]
270
271			+ matrix[2][2]*uZ0[currentMol];
272
273
274
275			normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
276
277			ub0x = ub0x/normalize;
278
279			ub0y = ub0y/normalize;
280
281			ub0z = ub0z/normalize;
282
283
284
285			// Theta is the dot product of the reference and new z-axes
286
287			theta = acos(ub0z);
288
289
290
291			return;
292
293			}
294
295
296
297			void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){
298
299			double zRotator[3][3];
300
301			double tempOmega;
302
303			double wholeTwoPis;
304
305			// Use the omega accumulated from the rotation propagation
306
307			omega = zAngle;
308
309
310
311			// translate the omega into a range between -PI and PI
312
313			if (omega < -PI){
314
315			tempOmega = omega / -TWO_PI;
316
317			wholeTwoPis = floor(tempOmega);
318
319			tempOmega = omega + TWO_PI*wholeTwoPis;
320
321			if (tempOmega < -PI)
322
323			omega = tempOmega + TWO_PI;
324
325			else
326
327			omega = tempOmega;
328
329			}
330
331			if (omega > PI){
332
333			tempOmega = omega / TWO_PI;
334
335			wholeTwoPis = floor(tempOmega);
336
337			tempOmega = omega - TWO_PI*wholeTwoPis;
338
339			if (tempOmega > PI)
340
341			omega = tempOmega - TWO_PI;
342
343			else
344
345			omega = tempOmega;
346
347			}
348
349
350
351			vb0x = sin(omega);
352
353			vb0y = cos(omega);
354
355			vb0z = 0.0;
356
357
358
359			normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
360
361			vb0x = vb0x/normalize;
362
363			vb0y = vb0y/normalize;
364
365			vb0z = vb0z/normalize;
366
367
368
369			return;
370
371			}
372
373
374
375			double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
376
377			double pos[3];
378
379			double A[3][3];
380
381			double tolerance;
382
383			double tempPotent;
384
385			double factor;
386
387			double spaceTrq[3];
388
389			double omegaPass;
390
391
392
393			tolerance = 5.72957795131e-7;
394
395
396
397			harmPotent = 0.0; // zero out the global harmonic potential variable
398
399
400
401			factor = 1 - pow(lambdaValue, lambdaK);
402
403
404
405			for (i=0; i<vecParticles.size(); i++){
406
407			if (vecParticles[i]->isDirectional()){
408
409			pos = vecParticles[i]->getPos();
410
411			vecParticles[i]->getA(A);
412
413			Calc_rVal( pos, i );
414
415			Calc_body_thetaVal( A, i );
416
417			omegaPass = vecParticles[i]->getZangle();
418
419			Calc_body_omegaVal( A, omegaPass );
420
421
422
423			// first we calculate the derivatives
424
425			dVdrx = -kDist*delRx;
426
427			dVdry = -kDist*delRy;
428
429			dVdrz = -kDist*delRz;
430
431
432
433			// uTx... and vTx... are the body-fixed z and y unit vectors
434
435			uTx = 0.0;
436
437			uTy = 0.0;
438
439			uTz = 1.0;
440
441			vTx = 0.0;
442
443			vTy = 1.0;
444
445			vTz = 0.0;
446
447
448
449			dVdux = 0;
450
451			dVduy = 0;
452
453			dVduz = 0;
454
455			dVdvx = 0;
456
457			dVdvy = 0;
458
459			dVdvz = 0;
460
461
462
463			if (fabs(theta) > tolerance) {
464
465			dVdux = -(kThetatheta/sin(theta))ub0x;
466
467			dVduy = -(kThetatheta/sin(theta))ub0y;
468
469			dVduz = -(kThetatheta/sin(theta))ub0z;
470
471			}
472
473
474
475			if (fabs(omega) > tolerance) {
476
477			dVdvx = -(kOmegaomega/sin(omega))vb0x;
478
479			dVdvy = -(kOmegaomega/sin(omega))vb0y;
480
481			dVdvz = -(kOmegaomega/sin(omega))vb0z;
482
483			}
484
485
486
487			// next we calculate the restraint forces and torques
488
489			restraintFrc[0] = dVdrx;
490
491			restraintFrc[1] = dVdry;
492
493			restraintFrc[2] = dVdrz;
494
495			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
496
497
498
499			restraintTrq[0] = 0.0;
500
501			restraintTrq[1] = 0.0;
502
503			restraintTrq[2] = 0.0;
504
505
506
507			if (fabs(omega) > tolerance) {
508
509			restraintTrq[0] += 0.0;
510
511			restraintTrq[1] += 0.0;
512
513			restraintTrq[2] += vTy*dVdvx;
514
515			tempPotent += 0.5(kOmegaomega*omega);
516
517			}
518
519			if (fabs(theta) > tolerance) {
520
521			restraintTrq[0] += (uTz*dVduy);
522
523			restraintTrq[1] += -(uTz*dVdux);
524
525			restraintTrq[2] += 0.0;
526
527			tempPotent += 0.5(kThetatheta*theta);
528
529			}
530
531
532
533			for (j = 0; j < 3; j++) {
534
535			restraintFrc[j] *= factor;
536
537			restraintTrq[j] *= factor;
538
539			}
540
541
542
543			harmPotent += tempPotent;
544
545
546
547			// now we need to convert from body-fixed torques to space-fixed torques
548
549			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
550
551			+ A[2][0]*restraintTrq[2];
552
553			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
554
555			+ A[2][1]*restraintTrq[2];
556
557			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
558
559			+ A[2][2]*restraintTrq[2];
560
561
562
563			// now it's time to pass these temporary forces and torques
564
565			// to the total forces and torques
566
567			vecParticles[i]->addFrc(restraintFrc);
568
569			vecParticles[i]->addTrq(spaceTrq);
570
571			}
572
573			}
574
575
576
577			// and we can return the appropriately scaled potential energy
578
579			tempPotent = harmPotent * factor;
580
581			return tempPotent;
582
583			}
584
585
586
587			void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
588
589			int idealSize;
590
591			double pos[3];
592
593			double A[3][3];
594
595			double RfromQ[3][3];
596
597			double quat0, quat1, quat2, quat3;
598
599			double dot;
600
601			std::vector<double> tempZangs;
602
603			const char *delimit = " \t\n;,";
604
605
606
607			//open the idealCrystal.in file and zAngle.ang file
608
609			strcpy(fileName, "idealCrystal.in");
610
611			strcpy(angleName, "zAngle.ang");
612
613
614
615			ifstream crystalIn(fileName);
616
617			ifstream angleIn(angleName);
618
619
620
621			// check to see if these files are present in the execution directory
622
623			if (!crystalIn) {
624
625			sprintf(painCave.errMsg,
626
627			"Restraints Error: Unable to open idealCrystal.in for reading.\n"
628
629			"\tMake sure a ref. crystal file is in the working directory.\n");
630
631			painCave.severity = OOPSE_ERROR;
632
633			painCave.isFatal = 1;
634
635			simError();
636
637			}
638
639
640
641			// it's not fatal to lack a zAngle.ang file, it just means you're starting
642
643			// from the ideal crystal state
644
645			if (!angleIn) {
646
647			sprintf(painCave.errMsg,
648
649			"Restraints Warning: The lack of a zAngle.ang file is mildly\n"
650
651			"\tunsettling... This means the simulation is starting from the\n"
652
653			"\tidealCrystal.in reference configuration, so the omega values\n"
654
655			"\twill all be set to zero. If this is not the case, the energy\n"
656
657			"\tcalculations will be wrong.\n");
658
659			painCave.severity = OOPSE_WARNING;
660
661			painCave.isFatal = 0;
662
663			simError();
664
665			}
666
667
668
669			// A rather specific reader for OOPSE .eor files...
670
671			// Let's read in the perfect crystal file
672
673			crystalIn.getline(inLine,999,'\n');
674
675			// check to see if the crystal file is the same length as starting config.
676
677			token = strtok(inLine,delimit);
678
679			strcpy(inValue,token);
680
681			idealSize = atoi(inValue);
682
683			if (idealSize != vecParticles.size()) {
684
685			sprintf(painCave.errMsg,
686
687			"Restraints Error: Reference crystal file is not valid.\n"
688
689			"\tMake sure the idealCrystal.in file is the same size as the\n"
690
691			"\tstarting configuration. Using an incompatable crystal will\n"
692
693			"\tlead to energy calculation failures.\n");
694
695			painCave.severity = OOPSE_ERROR;
696
697			painCave.isFatal = 1;
698
699			simError();
700
701			}
702
703			// else, the file is okay... let's continue
704
705			crystalIn.getline(inLine,999,'\n');
706
707
708
709			for (i=0; i<vecParticles.size(); i++) {
710
711			crystalIn.getline(inLine,999,'\n');
712
713			token = strtok(inLine,delimit);
714
715			token = strtok(NULL,delimit);
716
717			strcpy(inValue,token);
718
719			cofmPosX.push_back(atof(inValue));
720
721			token = strtok(NULL,delimit);
722
723			strcpy(inValue,token);
724
725			cofmPosY.push_back(atof(inValue));
726
727			token = strtok(NULL,delimit);
728
729			strcpy(inValue,token);
730
731			cofmPosZ.push_back(atof(inValue));
732
733			token = strtok(NULL,delimit);
734
735			token = strtok(NULL,delimit);
736
737			token = strtok(NULL,delimit);
738
739			token = strtok(NULL,delimit);
740
741			strcpy(inValue,token);
742
743			quat0 = atof(inValue);
744
745			token = strtok(NULL,delimit);
746
747			strcpy(inValue,token);
748
749			quat1 = atof(inValue);
750
751			token = strtok(NULL,delimit);
752
753			strcpy(inValue,token);
754
755			quat2 = atof(inValue);
756
757			token = strtok(NULL,delimit);
758
759			strcpy(inValue,token);
760
761			quat3 = atof(inValue);
762
763
764
765			// now build the rotation matrix and find the unit vectors
766
767			RfromQ[0][0] = quat0quat0 + quat1quat1 - quat2quat2 - quat3quat3;
768
769			RfromQ[0][1] = 2(quat1quat2 + quat0*quat3);
770
771			RfromQ[0][2] = 2(quat1quat3 - quat0*quat2);
772
773			RfromQ[1][0] = 2(quat1quat2 - quat0*quat3);
774
775			RfromQ[1][1] = quat0quat0 - quat1quat1 + quat2quat2 - quat3quat3;
776
777			RfromQ[1][2] = 2(quat2quat3 + quat0*quat1);
778
779			RfromQ[2][0] = 2(quat1quat3 + quat0*quat2);
780
781			RfromQ[2][1] = 2(quat2quat3 - quat0*quat1);
782
783			RfromQ[2][2] = quat0quat0 - quat1quat1 - quat2quat2 + quat3quat3;
784
785
786
787			normalize = sqrt(RfromQ[2][0]RfromQ[2][0] + RfromQ[2][1]RfromQ[2][1]
788
789			+ RfromQ[2][2]*RfromQ[2][2]);
790
791			uX0.push_back(RfromQ[2][0]/normalize);
792
793			uY0.push_back(RfromQ[2][1]/normalize);
794
795			uZ0.push_back(RfromQ[2][2]/normalize);
796
797
798
799			normalize = sqrt(RfromQ[1][0]RfromQ[1][0] + RfromQ[1][1]RfromQ[1][1]
800
801			+ RfromQ[1][2]*RfromQ[1][2]);
802
803			vX0.push_back(RfromQ[1][0]/normalize);
804
805			vY0.push_back(RfromQ[1][1]/normalize);
806
807			vZ0.push_back(RfromQ[1][2]/normalize);
808
809			}
810
811			crystalIn.close();
812
813
814
815			// now we read in the zAngle.ang file
816
817			if (angleIn){
818
819			angleIn.getline(inLine,999,'\n');
820
821			angleIn.getline(inLine,999,'\n');
822
823			while (!angleIn.eof()) {
824
825			token = strtok(inLine,delimit);
826
827			strcpy(inValue,token);
828
829			tempZangs.push_back(atof(inValue));
830
831			angleIn.getline(inLine,999,'\n');
832
833			}
834
835
836
837			// test to make sure the zAngle.ang file is the proper length
838
839			if (tempZangs.size() == vecParticles.size())
840
841			for (i=0; i<vecParticles.size(); i++)
842
843			vecParticles[i]->setZangle(tempZangs[i]);
844
845			else {
846
847			sprintf(painCave.errMsg,
848
849			"Restraints Error: the supplied zAngle file is not valid.\n"
850
851			"\tMake sure the zAngle.ang file matches with the initial\n"
852
853			"\tconfiguration (i.e. they're the same length). Using the wrong\n"
854
855			"\tzAngle file will lead to errors in the energy calculations.\n");
856
857			painCave.severity = OOPSE_ERROR;
858
859			painCave.isFatal = 1;
860
861			simError();
862
863			}
864
865			}
866
867			angleIn.close();
868
869
870
871			return;
872
873			}
874
875
876
877			void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
878
879
880
881			char zOutName[200];
882
883
884
885			strcpy(zOutName,"zAngle.ang");
886
887
888
889			ofstream angleOut(zOutName);
890
891			angleOut << "This file contains the omega values for the .eor file\n";
892
893			for (i=0; i<vecParticles.size(); i++) {
894
895			angleOut << vecParticles[i]->getZangle() << "\n";
896
897			}
898
899			return;
900
901			}
902
903
904
905			double Restraints::getVharm(){
906
907			return harmPotent;
908
909			}
910
911
912