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"
#include "io/basic_ifstrstream.hpp"

#ifdef IS_MPI
#include<mpi.h>
#include "brains/mpiSimulation.hpp"
#endif // is_mpi

#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], double refPosition[3]){
  delRx = position[0] - refPosition[0];
  delRy = position[1] - refPosition[1];
  delRz = position[2] - refPosition[2];

  return;
}

void Restraints::Calc_body_thetaVal(double matrix[3][3], double refUnit[3]){
  ub0x = matrix[0][0]*refUnit[0] + matrix[0][1]*refUnit[1]
    + matrix[0][2]*refUnit[2];
  ub0y = matrix[1][0]*refUnit[0] + matrix[1][1]*refUnit[1]
    + matrix[1][2]*refUnit[2];
  ub0z = matrix[2][0]*refUnit[0] + matrix[2][1]*refUnit[1]
    + matrix[2][2]*refUnit[2];

  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 refPos[3];
  double refVec[3];
  double tolerance;
  double tempPotent;
  double factor;
  double spaceTrq[3];
  double omegaPass;
  GenericData* data;
  DoubleGenericData* doubleData;

  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++){
    // obtain the current and reference positions
    vecParticles[i]->getPos(pos);

    data = vecParticles[i]->getProperty("refPosX");
    if (data){
      doubleData = dynamic_cast<DoubleGenericData*>(data);
      if (!doubleData){
        cerr << "Can't obtain refPosX from StuntDouble\n";
        return 0.0;
      }
      else refPos[0] = doubleData->getData();
    }
    data = vecParticles[i]->getProperty("refPosY");
    if (data){
      doubleData = dynamic_cast<DoubleGenericData*>(data);
      if (!doubleData){
        cerr << "Can't obtain refPosY from StuntDouble\n";
        return 0.0;
      }
      else refPos[1] = doubleData->getData();
    }
    data = vecParticles[i]->getProperty("refPosZ");
    if (data){
      doubleData = dynamic_cast<DoubleGenericData*>(data);
      if (!doubleData){
        cerr << "Can't obtain refPosZ from StuntDouble\n";
        return 0.0;
      }
      else refPos[2] = doubleData->getData();
    }

    // calculate the displacement
    Calc_rVal( pos, refPos );

    // calculate the derivatives
    dVdrx = -kDist*delRx;
    dVdry = -kDist*delRy;
    dVdrz = -kDist*delRz;
    
    // next we calculate the restraint forces
    restraintFrc[0] = dVdrx;
    restraintFrc[1] = dVdry;
    restraintFrc[2] = dVdrz;
    tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);

    // apply the lambda scaling factor to the forces
    for (j = 0; j < 3; j++) restraintFrc[j] *= factor;

    // and add the temporary force to the total force
    vecParticles[i]->addFrc(restraintFrc);

    // if the particle is directional, we accumulate the rot. restraints
    if (vecParticles[i]->isDirectional()){
 
      // get the current rotation matrix and reference vector
      vecParticles[i]->getA(A);
      
      data = vecParticles[i]->getProperty("refVectorX");
      if (data){
        doubleData = dynamic_cast<DoubleGenericData*>(data);
        if (!doubleData){
          cerr << "Can't obtain refVectorX from StuntDouble\n";
          return 0.0;
        }
        else refVec[0] = doubleData->getData();
      }
      data = vecParticles[i]->getProperty("refVectorY");
      if (data){
        doubleData = dynamic_cast<DoubleGenericData*>(data);
        if (!doubleData){
          cerr << "Can't obtain refVectorY from StuntDouble\n";
          return 0.0;
        }
        else refVec[1] = doubleData->getData();
      }
      data = vecParticles[i]->getProperty("refVectorZ");
      if (data){
        doubleData = dynamic_cast<DoubleGenericData*>(data);
        if (!doubleData){
          cerr << "Can't obtain refVectorZ from StuntDouble\n";
          return 0.0;
        }
        else refVec[2] = doubleData->getData();
      }
      
      // calculate the theta and omega displacements
      Calc_body_thetaVal( A, refVec );
      omegaPass = vecParticles[i]->getZangle();
      Calc_body_omegaVal( A, omegaPass );

      // 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.0;
      dVduy = 0.0;
      dVduz = 0.0;
      dVdvx = 0.0;
      dVdvy = 0.0;
      dVdvz = 0.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 torques
      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);
      }

      // apply the lambda scaling factor to these torques
      for (j = 0; j < 3; j++) restraintTrq[j] *= factor;

      // 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 pass this temporary torque vector to the total torque
      vecParticles[i]->addTrq(spaceTrq);
    }

    // update the total harmonic potential with this object's contribution
    harmPotent += tempPotent;
  }
  
  // we can finish by returning the appropriately scaled potential energy
  tempPotent = harmPotent * factor;
  return tempPotent;
}

