src/restraints/Restraints.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */

#include <stdlib.h>
#include <math.h>

using namespace std;

#include "restraints/Restraints.hpp"
#include "primitives/Molecule.hpp"
#include "utils/simError.h"

#define PI 3.14159265359
#define TWO_PI 6.28318530718

namespace oopse {
  
  Restraints::Restraints(SimInfo* info, double lambdaVal, double lambdaExp){
    info_ = info;
    Globals* simParam = info_->getSimParams();

    lambdaValue = lambdaVal;
    lambdaK = lambdaExp;
    
    if (simParam->getUseSolidThermInt()) {
      if (simParam->haveDistSpringConst()) {
        kDist = simParam->getDistSpringConst();
      }
      else{
        kDist = 6.0;
        sprintf(painCave.errMsg,
                "ThermoIntegration Warning: the spring constant for the\n"
                "\ttranslational restraint was not specified. OOPSE will use\n"
                "\ta default value of %f. To set it to something else, use\n"
                "\tthe thermIntDistSpringConst variable.\n",
                kDist);
        painCave.isFatal = 0;
        simError(); 
      }
      if (simParam->haveThetaSpringConst()) {
        kTheta = simParam->getThetaSpringConst();
      }
      else{
        kTheta = 7.5;
        sprintf(painCave.errMsg,
                "ThermoIntegration Warning: the spring constant for the\n"
                "\tdeflection orientational restraint was not specified.\n"
                "\tOOPSE will use a default value of %f. To set it to\n"
                "\tsomething else, use the thermIntThetaSpringConst variable.\n",
                kTheta);
        painCave.isFatal = 0;
        simError(); 
      }
      if (simParam->haveOmegaSpringConst()) {
        kOmega = simParam->getOmegaSpringConst();
      }
      else{
        kOmega = 13.5;
        sprintf(painCave.errMsg,
                "ThermoIntegration Warning: the spring constant for the\n"
                "\tspin orientational restraint was not specified. OOPSE\n"
                "\twill use a default value of %f. To set it to something\n"
                "\telse, use the thermIntOmegaSpringConst variable.\n",
                kOmega);
        painCave.isFatal = 0;
        simError(); 
      }
    }
    
    // build a RestReader and read in important information
    restRead_ = new RestReader(info_);
    restRead_->readIdealCrystal();
    restRead_->readZangle();
    
    delete restRead_;
    restRead_ = NULL;
    
  }
  
  Restraints::~Restraints(){
  }
  
  void Restraints::Calc_rVal(Vector3d &position, double refPosition[3]){
    delRx = position.x() - refPosition[0];
    delRy = position.y() - refPosition[1];
    delRz = position.z() - refPosition[2];
    
    return;
  }
  
