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:	507
Committed:	Fri Apr 15 22:04:00 2005 UTC (20 years, 6 months ago) by gezelter
File size:	12973 byte(s)
Log Message:	xemacs has been drafted to perform our indentation services
#	User	Rev	Content
1	gezelter	507	/*
2	chrisfen	417	* 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	chrisfen	43
42	gezelter	2	#include <stdlib.h>
43			#include <math.h>
44
45			using namespace std;
46
47	gezelter	46	#include "restraints/Restraints.hpp"
48	chrisfen	417	#include "primitives/Molecule.hpp"
49	gezelter	46	#include "utils/simError.h"
50	gezelter	2
51			#define PI 3.14159265359
52			#define TWO_PI 6.28318530718
53
54	chrisfen	417	namespace oopse {
55
56			Restraints::Restraints(SimInfo* info, double lambdaVal, double lambdaExp){
57			info_ = info;
58			Globals* simParam = info_->getSimParams();
59	gezelter	2
60	chrisfen	417	lambdaValue = lambdaVal;
61			lambdaK = lambdaExp;
62	gezelter	2
63	chrisfen	417	if (simParam->getUseSolidThermInt()) {
64			if (simParam->haveDistSpringConst()) {
65			kDist = simParam->getDistSpringConst();
66	chrisfen	43	}
67	chrisfen	417	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	chrisfen	43	}
78	chrisfen	417	if (simParam->haveThetaSpringConst()) {
79			kTheta = simParam->getThetaSpringConst();
80	chrisfen	43	}
81	chrisfen	417	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	chrisfen	43	}
92	chrisfen	417	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	gezelter	2	}
107	chrisfen	417
108			// build a RestReader and read in important information
109	chrisfen	431
110	chrisfen	417	restRead_ = new RestReader(info_);
111			restRead_->readIdealCrystal();
112			restRead_->readZangle();
113
114			delete restRead_;
115			restRead_ = NULL;
116
117	gezelter	2	}
118
119	chrisfen	417	Restraints::~Restraints(){
120			}
121	gezelter	2
122	chrisfen	417	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	gezelter	2
130	chrisfen	417	void Restraints::Calc_body_thetaVal(RotMat3x3d &matrix, double refUnit[3]){
131			ub0x = matrix(0,0)refUnit[0] + matrix(0,1)refUnit[1]
132	gezelter	507	+ matrix(0,2)*refUnit[2];
133	chrisfen	417	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	gezelter	2	}
148	chrisfen	417
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	chrisfen	221	}
166	chrisfen	417	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	chrisfen	221	}
175
176	chrisfen	417	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	chrisfen	221	}
370
371	gezelter	2	}
372	chrisfen	417
373			// we can finish by returning the appropriately scaled potential energy
374			tempPotent = harmPotent * factor;
375
376			return tempPotent;
377
378	gezelter	2	}
379	chrisfen	221
380	chrisfen	417	}// end namespace oopse