void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles,
                                   int currTime,
                                   int nIntObj){

  char zOutName[200];

  std::cerr << nIntObj << " is the number of integrable objects\n";

  //#ifndef IS_MPI
  
  strcpy(zOutName,"zAngle.ang");
  
  ofstream angleOut(zOutName);
  angleOut << currTime << ": omega values at this time\n";
  for (i=0; i<vecParticles.size(); i++) {
    angleOut << vecParticles[i]->getZangle() << "\n";
  }

  return;
}

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

Revision:	1772
Committed:	Tue Nov 23 22:48:31 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	11098 byte(s)
Log Message:	Improvements to restraints
#	User	Rev	Content
1	chrisfen	1534	// Thermodynamic integration is not multiprocessor friendly right now
2
3	gezelter	1490	#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	gezelter	1537	#include "restraints/Restraints.hpp"
15			#include "brains/SimInfo.hpp"
16			#include "utils/simError.h"
17	chrisfen	1731	#include "io/basic_ifstrstream.hpp"
18	gezelter	1490
19	chrisfen	1772	#ifdef IS_MPI
20			#include<mpi.h>
21			#include "brains/mpiSimulation.hpp"
22			#endif // is_mpi
23
24	gezelter	1490	#define PI 3.14159265359
25			#define TWO_PI 6.28318530718
26
27			Restraints::Restraints(double lambdaVal, double lambdaExp){
28			lambdaValue = lambdaVal;
29			lambdaK = lambdaExp;
30	chrisfen	1534	vector<double> resConsts;
31	gezelter	1490	const char *jolt = " \t\n;,";
32
33			#ifdef IS_MPI
34			if(worldRank == 0 ){
35			#endif // is_mpi
36
37			strcpy(springName, "HarmSpringConsts.txt");
38
39			ifstream springs(springName);
40
41			if (!springs) {
42			sprintf(painCave.errMsg,
43	chrisfen	1534	"Unable to open HarmSpringConsts.txt for reading, so the\n"
44			"\tdefault spring constants will be loaded. If you want\n"
45			"\tto specify spring constants, include a three line\n"
46			"\tHarmSpringConsts.txt file in the execution directory.\n");
47	gezelter	1490	painCave.severity = OOPSE_WARNING;
48			painCave.isFatal = 0;
49			simError();
50
51			// load default spring constants
52			kDist = 6; // spring constant in units of kcal/(mol*ang^2)
53			kTheta = 7.5; // in units of kcal/mol
54			kOmega = 13.5; // in units of kcal/mol
55			} else {
56
57			springs.getline(inLine,999,'\n');
58	chrisfen	1534	// the file is blank!
59			if (springs.eof()){
60			sprintf(painCave.errMsg,
61			"HarmSpringConsts.txt file is not valid.\n"
62			"\tThe file should contain four rows, the last three containing\n"
63			"\ta label and the spring constant value. They should be listed\n"
64			"\tin the following order: kDist (positional restrant), kTheta\n"
65			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
66			"\trestraint: rotation about the z-axis).\n");
67			painCave.severity = OOPSE_ERROR;
68			painCave.isFatal = 1;
69			simError();
70			}
71			// read in spring constants and check to make sure it is a valid file
72	gezelter	1490	springs.getline(inLine,999,'\n');
73	chrisfen	1534	while (!springs.eof()){
74			if (NULL != inLine){
75			token = strtok(inLine,jolt);
76			token = strtok(NULL,jolt);
77			if (NULL != token){
78			strcpy(inValue,token);
79			resConsts.push_back(atof(inValue));
80			}
81			}
82			springs.getline(inLine,999,'\n');
83			}
84			if (resConsts.size() == 3){
85			kDist = resConsts[0];
86			kTheta = resConsts[1];
87			kOmega = resConsts[2];
88			}
89			else {
90			sprintf(painCave.errMsg,
91			"HarmSpringConsts.txt file is not valid.\n"
92			"\tThe file should contain four rows, the last three containing\n"
93			"\ta label and the spring constant value. They should be listed\n"
94			"\tin the following order: kDist (positional restrant), kTheta\n"
95			"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
96			"\trestraint: rotation about the z-axis).\n");
97			painCave.severity = OOPSE_ERROR;
98			painCave.isFatal = 1;
99			simError();
100			}
101	gezelter	1490	}
102			#ifdef IS_MPI
103			}
104
105			MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
106			MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
107			MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
108
109			sprintf( checkPointMsg,
110			"Sucessfully opened and read spring file.\n");
111			MPIcheckPoint();
112
113			#endif // is_mpi
114
115			sprintf(painCave.errMsg,
116			"The spring constants for thermodynamic integration are:\n"
117			"\tkDist = %lf\n"
118			"\tkTheta = %lf\n"
119			"\tkOmega = %lf\n", kDist, kTheta, kOmega);
120			painCave.severity = OOPSE_INFO;
121			painCave.isFatal = 0;
122			simError();
123			}
124
125			Restraints::~Restraints(){
126			}
127
128	chrisfen	1772	void Restraints::Calc_rVal(double position[3], double refPosition[3]){
129			delRx = position[0] - refPosition[0];
130			delRy = position[1] - refPosition[1];
131			delRz = position[2] - refPosition[2];
132	gezelter	1490
133			return;
134			}
135
136	chrisfen	1772	void Restraints::Calc_body_thetaVal(double matrix[3][3], double refUnit[3]){
137			ub0x = matrix[0][0]refUnit[0] + matrix[0][1]refUnit[1]
138			+ matrix[0][2]*refUnit[2];
139			ub0y = matrix[1][0]refUnit[0] + matrix[1][1]refUnit[1]
140			+ matrix[1][2]*refUnit[2];
141			ub0z = matrix[2][0]refUnit[0] + matrix[2][1]refUnit[1]
142			+ matrix[2][2]*refUnit[2];
143	gezelter	1490
144			normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
145			ub0x = ub0x/normalize;
146			ub0y = ub0y/normalize;
147			ub0z = ub0z/normalize;
148
149			// Theta is the dot product of the reference and new z-axes
150			theta = acos(ub0z);
151
152			return;
153			}
154
155			void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){
156			double zRotator[3][3];
157			double tempOmega;
158			double wholeTwoPis;
159			// Use the omega accumulated from the rotation propagation
160			omega = zAngle;
161
162			// translate the omega into a range between -PI and PI
163			if (omega < -PI){
164			tempOmega = omega / -TWO_PI;
165			wholeTwoPis = floor(tempOmega);
166			tempOmega = omega + TWO_PI*wholeTwoPis;
167			if (tempOmega < -PI)
168			omega = tempOmega + TWO_PI;
169			else
170			omega = tempOmega;
171			}
172			if (omega > PI){
173			tempOmega = omega / TWO_PI;
174			wholeTwoPis = floor(tempOmega);
175			tempOmega = omega - TWO_PI*wholeTwoPis;
176			if (tempOmega > PI)
177			omega = tempOmega - TWO_PI;
178			else
179			omega = tempOmega;
180			}
181
182			vb0x = sin(omega);
183			vb0y = cos(omega);
184			vb0z = 0.0;
185
186			normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
187			vb0x = vb0x/normalize;
188			vb0y = vb0y/normalize;
189			vb0z = vb0z/normalize;
190
191			return;
192			}
193
194			double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
195			double pos[3];
196			double A[3][3];
197	chrisfen	1772	double refPos[3];
198			double refVec[3];
199	gezelter	1490	double tolerance;
200			double tempPotent;
201			double factor;
202			double spaceTrq[3];
203			double omegaPass;
204	chrisfen	1772	GenericData* data;
205			DoubleGenericData* doubleData;
206	gezelter	1490
207			tolerance = 5.72957795131e-7;
208
209			harmPotent = 0.0; // zero out the global harmonic potential variable
210
211			factor = 1 - pow(lambdaValue, lambdaK);
212
213			for (i=0; i<vecParticles.size(); i++){
214	chrisfen	1772	// obtain the current and reference positions
215			vecParticles[i]->getPos(pos);
216
217			data = vecParticles[i]->getProperty("refPosX");
218			if (data){
219			doubleData = dynamic_cast<DoubleGenericData*>(data);
220			if (!doubleData){
221			cerr << "Can't obtain refPosX from StuntDouble\n";
222			return 0.0;
223			}
224			else refPos[0] = doubleData->getData();
225			}
226			data = vecParticles[i]->getProperty("refPosY");
227			if (data){
228			doubleData = dynamic_cast<DoubleGenericData*>(data);
229			if (!doubleData){
230			cerr << "Can't obtain refPosY from StuntDouble\n";
231			return 0.0;
232			}
233			else refPos[1] = doubleData->getData();
234			}
235			data = vecParticles[i]->getProperty("refPosZ");
236			if (data){
237			doubleData = dynamic_cast<DoubleGenericData*>(data);
238			if (!doubleData){
239			cerr << "Can't obtain refPosZ from StuntDouble\n";
240			return 0.0;
241			}
242			else refPos[2] = doubleData->getData();
243			}
244
245			// calculate the displacement
246			Calc_rVal( pos, refPos );
247
248			// calculate the derivatives
249			dVdrx = -kDist*delRx;
250			dVdry = -kDist*delRy;
251			dVdrz = -kDist*delRz;
252
253			// next we calculate the restraint forces
254			restraintFrc[0] = dVdrx;
255			restraintFrc[1] = dVdry;
256			restraintFrc[2] = dVdrz;
257			tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
258
259			// apply the lambda scaling factor to the forces
260			for (j = 0; j < 3; j++) restraintFrc[j] *= factor;
261
262			// and add the temporary force to the total force
263			vecParticles[i]->addFrc(restraintFrc);
264
265			// if the particle is directional, we accumulate the rot. restraints
266	gezelter	1490	if (vecParticles[i]->isDirectional()){
267	chrisfen	1772
268			// get the current rotation matrix and reference vector
269	gezelter	1490	vecParticles[i]->getA(A);
270	chrisfen	1772
271			data = vecParticles[i]->getProperty("refVectorX");
272			if (data){
273			doubleData = dynamic_cast<DoubleGenericData*>(data);
274			if (!doubleData){
275			cerr << "Can't obtain refVectorX from StuntDouble\n";
276			return 0.0;
277			}
278			else refVec[0] = doubleData->getData();
279			}
280			data = vecParticles[i]->getProperty("refVectorY");
281			if (data){
282			doubleData = dynamic_cast<DoubleGenericData*>(data);
283			if (!doubleData){
284			cerr << "Can't obtain refVectorY from StuntDouble\n";
285			return 0.0;
286			}
287			else refVec[1] = doubleData->getData();
288			}
289			data = vecParticles[i]->getProperty("refVectorZ");
290			if (data){
291			doubleData = dynamic_cast<DoubleGenericData*>(data);
292			if (!doubleData){
293			cerr << "Can't obtain refVectorZ from StuntDouble\n";
294			return 0.0;
295			}
296			else refVec[2] = doubleData->getData();
297			}
298
299			// calculate the theta and omega displacements
300			Calc_body_thetaVal( A, refVec );
301	gezelter	1490	omegaPass = vecParticles[i]->getZangle();
302			Calc_body_omegaVal( A, omegaPass );
303
304			// uTx... and vTx... are the body-fixed z and y unit vectors
305			uTx = 0.0;
306			uTy = 0.0;
307			uTz = 1.0;
308			vTx = 0.0;
309			vTy = 1.0;
310			vTz = 0.0;
311
312	chrisfen	1772	dVdux = 0.0;
313			dVduy = 0.0;
314			dVduz = 0.0;
315			dVdvx = 0.0;
316			dVdvy = 0.0;
317			dVdvz = 0.0;
318	gezelter	1490
319			if (fabs(theta) > tolerance) {
320			dVdux = -(kThetatheta/sin(theta))ub0x;
321			dVduy = -(kThetatheta/sin(theta))ub0y;
322			dVduz = -(kThetatheta/sin(theta))ub0z;
323			}
324
325			if (fabs(omega) > tolerance) {
326			dVdvx = -(kOmegaomega/sin(omega))vb0x;
327			dVdvy = -(kOmegaomega/sin(omega))vb0y;
328			dVdvz = -(kOmegaomega/sin(omega))vb0z;
329			}
330
331	chrisfen	1772	// next we calculate the restraint torques
332	gezelter	1490	restraintTrq[0] = 0.0;
333			restraintTrq[1] = 0.0;
334			restraintTrq[2] = 0.0;
335
336			if (fabs(omega) > tolerance) {
337			restraintTrq[0] += 0.0;
338			restraintTrq[1] += 0.0;
339			restraintTrq[2] += vTy*dVdvx;
340			tempPotent += 0.5(kOmegaomega*omega);
341			}
342			if (fabs(theta) > tolerance) {
343			restraintTrq[0] += (uTz*dVduy);
344			restraintTrq[1] += -(uTz*dVdux);
345			restraintTrq[2] += 0.0;
346			tempPotent += 0.5(kThetatheta*theta);
347			}
348
349	chrisfen	1772	// apply the lambda scaling factor to these torques
350			for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
351	gezelter	1490
352			// now we need to convert from body-fixed torques to space-fixed torques
353			spaceTrq[0] = A[0][0]restraintTrq[0] + A[1][0]restraintTrq[1]
354			+ A[2][0]*restraintTrq[2];
355			spaceTrq[1] = A[0][1]restraintTrq[0] + A[1][1]restraintTrq[1]
356			+ A[2][1]*restraintTrq[2];
357			spaceTrq[2] = A[0][2]restraintTrq[0] + A[1][2]restraintTrq[1]
358			+ A[2][2]*restraintTrq[2];
359
360	chrisfen	1772	// now pass this temporary torque vector to the total torque
361	gezelter	1490	vecParticles[i]->addTrq(spaceTrq);
362			}
363	chrisfen	1772
364			// update the total harmonic potential with this object's contribution
365			harmPotent += tempPotent;
366	gezelter	1490	}
367	chrisfen	1772
368			// we can finish by returning the appropriately scaled potential energy
369	gezelter	1490	tempPotent = harmPotent * factor;
370			return tempPotent;
371			}
372
373	chrisfen	1772	void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles,
374			int currTime,
375			int nIntObj){
376	gezelter	1490
377	chrisfen	1772	char zOutName[200];
378	gezelter	1490
379	chrisfen	1772	std::cerr << nIntObj << " is the number of integrable objects\n";
380	gezelter	1490
381	chrisfen	1772	//#ifndef IS_MPI
382	gezelter	1490
383			strcpy(zOutName,"zAngle.ang");
384	chrisfen	1772
385	gezelter	1490	ofstream angleOut(zOutName);
386	chrisfen	1772	angleOut << currTime << ": omega values at this time\n";
387	gezelter	1490	for (i=0; i<vecParticles.size(); i++) {
388			angleOut << vecParticles[i]->getZangle() << "\n";
389			}
390	chrisfen	1772
391	gezelter	1490	return;
392			}
393
394			double Restraints::getVharm(){
395			return harmPotent;
396			}
397