From da9faf032d61faf8591b6f0ff4431f501d1ce1dd Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sun, 8 Jun 2014 14:07:06 +0100
Subject: [PATCH] Made it so that my RMSD uses code from old RMSD
---
src/colvar/PathMSDBase.cpp | 2 +-
src/reference/OptimalRMSD.cpp | 202 +--------------------------------
src/reference/ReferenceAtoms.h | 7 --
src/reference/SimpleRMSD.cpp | 66 +----------
src/tools/RMSD.cpp | 7 +-
src/tools/RMSD.h | 4 +-
6 files changed, 13 insertions(+), 275 deletions(-)
diff --git a/src/colvar/PathMSDBase.cpp b/src/colvar/PathMSDBase.cpp
index 7537ee824..c1f8f9411 100644
--- a/src/colvar/PathMSDBase.cpp
+++ b/src/colvar/PathMSDBase.cpp
@@ -62,7 +62,7 @@ nframes(0)
bool do_read=true;
while (do_read){
PDB mypdb;
- RMSD mymsd(log);
+ RMSD mymsd;
do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
if(do_read){
unsigned nat=0;
diff --git a/src/reference/OptimalRMSD.cpp b/src/reference/OptimalRMSD.cpp
index d73fe473b..f9a0ef336 100644
--- a/src/reference/OptimalRMSD.cpp
+++ b/src/reference/OptimalRMSD.cpp
@@ -22,22 +22,18 @@
#include "MetricRegister.h"
#include "RMSDBase.h"
#include "tools/Matrix.h"
+#include "tools/RMSD.h"
namespace PLMD{
class OptimalRMSD : public RMSDBase {
private:
bool fast;
+ RMSD myrmsd;
public:
OptimalRMSD(const ReferenceConfigurationOptions& ro);
void read( const PDB& );
double calc( const std::vector<Vector>& pos, const bool& squared );
-
- template <bool safe,bool alEqDis>
- double optimalAlignment(const std::vector<double> & align,
- const std::vector<double> & displace,
- const std::vector<Vector> & positions,
- bool squared=false);
};
PLUMED_REGISTER_METRIC(OptimalRMSD,"OPTIMAL")
@@ -55,198 +51,12 @@ void OptimalRMSD::read( const PDB& pdb ){
double OptimalRMSD::calc( const std::vector<Vector>& pos, const bool& squared ){
if( fast ){
- if( getAlign()==getDisplace() ) return optimalAlignment<false,true>(getAlign(),getDisplace(),pos,squared);
- return optimalAlignment<false,false>(getAlign(),getDisplace(),pos,squared);
+ if( getAlign()==getDisplace() ) return myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ return myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
} else {
- if( getAlign()==getDisplace() ) return optimalAlignment<true,true>(getAlign(),getDisplace(),pos,squared);
- return optimalAlignment<true,false>(getAlign(),getDisplace(),pos,squared);
- }
-}
-
-template <bool safe,bool alEqDis>
-double OptimalRMSD::optimalAlignment(const std::vector<double> & align,
- const std::vector<double> & displace,
- const std::vector<Vector> & positions,
- bool squared){
- double dist(0);
- unsigned n=getNumberOfReferencePositions();
-// This is the trace of positions*positions + reference*reference
- double rr00(0);
- double rr11(0);
-// This is positions*reference
- Tensor rr01; Vector rpos;
-
- Vector cpositions;
-
-// first expensive loop: compute centers
- for(unsigned iat=0;iat<n;iat++){
- double w=align[iat];
- cpositions+=positions[iat]*w;
- }
-
-// second expensive loop: compute second moments wrt centers
- for(unsigned iat=0;iat<n;iat++){
- double w=align[iat]; rpos=getReferencePosition(iat);
- rr00+=dotProduct(positions[iat]-cpositions,positions[iat]-cpositions)*w;
- rr11+=dotProduct(rpos,rpos)*w;
- rr01+=Tensor(positions[iat]-cpositions,rpos)*w;
- }
-
- Matrix<double> m=Matrix<double>(4,4);
- m[0][0]=2.0*(-rr01[0][0]-rr01[1][1]-rr01[2][2]);
- m[1][1]=2.0*(-rr01[0][0]+rr01[1][1]+rr01[2][2]);
- m[2][2]=2.0*(+rr01[0][0]-rr01[1][1]+rr01[2][2]);
- m[3][3]=2.0*(+rr01[0][0]+rr01[1][1]-rr01[2][2]);
- m[0][1]=2.0*(-rr01[1][2]+rr01[2][1]);
- m[0][2]=2.0*(+rr01[0][2]-rr01[2][0]);
- m[0][3]=2.0*(-rr01[0][1]+rr01[1][0]);
- m[1][2]=2.0*(-rr01[0][1]-rr01[1][0]);
- m[1][3]=2.0*(-rr01[0][2]-rr01[2][0]);
- m[2][3]=2.0*(-rr01[1][2]-rr01[2][1]);
- m[1][0] = m[0][1];
- m[2][0] = m[0][2];
- m[2][1] = m[1][2];
- m[3][0] = m[0][3];
- m[3][1] = m[1][3];
- m[3][2] = m[2][3];
-
- Tensor dm_drr01[4][4];
- if(!alEqDis){
- dm_drr01[0][0] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,-1.0);
- dm_drr01[1][1] = 2.0*Tensor(-1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,+1.0);
- dm_drr01[2][2] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,-1.0, 0.0, 0.0, 0.0,+1.0);
- dm_drr01[3][3] = 2.0*Tensor(+1.0, 0.0, 0.0, 0.0,+1.0, 0.0, 0.0, 0.0,-1.0);
- dm_drr01[0][1] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,+1.0, 0.0);
- dm_drr01[0][2] = 2.0*Tensor( 0.0, 0.0,+1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
- dm_drr01[0][3] = 2.0*Tensor( 0.0,-1.0, 0.0, +1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
- dm_drr01[1][2] = 2.0*Tensor( 0.0,-1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
- dm_drr01[1][3] = 2.0*Tensor( 0.0, 0.0,-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0);
- dm_drr01[2][3] = 2.0*Tensor( 0.0, 0.0, 0.0, 0.0, 0.0,-1.0, 0.0,-1.0, 0.0);
- dm_drr01[1][0] = dm_drr01[0][1];
- dm_drr01[2][0] = dm_drr01[0][2];
- dm_drr01[2][1] = dm_drr01[1][2];
- dm_drr01[3][0] = dm_drr01[0][3];
- dm_drr01[3][1] = dm_drr01[1][3];
- dm_drr01[3][2] = dm_drr01[2][3];
- }
-
- std::vector<double> eigenvals;
- Matrix<double> eigenvecs;
- int diagerror=diagMat(m, eigenvals, eigenvecs );
-
- if (diagerror!=0){
- std::string sdiagerror;
- Tools::convert(diagerror,sdiagerror);
- std::string msg="DIAGONALIZATION FAILED WITH ERROR CODE "+sdiagerror;
- plumed_merror(msg);
+ if( getAlign()==getDisplace() ) return myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ return myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
}
-
- dist=eigenvals[0]+rr00+rr11;
-
- Matrix<double> ddist_dm(4,4);
-
- Vector4d q(eigenvecs[0][0],eigenvecs[0][1],eigenvecs[0][2],eigenvecs[0][3]);
-
- Tensor dq_drr01[4];
- if(!alEqDis){
- double dq_dm[4][4][4];
- for(unsigned i=0;i<4;i++) for(unsigned j=0;j<4;j++) for(unsigned k=0;k<4;k++){
- double tmp=0.0;
-// perturbation theory for matrix m
- for(unsigned l=1;l<4;l++) tmp+=eigenvecs[l][j]*eigenvecs[l][i]/(eigenvals[0]-eigenvals[l])*eigenvecs[0][k];
- dq_dm[i][j][k]=tmp;
- }
-// propagation to _drr01
- for(unsigned i=0;i<4;i++){
- Tensor tmp;
- for(unsigned j=0;j<4;j++) for(unsigned k=0;k<4;k++) {
- tmp+=dq_dm[i][j][k]*dm_drr01[j][k];
- }
- dq_drr01[i]=tmp;
- }
- }
-
-// This is the rotation matrix that brings reference to positions
-// i.e. matmul(rotation,reference[iat])+shift is fitted to positions[iat]
-
- Tensor rotation;
- rotation[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3];
- rotation[1][1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3];
- rotation[2][2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3];
- rotation[0][1]=2*(+q[0]*q[3]+q[1]*q[2]);
- rotation[0][2]=2*(-q[0]*q[2]+q[1]*q[3]);
- rotation[1][2]=2*(+q[0]*q[1]+q[2]*q[3]);
- rotation[1][0]=2*(-q[0]*q[3]+q[1]*q[2]);
- rotation[2][0]=2*(+q[0]*q[2]+q[1]*q[3]);
- rotation[2][1]=2*(-q[0]*q[1]+q[2]*q[3]);
-
-
- Tensor drotation_drr01[3][3];
- if(!alEqDis){
- drotation_drr01[0][0]=2*q[0]*dq_drr01[0]+2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
- drotation_drr01[1][1]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]+2*q[2]*dq_drr01[2]-2*q[3]*dq_drr01[3];
- drotation_drr01[2][2]=2*q[0]*dq_drr01[0]-2*q[1]*dq_drr01[1]-2*q[2]*dq_drr01[2]+2*q[3]*dq_drr01[3];
- drotation_drr01[0][1]=2*(+(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
- drotation_drr01[0][2]=2*(-(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
- drotation_drr01[1][2]=2*(+(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
- drotation_drr01[1][0]=2*(-(q[0]*dq_drr01[3]+dq_drr01[0]*q[3])+(q[1]*dq_drr01[2]+dq_drr01[1]*q[2]));
- drotation_drr01[2][0]=2*(+(q[0]*dq_drr01[2]+dq_drr01[0]*q[2])+(q[1]*dq_drr01[3]+dq_drr01[1]*q[3]));
- drotation_drr01[2][1]=2*(-(q[0]*dq_drr01[1]+dq_drr01[0]*q[1])+(q[2]*dq_drr01[3]+dq_drr01[2]*q[3]));
- }
-
- double prefactor=2.0;
-
- if(!squared && alEqDis) prefactor*=0.5/sqrt(dist);
-
-// if "safe", recompute dist here to a better accuracy
- if(safe || !alEqDis) dist=0.0;
-
-// If safe is set to "false", MSD is taken from the eigenvalue of the M matrix
-// If safe is set to "true", MSD is recomputed from the rotational matrix
-// For some reason, this last approach leads to less numerical noise but adds an overhead
-
- Tensor ddist_drotation;
- Vector ddist_dcpositions;
-
-// third expensive loop: derivatives
- for(unsigned iat=0;iat<n;iat++){
- Vector d(positions[iat]-cpositions - matmul(rotation,getReferencePosition(iat)));
- if(alEqDis){
-// there is no need for derivatives of rotation and shift here as it is by construction zero
-// (similar to Hellman-Feynman forces)
- addAtomicDerivatives( iat, prefactor*align[iat]*d );
- if(safe) dist+=align[iat]*modulo2(d);
- } else {
-// the case for align != displace is different, sob:
- dist+=displace[iat]*modulo2(d);
-// these are the derivatives assuming the roto-translation as frozen
- atom_ders[iat]=2*displace[iat]*d;
-// here I accumulate derivatives wrt rotation matrix ..
- ddist_drotation+=-2*displace[iat]*extProduct(d,getReferencePosition(iat));
-// .. and cpositions
- ddist_dcpositions+=-2*displace[iat]*d;
- }
- }
-
- if(!alEqDis){
- Tensor ddist_drr01;
- for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) ddist_drr01+=ddist_drotation[i][j]*drotation_drr01[i][j];
- for(unsigned iat=0;iat<n;iat++){
-// this is propagating to positions.
-// I am implicitly using the derivative of rr01 wrt positions here
- addAtomicDerivatives( iat, matmul(ddist_drr01,getReferencePosition(iat))*align[iat] );
- addAtomicDerivatives( iat, ddist_dcpositions*align[iat] );
- }
- }
- if(!squared){
- dist=sqrt(dist);
- if(!alEqDis){
- double xx=0.5/dist;
- for(unsigned iat=0;iat<atom_ders.size();iat++) atom_ders[iat]*=xx;
- }
- }
-
- return dist;
}
}
diff --git a/src/reference/ReferenceAtoms.h b/src/reference/ReferenceAtoms.h
index 2ea7ca1f2..1772c62d1 100644
--- a/src/reference/ReferenceAtoms.h
+++ b/src/reference/ReferenceAtoms.h
@@ -100,8 +100,6 @@ public:
virtual void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in )=0;
/// Print the atomic positions
void printAtoms( OFile& ofile ) const ;
-/// Return all the reference positions
- const std::vector<Vector>& getReferencePositions();
/// Return all atom indexes
const std::vector<AtomNumber>& getAbsoluteIndexes();
};
@@ -160,11 +158,6 @@ void ReferenceAtoms::addBoxDerivatives( const Tensor& vir ){
virialWasSet=true; virial+=vir;
}
-inline
-const std::vector<Vector>& ReferenceAtoms::getReferencePositions(){
- return reference_atoms;
-}
-
inline
const std::vector<AtomNumber>& ReferenceAtoms::getAbsoluteIndexes(){
return indices;
diff --git a/src/reference/SimpleRMSD.cpp b/src/reference/SimpleRMSD.cpp
index 06ec688f3..a40850c13 100644
--- a/src/reference/SimpleRMSD.cpp
+++ b/src/reference/SimpleRMSD.cpp
@@ -21,20 +21,17 @@
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "RMSDBase.h"
#include "MetricRegister.h"
+#include "tools/RMSD.h"
namespace PLMD{
class SimpleRMSD : public RMSDBase {
private:
- std::vector<double> weights;
+ RMSD myrmsd;
public:
SimpleRMSD( const ReferenceConfigurationOptions& ro );
void read( const PDB& );
double calc( const std::vector<Vector>& pos, const bool& squared );
- double simpleAlignment(const std::vector<double> & align,
- const std::vector<double> & displace,
- const std::vector<Vector> & positions,
- bool squared);
};
PLUMED_REGISTER_METRIC(SimpleRMSD,"SIMPLE")
@@ -46,66 +43,11 @@ RMSDBase( ro )
}
void SimpleRMSD::read( const PDB& pdb ){
- readAtomsFromPDB( pdb ); weights.resize( pdb.getAtomNumbers().size() );
- for(unsigned i=0;i<weights.size();++i) weights[i]=getDisplace()[i];
+ readReference( pdb );
}
double SimpleRMSD::calc( const std::vector<Vector>& pos, const bool& squared ){
- return simpleAlignment( getAlign(), getDisplace(), pos, squared );
-}
-
-double SimpleRMSD::simpleAlignment(const std::vector<double> & align,
- const std::vector<double> & displace,
- const std::vector<Vector> & positions,
- bool squared) {
-
- double dist(0), anorm(0), dnorm(0);
- unsigned n=getNumberOfReferencePositions();
- // Number of members of align and displace is the number of reference positions
- plumed_dbg_assert( n==align.size() && n==displace.size() );
- // Positions array might contain vectors that are not particularly interesting
- plumed_dbg_assert( positions.size()==getNumberOfAtoms() );
-
- Vector iref, apositions, areference;
- Vector dpositions, dreference;
-
- for(unsigned i=0;i<n;i++){
- unsigned iatom=getAtomIndex(i);
- double aw=align[i];
- double dw=displace[i];
- anorm+=aw;
- dnorm+=dw;
- iref=getReferencePosition(i);
- apositions+=positions[iatom]*aw;
- areference+=iref*aw;
- dpositions+=positions[iatom]*dw;
- dreference+=iref*dw;
- }
-
- double invdnorm=1.0/dnorm;
- double invanorm=1.0/anorm;
-
- apositions*=invanorm;
- areference*=invanorm;
- dpositions*=invdnorm;
- dreference*=invdnorm;
-
- Vector shift=((apositions-areference)-(dpositions-dreference));
- for(unsigned i=0;i<n;i++){
- unsigned iatom=getAtomIndex(i);
- Vector d=(positions[iatom]-apositions)-(getReferencePosition(i)-areference);
- dist+=displace[i]*d.modulo2();
- addAtomicDerivatives( i,2*(invdnorm*displace[i]*d+invanorm*align[i]*shift) );
- }
- dist*=invdnorm;
-
- if(!squared){
- // sqrt and normalization
- dist=sqrt(dist);
- ///// sqrt and normalization on derivatives
- for(unsigned i=0;i<getNumberOfAtoms();i++){atom_ders[i]*=(0.5/dist);}
- }
- return dist;
+ return myrmsd.simpleAlignment( getAlign(), getDisplace(), pos, getReferencePositions(), atom_ders, squared );
}
}
diff --git a/src/tools/RMSD.cpp b/src/tools/RMSD.cpp
index 9c3e3a9d6..9a400b83f 100644
--- a/src/tools/RMSD.cpp
+++ b/src/tools/RMSD.cpp
@@ -32,9 +32,7 @@
using namespace std;
namespace PLMD{
-RMSD::RMSD(Log & log ):
- alignmentMethod(SIMPLE),
- log(&log){}
+RMSD::RMSD() : alignmentMethod(SIMPLE) {}
void RMSD::set(const PDB&pdb, string mytype ){
@@ -49,15 +47,12 @@ void RMSD::setType(string mytype){
alignmentMethod=SIMPLE; // initialize with the simplest case: no rotation
if (mytype=="SIMPLE"){
alignmentMethod=SIMPLE;
- log->printf("RMSD IS DONE WITH SIMPLE METHOD(NO ROTATION)\n");
}
else if (mytype=="OPTIMAL"){
alignmentMethod=OPTIMAL;
- log->printf("RMSD IS DONE WITH OPTIMAL ALIGNMENT METHOD\n");
}
else if (mytype=="OPTIMAL-FAST"){
alignmentMethod=OPTIMAL_FAST;
- log->printf("RMSD IS DONE WITH OPTIMAL-FAST ALIGNMENT METHOD (fast version, numerically less stable, only valid with align==displace)\n");
}
else plumed_merror("unknown RMSD type" + mytype);
diff --git a/src/tools/RMSD.h b/src/tools/RMSD.h
index 5d1636c3c..2d79a4c30 100644
--- a/src/tools/RMSD.h
+++ b/src/tools/RMSD.h
@@ -68,11 +68,9 @@ class RMSD
std::vector<double> align;
// Weights for deviation
std::vector<double> displace;
-// Logfile
- Log *log;
public:
/// Constructor
- RMSD(Log & log );
+ RMSD();
/// clear the structure
void clear();
/// set reference, align and displace from input pdb structure
--
GitLab