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:	417
Committed:	Thu Mar 10 15:10:24 2005 UTC (20 years, 1 month ago) by chrisfen
File size:	12967 byte(s)
Log Message:	First commit of the new restraints code
#	User	Rev	Content
1	chrisfen	417	/*
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	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			restRead_ = new RestReader(info_);
110			restRead_->readIdealCrystal();
111			restRead_->readZangle();
112
113			delete restRead_;
114			restRead_ = NULL;
115
116	gezelter	2	}
117
118	chrisfen	417	Restraints::~Restraints(){
119			}
120	gezelter	2
121	chrisfen	417	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	gezelter	2
129	chrisfen	417	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	gezelter	2	}
147	chrisfen	417
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	chrisfen	221	}
165	chrisfen	417	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	chrisfen	221	}
174
175	chrisfen	417	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	chrisfen	221	}
369
370	gezelter	2	}
371	chrisfen	417
372			// we can finish by returning the appropriately scaled potential energy
373			tempPotent = harmPotent * factor;
374
375			return tempPotent;
376
377	gezelter	2	}
378	chrisfen	221
379	chrisfen	417	}// end namespace oopse