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:	2115
Committed:	Fri Mar 11 00:43:19 2005 UTC (19 years, 4 months ago) by chrisfen
File size:	12972 byte(s)
Log Message:	fixed a bug in MPI restraints
#	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
110	restRead_ = new RestReader(info_);
111	restRead_->readIdealCrystal();
112	restRead_->readZangle();
113
114	delete restRead_;
115	restRead_ = NULL;
116
117	}
118
119	Restraints::~Restraints(){
120	}
121
122	void Restraints::Calc_rVal(Vector3d &position, double refPosition[3]){
123	delRx = position.x() - refPosition[0];
124	delRy = position.y() - refPosition[1];
125	delRz = position.z() - refPosition[2];
126
127	return;
128	}
129
130	void Restraints::Calc_body_thetaVal(RotMat3x3d &matrix, double refUnit[3]){
131	ub0x = matrix(0,0)refUnit[0] + matrix(0,1)refUnit[1]
132	+ matrix(0,2)*refUnit[2];
133	ub0y = matrix(1,0)refUnit[0] + matrix(1,1)refUnit[1]
134	+ matrix(1,2)*refUnit[2];
135	ub0z = matrix(2,0)refUnit[0] + matrix(2,1)refUnit[1]
136	+ matrix(2,2)*refUnit[2];
137
138	normalize = sqrt(ub0xub0x + ub0yub0y + ub0z*ub0z);
139	ub0x = ub0x/normalize;
140	ub0y = ub0y/normalize;
141	ub0z = ub0z/normalize;
142
143	// Theta is the dot product of the reference and new z-axes
144	theta = acos(ub0z);
145
146	return;
147	}
148
149	void Restraints::Calc_body_omegaVal(double zAngle){
150	double zRotator[3][3];
151	double tempOmega;
152	double wholeTwoPis;
153	// Use the omega accumulated from the rotation propagation
154	omega = zAngle;
155
156	// translate the omega into a range between -PI and PI
157	if (omega < -PI){
158	tempOmega = omega / -TWO_PI;
159	wholeTwoPis = floor(tempOmega);
160	tempOmega = omega + TWO_PI*wholeTwoPis;
161	if (tempOmega < -PI)
162	omega = tempOmega + TWO_PI;
163	else
164	omega = tempOmega;
165	}
166	if (omega > PI){
167	tempOmega = omega / TWO_PI;
168	wholeTwoPis = floor(tempOmega);
169	tempOmega = omega - TWO_PI*wholeTwoPis;
170	if (tempOmega > PI)
171	omega = tempOmega - TWO_PI;
172	else
173	omega = tempOmega;
174	}
175
176	vb0x = sin(omega);
177	vb0y = cos(omega);
178	vb0z = 0.0;
179
180	normalize = sqrt(vb0xvb0x + vb0yvb0y + vb0z*vb0z);
181	vb0x = vb0x/normalize;
182	vb0y = vb0y/normalize;
183	vb0z = vb0z/normalize;
184
185	return;
186	}
187
188	double Restraints::Calc_Restraint_Forces(){
189	SimInfo::MoleculeIterator mi;
190	Molecule* mol;
191	Molecule::IntegrableObjectIterator ii;
192	StuntDouble* integrableObject;
193	Vector3d pos;
194	RotMat3x3d A;
195	double refPos[3];
196	double refVec[3];
197	double tolerance;
198	double tempPotent;
199	double factor;
200	double spaceTrq[3];
201	double omegaPass;
202	GenericData* data;
203	DoubleGenericData* doubleData;
204
205	tolerance = 5.72957795131e-7;
206
207	harmPotent = 0.0; // zero out the global harmonic potential variable
208
209	factor = 1 - pow(lambdaValue, lambdaK);
210
211	for (mol = info_->beginMolecule(mi); mol != NULL;
212	mol = info_->nextMolecule(mi)) {
213	for (integrableObject = mol->beginIntegrableObject(ii);
214	integrableObject != NULL;
215	integrableObject = mol->nextIntegrableObject(ii)) {
216
217	// obtain the current and reference positions
218	pos = integrableObject->getPos();
219
220	data = integrableObject->getPropertyByName("refPosX");
221	if (data){
222	doubleData = dynamic_cast<DoubleGenericData*>(data);
223	if (!doubleData){
224	cerr << "Can't obtain refPosX from StuntDouble\n";
225	return 0.0;
226	}
227	else refPos[0] = doubleData->getData();
228	}
229	data = integrableObject->getPropertyByName("refPosY");
230	if (data){
231	doubleData = dynamic_cast<DoubleGenericData*>(data);
232	if (!doubleData){
233	cerr << "Can't obtain refPosY from StuntDouble\n";
234	return 0.0;
235	}
236	else refPos[1] = doubleData->getData();
237	}
238	data = integrableObject->getPropertyByName("refPosZ");
239	if (data){
240	doubleData = dynamic_cast<DoubleGenericData*>(data);
241	if (!doubleData){
242	cerr << "Can't obtain refPosZ from StuntDouble\n";
243	return 0.0;
244	}
245	else refPos[2] = doubleData->getData();
246	}
247
248	// calculate the displacement
249	Calc_rVal( pos, refPos );
250
251	// calculate the derivatives
252	dVdrx = -kDist*delRx;
253	dVdry = -kDist*delRy;
254	dVdrz = -kDist*delRz;
255
256	// next we calculate the restraint forces
257	restraintFrc[0] = dVdrx;
258	restraintFrc[1] = dVdry;
259	restraintFrc[2] = dVdrz;
260	tempPotent = 0.5kDist(delRxdelRx + delRydelRy + delRz*delRz);
261
262	// apply the lambda scaling factor to the forces
263	for (j = 0; j < 3; j++) restraintFrc[j] *= factor;
264
265	// and add the temporary force to the total force
266	integrableObject->addFrc(restraintFrc);
267
268	// if the particle is directional, we accumulate the rot. restraints
269	if (integrableObject->isDirectional()){
270
271	// get the current rotation matrix and reference vector
272	A = integrableObject->getA();
273
274	data = integrableObject->getPropertyByName("refVectorX");
275	if (data){
276	doubleData = dynamic_cast<DoubleGenericData*>(data);
277	if (!doubleData){
278	cerr << "Can't obtain refVectorX from StuntDouble\n";
279	return 0.0;
280	}
281	else refVec[0] = doubleData->getData();
282	}
283	data = integrableObject->getPropertyByName("refVectorY");
284	if (data){
285	doubleData = dynamic_cast<DoubleGenericData*>(data);
286	if (!doubleData){
287	cerr << "Can't obtain refVectorY from StuntDouble\n";
288	return 0.0;
289	}
290	else refVec[1] = doubleData->getData();
291	}
292	data = integrableObject->getPropertyByName("refVectorZ");
293	if (data){
294	doubleData = dynamic_cast<DoubleGenericData*>(data);
295	if (!doubleData){
296	cerr << "Can't obtain refVectorZ from StuntDouble\n";
297	return 0.0;
298	}
299	else refVec[2] = doubleData->getData();
300	}
301
302	// calculate the theta and omega displacements
303	Calc_body_thetaVal( A, refVec );
304	omegaPass = integrableObject->getZangle();
305	Calc_body_omegaVal( omegaPass );
306
307	// uTx... and vTx... are the body-fixed z and y unit vectors
308	uTx = 0.0;
309	uTy = 0.0;
310	uTz = 1.0;
311	vTx = 0.0;
312	vTy = 1.0;
313	vTz = 0.0;
314
315	dVdux = 0.0;
316	dVduy = 0.0;
317	dVduz = 0.0;
318	dVdvx = 0.0;
319	dVdvy = 0.0;
320	dVdvz = 0.0;
321
322	if (fabs(theta) > tolerance) {
323	dVdux = -(kThetatheta/sin(theta))ub0x;
324	dVduy = -(kThetatheta/sin(theta))ub0y;
325	dVduz = -(kThetatheta/sin(theta))ub0z;
326	}
327
328	if (fabs(omega) > tolerance) {
329	dVdvx = -(kOmegaomega/sin(omega))vb0x;
330	dVdvy = -(kOmegaomega/sin(omega))vb0y;
331	dVdvz = -(kOmegaomega/sin(omega))vb0z;
332	}
333
334	// next we calculate the restraint torques
335	restraintTrq[0] = 0.0;
336	restraintTrq[1] = 0.0;
337	restraintTrq[2] = 0.0;
338
339	if (fabs(omega) > tolerance) {
340	restraintTrq[0] += 0.0;
341	restraintTrq[1] += 0.0;
342	restraintTrq[2] += vTy*dVdvx;
343	tempPotent += 0.5(kOmegaomega*omega);
344	}
345	if (fabs(theta) > tolerance) {
346	restraintTrq[0] += (uTz*dVduy);
347	restraintTrq[1] += -(uTz*dVdux);
348	restraintTrq[2] += 0.0;
349	tempPotent += 0.5(kThetatheta*theta);
350	}
351
352	// apply the lambda scaling factor to these torques
353	for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
354
355	// now we need to convert from body-fixed to space-fixed torques
356	spaceTrq[0] = A(0,0)restraintTrq[0] + A(1,0)restraintTrq[1]
357	+ A(2,0)*restraintTrq[2];
358	spaceTrq[1] = A(0,1)restraintTrq[0] + A(1,1)restraintTrq[1]
359	+ A(2,1)*restraintTrq[2];
360	spaceTrq[2] = A(0,2)restraintTrq[0] + A(1,2)restraintTrq[1]
361	+ A(2,2)*restraintTrq[2];
362
363	// now pass this temporary torque vector to the total torque
364	integrableObject->addTrq(spaceTrq);
365	}
366
367	// update the total harmonic potential with this object's contribution
368	harmPotent += tempPotent;
369	}
370
371	}
372
373	// we can finish by returning the appropriately scaled potential energy
374	tempPotent = harmPotent * factor;
375
376	return tempPotent;
377
378	}
379
380	}// end namespace oopse