So, I've made a piece of code that (based off the torsion code in
OBMol::SetTorsion) changes the angle between three atoms.
I'm using this to update the angles of an added fragment, based on a
referenced (complete) fragment. I'm having some trouble getting it to work
properly (the angles it updates to are incorrect), and I must have an error
somewhere but I can't identify where it is.
Any help would be appreciated.
void setAngle(OBMol &mol, OBAtom * vertex, OBAtom * b, OBAtom * c, double
ang) {
vector<int> atoms;
double sn, cs, t, x, y, z, tx, ty, tz;
double m[9];
//! locates all atoms for which there exists a path to 'end'
//! without going through 'bgn'
//! children must not include 'end'
// void OBMol::FindChildren(vector<OBAtom*> &children,OBAtom *bgn,OBAtom
*end)
//need to update location of c, also
mol.FindChildren(atoms, vertex->GetIdx(), c->GetIdx());
atoms.push_back(c->GetIdx());
cout << "Updating: " << atoms.size() << " atoms: ";
for (int a = 0; a < atoms.size(); a++) {
cout << describe(mol.GetAtom(atoms.at(a))) << " ";
}
cout << endl;
//used later to convert system to new frame of reference when we apply
rotation matrix
vector3 vertexvec = vertex->GetVector();
vector3 bvec = b->GetVector();
vector3 cvec = c->GetVector();
//need normalized rotation axis to perform rotation.
vector3 rotationAxisU = OBAPI::cross(bvec, cvec);
vector3 rotationAxisN = rotationAxisU.normalize();
// b1/c1 are "true" vectors, used for calculating angle between them.
vector3 b1 = bvec - vertexvec;
vector3 c1 = cvec - vertexvec;
double angle = acos( OBAPI::dot(b1.normalize(), c1.normalize()));
cout << "DOT: " << OBAPI::dot(b1.normalize(), c1.normalize()) << endl;
cout << "current angle: " << angle << " requested: " << ang << endl;
double diffAng = ang - angle;
cout << "diff: " << diffAng << endl;
//set up rotation matrix
sn = sin(diffAng);
cs = cos(diffAng);
t = 1 - cs;
cout << "Debug: " << sn << " " << cs << " " << t << endl;
x = rotationAxisN.x();
y = rotationAxisN.y();
z = rotationAxisN.z();
m[0] = t * x * x + cs;
m[1] = t * x * y + sn * z;
m[2] = t * x * z - sn * y;
m[3] = t * x * y - sn * z;
m[4] = t * y * y + cs;
m[5] = t * y * z + sn * x;
m[6] = t * x * z + sn * y;
m[7] = t * y * z - sn * x;
m[8] = t * z * z + cs;
//rotation matrix set
tx = vertexvec.x();
ty = vertexvec.y();
tz = vertexvec.z();
//move atoms
vector<int>::iterator i;
int j;
for (i = atoms.begin(); i != atoms.end(); ++i) {
j = *i;
OBAtom * atom = mol.GetAtom(j);
vector3 tempvec = atom->GetVector();
vector3 tempvecN = atom->GetVector();
double xtemp = tempvec.x();
double ytemp = tempvec.y();
double ztemp = tempvec.z();
xtemp -= tx;
ytemp -= ty;
ztemp -= tz;
x = xtemp * m[0] + ytemp * m[1] + ztemp * m[2];
y = xtemp * m[3] + ytemp * m[4] + ztemp * m[5];
z = xtemp * m[6] + ytemp * m[7] + ztemp * m[8];
xtemp = x;
ytemp = y;
ztemp = z;
xtemp += tx;
ytemp += ty;
ztemp += tz;
vector3 newvec(xtemp, ytemp, ztemp);
cout << "->updating atom : " << describe(atom) << " from pos : "
<< atom->GetVector() << " to pos: " << newvec
<< endl;
atom->SetVector(newvec);
}
}
--
View this message in context:
http://forums.openbabel.org/OBMol-SetAngle-code-problem-c-tp4657153.html
Sent from the openbabel-devel mailing list archive at Nabble.com.
------------------------------------------------------------------------------
CenturyLink Cloud: The Leader in Enterprise Cloud Services.
Learn Why More Businesses Are Choosing CenturyLink Cloud For
Critical Workloads, Development Environments & Everything In Between.
Get a Quote or Start a Free Trial Today.
http://pubads.g.doubleclick.net/gampad/clk?id=119420431&iu=/4140/ostg.clktrk
_______________________________________________
OpenBabel-Devel mailing list
OpenBabel-Devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-devel