  void Restraints::Calc_body_thetaVal(RotMat3x3d &matrix, 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 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(){
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::IntegrableObjectIterator ii;
    StuntDouble* integrableObject;
    Vector3d pos;
    RotMat3x3d A;
    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 (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
      for (integrableObject = mol->beginIntegrableObject(ii); 
           integrableObject != NULL; 
           integrableObject = mol->nextIntegrableObject(ii)) {
        
        // obtain the current and reference positions
        pos = integrableObject->getPos();
        
        data = integrableObject->getPropertyByName("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 = integrableObject->getPropertyByName("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 = integrableObject->getPropertyByName("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
        integrableObject->addFrc(restraintFrc);
        
        // if the particle is directional, we accumulate the rot. restraints
        if (integrableObject->isDirectional()){
          
          // get the current rotation matrix and reference vector
          A = integrableObject->getA();
          
          data = integrableObject->getPropertyByName("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 = integrableObject->getPropertyByName("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 = integrableObject->getPropertyByName("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 = integrableObject->getZangle();
          Calc_body_omegaVal( 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 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
          integrableObject->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;
    
  }
  
}// end namespace oopse
Revision:	2101
Committed:	Thu Mar 10 15:10:24 2005 UTC (19 years, 4 months ago) by chrisfen
File size:	12967 byte(s)
Log Message:	First commit of the new restraints code
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42	#include <stdlib.h>
43	#include <math.h>
44
45	using namespace std;
46
47	#include "restraints/Restraints.hpp"
48	#include "primitives/Molecule.hpp"
49	#include "utils/simError.h"
50
51	#define PI 3.14159265359
52	#define TWO_PI 6.28318530718
53
54	namespace oopse {
55
56	Restraints::Restraints(SimInfo* info, double lambdaVal, double lambdaExp){
57	info_ = info;
58	Globals* simParam = info_->getSimParams();
59
60	lambdaValue = lambdaVal;
61	lambdaK = lambdaExp;
62
63	if (simParam->getUseSolidThermInt()) {
64	if (simParam->haveDistSpringConst()) {
65	kDist = simParam->getDistSpringConst();
66	}
67	else{
68	kDist = 6.0;
69	sprintf(painCave.errMsg,
70	"ThermoIntegration Warning: the spring constant for the\n"
71	"\ttranslational restraint was not specified. OOPSE will use\n"
72	"\ta default value of %f. To set it to something else, use\n"
73	"\tthe thermIntDistSpringConst variable.\n",
74	kDist);
75	painCave.isFatal = 0;
76	simError();
77	}
78	if (simParam->haveThetaSpringConst()) {
79	kTheta = simParam->getThetaSpringConst();
80	}
81	else{
82	kTheta = 7.5;
83	sprintf(painCave.errMsg,
84	"ThermoIntegration Warning: the spring constant for the\n"
85	"\tdeflection orientational restraint was not specified.\n"
86	"\tOOPSE will use a default value of %f. To set it to\n"
87	"\tsomething else, use the thermIntThetaSpringConst variable.\n",
88	kTheta);
89	painCave.isFatal = 0;
90	simError();
91	}
92	if (simParam->haveOmegaSpringConst()) {
93	kOmega = simParam->getOmegaSpringConst();
94	}
95	else{
96	kOmega = 13.5;
97	sprintf(painCave.errMsg,
98	"ThermoIntegration Warning: the spring constant for the\n"
99	"\tspin orientational restraint was not specified. OOPSE\n"
100	"\twill use a default value of %f. To set it to something\n"
101	"\telse, use the thermIntOmegaSpringConst variable.\n",
102	kOmega);
103	painCave.isFatal = 0;
104	simError();
105	}
106	}
107
108	// build a RestReader and read in important information
109	restRead_ = new RestReader(info_);
110	restRead_->readIdealCrystal();
111	restRead_->readZangle();
112
113	delete restRead_;
114	restRead_ = NULL;
115
116	}
117
118	Restraints::~Restraints(){
119	}
120
121	void Restraints::Calc_rVal(Vector3d &position, double refPosition[3]){
122	delRx = position.x() - refPosition[0];
123	delRy = position.y() - refPosition[1];
124	delRz = position.z() - refPosition[2];
125
126	return;
127	}
128
129	void Restraints::Calc_body_thetaVal(RotMat3x3d &matrix, double refUnit[3]){
130	ub0x = matrix(0,0)refUnit[0] + matrix(0,1)refUnit[1]
131	+ matrix(0,2)*refUnit[2];
132	ub0y = matrix(1,0)refUnit[0] + matrix(1,1)refUnit[1]
133	+ matrix(1,2)*refUnit[2];
134	ub0z = matrix(2,0)refUnit[0] + matrix(2,1)refUnit[1]
135	+ matrix(2,2)*refUnit[2];
136
137	normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
138	ub0x = ub0x/normalize;
139	ub0y = ub0y/normalize;
140	ub0z = ub0z/normalize;
141
142	// Theta is the dot product of the reference and new z-axes
143	theta = acos(ub0z);
144
145	return;
146	}
147
148	void Restraints::Calc_body_omegaVal(double zAngle){
149	double zRotator[3][3];
150	double tempOmega;
151	double wholeTwoPis;
152	// Use the omega accumulated from the rotation propagation
153	omega = zAngle;
154
155	// translate the omega into a range between -PI and PI
156	if (omega < -PI){
157	tempOmega = omega / -TWO_PI;
158	wholeTwoPis = floor(tempOmega);
159	tempOmega = omega + TWO_PI*wholeTwoPis;
160	if (tempOmega < -PI)
161	omega = tempOmega + TWO_PI;
162	else
163	omega = tempOmega;
164	}
165	if (omega > PI){
166	tempOmega = omega / TWO_PI;
167	wholeTwoPis = floor(tempOmega);
168	tempOmega = omega - TWO_PI*wholeTwoPis;
169	if (tempOmega > PI)
170	omega = tempOmega - TWO_PI;
171	else
172	omega = tempOmega;
173	}
174
175	vb0x = sin(omega);
176	vb0y = cos(omega);
177	vb0z = 0.0;
178
179	normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
180	vb0x = vb0x/normalize;
181	vb0y = vb0y/normalize;
182	vb0z = vb0z/normalize;
183
184	return;
185	}
186
187	double Restraints::Calc_Restraint_Forces(){
188	SimInfo::MoleculeIterator mi;
189	Molecule* mol;
190	Molecule::IntegrableObjectIterator ii;
191	StuntDouble* integrableObject;
192	Vector3d pos;
193	RotMat3x3d A;
194	double refPos[3];
195	double refVec[3];
196	double tolerance;
197	double tempPotent;
198	double factor;
199	double spaceTrq[3];
200	double omegaPass;
201	GenericData* data;
202	DoubleGenericData* doubleData;
203
204	tolerance = 5.72957795131e-7;
205
206	harmPotent = 0.0; // zero out the global harmonic potential variable
207
208	factor = 1 - pow(lambdaValue, lambdaK);
209
210	for (mol = info_->beginMolecule(mi); mol != NULL;
211	mol = info_->nextMolecule(mi)) {
212	for (integrableObject = mol->beginIntegrableObject(ii);
213	integrableObject != NULL;
214	integrableObject = mol->nextIntegrableObject(ii)) {
215
216	// obtain the current and reference positions
217	pos = integrableObject->getPos();
218
219	data = integrableObject->getPropertyByName("refPosX");
220	if (data){
221	doubleData = dynamic_cast<DoubleGenericData*>(data);
222	if (!doubleData){
223	cerr << "Can't obtain refPosX from StuntDouble\n";
224	return 0.0;
225	}
226	else refPos[0] = doubleData->getData();
227	}
228	data = integrableObject->getPropertyByName("refPosY");
229	if (data){
230	doubleData = dynamic_cast<DoubleGenericData*>(data);
231	if (!doubleData){
232	cerr << "Can't obtain refPosY from StuntDouble\n";
233	return 0.0;
234	}
235	else refPos[1] = doubleData->getData();
236	}
237	data = integrableObject->getPropertyByName("refPosZ");
238	if (data){
239	doubleData = dynamic_cast<DoubleGenericData*>(data);
240	if (!doubleData){
241	cerr << "Can't obtain refPosZ from StuntDouble\n";
242	return 0.0;
243	}
244	else refPos[2] = doubleData->getData();
245	}
246
247	// calculate the displacement
248	Calc_rVal( pos, refPos );
249
250	// calculate the derivatives
251	dVdrx = -kDist*delRx;
252	dVdry = -kDist*delRy;
253	dVdrz = -kDist*delRz;
254
255	// next we calculate the restraint forces
256	restraintFrc[0] = dVdrx;
257	restraintFrc[1] = dVdry;
258	restraintFrc[2] = dVdrz;
259	tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
260
261	// apply the lambda scaling factor to the forces
262	for (j = 0; j < 3; j++) restraintFrc[j] *= factor;
263
264	// and add the temporary force to the total force
265	integrableObject->addFrc(restraintFrc);
266
267	// if the particle is directional, we accumulate the rot. restraints
268	if (integrableObject->isDirectional()){
269
270	// get the current rotation matrix and reference vector
271	A = integrableObject->getA();
272
273	data = integrableObject->getPropertyByName("refVectorX");
274	if (data){
275	doubleData = dynamic_cast<DoubleGenericData*>(data);
276	if (!doubleData){
277	cerr << "Can't obtain refVectorX from StuntDouble\n";
278	return 0.0;
279	}
280	else refVec[0] = doubleData->getData();
281	}
282	data = integrableObject->getPropertyByName("refVectorY");
283	if (data){
284	doubleData = dynamic_cast<DoubleGenericData*>(data);
285	if (!doubleData){
286	cerr << "Can't obtain refVectorY from StuntDouble\n";
287	return 0.0;
288	}
289	else refVec[1] = doubleData->getData();
290	}
291	data = integrableObject->getPropertyByName("refVectorZ");
292	if (data){
293	doubleData = dynamic_cast<DoubleGenericData*>(data);
294	if (!doubleData){
295	cerr << "Can't obtain refVectorZ from StuntDouble\n";
296	return 0.0;
297	}
298	else refVec[2] = doubleData->getData();
299	}
300
301	// calculate the theta and omega displacements
302	Calc_body_thetaVal( A, refVec );
303	omegaPass = integrableObject->getZangle();
304	Calc_body_omegaVal( omegaPass );
305
306	// uTx... and vTx... are the body-fixed z and y unit vectors
307	uTx = 0.0;
308	uTy = 0.0;
309	uTz = 1.0;
310	vTx = 0.0;
311	vTy = 1.0;
312	vTz = 0.0;
313
314	dVdux = 0.0;
315	dVduy = 0.0;
316	dVduz = 0.0;
317	dVdvx = 0.0;
318	dVdvy = 0.0;
319	dVdvz = 0.0;
320
321	if (fabs(theta) > tolerance) {
322	dVdux = -(kThetatheta/sin(theta))ub0x;
323	dVduy = -(kThetatheta/sin(theta))ub0y;
324	dVduz = -(kThetatheta/sin(theta))ub0z;
325	}
326
327	if (fabs(omega) > tolerance) {
328	dVdvx = -(kOmegaomega/sin(omega))vb0x;
329	dVdvy = -(kOmegaomega/sin(omega))vb0y;
330	dVdvz = -(kOmegaomega/sin(omega))vb0z;
331	}
332
333	// next we calculate the restraint torques
334	restraintTrq[0] = 0.0;
335	restraintTrq[1] = 0.0;
336	restraintTrq[2] = 0.0;
337
338	if (fabs(omega) > tolerance) {
339	restraintTrq[0] += 0.0;
340	restraintTrq[1] += 0.0;
341	restraintTrq[2] += vTy*dVdvx;
342	tempPotent += 0.5(kOmegaomega*omega);
343	}
344	if (fabs(theta) > tolerance) {
345	restraintTrq[0] += (uTz*dVduy);
346	restraintTrq[1] += -(uTz*dVdux);
347	restraintTrq[2] += 0.0;
348	tempPotent += 0.5(kThetatheta*theta);
349	}
350
351	// apply the lambda scaling factor to these torques
352	for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
353
354	// now we need to convert from body-fixed to space-fixed torques
355	spaceTrq[0] = A(0,0)restraintTrq[0] + A(1,0)restraintTrq[1]
356	+ A(2,0)*restraintTrq[2];
357	spaceTrq[1] = A(0,1)restraintTrq[0] + A(1,1)restraintTrq[1]
358	+ A(2,1)*restraintTrq[2];
359	spaceTrq[2] = A(0,2)restraintTrq[0] + A(1,2)restraintTrq[1]
360	+ A(2,2)*restraintTrq[2];
361
362	// now pass this temporary torque vector to the total torque
363	integrableObject->addTrq(spaceTrq);
364	}
365
366	// update the total harmonic potential with this object's contribution
367	harmPotent += tempPotent;
368	}
369
370	}
371
372	// we can finish by returning the appropriately scaled potential energy
373	tempPotent = harmPotent * factor;
374
375	return tempPotent;
376
377	}
378
379	}// end namespace oopse