Re: custom forces in Molmodel
Posted: Tue Jun 09, 2009 4:24 pm
You should create a new Custom Force which queries the CompoundSystem for each of its Compounds, then queries each compound for each of its atom indices. Once you have that you can get each atom's body index and its cartesian coordinates.
When the Force's calcForce method is called it should apply the appropriate force to each atom's body. If the kinematics are coarse grained the force will be applied to the origin of the body frame, and so you would have to add a torque to compensate. But I'm assuming for a first cut you will be doing an atomistic simulation.
RNABuilder's TwoTransformForces.h does something much like this, except it queries a list of interactions rather than atoms. I am appending it at the very bottom of this email. It's an overcomplicated example though. You really ought to start from scratch. Peter documented the Custom Force class pretty well if I recall correctly.
Are you working on protein or RNA?
I submitted this to the forum.
Sam
On Jun 9, 2009, at 2:36 PM, Balaraman, Gouthaman wrote:
Hey Sherm,
The force module is initialized using a pdb, psf, and parameter file. This code is part of a normal Cartesian MD code. I have written a wrapper such that, after the initialization, I can pass a set of coordinates to the function, and get back the cartesian forces.
I believe if I use the Molmodels objects, that would make the clustering part very simple. So I would definitely like to use it, if it is possible. In the future, I would also like to experiment with different types of coarse graining. In that case, can I still use Molmodel's objects to extend to arbitrary clustering?
Thanks,
Goutham
-----Original Message-----
From: Michael Sherman [mailto:msherman@stanford.edu]
Sent: Tue 6/9/2009 2:24 PM
To: Balaraman, Gouthaman
Cc: 'Christopher Bruns'; 'Samuel Coulbourn Flores'
Subject: RE: custom forces in Molmodel
Hi, Goutham.
How does your force field know the bonded structure of the molecule, and how
does that relate to the model built by Molmodel? Are you planning to use
Molmodel to construct the mechanical model of the molecules (via Molmodel's
Protein() or RNA() objects, for example)?
Regards,
Sherm
From: Balaraman, Gouthaman [mailto:GBalaraman@coh.org]
Sent: Tuesday, June 09, 2009 12:40 PM
To: Michael Sherman
Cc: Christopher Bruns; Samuel Coulbourn Flores
Subject: RE: custom forces in Molmodel
Hey All,
I have an existing force field that I want to use as is. This force module
will take cartesian positions as input, and return a force vector in
cartesian coordinates. The force module will take care of reading force
parameters, and evaluating the forces. I want to incorporate this force
engine into molmodel. I want to start with all torsions for the clustering.
In the longer run, I would like to define my own coarse graining as well.
I really appreciate your help on this aspect.
Thank you,
Goutham
-----Original Message-----
From: Michael Sherman [mailto:msherman@stanford.edu]
Sent: Mon 6/8/2009 8:12 PM
To: Balaraman, Gouthaman
Cc: 'Christopher Bruns'; Samuel Coulbourn Flores
Subject: RE: custom forces in Molmodel
Hi, Goutham.
Nice to hear from you again and I'm very glad you're interested in using
Molmodel. I'm including Chris Bruns and Sam Flores on this email - Chris is
the primary author of Molmodel, and Sam is the most advanced Molmodel user
at Stanford. We will need some more information in order to figure out the
best way for you to get your force field working on a system that has been
modeled using Molmodel.
We are about to release version 1.6 of Molmodel in which DuMM has been made
more flexible (and it now uses OpenMM for acceleration when it can also).
Depending on the functional form of your force field, DuMM may be able to
accommodate it (DuMM is neutral about the parameterization). Otherwise there
are a variety of ways to apply different forces - Sam has used Molmodel
successfully with several special purpose RNA force fields. Any combination
is plausible - you can avoid DuMM altogether, you can use it for part of the
force field and something else for other parts, or you can use DuMM directly
(maybe with some hacks).
Please tell us what you have in mind. Do you have an existing force field
that you want to use as-is? What are the functional forms involved? What
mechanical model are you using (all torsions, or something coarser?)
Chris, should we move this discussion to one of the Molmodel forums on
simtk.org?
Best regards,
Sherm
From: Balaraman, Gouthaman [mailto:GBalaraman@coh.org]
Sent: Monday, June 08, 2009 5:42 PM
To: msherman@stanford.edu
Subject: custom forces in Molmodel
Hey Mike,
Thanks for getting back on my question. I am a postdoc at Beckman Institute
of City of Hope in Prof Nagarajan Vaidehi's group. I had a chance to meet
you in an OpenMM session.
I want to use my own force engine, in Molmodel. Can you give me any pointers
on how I can incorporate that within Molmodel?
Thanks,
Goutham
You have been contacted by a member of Simtk.org.
From: Michael Sherman
Email: msherman@stanford.edu
The e-mail was sent via the "contact" link on Simtk.org by the above member.
It is not screened by Simtk.org or its staff.
*** MESSAGE BEGINS ***
Hi, Goutham.
Ayman Habib forwarded me your post to the OpenSim forum. Please send me your
actual email address so I can connect you with the right people here to talk
about custom forces in Molmodel.
Best regards,
Michael Sherman (Sherm)
---------------------------------------------------------------------
SECURITY/CONFIDENTIALITY WARNING: This message and any attachments are
intended solely for the individual or entity to which they are addressed.
This communication may contain information that is privileged, confidential,
or exempt from disclosure under applicable law (e.g., personal health
information, research data, financial information). Because this e-mail has
been sent without encryption, individuals other than the intended recipient
may be able to view the information, forward it to others or tamper with the
information without the knowledge or consent of the sender. If you are not
the intended recipient, or the employee or person responsible for delivering
the message to the intended recipient, any dissemination, distribution or
copying of the communication is strictly prohibited. If you received the
communication in error, please notify the sender immediately by replying to
this message and deleting the message and any accompanying files from your
system. If, due to the security risks, you do not wish to receive further
communications via e-mail, please reply to this message and inform the
sender that you do not wish to receive further e-mail from the sender.
---------------------------------------------------------------------
Samuel Coulbourn Flores, PhD
Altman Lab
650.644.8416
花山
科学者
生物工学部
スタンフォ一ド大学
スタンフォ一ド、カリフォルニア、米国
using namespace SimTK;
using namespace std;
std::ostream& operator<<(std::ostream& o, const ParameterReader&) {
assert(false);
return o;
}
class AllTwoTransformLinearSprings : public Force::Custom::Implementation {
private:
SimbodyMatterSubsystem& matter;
ParameterReader& myParameterReader;
LeontisWesthofClass& myLeontisWesthofClass;
Biopolymer * myChain;
mutable int parameterReaderIndex;
public:
AllTwoTransformLinearSprings (SimbodyMatterSubsystem& matter,ParameterReader& myParameterReader, LeontisWesthofClass& myLeontisWesthofClass, Biopolymer * myChain ) : matter(matter),myParameterReader(myParameterReader), myLeontisWesthofClass (myLeontisWesthofClass), myChain(myChain)
{
}
void realizeTopology(State& state) const {
parameterReaderIndex = state.allocateDiscreteVariable(matter.getMySubsystemIndex(),
Stage::Dynamics, new Value<ParameterReader>(myParameterReader));
//return 0;
//implementation->realizeTopology(state);
}
void setParameterReader(State& state, const ParameterReader& myParameterReader) const {
static_cast<Value<ParameterReader>&>(state.updDiscreteVariable( matter.getMySubsystemIndex(),
SimTK::DiscreteVariableIndex(parameterReaderIndex))).upd() = myParameterReader;
}
const ParameterReader & getParameterReader(const State& state) const {
return static_cast<const Value<ParameterReader>&>(state.getDiscreteVariable(matter.getMySubsystemIndex(),
SimTK::DiscreteVariableIndex(parameterReaderIndex))).get();
}
void calcAxes (const State& state,LeontisWesthofBondRow myLeontisWesthofBondRow,int residueNumber1,int residueNumber2,int i,int j,Vec3 & xAxisVector1,Vec3 & yAxisVector1, Vec3 & zAxisVector1,Vec3 & xAxisVector2,Vec3 & yAxisVector2 , Vec3 & zAxisVector2,Vec3 & glycosidicNitrogenAtom1LocationInGround,Vec3 & glycosidicNitrogenAtom2LocationInGround) const {
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
stringstream ss1first;
ss1first<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[1];
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myLeontisWesthofBondRow ="<<myLeontisWesthofBondRow.residue1Atom[0]<< ","<<myLeontisWesthofBondRow.residue1Atom[1]<<endl;//ss1first.str() = "<<ss1first.str()<<endl;
stringstream ss1second;
ss1second<< residueNumber1<<"/"<<(myLeontisWesthofBondRow.residue1Atom[2]);
stringstream ss1c1p ;
ss1c1p<< residueNumber1<<"/C1*";//<<(myLeontisWesthofBondRow.residue1Atom[2]);
stringstream ss2first;
ss2first<< residueNumber2<<"/"<<myLeontisWesthofBondRow.residue2Atom[1];
stringstream ss2second;
ss2second<< residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[2]);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] ss3.str()="<<ss3.str()<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myChain.getAtomIndex(ss3.str()) ="<<myChain.getAtomIndex(ss3.str())<<endl;
glycosidicNitrogenAtom1LocationInGround = myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss3.str()));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] glycosidicNitrogenAtom1LocationInGround ="<<glycosidicNitrogenAtom1LocationInGround<<endl;
glycosidicNitrogenAtom2LocationInGround = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1first.str() = "<<ss1first.str()<<endl;
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1second.str() = "<<ss1second.str()<<endl;
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1c1p.str() = "<<ss1c1p.str()<<endl;
//Vec3 firstRingAtomvector1 = myChain.getAtomLocationInMobilizedBodyFrame(myChain.getAtomIndex(ss1first.str())) - station1;
Vec3 firstRingAtomvector1 = myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss1first.str())) - glycosidicNitrogenAtom1LocationInGround;// a myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss3.str()));
Vec3 c1pRingAtomvector1 = myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss1c1p.str())) - glycosidicNitrogenAtom1LocationInGround;// myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss3.str()));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] c1p RingAtomvector1 =" << c1pRingAtomvector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] firstRingAtomvector1 =" <<firstRingAtomvector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] secondRingAtomvector1 goes from atom =" <<ss1second.str()<<" to "<<ss3.str()<<endl;//secondRingAtomvector1<<endl;
Vec3 secondRingAtomvector1 = myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss1second.str())) - glycosidicNitrogenAtom1LocationInGround;//myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss3.str()));
Vec3 c1pRingAtomvector1norm = c1pRingAtomvector1/c1pRingAtomvector1.norm();
double alpha= asin((c1pRingAtomvector1norm%firstRingAtomvector1).norm()/c1pRingAtomvector1norm.norm()/firstRingAtomvector1.norm());
double beta= asin((c1pRingAtomvector1norm%secondRingAtomvector1).norm()/c1pRingAtomvector1norm.norm()/secondRingAtomvector1.norm());
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] alpha = "<<alpha<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] beta = "<<beta<<endl;
double A = sin(beta)*c1pRingAtomvector1norm.norm()/sin(180/Rad2Deg-alpha-beta)/firstRingAtomvector1.norm();
double B = sin(alpha)*c1pRingAtomvector1norm.norm()/sin(180/Rad2Deg-alpha-beta)/secondRingAtomvector1.norm();
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] A = "<<A <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] B = "<<B <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] secondRingAtomvector1 =" <<secondRingAtomvector1<<endl;
Vec3 firstRingAtomvector2 = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss2first.str())) -glycosidicNitrogenAtom2LocationInGround;// myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
Vec3 secondRingAtomvector2 = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss2second.str())) - glycosidicNitrogenAtom2LocationInGround;//myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
//Vec3 xAxisVector1
//Vec3 xAxisVector2;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] Residue 1 type = "<< myLeontisWesthofBondRow.pdbResidueName1 <<endl;
if ((myLeontisWesthofBondRow.pdbResidueName1.compare("A ") == 0) || (myLeontisWesthofBondRow.pdbResidueName1.compare("G ") == 0) ) { //if purine
xAxisVector1 = -5.88327 * firstRingAtomvector1 - 6.13617 * secondRingAtomvector1;
//xAxisVector1 = -5.90016805717817 * firstRingAtomvector1 - 6.13958334585491 * secondRingAtomvector1;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] Residue type should be = "<< myLeontisWesthofBondRow.pdbResidueName1 <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] xAxisVector1 = "<<xAxisVector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] should match = "<< i <<endl;
}
else if ((myLeontisWesthofBondRow.pdbResidueName1.compare("C ") == 0))
xAxisVector1 = -7.83435 * firstRingAtomvector1 -6.99265 *secondRingAtomvector1;
//xAxisVector1 = -7.85229031618972 * firstRingAtomvector1 -6.99951315477877 *secondRingAtomvector1;
else if ((myLeontisWesthofBondRow.pdbResidueName1.compare("U ")) == 0)
xAxisVector1 = -7.3491 * firstRingAtomvector1 -6.47606 *secondRingAtomvector1;
else { cout <<"[TwoTransformForces.h] Unrecognized residue type"<<endl; assert(0);} // trap errors
if ((myLeontisWesthofBondRow.pdbResidueName2.compare("A ") == 0) || (myLeontisWesthofBondRow.pdbResidueName2.compare("G ") == 0)) //if purine
xAxisVector2 = -5.88327 * firstRingAtomvector2 -6.13617 *secondRingAtomvector2;
else if ((myLeontisWesthofBondRow.pdbResidueName2.compare("C ") == 0))
xAxisVector2 = -7.83435 * firstRingAtomvector2 -6.99265 *secondRingAtomvector2;
else if ((myLeontisWesthofBondRow.pdbResidueName2.compare("U ")) == 0)
xAxisVector2 = -7.3491 * firstRingAtomvector2 -6.47606 *secondRingAtomvector2;
else { cout <<"[TwoTransformForces.h] Unrecognized residue type"<<endl; assert(0);} // trap errors
//xAxisVector1 = xAxisVector1/xAxisVector1.norm();
//xAxisVector2 = xAxisVector2/xAxisVector2.norm();
zAxisVector1 = (firstRingAtomvector1%secondRingAtomvector1);
zAxisVector1 = zAxisVector1/zAxisVector1.norm();
zAxisVector2 = (firstRingAtomvector2%secondRingAtomvector2);
zAxisVector2 = zAxisVector2/zAxisVector2.norm();
yAxisVector1 = zAxisVector1%xAxisVector1;
yAxisVector1= yAxisVector1/yAxisVector1.norm();
yAxisVector2 = zAxisVector2%xAxisVector2;
yAxisVector2= yAxisVector2/yAxisVector2.norm();
}
int isThisATwoTransformForce(string myBPEdge) const { // string const & myBPEdge) {
if (((myBPEdge).compare("WatsonCrick") == 0) ||
((myBPEdge).compare("Hoogsteen" ) == 0) ||
((myBPEdge).compare("SugarEdge" ) == 0) ||
((myBPEdge).compare("Stacking" ) == 0) ||
((myBPEdge).compare("AntiParallelStacking") == 0) ||
((myBPEdge).compare("HelicalStacking") == 0) ||
((myBPEdge).compare("Custom" ) == 0) )
return 1;
else
return 0;
}
void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,
Vector_<Vec3>& particleForces, Vector& mobilityForces) const
{
//cout<<"[TwoTransformForces.h] about to call getParameterReader(state)"<<endl;
const ParameterReader & myParameterReader = getParameterReader(state);
//cout<<"[TwoTransformForces.h] done with call getParameterReader(state)"<<endl;
MobilizedBody body1;
MobilizedBody body2;
Transform transform1;
Transform transform2;
float forceConstant;
float torqueConstant;
float dutyCycle; //must be between 0 and 1. at 1, force is applied all the time. at 0, basically never applied.
double scrubberPeriod;
float cutoffRadius;
for (int r=0;r<myParameterReader.numBasePairs;r++)
/*if (((myParameterReader.firstBPEdge[r]).compare("Parallelness") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("ChiBond") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("Rigid") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphere") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphereSmall") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("PointToPlane") != 0)) */
/*if (((myParameterReader.firstBPEdge[r]).compare("WatsonCrick") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Hoogsteen" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("SugarEdge" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Stacking" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("HelicalStacking") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Custom" ) == 0) )*/
//string temp ; temp = myParameterReader.firstBPEdge[r];
if (isThisATwoTransformForce(myParameterReader.firstBPEdge[r]) )
//if (isThisATwoTransformForce(temp))//myParameterReader.firstBPEdge[r]) )
{
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] doing base pair #"<<r<<endl;
string chainId1=myParameterReader.firstBPChain[r];
string chainId2=myParameterReader.secondBPChain[r];
string bondingEdge1=myParameterReader.firstBPEdge[r];
string bondingEdge2=myParameterReader.secondBPEdge[r];
string glycosidicBondOrientation=myParameterReader.orientationBP[r] ;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] bondingEdge1, bondingEdge2: "<<bondingEdge1<<","<< bondingEdge2<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] chainId1, chainId2: "<<chainId1<<","<< chainId2<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.firstBPResidue[r], myParameterReader.secondBPResidue[r]: "<<myParameterReader.firstBPResidue[r]<<","<< myParameterReader.secondBPResidue[r]<<endl;
//cout<<"[TwoTransformForces.h] check 1 "<<r<<endl;
int i = -11111;
int j = -11111;
for (int k = 0; k < myParameterReader.numChains; k++) {
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId1.c_str()) == 0) i = k;
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId2.c_str()) == 0) j = k;
}
assert(i>=0);
assert(j>=0);
int residueNumber1=myParameterReader.firstBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str());
int residueNumber2=myParameterReader.secondBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str());
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] chainId1 myParameterReader.firstBPResidue[r] residueNumber1 chainId2 myParameterReader.secondBPResidue[r] residueNumber2 bondingEdge1 bondingEdge2 "<<chainId1<<","<< myParameterReader.firstBPResidue[r]<<","<< residueNumber1 <<","<< chainId2<<","<< myParameterReader.secondBPResidue[r]<<","<< residueNumber2<<","<< bondingEdge1 <<","<< bondingEdge2 <<endl;
string myPdbResidueName1 = (myChain[i].updResidue(Compound::Index(residueNumber1))).getPdbResidueName();
string myPdbResidueName2 = (myChain[j].updResidue(Compound::Index(residueNumber2))).getPdbResidueName();
//cout<<"[TwoTransformForces.h] check 1.5 "<<r<<endl;
LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myLeontisWesthofBondRow = "<<myLeontisWesthofBondRow.residue1Atom[1]<<"," << myLeontisWesthofBondRow.residue1Atom[2]<<","<<myLeontisWesthofBondRow.residue2Atom[1] <<endl;
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
//body1 = matter.updMobilizedBody((MobilizedBodyIndex (1)));
//body1 = matter.updMobilizedBody((myChain[i]).getAtomMobilizedBodyIndex(Compound::AtomIndex(1)));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] getting body1 from string :"<<ss3.str()<<endl;
body1 = matter.updMobilizedBody(myChain[i].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[i].getAtomIndex(ss3.str()))));
body2 = matter.updMobilizedBody(myChain[j].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[j].getAtomIndex(ss4.str()))));
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = Vec3(0);
//cout<<"[TwoTransformForces.h] station2 ="<<station2<<endl;
Rotation rotation1 = Rotation(myLeontisWesthofBondRow.rotationAngle,myLeontisWesthofBondRow.rotationAxis);
Rotation rotation2;
rotation2.setRotationToIdentityMatrix();
forceConstant = myLeontisWesthofBondRow.springConstant[0];
torqueConstant = myLeontisWesthofBondRow.torqueConstant;
dutyCycle = myParameterReader.dutyCycle;
scrubberPeriod = myParameterReader.scrubberPeriod;
cutoffRadius = myParameterReader.cutoffRadius;
transform1 = Transform(rotation1, station1);
transform2 = Transform(rotation2, station2);
if ((int)(fmod(state.getTime(),scrubberPeriod)/scrubberPeriod+dutyCycle))
{
// put in conditionals to set these dependent on whetehr they are purines or pyrimidines..
//purines
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] getting XGB1 from string :"<<ss3.str()<<endl;
//cout<<"TwoTransformForces.h] X_GB1.R() ="<<X_GB1.R()<<endl;
//cout<<"TwoTransformForces.h] X_GB1.T() ="<<X_GB1.T()<<endl;
Vec3 xAxisVector1 ;
Vec3 yAxisVector1;
Vec3 zAxisVector1;
Vec3 xAxisVector2;
Vec3 yAxisVector2;
Vec3 zAxisVector2;
Vec3 glycosidicNitrogenAtom1LocationInGround;
Vec3 glycosidicNitrogenAtom2LocationInGround;
calcAxes(state,myLeontisWesthofBondRow,residueNumber1,residueNumber2,i,j,xAxisVector1,yAxisVector1,zAxisVector1,xAxisVector2,yAxisVector2,zAxisVector2,glycosidicNitrogenAtom1LocationInGround,glycosidicNitrogenAtom2LocationInGround);
Rotation rotation1FromRingAtoms(Mat33(xAxisVector1,yAxisVector1,zAxisVector1));
Rotation rotation2FromRingAtoms(Mat33(xAxisVector2,yAxisVector2,zAxisVector2));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] computed rotation for residue 1 based on ring atoms to Mat33 = "<<rotation1FromRingAtoms<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] computed rotation for residue 1 based on two-vector constructor = "<<Rotation(SimTK::UnitVec3(xAxisVector1),CoordinateAxis::XCoordinateAxis(),zAxisVector1, CoordinateAxis::ZCoordinateAxis())<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] check 1 old version of X_GB1 = "<< matter.getMobilizedBody(body1).getBodyTransform(state)<<endl;;
//const Transform X_GB1 = matter.getMobilizedBody(body1).getBodyTransform(state);
const Transform X_GB1(rotation1FromRingAtoms,glycosidicNitrogenAtom1LocationInGround);
const Transform X_GB2(rotation2FromRingAtoms,glycosidicNitrogenAtom2LocationInGround);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] check 4 new version of X_GB2 = "<< X_GB2<<endl;
//const Transform& X_GB2 = matter.getMobilizedBody(body2).getBodyTransform(state);
/*const Transform& X_GB2 = (
Rotation(SimTK::UnitVec3(xAxisVector2),CoordinateAxis::XCoordinateAxis(),zAxisVector2, CoordinateAxis::ZCoordinateAxis()),
(matter.getMobilizedBody(body2).getBodyTransform(state)).T()
);
*/
//Vec3 station1 = transform1.T();
//Vec3 station2 = transform2.T();
//cout<<"[TwoTransformLinearSpring] transform1, transform2 :"<<transform1<<","<<transform2<<endl;
//Rotation rotation1 = transform1.R();
//Rotation rotation2 = transform2.R();
//cout<<"[TwoTransformForces.h] station1,station2:"<<station1<<","<<station2<<endl;
const Vec3 s1_G = X_GB1.R() * station1;
const Vec3 s2_G = X_GB2.R() * station2;
const Rotation rot1_G = X_GB1.R() * rotation1;
const Rotation rot2_G = X_GB2.R() * rotation2;
const Vec3 p1_G = X_GB1.T() + s1_G; // station measured from ground origin
const Vec3 p2_G = X_GB2.T() + s2_G;
const Vec3 r_G = p2_G - p1_G; // vector from point1 to point2
const Real d = r_G.norm(); // distance between the points
Real stretch = d ; //d - x0; // + -> tension, - -> compression
//if (stretch > cutoffRadius) stretch = 0;//cutoffRadius ; //limits the force to a constant after given distance.
//if (stretch < cutoffRadius) stretch = 1000000;//cutoffRadius ; //limits the force to a constant after given distance.
//const Real frcScalar = stretch*forceConstant;
const Vec4 rotationAngleAxis = (rot1_G*(~rot2_G)).convertRotationToAngleAxis();
Vec3 torque;
float scale = -2* 1000;
float xoffset = -1.2;
float xmultiplier = .4;
float yoffset = 3.5;
float myFrcScalar;
float A, B, C;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.potentialType :"<<myParameterReader.potentialType<<endl;
for (int i = 0; i<3; i++) torque[i] = rotationAngleAxis[i+1];
torque *= -rotationAngleAxis[0] * torqueConstant ;
if ((myParameterReader.potentialType).compare("Harmonic" ) == 0)
{if (stretch > cutoffRadius) stretch = 0;myFrcScalar = forceConstant*stretch; }
else if ((myParameterReader.potentialType).compare("HarmonicLinear" ) == 0)
{if (stretch > cutoffRadius) stretch = cutoffRadius;myFrcScalar = forceConstant*stretch; }
else if ((myParameterReader.potentialType).compare("Inverse") == 0)
{if (stretch < cutoffRadius) {myFrcScalar = 0;}
else {float C = forceConstant*cutoffRadius; myFrcScalar = -C/stretch/stretch; } // C should be negatve
}
else if ((myParameterReader.potentialType).compare("HarmonicInverse") == 0)
{if (stretch < cutoffRadius)
{
C = forceConstant*cutoffRadius;
A = -forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
myFrcScalar = 2*A*stretch;
torque *= (A*stretch*stretch + B) /forceConstant; //*(-torqueConstant)
}
else {
C = forceConstant*cutoffRadius; myFrcScalar = -C/stretch/stretch;
A = -forceConstant/2/cutoffRadius/cutoffRadius;
B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
torque *= (C/stretch)/forceConstant;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" stretch, myFrcScalar: "<<stretch<<","<<myFrcScalar<<endl;
}
else if ((myParameterReader.potentialType).compare("HarmonicInverseCube") == 0)
{if (stretch < cutoffRadius)
{
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius;
A = -3*forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
//float B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
myFrcScalar = 2*A*stretch;
}
else {
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius; myFrcScalar = -3*C/stretch/stretch/stretch/stretch;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" stretch, myFrcScalar: "<<stretch<<","<<myFrcScalar<<endl;
} else {assert (0);}
myFrcScalar *= myParameterReader.twoTransformForceMultiplier;
torque *= myParameterReader.twoTransformForceMultiplier;
const Real frcScalar = myFrcScalar;
const Vec3 f1_G = (frcScalar/d) * r_G;
if (myParameterReader.verbose) cout<<__FILE__<<":"<<__LINE__<<" : frcScalar = "<<frcScalar<<endl;
//const Vec4 rotationAngleAxis = (transform1.R()*(~(transform2.R()))).convertRotationToAngleAxis();
//const Vec4 rotationAngleAxis = (transform2.R()*(~(transform1.R()))).convertRotationToAngleAxis();
//if (( stretch > cutoffRadius ) || (fabs(rotationAngleAxis[0] ) > (90 /Rad2Deg))) { torque =0; }
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] d, forceConstant, stretch, f1_G, r_G,frcScalar ,torque :"<<d<<","<<forceConstant<<","<<stretch<<","<<f1_G<<","<<r_G<<","<<frcScalar<<","<< torque<<endl;
//cout<<"[TwoTransformForces.h] s1_G , s2_G :"<<s1_G<<","<<s2_G <<endl;
//cout<<"[TwoTransformForces.h] check 5 "<<r<<endl;
if (myParameterReader.verbose) {cout<<"[TwoTransformForces.h] torque on body 1 ="<<torque + s1_G % f1_G<<endl;}
if (myParameterReader.verbose) {cout<<"[TwoTransformForces.h] torque on body 2 ="<<-torque - s2_G % f1_G<<endl;}
bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + (-(matter.getMobilizedBody(body1).getBodyTransform(state)).T()+p1_G) % f1_G, f1_G);
//bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + s1_G % f1_G, f1_G);
bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + (-(matter.getMobilizedBody(body2).getBodyTransform(state)).T()+p2_G) % f1_G, f1_G);
//bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + s2_G % f1_G, f1_G);
//bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + 0*s1_G % f1_G, f1_G);
//bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + 0*s2_G % f1_G, f1_G);
//cout<<"[TwoTransformForces.h] new body1.getMobilizedBodyIndex() bodyForces [body1.getMobilizedBodyIndex()] , r "<<body1.getMobilizedBodyIndex()<<" , " <<bodyForces[body1.getMobilizedBodyIndex()]<<" , "<<r<<endl;
//cout<<"[TwoTransformForces.h] new body2.getMobilizedBodyIndex() bodyForces [body2.getMobilizedBodyIndex()] , r "<<body2.getMobilizedBodyIndex()<<" , " <<bodyForces[body2.getMobilizedBodyIndex()]<<" , "<<r<<endl;
}
//cout<<"[TwoTransformForces.h] check 7 "<<r<<endl;
}
}
Real calcPotentialEnergy(const State& state) const {
double energy = 0.0;
MobilizedBody body1;
MobilizedBody body2;
Transform transform1;
Transform transform2;
float forceConstant;
float torqueConstant;
float dutyCycle; //must be between 0 and 1. at 1, force is applied all the time. at 0, basically never applied.
double scrubberPeriod;
float cutoffRadius;
for (int r=0;r<myParameterReader.numBasePairs;r++)
if (isThisATwoTransformForce(myParameterReader.firstBPEdge[r]))
/*if (((myParameterReader.firstBPEdge[r]).compare("WatsonCrick") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Hoogsteen" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("SugarEdge" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Stacking" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("HelicalStacking") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Custom" ) == 0) )*/
/*if (((myParameterReader.firstBPEdge[r]).compare("Parallelness") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("ChiBond") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("Rigid") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphere") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphereSmall") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("PointToPlane") != 0)) */
//for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++)
{
//const ParameterReader & myParameterReader = getParameterReader(state);
string chainId1=myParameterReader.firstBPChain[r];
string chainId2=myParameterReader.secondBPChain[r];
string bondingEdge1=myParameterReader.firstBPEdge[r];
string bondingEdge2=myParameterReader.secondBPEdge[r];
string glycosidicBondOrientation=myParameterReader.orientationBP[r] ;
int i = -11111;
int j = -11111;
for (int k = 0; k < myParameterReader.numChains; k++) {
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId1.c_str()) == 0) i = k;
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId2.c_str()) == 0) j = k;
}
assert (i>=0);
//SimTK_ASSERT_ALWAYS(i>=0);
assert(j>=0);
int residueNumber1=myParameterReader.firstBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str());
int residueNumber2=myParameterReader.secondBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str());
string myPdbResidueName1 = (myChain[i].updResidue(Compound::Index(residueNumber1))).getPdbResidueName();
string myPdbResidueName2 = (myChain[j].updResidue(Compound::Index(residueNumber2))).getPdbResidueName();
LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
//LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
//cout<<matter.updMobilizedBody(MobilizedBodyIndex(1))<<endl;
body1 = matter.updMobilizedBody((myChain[i]).getAtomMobilizedBodyIndex(Compound::AtomIndex((myChain[i]).getAtomIndex(ss3.str()))));
body2 = matter.updMobilizedBody(myChain[j].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[j].getAtomIndex(ss4.str()))));
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = myChain[j].getAtomLocationInMobilizedBodyFrame(myChain[j].getAtomIndex(ss4.str()));
Rotation rotation1 = Rotation(myLeontisWesthofBondRow.rotationAngle,myLeontisWesthofBondRow.rotationAxis);
Rotation rotation2;
rotation2.setRotationToIdentityMatrix();
forceConstant = myLeontisWesthofBondRow.springConstant[0];
torqueConstant = myLeontisWesthofBondRow.torqueConstant;
dutyCycle = myParameterReader.dutyCycle;
scrubberPeriod = myParameterReader.scrubberPeriod;
cutoffRadius = myParameterReader.cutoffRadius;
transform1 = Transform(rotation1, station1);
transform2 = Transform(rotation2, station2);
//const MobilizedBody& body1 = matter.getMobilizedBody(i);
//for (MobilizedBodyIndex j(0); j < i; j++)
{
//const MobilizedBody& body2 = matter.getMobilizedBody(j);
Vec3 r = body1.getBodyOriginLocation(state)-
body2.getBodyOriginLocation(state);
// still haven't added the torsional part of the energy.
float scale = -2;
float xoffset = -1.2;
float xmultiplier = .4;
float yoffset = 3.5;
//energy += scale/((r.norm()*10-xoffset)*xmultiplier)+yoffset;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r.norm(): "<<r.norm()<<","<<endl;
//scf
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = Vec3(0);
Vec3 xAxisVector1 ;
Vec3 yAxisVector1;
Vec3 zAxisVector1;
Vec3 xAxisVector2;
Vec3 yAxisVector2;
Vec3 zAxisVector2;
Vec3 glycosidicNitrogenAtom1LocationInGround;
Vec3 glycosidicNitrogenAtom2LocationInGround;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" calling calcAxes for residueNumber1,residueNumber2,i,j: "<<residueNumber1<<","<<residueNumber2<<","<<i<<","<< j<<endl;
calcAxes(state,myLeontisWesthofBondRow,residueNumber1,residueNumber2,i,j,xAxisVector1,yAxisVector1,zAxisVector1,xAxisVector2,yAxisVector2,zAxisVector2,glycosidicNitrogenAtom1LocationInGround,glycosidicNitrogenAtom2LocationInGround);
Rotation rotation1FromRingAtoms(Mat33(xAxisVector1,yAxisVector1,zAxisVector1));
Rotation rotation2FromRingAtoms(Mat33(xAxisVector2,yAxisVector2,zAxisVector2));
const Transform X_GB1(rotation1FromRingAtoms,glycosidicNitrogenAtom1LocationInGround);
const Transform X_GB2(rotation2FromRingAtoms,glycosidicNitrogenAtom2LocationInGround);
const Vec3 s1_G = X_GB1.R() * station1;
const Vec3 s2_G = X_GB2.R() * station2;
const Vec3 p1_G = X_GB1.T() + s1_G; // station measured from ground origin
const Vec3 p2_G = X_GB2.T() + s2_G;
const Vec3 r_G = p2_G - p1_G; // vector from point1 to point2
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy :: "<<energy <<","<<endl;
float A, B, C;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.potentialType :"<<myParameterReader.potentialType<<endl;
if ((myParameterReader.potentialType).compare("Harmonic" ) == 0)
{if (r_G.norm() > cutoffRadius) energy += 0;
else energy += forceConstant*(-cutoffRadius*cutoffRadius + r_G.norm()*r_G.norm())/2*myParameterReader.twoTransformForceMultiplier;}
else if ((myParameterReader.potentialType).compare("HarmonicLinear" ) == 0)
{if (r_G.norm() > cutoffRadius) energy += (r_G.norm()-cutoffRadius)*forceConstant*myParameterReader.twoTransformForceMultiplier;
else energy += forceConstant*(-cutoffRadius*cutoffRadius + r_G.norm()*r_G.norm())/2*myParameterReader.twoTransformForceMultiplier;}
else if ((myParameterReader.potentialType).compare("Inverse") == 0)
{if (r_G.norm() <= cutoffRadius) energy += C/cutoffRadius*myParameterReader.twoTransformForceMultiplier;
else {float C = forceConstant*cutoffRadius; energy += C/r_G.norm()*myParameterReader.twoTransformForceMultiplier; } // C should be negatve
}
else if ((myParameterReader.potentialType).compare("HarmonicInverse") == 0)
{if (r_G.norm() <= cutoffRadius)
{
C = forceConstant*cutoffRadius;
A = -forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
float B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
energy += (B+ A*r_G.norm()*r_G.norm())*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" small separation, just added this much to energy:"<< B+ A*r_G.norm()*r_G.norm() <<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
}
else {
C = forceConstant*cutoffRadius; energy += C/r_G.norm()*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" large separation, energy is now:"<< energy <<endl;
//if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" large separation, just added this much to energy:"<< C/r_G.norm()<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
} //else assert (0);//SimTK_ERRCHK1_ALWAYS(0,"[TwoTransformForces.h]","unrecognized potentialType: %s",myParameterReader.potentialType);
else if ((myParameterReader.potentialType).compare("HarmonicInverseCube") == 0)
{if (r_G.norm() < cutoffRadius)
{
A = -3*forceConstant/2/cutoffRadius/cutoffRadius;
float B = 5/2*forceConstant;
energy += (B+ A*r_G.norm()*r_G.norm())*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
}
else {
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius; // myFrcScalar = -3*C/r_G.norm()/r_G.norm()/r_G.norm()/r_G.norm();
energy += C/r_G.norm()/r_G.norm()/r_G.norm() *myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
}
//if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(), myFrcScalar: "<<r_G.norm()<<","<<myFrcScalar<<endl;
} else assert(0);
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy : "<<energy <<","<<endl;
/*if (r.norm() > cutoffRadius)
energy += forceConstant/r.norm(); //scale/((r.norm()*10-xoffset)*xmultiplier)+yoffset;
else
energy += forceConstant/cutoffRadius;*/
//cout <<"[TwoTransformForces.h] distance ="<<r.norm()<<endl;
//energy -= forceConstant*.5*cutoffRadius*cutoffRadius;
//cout <<"[TwoTransformForces.h] energy after translational part="<<energy<<endl;
// torsional part
Rotation rotation1 = transform1.R();
Rotation rotation2 = transform2.R();
const Rotation rot1_G = X_GB1.R() * rotation1;
const Rotation rot2_G = X_GB2.R() * rotation2;
const Vec4 rotationAngleAxis = (rot1_G*(~rot2_G)).convertRotationToAngleAxis();
/*Vec3 torque;
for (int i = 0; i<3; i++) torque[i] = rotationAngleAxis[i+1];
torque = -rotationAngleAxis[0] * torqueConstant * torque;*/
if (myParameterReader.verbose) printf ("r.norm(),energy:%f,%f ", r.norm(),energy);
//scf
if (r_G.norm() <= cutoffRadius) {energy += 0.5*((rotationAngleAxis[0]*rotationAngleAxis[0]-180*180/Rad2Deg/Rad2Deg)*torqueConstant)*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" rotationAngleAxis[0] = "<<rotationAngleAxis[0]<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" calculating torsional energy. r_G.norm(),cutoffRadius: "<<r_G.norm()<<","<<cutoffRadius<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy : "<<energy <<","<<endl;
} // add nothing if r.norm()>cutoffRadius
else
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" Not Calculating torsional energy. r_G.norm(),cutoffRadius: "<<r_G.norm()<<","<<cutoffRadius<<endl;
if (myParameterReader.verbose) printf ("r.norm(),energy:%f,%f ", r.norm(),energy);
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r.norm(),energy: "<<r.norm()<<","<<energy<<endl;
//
//cout <<"[TwoTransformForces.h] energy after torsional part ="<<energy<<endl;
}
}
return energy;
}
bool dependsOnlyOnPositions() const {
return true;
}
} ;
/*
*/
When the Force's calcForce method is called it should apply the appropriate force to each atom's body. If the kinematics are coarse grained the force will be applied to the origin of the body frame, and so you would have to add a torque to compensate. But I'm assuming for a first cut you will be doing an atomistic simulation.
RNABuilder's TwoTransformForces.h does something much like this, except it queries a list of interactions rather than atoms. I am appending it at the very bottom of this email. It's an overcomplicated example though. You really ought to start from scratch. Peter documented the Custom Force class pretty well if I recall correctly.
Are you working on protein or RNA?
I submitted this to the forum.
Sam
On Jun 9, 2009, at 2:36 PM, Balaraman, Gouthaman wrote:
Hey Sherm,
The force module is initialized using a pdb, psf, and parameter file. This code is part of a normal Cartesian MD code. I have written a wrapper such that, after the initialization, I can pass a set of coordinates to the function, and get back the cartesian forces.
I believe if I use the Molmodels objects, that would make the clustering part very simple. So I would definitely like to use it, if it is possible. In the future, I would also like to experiment with different types of coarse graining. In that case, can I still use Molmodel's objects to extend to arbitrary clustering?
Thanks,
Goutham
-----Original Message-----
From: Michael Sherman [mailto:msherman@stanford.edu]
Sent: Tue 6/9/2009 2:24 PM
To: Balaraman, Gouthaman
Cc: 'Christopher Bruns'; 'Samuel Coulbourn Flores'
Subject: RE: custom forces in Molmodel
Hi, Goutham.
How does your force field know the bonded structure of the molecule, and how
does that relate to the model built by Molmodel? Are you planning to use
Molmodel to construct the mechanical model of the molecules (via Molmodel's
Protein() or RNA() objects, for example)?
Regards,
Sherm
From: Balaraman, Gouthaman [mailto:GBalaraman@coh.org]
Sent: Tuesday, June 09, 2009 12:40 PM
To: Michael Sherman
Cc: Christopher Bruns; Samuel Coulbourn Flores
Subject: RE: custom forces in Molmodel
Hey All,
I have an existing force field that I want to use as is. This force module
will take cartesian positions as input, and return a force vector in
cartesian coordinates. The force module will take care of reading force
parameters, and evaluating the forces. I want to incorporate this force
engine into molmodel. I want to start with all torsions for the clustering.
In the longer run, I would like to define my own coarse graining as well.
I really appreciate your help on this aspect.
Thank you,
Goutham
-----Original Message-----
From: Michael Sherman [mailto:msherman@stanford.edu]
Sent: Mon 6/8/2009 8:12 PM
To: Balaraman, Gouthaman
Cc: 'Christopher Bruns'; Samuel Coulbourn Flores
Subject: RE: custom forces in Molmodel
Hi, Goutham.
Nice to hear from you again and I'm very glad you're interested in using
Molmodel. I'm including Chris Bruns and Sam Flores on this email - Chris is
the primary author of Molmodel, and Sam is the most advanced Molmodel user
at Stanford. We will need some more information in order to figure out the
best way for you to get your force field working on a system that has been
modeled using Molmodel.
We are about to release version 1.6 of Molmodel in which DuMM has been made
more flexible (and it now uses OpenMM for acceleration when it can also).
Depending on the functional form of your force field, DuMM may be able to
accommodate it (DuMM is neutral about the parameterization). Otherwise there
are a variety of ways to apply different forces - Sam has used Molmodel
successfully with several special purpose RNA force fields. Any combination
is plausible - you can avoid DuMM altogether, you can use it for part of the
force field and something else for other parts, or you can use DuMM directly
(maybe with some hacks).
Please tell us what you have in mind. Do you have an existing force field
that you want to use as-is? What are the functional forms involved? What
mechanical model are you using (all torsions, or something coarser?)
Chris, should we move this discussion to one of the Molmodel forums on
simtk.org?
Best regards,
Sherm
From: Balaraman, Gouthaman [mailto:GBalaraman@coh.org]
Sent: Monday, June 08, 2009 5:42 PM
To: msherman@stanford.edu
Subject: custom forces in Molmodel
Hey Mike,
Thanks for getting back on my question. I am a postdoc at Beckman Institute
of City of Hope in Prof Nagarajan Vaidehi's group. I had a chance to meet
you in an OpenMM session.
I want to use my own force engine, in Molmodel. Can you give me any pointers
on how I can incorporate that within Molmodel?
Thanks,
Goutham
You have been contacted by a member of Simtk.org.
From: Michael Sherman
Email: msherman@stanford.edu
The e-mail was sent via the "contact" link on Simtk.org by the above member.
It is not screened by Simtk.org or its staff.
*** MESSAGE BEGINS ***
Hi, Goutham.
Ayman Habib forwarded me your post to the OpenSim forum. Please send me your
actual email address so I can connect you with the right people here to talk
about custom forces in Molmodel.
Best regards,
Michael Sherman (Sherm)
---------------------------------------------------------------------
SECURITY/CONFIDENTIALITY WARNING: This message and any attachments are
intended solely for the individual or entity to which they are addressed.
This communication may contain information that is privileged, confidential,
or exempt from disclosure under applicable law (e.g., personal health
information, research data, financial information). Because this e-mail has
been sent without encryption, individuals other than the intended recipient
may be able to view the information, forward it to others or tamper with the
information without the knowledge or consent of the sender. If you are not
the intended recipient, or the employee or person responsible for delivering
the message to the intended recipient, any dissemination, distribution or
copying of the communication is strictly prohibited. If you received the
communication in error, please notify the sender immediately by replying to
this message and deleting the message and any accompanying files from your
system. If, due to the security risks, you do not wish to receive further
communications via e-mail, please reply to this message and inform the
sender that you do not wish to receive further e-mail from the sender.
---------------------------------------------------------------------
Samuel Coulbourn Flores, PhD
Altman Lab
650.644.8416
花山
科学者
生物工学部
スタンフォ一ド大学
スタンフォ一ド、カリフォルニア、米国
using namespace SimTK;
using namespace std;
std::ostream& operator<<(std::ostream& o, const ParameterReader&) {
assert(false);
return o;
}
class AllTwoTransformLinearSprings : public Force::Custom::Implementation {
private:
SimbodyMatterSubsystem& matter;
ParameterReader& myParameterReader;
LeontisWesthofClass& myLeontisWesthofClass;
Biopolymer * myChain;
mutable int parameterReaderIndex;
public:
AllTwoTransformLinearSprings (SimbodyMatterSubsystem& matter,ParameterReader& myParameterReader, LeontisWesthofClass& myLeontisWesthofClass, Biopolymer * myChain ) : matter(matter),myParameterReader(myParameterReader), myLeontisWesthofClass (myLeontisWesthofClass), myChain(myChain)
{
}
void realizeTopology(State& state) const {
parameterReaderIndex = state.allocateDiscreteVariable(matter.getMySubsystemIndex(),
Stage::Dynamics, new Value<ParameterReader>(myParameterReader));
//return 0;
//implementation->realizeTopology(state);
}
void setParameterReader(State& state, const ParameterReader& myParameterReader) const {
static_cast<Value<ParameterReader>&>(state.updDiscreteVariable( matter.getMySubsystemIndex(),
SimTK::DiscreteVariableIndex(parameterReaderIndex))).upd() = myParameterReader;
}
const ParameterReader & getParameterReader(const State& state) const {
return static_cast<const Value<ParameterReader>&>(state.getDiscreteVariable(matter.getMySubsystemIndex(),
SimTK::DiscreteVariableIndex(parameterReaderIndex))).get();
}
void calcAxes (const State& state,LeontisWesthofBondRow myLeontisWesthofBondRow,int residueNumber1,int residueNumber2,int i,int j,Vec3 & xAxisVector1,Vec3 & yAxisVector1, Vec3 & zAxisVector1,Vec3 & xAxisVector2,Vec3 & yAxisVector2 , Vec3 & zAxisVector2,Vec3 & glycosidicNitrogenAtom1LocationInGround,Vec3 & glycosidicNitrogenAtom2LocationInGround) const {
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
stringstream ss1first;
ss1first<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[1];
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myLeontisWesthofBondRow ="<<myLeontisWesthofBondRow.residue1Atom[0]<< ","<<myLeontisWesthofBondRow.residue1Atom[1]<<endl;//ss1first.str() = "<<ss1first.str()<<endl;
stringstream ss1second;
ss1second<< residueNumber1<<"/"<<(myLeontisWesthofBondRow.residue1Atom[2]);
stringstream ss1c1p ;
ss1c1p<< residueNumber1<<"/C1*";//<<(myLeontisWesthofBondRow.residue1Atom[2]);
stringstream ss2first;
ss2first<< residueNumber2<<"/"<<myLeontisWesthofBondRow.residue2Atom[1];
stringstream ss2second;
ss2second<< residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[2]);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] ss3.str()="<<ss3.str()<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myChain.getAtomIndex(ss3.str()) ="<<myChain.getAtomIndex(ss3.str())<<endl;
glycosidicNitrogenAtom1LocationInGround = myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss3.str()));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] glycosidicNitrogenAtom1LocationInGround ="<<glycosidicNitrogenAtom1LocationInGround<<endl;
glycosidicNitrogenAtom2LocationInGround = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1first.str() = "<<ss1first.str()<<endl;
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1second.str() = "<<ss1second.str()<<endl;
//if (myParameterReader.verbose) cout<<"TwoTransformForces.h] ss1c1p.str() = "<<ss1c1p.str()<<endl;
//Vec3 firstRingAtomvector1 = myChain.getAtomLocationInMobilizedBodyFrame(myChain.getAtomIndex(ss1first.str())) - station1;
Vec3 firstRingAtomvector1 = myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss1first.str())) - glycosidicNitrogenAtom1LocationInGround;// a myChain.calcAtomLocationInGroundFrame(state,myChain.getAtomIndex(ss3.str()));
Vec3 c1pRingAtomvector1 = myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss1c1p.str())) - glycosidicNitrogenAtom1LocationInGround;// myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss3.str()));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] c1p RingAtomvector1 =" << c1pRingAtomvector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] firstRingAtomvector1 =" <<firstRingAtomvector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] secondRingAtomvector1 goes from atom =" <<ss1second.str()<<" to "<<ss3.str()<<endl;//secondRingAtomvector1<<endl;
Vec3 secondRingAtomvector1 = myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss1second.str())) - glycosidicNitrogenAtom1LocationInGround;//myChain[i].calcAtomLocationInGroundFrame(state,myChain[i].getAtomIndex(ss3.str()));
Vec3 c1pRingAtomvector1norm = c1pRingAtomvector1/c1pRingAtomvector1.norm();
double alpha= asin((c1pRingAtomvector1norm%firstRingAtomvector1).norm()/c1pRingAtomvector1norm.norm()/firstRingAtomvector1.norm());
double beta= asin((c1pRingAtomvector1norm%secondRingAtomvector1).norm()/c1pRingAtomvector1norm.norm()/secondRingAtomvector1.norm());
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] alpha = "<<alpha<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] beta = "<<beta<<endl;
double A = sin(beta)*c1pRingAtomvector1norm.norm()/sin(180/Rad2Deg-alpha-beta)/firstRingAtomvector1.norm();
double B = sin(alpha)*c1pRingAtomvector1norm.norm()/sin(180/Rad2Deg-alpha-beta)/secondRingAtomvector1.norm();
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] A = "<<A <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] B = "<<B <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] secondRingAtomvector1 =" <<secondRingAtomvector1<<endl;
Vec3 firstRingAtomvector2 = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss2first.str())) -glycosidicNitrogenAtom2LocationInGround;// myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
Vec3 secondRingAtomvector2 = myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss2second.str())) - glycosidicNitrogenAtom2LocationInGround;//myChain[j].calcAtomLocationInGroundFrame(state,myChain[j].getAtomIndex(ss4.str()));
//Vec3 xAxisVector1
//Vec3 xAxisVector2;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] Residue 1 type = "<< myLeontisWesthofBondRow.pdbResidueName1 <<endl;
if ((myLeontisWesthofBondRow.pdbResidueName1.compare("A ") == 0) || (myLeontisWesthofBondRow.pdbResidueName1.compare("G ") == 0) ) { //if purine
xAxisVector1 = -5.88327 * firstRingAtomvector1 - 6.13617 * secondRingAtomvector1;
//xAxisVector1 = -5.90016805717817 * firstRingAtomvector1 - 6.13958334585491 * secondRingAtomvector1;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] Residue type should be = "<< myLeontisWesthofBondRow.pdbResidueName1 <<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] xAxisVector1 = "<<xAxisVector1<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] should match = "<< i <<endl;
}
else if ((myLeontisWesthofBondRow.pdbResidueName1.compare("C ") == 0))
xAxisVector1 = -7.83435 * firstRingAtomvector1 -6.99265 *secondRingAtomvector1;
//xAxisVector1 = -7.85229031618972 * firstRingAtomvector1 -6.99951315477877 *secondRingAtomvector1;
else if ((myLeontisWesthofBondRow.pdbResidueName1.compare("U ")) == 0)
xAxisVector1 = -7.3491 * firstRingAtomvector1 -6.47606 *secondRingAtomvector1;
else { cout <<"[TwoTransformForces.h] Unrecognized residue type"<<endl; assert(0);} // trap errors
if ((myLeontisWesthofBondRow.pdbResidueName2.compare("A ") == 0) || (myLeontisWesthofBondRow.pdbResidueName2.compare("G ") == 0)) //if purine
xAxisVector2 = -5.88327 * firstRingAtomvector2 -6.13617 *secondRingAtomvector2;
else if ((myLeontisWesthofBondRow.pdbResidueName2.compare("C ") == 0))
xAxisVector2 = -7.83435 * firstRingAtomvector2 -6.99265 *secondRingAtomvector2;
else if ((myLeontisWesthofBondRow.pdbResidueName2.compare("U ")) == 0)
xAxisVector2 = -7.3491 * firstRingAtomvector2 -6.47606 *secondRingAtomvector2;
else { cout <<"[TwoTransformForces.h] Unrecognized residue type"<<endl; assert(0);} // trap errors
//xAxisVector1 = xAxisVector1/xAxisVector1.norm();
//xAxisVector2 = xAxisVector2/xAxisVector2.norm();
zAxisVector1 = (firstRingAtomvector1%secondRingAtomvector1);
zAxisVector1 = zAxisVector1/zAxisVector1.norm();
zAxisVector2 = (firstRingAtomvector2%secondRingAtomvector2);
zAxisVector2 = zAxisVector2/zAxisVector2.norm();
yAxisVector1 = zAxisVector1%xAxisVector1;
yAxisVector1= yAxisVector1/yAxisVector1.norm();
yAxisVector2 = zAxisVector2%xAxisVector2;
yAxisVector2= yAxisVector2/yAxisVector2.norm();
}
int isThisATwoTransformForce(string myBPEdge) const { // string const & myBPEdge) {
if (((myBPEdge).compare("WatsonCrick") == 0) ||
((myBPEdge).compare("Hoogsteen" ) == 0) ||
((myBPEdge).compare("SugarEdge" ) == 0) ||
((myBPEdge).compare("Stacking" ) == 0) ||
((myBPEdge).compare("AntiParallelStacking") == 0) ||
((myBPEdge).compare("HelicalStacking") == 0) ||
((myBPEdge).compare("Custom" ) == 0) )
return 1;
else
return 0;
}
void calcForce(const State& state, Vector_<SpatialVec>& bodyForces,
Vector_<Vec3>& particleForces, Vector& mobilityForces) const
{
//cout<<"[TwoTransformForces.h] about to call getParameterReader(state)"<<endl;
const ParameterReader & myParameterReader = getParameterReader(state);
//cout<<"[TwoTransformForces.h] done with call getParameterReader(state)"<<endl;
MobilizedBody body1;
MobilizedBody body2;
Transform transform1;
Transform transform2;
float forceConstant;
float torqueConstant;
float dutyCycle; //must be between 0 and 1. at 1, force is applied all the time. at 0, basically never applied.
double scrubberPeriod;
float cutoffRadius;
for (int r=0;r<myParameterReader.numBasePairs;r++)
/*if (((myParameterReader.firstBPEdge[r]).compare("Parallelness") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("ChiBond") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("Rigid") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphere") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphereSmall") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("PointToPlane") != 0)) */
/*if (((myParameterReader.firstBPEdge[r]).compare("WatsonCrick") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Hoogsteen" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("SugarEdge" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Stacking" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("HelicalStacking") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Custom" ) == 0) )*/
//string temp ; temp = myParameterReader.firstBPEdge[r];
if (isThisATwoTransformForce(myParameterReader.firstBPEdge[r]) )
//if (isThisATwoTransformForce(temp))//myParameterReader.firstBPEdge[r]) )
{
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] doing base pair #"<<r<<endl;
string chainId1=myParameterReader.firstBPChain[r];
string chainId2=myParameterReader.secondBPChain[r];
string bondingEdge1=myParameterReader.firstBPEdge[r];
string bondingEdge2=myParameterReader.secondBPEdge[r];
string glycosidicBondOrientation=myParameterReader.orientationBP[r] ;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] bondingEdge1, bondingEdge2: "<<bondingEdge1<<","<< bondingEdge2<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] chainId1, chainId2: "<<chainId1<<","<< chainId2<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.firstBPResidue[r], myParameterReader.secondBPResidue[r]: "<<myParameterReader.firstBPResidue[r]<<","<< myParameterReader.secondBPResidue[r]<<endl;
//cout<<"[TwoTransformForces.h] check 1 "<<r<<endl;
int i = -11111;
int j = -11111;
for (int k = 0; k < myParameterReader.numChains; k++) {
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId1.c_str()) == 0) i = k;
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId2.c_str()) == 0) j = k;
}
assert(i>=0);
assert(j>=0);
int residueNumber1=myParameterReader.firstBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str());
int residueNumber2=myParameterReader.secondBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str());
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] chainId1 myParameterReader.firstBPResidue[r] residueNumber1 chainId2 myParameterReader.secondBPResidue[r] residueNumber2 bondingEdge1 bondingEdge2 "<<chainId1<<","<< myParameterReader.firstBPResidue[r]<<","<< residueNumber1 <<","<< chainId2<<","<< myParameterReader.secondBPResidue[r]<<","<< residueNumber2<<","<< bondingEdge1 <<","<< bondingEdge2 <<endl;
string myPdbResidueName1 = (myChain[i].updResidue(Compound::Index(residueNumber1))).getPdbResidueName();
string myPdbResidueName2 = (myChain[j].updResidue(Compound::Index(residueNumber2))).getPdbResidueName();
//cout<<"[TwoTransformForces.h] check 1.5 "<<r<<endl;
LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myLeontisWesthofBondRow = "<<myLeontisWesthofBondRow.residue1Atom[1]<<"," << myLeontisWesthofBondRow.residue1Atom[2]<<","<<myLeontisWesthofBondRow.residue2Atom[1] <<endl;
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
//body1 = matter.updMobilizedBody((MobilizedBodyIndex (1)));
//body1 = matter.updMobilizedBody((myChain[i]).getAtomMobilizedBodyIndex(Compound::AtomIndex(1)));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] getting body1 from string :"<<ss3.str()<<endl;
body1 = matter.updMobilizedBody(myChain[i].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[i].getAtomIndex(ss3.str()))));
body2 = matter.updMobilizedBody(myChain[j].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[j].getAtomIndex(ss4.str()))));
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = Vec3(0);
//cout<<"[TwoTransformForces.h] station2 ="<<station2<<endl;
Rotation rotation1 = Rotation(myLeontisWesthofBondRow.rotationAngle,myLeontisWesthofBondRow.rotationAxis);
Rotation rotation2;
rotation2.setRotationToIdentityMatrix();
forceConstant = myLeontisWesthofBondRow.springConstant[0];
torqueConstant = myLeontisWesthofBondRow.torqueConstant;
dutyCycle = myParameterReader.dutyCycle;
scrubberPeriod = myParameterReader.scrubberPeriod;
cutoffRadius = myParameterReader.cutoffRadius;
transform1 = Transform(rotation1, station1);
transform2 = Transform(rotation2, station2);
if ((int)(fmod(state.getTime(),scrubberPeriod)/scrubberPeriod+dutyCycle))
{
// put in conditionals to set these dependent on whetehr they are purines or pyrimidines..
//purines
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] getting XGB1 from string :"<<ss3.str()<<endl;
//cout<<"TwoTransformForces.h] X_GB1.R() ="<<X_GB1.R()<<endl;
//cout<<"TwoTransformForces.h] X_GB1.T() ="<<X_GB1.T()<<endl;
Vec3 xAxisVector1 ;
Vec3 yAxisVector1;
Vec3 zAxisVector1;
Vec3 xAxisVector2;
Vec3 yAxisVector2;
Vec3 zAxisVector2;
Vec3 glycosidicNitrogenAtom1LocationInGround;
Vec3 glycosidicNitrogenAtom2LocationInGround;
calcAxes(state,myLeontisWesthofBondRow,residueNumber1,residueNumber2,i,j,xAxisVector1,yAxisVector1,zAxisVector1,xAxisVector2,yAxisVector2,zAxisVector2,glycosidicNitrogenAtom1LocationInGround,glycosidicNitrogenAtom2LocationInGround);
Rotation rotation1FromRingAtoms(Mat33(xAxisVector1,yAxisVector1,zAxisVector1));
Rotation rotation2FromRingAtoms(Mat33(xAxisVector2,yAxisVector2,zAxisVector2));
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] computed rotation for residue 1 based on ring atoms to Mat33 = "<<rotation1FromRingAtoms<<endl;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] computed rotation for residue 1 based on two-vector constructor = "<<Rotation(SimTK::UnitVec3(xAxisVector1),CoordinateAxis::XCoordinateAxis(),zAxisVector1, CoordinateAxis::ZCoordinateAxis())<<endl;
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] check 1 old version of X_GB1 = "<< matter.getMobilizedBody(body1).getBodyTransform(state)<<endl;;
//const Transform X_GB1 = matter.getMobilizedBody(body1).getBodyTransform(state);
const Transform X_GB1(rotation1FromRingAtoms,glycosidicNitrogenAtom1LocationInGround);
const Transform X_GB2(rotation2FromRingAtoms,glycosidicNitrogenAtom2LocationInGround);
//if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] check 4 new version of X_GB2 = "<< X_GB2<<endl;
//const Transform& X_GB2 = matter.getMobilizedBody(body2).getBodyTransform(state);
/*const Transform& X_GB2 = (
Rotation(SimTK::UnitVec3(xAxisVector2),CoordinateAxis::XCoordinateAxis(),zAxisVector2, CoordinateAxis::ZCoordinateAxis()),
(matter.getMobilizedBody(body2).getBodyTransform(state)).T()
);
*/
//Vec3 station1 = transform1.T();
//Vec3 station2 = transform2.T();
//cout<<"[TwoTransformLinearSpring] transform1, transform2 :"<<transform1<<","<<transform2<<endl;
//Rotation rotation1 = transform1.R();
//Rotation rotation2 = transform2.R();
//cout<<"[TwoTransformForces.h] station1,station2:"<<station1<<","<<station2<<endl;
const Vec3 s1_G = X_GB1.R() * station1;
const Vec3 s2_G = X_GB2.R() * station2;
const Rotation rot1_G = X_GB1.R() * rotation1;
const Rotation rot2_G = X_GB2.R() * rotation2;
const Vec3 p1_G = X_GB1.T() + s1_G; // station measured from ground origin
const Vec3 p2_G = X_GB2.T() + s2_G;
const Vec3 r_G = p2_G - p1_G; // vector from point1 to point2
const Real d = r_G.norm(); // distance between the points
Real stretch = d ; //d - x0; // + -> tension, - -> compression
//if (stretch > cutoffRadius) stretch = 0;//cutoffRadius ; //limits the force to a constant after given distance.
//if (stretch < cutoffRadius) stretch = 1000000;//cutoffRadius ; //limits the force to a constant after given distance.
//const Real frcScalar = stretch*forceConstant;
const Vec4 rotationAngleAxis = (rot1_G*(~rot2_G)).convertRotationToAngleAxis();
Vec3 torque;
float scale = -2* 1000;
float xoffset = -1.2;
float xmultiplier = .4;
float yoffset = 3.5;
float myFrcScalar;
float A, B, C;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.potentialType :"<<myParameterReader.potentialType<<endl;
for (int i = 0; i<3; i++) torque[i] = rotationAngleAxis[i+1];
torque *= -rotationAngleAxis[0] * torqueConstant ;
if ((myParameterReader.potentialType).compare("Harmonic" ) == 0)
{if (stretch > cutoffRadius) stretch = 0;myFrcScalar = forceConstant*stretch; }
else if ((myParameterReader.potentialType).compare("HarmonicLinear" ) == 0)
{if (stretch > cutoffRadius) stretch = cutoffRadius;myFrcScalar = forceConstant*stretch; }
else if ((myParameterReader.potentialType).compare("Inverse") == 0)
{if (stretch < cutoffRadius) {myFrcScalar = 0;}
else {float C = forceConstant*cutoffRadius; myFrcScalar = -C/stretch/stretch; } // C should be negatve
}
else if ((myParameterReader.potentialType).compare("HarmonicInverse") == 0)
{if (stretch < cutoffRadius)
{
C = forceConstant*cutoffRadius;
A = -forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
myFrcScalar = 2*A*stretch;
torque *= (A*stretch*stretch + B) /forceConstant; //*(-torqueConstant)
}
else {
C = forceConstant*cutoffRadius; myFrcScalar = -C/stretch/stretch;
A = -forceConstant/2/cutoffRadius/cutoffRadius;
B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
torque *= (C/stretch)/forceConstant;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" stretch, myFrcScalar: "<<stretch<<","<<myFrcScalar<<endl;
}
else if ((myParameterReader.potentialType).compare("HarmonicInverseCube") == 0)
{if (stretch < cutoffRadius)
{
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius;
A = -3*forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
//float B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
myFrcScalar = 2*A*stretch;
}
else {
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius; myFrcScalar = -3*C/stretch/stretch/stretch/stretch;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" stretch, myFrcScalar: "<<stretch<<","<<myFrcScalar<<endl;
} else {assert (0);}
myFrcScalar *= myParameterReader.twoTransformForceMultiplier;
torque *= myParameterReader.twoTransformForceMultiplier;
const Real frcScalar = myFrcScalar;
const Vec3 f1_G = (frcScalar/d) * r_G;
if (myParameterReader.verbose) cout<<__FILE__<<":"<<__LINE__<<" : frcScalar = "<<frcScalar<<endl;
//const Vec4 rotationAngleAxis = (transform1.R()*(~(transform2.R()))).convertRotationToAngleAxis();
//const Vec4 rotationAngleAxis = (transform2.R()*(~(transform1.R()))).convertRotationToAngleAxis();
//if (( stretch > cutoffRadius ) || (fabs(rotationAngleAxis[0] ) > (90 /Rad2Deg))) { torque =0; }
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] d, forceConstant, stretch, f1_G, r_G,frcScalar ,torque :"<<d<<","<<forceConstant<<","<<stretch<<","<<f1_G<<","<<r_G<<","<<frcScalar<<","<< torque<<endl;
//cout<<"[TwoTransformForces.h] s1_G , s2_G :"<<s1_G<<","<<s2_G <<endl;
//cout<<"[TwoTransformForces.h] check 5 "<<r<<endl;
if (myParameterReader.verbose) {cout<<"[TwoTransformForces.h] torque on body 1 ="<<torque + s1_G % f1_G<<endl;}
if (myParameterReader.verbose) {cout<<"[TwoTransformForces.h] torque on body 2 ="<<-torque - s2_G % f1_G<<endl;}
bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + (-(matter.getMobilizedBody(body1).getBodyTransform(state)).T()+p1_G) % f1_G, f1_G);
//bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + s1_G % f1_G, f1_G);
bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + (-(matter.getMobilizedBody(body2).getBodyTransform(state)).T()+p2_G) % f1_G, f1_G);
//bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + s2_G % f1_G, f1_G);
//bodyForces[body1.getMobilizedBodyIndex()] += SpatialVec(torque + 0*s1_G % f1_G, f1_G);
//bodyForces[body2.getMobilizedBodyIndex()] -= SpatialVec(torque + 0*s2_G % f1_G, f1_G);
//cout<<"[TwoTransformForces.h] new body1.getMobilizedBodyIndex() bodyForces [body1.getMobilizedBodyIndex()] , r "<<body1.getMobilizedBodyIndex()<<" , " <<bodyForces[body1.getMobilizedBodyIndex()]<<" , "<<r<<endl;
//cout<<"[TwoTransformForces.h] new body2.getMobilizedBodyIndex() bodyForces [body2.getMobilizedBodyIndex()] , r "<<body2.getMobilizedBodyIndex()<<" , " <<bodyForces[body2.getMobilizedBodyIndex()]<<" , "<<r<<endl;
}
//cout<<"[TwoTransformForces.h] check 7 "<<r<<endl;
}
}
Real calcPotentialEnergy(const State& state) const {
double energy = 0.0;
MobilizedBody body1;
MobilizedBody body2;
Transform transform1;
Transform transform2;
float forceConstant;
float torqueConstant;
float dutyCycle; //must be between 0 and 1. at 1, force is applied all the time. at 0, basically never applied.
double scrubberPeriod;
float cutoffRadius;
for (int r=0;r<myParameterReader.numBasePairs;r++)
if (isThisATwoTransformForce(myParameterReader.firstBPEdge[r]))
/*if (((myParameterReader.firstBPEdge[r]).compare("WatsonCrick") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Hoogsteen" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("SugarEdge" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Stacking" ) == 0) ||
((myParameterReader.firstBPEdge[r]).compare("HelicalStacking") == 0) ||
((myParameterReader.firstBPEdge[r]).compare("Custom" ) == 0) )*/
/*if (((myParameterReader.firstBPEdge[r]).compare("Parallelness") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("ChiBond") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("Rigid") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphere") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("HardSphereSmall") != 0) &&
((myParameterReader.firstBPEdge[r]).compare("PointToPlane") != 0)) */
//for (MobilizedBodyIndex i(0); i < matter.getNBodies(); i++)
{
//const ParameterReader & myParameterReader = getParameterReader(state);
string chainId1=myParameterReader.firstBPChain[r];
string chainId2=myParameterReader.secondBPChain[r];
string bondingEdge1=myParameterReader.firstBPEdge[r];
string bondingEdge2=myParameterReader.secondBPEdge[r];
string glycosidicBondOrientation=myParameterReader.orientationBP[r] ;
int i = -11111;
int j = -11111;
for (int k = 0; k < myParameterReader.numChains; k++) {
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId1.c_str()) == 0) i = k;
if (strcmp((myParameterReader.chainId[k]).c_str(),chainId2.c_str()) == 0) j = k;
}
assert (i>=0);
//SimTK_ASSERT_ALWAYS(i>=0);
assert(j>=0);
int residueNumber1=myParameterReader.firstBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[i]).c_str());
int residueNumber2=myParameterReader.secondBPResidue[r]-myParameterReader.getFirstResidueNumbers((myParameterReader.chainId[j]).c_str());
string myPdbResidueName1 = (myChain[i].updResidue(Compound::Index(residueNumber1))).getPdbResidueName();
string myPdbResidueName2 = (myChain[j].updResidue(Compound::Index(residueNumber2))).getPdbResidueName();
LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
stringstream ss3;
ss3<<residueNumber1<<"/"<<myLeontisWesthofBondRow.residue1Atom[0];
stringstream ss4;
ss4<<residueNumber2<<"/"<<(myLeontisWesthofBondRow.residue2Atom[0]);
//LeontisWesthofBondRow myLeontisWesthofBondRow = myLeontisWesthofClass.getLeontisWesthofBondRow(myPdbResidueName1,bondingEdge1,myPdbResidueName2,bondingEdge2,glycosidicBondOrientation);
//cout<<matter.updMobilizedBody(MobilizedBodyIndex(1))<<endl;
body1 = matter.updMobilizedBody((myChain[i]).getAtomMobilizedBodyIndex(Compound::AtomIndex((myChain[i]).getAtomIndex(ss3.str()))));
body2 = matter.updMobilizedBody(myChain[j].getAtomMobilizedBodyIndex(Compound::AtomIndex(myChain[j].getAtomIndex(ss4.str()))));
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = myChain[j].getAtomLocationInMobilizedBodyFrame(myChain[j].getAtomIndex(ss4.str()));
Rotation rotation1 = Rotation(myLeontisWesthofBondRow.rotationAngle,myLeontisWesthofBondRow.rotationAxis);
Rotation rotation2;
rotation2.setRotationToIdentityMatrix();
forceConstant = myLeontisWesthofBondRow.springConstant[0];
torqueConstant = myLeontisWesthofBondRow.torqueConstant;
dutyCycle = myParameterReader.dutyCycle;
scrubberPeriod = myParameterReader.scrubberPeriod;
cutoffRadius = myParameterReader.cutoffRadius;
transform1 = Transform(rotation1, station1);
transform2 = Transform(rotation2, station2);
//const MobilizedBody& body1 = matter.getMobilizedBody(i);
//for (MobilizedBodyIndex j(0); j < i; j++)
{
//const MobilizedBody& body2 = matter.getMobilizedBody(j);
Vec3 r = body1.getBodyOriginLocation(state)-
body2.getBodyOriginLocation(state);
// still haven't added the torsional part of the energy.
float scale = -2;
float xoffset = -1.2;
float xmultiplier = .4;
float yoffset = 3.5;
//energy += scale/((r.norm()*10-xoffset)*xmultiplier)+yoffset;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r.norm(): "<<r.norm()<<","<<endl;
//scf
Vec3 station1 = myLeontisWesthofBondRow.attachmentPoint;
Vec3 station2 = Vec3(0);
Vec3 xAxisVector1 ;
Vec3 yAxisVector1;
Vec3 zAxisVector1;
Vec3 xAxisVector2;
Vec3 yAxisVector2;
Vec3 zAxisVector2;
Vec3 glycosidicNitrogenAtom1LocationInGround;
Vec3 glycosidicNitrogenAtom2LocationInGround;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" calling calcAxes for residueNumber1,residueNumber2,i,j: "<<residueNumber1<<","<<residueNumber2<<","<<i<<","<< j<<endl;
calcAxes(state,myLeontisWesthofBondRow,residueNumber1,residueNumber2,i,j,xAxisVector1,yAxisVector1,zAxisVector1,xAxisVector2,yAxisVector2,zAxisVector2,glycosidicNitrogenAtom1LocationInGround,glycosidicNitrogenAtom2LocationInGround);
Rotation rotation1FromRingAtoms(Mat33(xAxisVector1,yAxisVector1,zAxisVector1));
Rotation rotation2FromRingAtoms(Mat33(xAxisVector2,yAxisVector2,zAxisVector2));
const Transform X_GB1(rotation1FromRingAtoms,glycosidicNitrogenAtom1LocationInGround);
const Transform X_GB2(rotation2FromRingAtoms,glycosidicNitrogenAtom2LocationInGround);
const Vec3 s1_G = X_GB1.R() * station1;
const Vec3 s2_G = X_GB2.R() * station2;
const Vec3 p1_G = X_GB1.T() + s1_G; // station measured from ground origin
const Vec3 p2_G = X_GB2.T() + s2_G;
const Vec3 r_G = p2_G - p1_G; // vector from point1 to point2
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy :: "<<energy <<","<<endl;
float A, B, C;
if (myParameterReader.verbose) cout<<"[TwoTransformForces.h] myParameterReader.potentialType :"<<myParameterReader.potentialType<<endl;
if ((myParameterReader.potentialType).compare("Harmonic" ) == 0)
{if (r_G.norm() > cutoffRadius) energy += 0;
else energy += forceConstant*(-cutoffRadius*cutoffRadius + r_G.norm()*r_G.norm())/2*myParameterReader.twoTransformForceMultiplier;}
else if ((myParameterReader.potentialType).compare("HarmonicLinear" ) == 0)
{if (r_G.norm() > cutoffRadius) energy += (r_G.norm()-cutoffRadius)*forceConstant*myParameterReader.twoTransformForceMultiplier;
else energy += forceConstant*(-cutoffRadius*cutoffRadius + r_G.norm()*r_G.norm())/2*myParameterReader.twoTransformForceMultiplier;}
else if ((myParameterReader.potentialType).compare("Inverse") == 0)
{if (r_G.norm() <= cutoffRadius) energy += C/cutoffRadius*myParameterReader.twoTransformForceMultiplier;
else {float C = forceConstant*cutoffRadius; energy += C/r_G.norm()*myParameterReader.twoTransformForceMultiplier; } // C should be negatve
}
else if ((myParameterReader.potentialType).compare("HarmonicInverse") == 0)
{if (r_G.norm() <= cutoffRadius)
{
C = forceConstant*cutoffRadius;
A = -forceConstant/2/cutoffRadius/cutoffRadius;// -C/2/cutoffRadius/cutoffRadius/cutoffRadius;
float B = +C/cutoffRadius-A*cutoffRadius*cutoffRadius;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
energy += (B+ A*r_G.norm()*r_G.norm())*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" small separation, just added this much to energy:"<< B+ A*r_G.norm()*r_G.norm() <<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
}
else {
C = forceConstant*cutoffRadius; energy += C/r_G.norm()*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" large separation, energy is now:"<< energy <<endl;
//if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" large separation, just added this much to energy:"<< C/r_G.norm()<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
}
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(): "<<r_G.norm()<<","<<endl;
} //else assert (0);//SimTK_ERRCHK1_ALWAYS(0,"[TwoTransformForces.h]","unrecognized potentialType: %s",myParameterReader.potentialType);
else if ((myParameterReader.potentialType).compare("HarmonicInverseCube") == 0)
{if (r_G.norm() < cutoffRadius)
{
A = -3*forceConstant/2/cutoffRadius/cutoffRadius;
float B = 5/2*forceConstant;
energy += (B+ A*r_G.norm()*r_G.norm())*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" A,C :"<<A<<","<<C<<endl;
}
else {
C = forceConstant*cutoffRadius*cutoffRadius*cutoffRadius; // myFrcScalar = -3*C/r_G.norm()/r_G.norm()/r_G.norm()/r_G.norm();
energy += C/r_G.norm()/r_G.norm()/r_G.norm() *myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" C :"<<C<<endl;
}
//if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r_G.norm(), myFrcScalar: "<<r_G.norm()<<","<<myFrcScalar<<endl;
} else assert(0);
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy : "<<energy <<","<<endl;
/*if (r.norm() > cutoffRadius)
energy += forceConstant/r.norm(); //scale/((r.norm()*10-xoffset)*xmultiplier)+yoffset;
else
energy += forceConstant/cutoffRadius;*/
//cout <<"[TwoTransformForces.h] distance ="<<r.norm()<<endl;
//energy -= forceConstant*.5*cutoffRadius*cutoffRadius;
//cout <<"[TwoTransformForces.h] energy after translational part="<<energy<<endl;
// torsional part
Rotation rotation1 = transform1.R();
Rotation rotation2 = transform2.R();
const Rotation rot1_G = X_GB1.R() * rotation1;
const Rotation rot2_G = X_GB2.R() * rotation2;
const Vec4 rotationAngleAxis = (rot1_G*(~rot2_G)).convertRotationToAngleAxis();
/*Vec3 torque;
for (int i = 0; i<3; i++) torque[i] = rotationAngleAxis[i+1];
torque = -rotationAngleAxis[0] * torqueConstant * torque;*/
if (myParameterReader.verbose) printf ("r.norm(),energy:%f,%f ", r.norm(),energy);
//scf
if (r_G.norm() <= cutoffRadius) {energy += 0.5*((rotationAngleAxis[0]*rotationAngleAxis[0]-180*180/Rad2Deg/Rad2Deg)*torqueConstant)*myParameterReader.twoTransformForceMultiplier;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" rotationAngleAxis[0] = "<<rotationAngleAxis[0]<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" calculating torsional energy. r_G.norm(),cutoffRadius: "<<r_G.norm()<<","<<cutoffRadius<<endl;
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" energy : "<<energy <<","<<endl;
} // add nothing if r.norm()>cutoffRadius
else
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" Not Calculating torsional energy. r_G.norm(),cutoffRadius: "<<r_G.norm()<<","<<cutoffRadius<<endl;
if (myParameterReader.verbose) printf ("r.norm(),energy:%f,%f ", r.norm(),energy);
if (myParameterReader.verbose) cout<< __FILE__<<":"<<__LINE__ <<" r.norm(),energy: "<<r.norm()<<","<<energy<<endl;
//
//cout <<"[TwoTransformForces.h] energy after torsional part ="<<energy<<endl;
}
}
return energy;
}
bool dependsOnlyOnPositions() const {
return true;
}
} ;
/*
*/