1 /**********************************************************************
2 align.cpp - Align two molecules or vectors of vector3
3 
4 Copyright (C) 2010 by Noel M. O'Boyle
5 
6 This file is part of the Open Babel project.
7 For more information, see <http://openbabel.org/>
8 
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation version 2 of the License.
12 
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 GNU General Public License for more details.
17 ***********************************************************************/
18 
19 #include <openbabel/babelconfig.h>
20 
21 #include <vector>
22 #include <climits> // UINT_MAX
23 
24 #include <openbabel/math/align.h>
25 #include <openbabel/atom.h>
26 #include <openbabel/oberror.h>
27 #include <openbabel/obiter.h>
28 #include <openbabel/graphsym.h>
29 #include <openbabel/math/vector3.h>
30 #include <openbabel/elements.h>
31 
32 #include <Eigen/Core>
33 #include <Eigen/Geometry>
34 #include <Eigen/SVD>
35 #include <Eigen/LU>
36 
37 using namespace std;
38 
39 namespace OpenBabel
40 {
41   extern OBMessageHandler obErrorLog;
42 
OBAlign(bool includeH,bool symmetry)43   OBAlign::OBAlign(bool includeH, bool symmetry) : _method(OBAlign::Kabsch)
44   {
45     _ready = false;
46     _symmetry = symmetry;
47     _includeH = includeH;
48     _prefmol = nullptr;
49   }
50 
OBAlign(const vector<vector3> & ref,const vector<vector3> & target)51   OBAlign::OBAlign(const vector<vector3> &ref, const vector<vector3> &target) : _method(OBAlign::Kabsch)
52   {
53     SetRef(ref);
54     SetTarget(target);
55     _symmetry = false;
56     _prefmol = nullptr;
57   }
58 
OBAlign(const OBMol & refmol,const OBMol & targetmol,bool includeH,bool symmetry)59   OBAlign::OBAlign(const OBMol &refmol, const OBMol &targetmol, bool includeH, bool symmetry) : _method(OBAlign::Kabsch)
60   {
61     _symmetry = symmetry;
62     _includeH = includeH;
63     SetRefMol(refmol);
64     SetTargetMol(targetmol);
65   }
66 
VectorsToMatrix(const vector<vector3> * pcoords,Eigen::MatrixXd & coords)67   void OBAlign::VectorsToMatrix(const vector<vector3> *pcoords, Eigen::MatrixXd &coords) {
68 
69     vector<vector3>::size_type N = pcoords->size();
70     coords.resize(3, N);
71 
72     // Create a 3xN matrix of the coords
73     vector<vector3>::const_iterator it;
74     vector<vector3>::size_type colm;
75     for (colm=0,it=pcoords->begin();colm<N;++colm,++it)
76       coords.col(colm) = Eigen::Vector3d( it->AsArray() );
77   }
78 
MoveToOrigin(Eigen::MatrixXd & coords)79   Eigen::Vector3d OBAlign::MoveToOrigin(Eigen::MatrixXd &coords) {
80 
81     vector<vector3>::size_type N = coords.cols();
82 
83     // Find the centroid
84     Eigen::Vector3d centroid;
85     centroid = coords.rowwise().sum() / N;
86 
87     // Subtract the centroids
88     for (vector<vector3>::size_type i=0; i<N; ++i)
89       coords.col(i) -= centroid;
90     return centroid;
91   }
92 
SetRef(const vector<vector3> & ref)93   void OBAlign::SetRef(const vector<vector3> &ref) {
94     _pref = &ref;
95     VectorsToMatrix(_pref, _mref);
96     _ref_centr = MoveToOrigin(_mref);
97 
98     _ready = false;
99   }
100 
SetTarget(const vector<vector3> & target)101   void OBAlign::SetTarget(const vector<vector3> &target) {
102     _ptarget = &target;
103     VectorsToMatrix(_ptarget, _mtarget);
104     _target_centr = MoveToOrigin(_mtarget);
105 
106     _ready = false;
107   }
108 
SetRefMol(const OBMol & refmol)109   void OBAlign::SetRefMol(const OBMol &refmol) {
110     _prefmol = &refmol;
111 
112     // Set up the BitVec for the hydrogens and store the refmol coords
113     _frag_atoms.Clear();
114     _frag_atoms.Resize(refmol.NumAtoms() + 1);
115     _refmol_coords.resize(0);
116     OBAtom* atom;
117     int delta = 1;
118     _newidx.resize(0);
119 
120     for (unsigned int i=1; i<=refmol.NumAtoms(); ++i) {
121       atom = refmol.GetAtom(i);
122       if (_includeH || atom->GetAtomicNum() != OBElements::Hydrogen) {
123         _frag_atoms.SetBitOn(i);
124         _newidx.push_back(i - delta);
125         _refmol_coords.push_back(atom->GetVector());
126       }
127       else {
128         delta++;
129         _newidx.push_back(UINT_MAX);
130       }
131     }
132     SetRef(_refmol_coords);
133 
134     if (_symmetry) {
135       FindAutomorphisms((OBMol*)&refmol, _aut, _frag_atoms);
136     }
137   }
138 
SetTargetMol(const OBMol & targetmol)139   void OBAlign::SetTargetMol(const OBMol &targetmol) {
140     _ptargetmol = &targetmol;
141     _targetmol_coords.resize(0);
142     OBAtom const *atom;
143     for (unsigned int i=1; i<=targetmol.NumAtoms(); ++i) {
144       atom = targetmol.GetAtom(i);
145       if (_includeH || atom->GetAtomicNum() != OBElements::Hydrogen)
146         _targetmol_coords.push_back(atom->GetVector());
147     }
148     SetTarget(_targetmol_coords);
149   }
150 
SetMethod(OBAlign::AlignMethod method)151   void OBAlign::SetMethod(OBAlign::AlignMethod method) {
152     _method = method;
153   }
154 
155 /* Evaluates the Newton-Raphson correction for the Horn quartic.
156    only 11 FLOPs */
eval_horn_NR_corrxn(const vector<double> & c,const double x)157   static double eval_horn_NR_corrxn(const vector<double> &c, const double x)
158   {
159     double x2 = x*x;
160     double b = (x2 + c[2])*x;
161     double a = b + c[1];
162 
163     return((a*x + c[0])/(2.0*x2*x + b + a));
164   }
165 
166   /* Newton-Raphson root finding */
QCProot(const vector<double> & coeff,double guess,const double delta)167   static double QCProot(const vector<double> &coeff, double guess, const double delta)
168   {
169     int             i;
170     double          oldg;
171     double initialg = guess;
172 
173     for (i = 0; i < 50; ++i)
174     {
175         oldg = guess;
176         /* guess -= (eval_horn_quart(coeff, guess) / eval_horn_quart_deriv(coeff, guess)); */
177         guess -= eval_horn_NR_corrxn(coeff, guess);
178 
179         if (fabs(guess - oldg) < fabs(delta*guess))
180             return(guess);
181     }
182 
183     return initialg + 1.0; // Failed to converge!
184   }
185 
CalcQuarticCoeffs(const Eigen::Matrix3d & M)186   vector<double> CalcQuarticCoeffs(const Eigen::Matrix3d &M)
187   {
188     vector<double> coeff(4);
189 
190     double          Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz;
191     double          Szz2, Syy2, Sxx2, Sxy2, Syz2, Sxz2, Syx2, Szy2, Szx2,
192                     SyzSzymSyySzz2, Sxx2Syy2Szz2Syz2Szy2, Sxy2Sxz2Syx2Szx2,
193                     SxzpSzx, SyzpSzy, SxypSyx, SyzmSzy,
194                     SxzmSzx, SxymSyx, SxxpSyy, SxxmSyy;
195 
196 #ifdef HAVE_EIGEN3
197     Eigen::MatrixXd M_sqr = M.array().square();
198 #else
199     Eigen::MatrixXd M_sqr = M.cwise().square();
200 #endif
201 
202     Sxx = M(0, 0);
203     Sxy = M(1, 0);
204     Sxz = M(2, 0);
205     Syx = M(0, 1);
206     Syy = M(1, 1);
207     Syz = M(2, 1);
208     Szx = M(0, 2);
209     Szy = M(1, 2);
210     Szz = M(2, 2);
211 
212     Sxx2 = Sxx * Sxx;
213     Syy2 = Syy * Syy;
214     Szz2 = Szz * Szz;
215 
216     Sxy2 = Sxy * Sxy;
217     Syz2 = Syz * Syz;
218     Sxz2 = Sxz * Sxz;
219 
220     Syx2 = Syx * Syx;
221     Szy2 = Szy * Szy;
222     Szx2 = Szx * Szx;
223 
224     SyzSzymSyySzz2 = 2.0*(Syz*Szy - Syy*Szz);
225     Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2;
226 
227     /* coeff[4] = 1.0; */
228     /* coeff[3] = 0.0; */
229     // coeff[2] = -2.0 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2);
230     coeff[2] = -2.0 * M_sqr.sum();
231     coeff[1] = 8.0 * (Sxx*Syz*Szy + Syy*Szx*Sxz + Szz*Sxy*Syx - Sxx*Syy*Szz - Syz*Szx*Sxy - Szy*Syx*Sxz);
232 
233     SxzpSzx = Sxz+Szx;
234     SyzpSzy = Syz+Szy;
235     SxypSyx = Sxy+Syx;
236     SyzmSzy = Syz-Szy;
237     SxzmSzx = Sxz-Szx;
238     SxymSyx = Sxy-Syx;
239     SxxpSyy = Sxx+Syy;
240     SxxmSyy = Sxx-Syy;
241     Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2;
242 
243     coeff[0] = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2
244              + (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2)
245              + (-(SxzpSzx)*(SyzmSzy)+(SxymSyx)*(SxxmSyy-Szz)) * (-(SxzmSzx)*(SyzpSzy)+(SxymSyx)*(SxxmSyy+Szz))
246              + (-(SxzpSzx)*(SyzpSzy)-(SxypSyx)*(SxxpSyy-Szz)) * (-(SxzmSzx)*(SyzmSzy)-(SxypSyx)*(SxxpSyy+Szz))
247              + (+(SxypSyx)*(SyzpSzy)+(SxzpSzx)*(SxxmSyy+Szz)) * (-(SxymSyx)*(SyzmSzy)+(SxzpSzx)*(SxxpSyy+Szz))
248              + (+(SxypSyx)*(SyzmSzy)+(SxzmSzx)*(SxxmSyy-Szz)) * (-(SxymSyx)*(SyzpSzy)+(SxzmSzx)*(SxxpSyy-Szz));
249 
250     return coeff;
251   }
252 
TheobaldAlign(const Eigen::MatrixXd & mtarget)253   void OBAlign::TheobaldAlign(const Eigen::MatrixXd &mtarget)
254   {
255     // M = B(t) times A (where A, B are N x 3 matrices)
256     Eigen::Matrix3d M = mtarget * _mref.transpose();
257 
258     // Maximum value for lambda is (Ga + Gb) / 2
259     double innerprod = mtarget.squaredNorm() + _mref.squaredNorm();
260 
261     vector<double> coeffs = CalcQuarticCoeffs(M);
262     double lambdamax = QCProot(coeffs, 0.5 * innerprod, 1e-6);
263     if (lambdamax > (0.5 * innerprod))
264       _fail = true;
265     else {
266       double sqrdev = innerprod - (2.0 * lambdamax);
267       _rmsd = sqrt(sqrdev / mtarget.cols());
268     }
269   }
270 
SimpleAlign(const Eigen::MatrixXd & mtarget)271   void OBAlign::SimpleAlign(const Eigen::MatrixXd &mtarget)
272   {
273     // Covariance matrix C = X times Y(t)
274     Eigen::Matrix3d C = _mref * mtarget.transpose();
275 
276     // Singular Value Decomposition of C into USV(t)
277 #ifdef HAVE_EIGEN3
278     Eigen::JacobiSVD<Eigen::Matrix3d> svd(C, Eigen::ComputeFullU | Eigen::ComputeFullV);
279 #else
280     Eigen::SVD<Eigen::Matrix3d> svd(C);
281 #endif
282 
283     // Prepare matrix T
284     double sign = (C.determinant() > 0) ? 1. : -1.; // Sign of determinant
285     Eigen::Matrix3d T = Eigen::Matrix3d::Identity();
286     T(2,2) = sign;
287 
288     // Optimal rotation matrix, U, is V T U(t)
289     _rotMatrix = svd.matrixV() * T * svd.matrixU().transpose();
290 
291     // Rotate target using rotMatrix
292     _result = _rotMatrix.transpose() * mtarget;
293 
294     Eigen::MatrixXd deviation = _result - _mref;
295 #ifdef HAVE_EIGEN3
296     Eigen::MatrixXd sqr = deviation.array().square();
297 #else
298     Eigen::MatrixXd sqr = deviation.cwise().square();
299 #endif
300     double sum = sqr.sum();
301     _rmsd = sqrt( sum / sqr.cols() );
302 
303   }
304 
Align()305   bool OBAlign::Align()
306   {
307     vector<vector3>::size_type N = _ptarget->size();
308 
309     if (_pref->size() != N) {
310       obErrorLog.ThrowError(__FUNCTION__, "Cannot align the reference and target as they are of different size" , obError);
311       return false;
312     }
313 
314     if (!_symmetry || _aut.size() == 1) {
315       if (_method == OBAlign::Kabsch)
316         SimpleAlign(_mtarget);
317       else
318         TheobaldAlign(_mtarget);
319     }
320     else {  // Iterate over the automorphisms
321 
322       // ...for storing the results from the lowest rmsd to date
323       double min_rmsd = DBL_MAX;
324       Eigen::MatrixXd result, rotMatrix;
325 
326       // Try all of the symmetry-allowed permutations
327       OBIsomorphismMapper::Mappings::const_iterator cit;
328       Eigen::MatrixXd mtarget(_mtarget.rows(), _mtarget.cols());
329 
330       for (unsigned int k = 0; k < _aut.size(); ++k) {
331         // Rearrange columns of _mtarget for this permutation
332         unsigned int i=0;
333         for (unsigned int j=1; j<=_prefmol->NumAtoms(); ++j) {
334           if (_frag_atoms.BitIsSet(j)) {
335             for (std::size_t l = 0; l < _aut[k].size(); ++l)
336               if (_aut[k][l].first == j - 1) {
337                 mtarget.col(i) = _mtarget.col(_newidx[_aut[k][l].second]);
338                 break;
339               }
340             i++;
341           }
342         }
343         if (_method == OBAlign::Kabsch)
344           SimpleAlign(mtarget);
345         else
346           TheobaldAlign(mtarget);
347         if (_rmsd < min_rmsd) {
348           min_rmsd = _rmsd;
349           result = _result;
350           rotMatrix = _rotMatrix;
351         }
352       }
353 
354       // Restore the best answer from memory
355       _rmsd = min_rmsd;
356       _result = result;
357       _rotMatrix = rotMatrix;
358     }
359 
360     _ready = true;
361     return true;
362   }
363 
GetRotMatrix()364   matrix3x3 OBAlign::GetRotMatrix()
365   {
366     if (!_ready) {
367       obErrorLog.ThrowError(__FUNCTION__, "Rotation matrix not available until you call Align()" , obError);
368       return matrix3x3();
369     }
370 
371     // Convert Eigen::Matrix to matrix3x3
372     double rot[3][3];
373     for (int row=0; row<3; ++row)
374        for (int col=0; col<3; ++col)
375          rot[col][row] = _rotMatrix(row, col); // Return in form suitable for use in expressions like "result *= rotMatrix";
376     matrix3x3 rotmat = matrix3x3(rot);
377 
378     return rotmat;
379   }
380 
GetAlignment()381   vector<vector3> OBAlign::GetAlignment() {
382     vector<vector3> aligned_coords;
383     if (!_ready) {
384       obErrorLog.ThrowError(__FUNCTION__, "Alignment not available until you call Align()" , obError);
385       return aligned_coords;
386     }
387 
388     if (!_prefmol || _includeH) {
389       // Add back the centroid of the reference and convert to vv3
390       Eigen::Vector3d tmp;
391       aligned_coords.reserve(_result.cols());
392       for (int i=0; i<_result.cols(); ++i) {
393         tmp = _result.col(i) + _ref_centr;
394         aligned_coords.push_back(vector3(tmp(0), tmp(1), tmp(2)));
395       }
396     }
397     else { // Need to deal with the case where hydrogens were excluded
398       vector<vector3> target_coords;
399       for (unsigned int i=1; i<=_ptargetmol->NumAtoms(); ++i)
400         target_coords.push_back(_ptargetmol->GetAtom(i)->GetVector());
401       Eigen::MatrixXd mtarget;
402       VectorsToMatrix(&target_coords, mtarget);
403 
404       // Subtract the centroid of the non-H atoms
405       for (unsigned int i=0; i<mtarget.cols(); ++i)
406         mtarget.col(i) -= _target_centr;
407 
408       // Rotate
409       Eigen::MatrixXd result = mtarget.transpose() * _rotMatrix;
410       result.transposeInPlace();
411 
412       // Add back the centroid of the reference and convert to vv3
413       Eigen::Vector3d tmp;
414       aligned_coords.reserve(_result.cols());
415       for (int i=0; i<result.cols(); ++i) {
416         tmp = result.col(i) + _ref_centr;
417         aligned_coords.push_back(vector3(tmp(0), tmp(1), tmp(2)));
418       }
419     }
420 
421     return aligned_coords;
422   }
423 
UpdateCoords(OBMol * target)424   bool OBAlign::UpdateCoords(OBMol* target) {
425     if (!_ready) {
426       obErrorLog.ThrowError(__FUNCTION__, "Alignment not available until you call Align()" , obError);
427       return false;
428     }
429 
430     vector<vector3> newcoords = GetAlignment();
431     if (newcoords.size() != target->NumAtoms()) {
432       obErrorLog.ThrowError(__FUNCTION__, "Cannot update the target molecule with the alignment coordinates as they are of different size" , obError);
433       return false;
434     }
435 
436     int i = 0;
437     FOR_ATOMS_OF_MOL(a, *target) {
438       a->SetVector(newcoords.at(i));
439       i++;
440     }
441 
442     return true;
443   }
444 
GetRMSD()445   double OBAlign::GetRMSD() {
446     if (!_ready) {
447       obErrorLog.ThrowError(__FUNCTION__, "RMSD not available until you call Align()" , obError);
448       return nan("");
449     }
450 
451     return _rmsd;
452   }
453 
454 } // namespace OpenBabel
455 
456 //! \file align.cpp
457 //! \brief Handle 3D coordinates.
458