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
#	Content
1	// Thermodynamic integration is not multiprocessor friendly right now
2
3	#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	#include "restraints/Restraints.hpp"
15	#include "brains/SimInfo.hpp"
16	#include "utils/simError.h"
17	#include "io/basic_ifstrstream.hpp"
18
19	#ifdef IS_MPI
20	#include<mpi.h>
21	#include "brains/mpiSimulation.hpp"
22	#endif // is_mpi
23
24	#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	vector<double> resConsts;
31	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	"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	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	// 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	springs.getline(inLine,999,'\n');
73	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	}
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	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
133	return;
134	}
135
136	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
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	double refPos[3];
198	double refVec[3];
199	double tolerance;
200	double tempPotent;
201	double factor;
202	double spaceTrq[3];
203	double omegaPass;
204	GenericData* data;
205	DoubleGenericData* doubleData;
206
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	// 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	if (vecParticles[i]->isDirectional()){
267
268	// get the current rotation matrix and reference vector
269	vecParticles[i]->getA(A);
270
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	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	dVdux = 0.0;
313	dVduy = 0.0;
314	dVduz = 0.0;
315	dVdvx = 0.0;
316	dVdvy = 0.0;
317	dVdvz = 0.0;
318
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	// next we calculate the restraint torques
332	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	// apply the lambda scaling factor to these torques
350	for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
351
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	// now pass this temporary torque vector to the total torque
361	vecParticles[i]->addTrq(spaceTrq);
362	}
363
364	// update the total harmonic potential with this object's contribution
365	harmPotent += tempPotent;
366	}
367
368	// we can finish by returning the appropriately scaled potential energy
369	tempPotent = harmPotent * factor;
370	return tempPotent;
371	}
372
373	void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles,
374	int currTime,
375	int nIntObj){
376
377	char zOutName[200];
378
379	std::cerr << nIntObj << " is the number of integrable objects\n";
380
381	//#ifndef IS_MPI
382
383	strcpy(zOutName,"zAngle.ang");
384
385	ofstream angleOut(zOutName);
386	angleOut << currTime << ": omega values at this time\n";
387	for (i=0; i<vecParticles.size(); i++) {
388	angleOut << vecParticles[i]->getZangle() << "\n";
389	}
390
391	return;
392	}
393
394	double Restraints::getVharm(){
395	return harmPotent;
396	}
397