From 4cf58b08d6b6a1136dee029136081f97dc4a4166 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Mon, 27 Feb 2017 16:59:20 +0100
Subject: [PATCH] Removed unused part of the code
These files were related to the original implementation of RMSD.
I think they have not been used for years, so I remove them.
[makedoc]
---
src/tools/Kearsley.cpp | 874 ---------------------------------
src/tools/Kearsley.h | 239 ---------
src/tools/OptimalAlignment.cpp | 399 ---------------
src/tools/OptimalAlignment.h | 82 ----
4 files changed, 1594 deletions(-)
delete mode 100644 src/tools/Kearsley.cpp
delete mode 100644 src/tools/Kearsley.h
delete mode 100644 src/tools/OptimalAlignment.cpp
delete mode 100644 src/tools/OptimalAlignment.h
diff --git a/src/tools/Kearsley.cpp b/src/tools/Kearsley.cpp
deleted file mode 100644
index 0061b5e9a..000000000
--- a/src/tools/Kearsley.cpp
+++ /dev/null
@@ -1,874 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- Copyright (c) 2011-2016 The plumed team
- (see the PEOPLE file at the root of the distribution for a list of names)
-
- See http://www.plumed.org for more information.
-
- This file is part of plumed, version 2.
-
- plumed is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- plumed is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with plumed. If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "Kearsley.h"
-#include <cmath>
-#include <iostream>
-#include <cstdlib>
-#include "Matrix.h"
-#include "Tensor.h"
-#include "Log.h"
-#include "Matrix.h"
-#include "Random.h"
-
-using namespace std;
-namespace PLMD{
-
-// put some notes
-
-Kearsley::Kearsley(const vector<Vector> &p0, const vector<Vector> &p1, const vector<double> &align, Log* &log):
- log(log),
- p0(p0),
- p1(p1),
- align(align),
- com0_is_removed(false),
- com1_is_removed(false),
- err(0.0)
-{
- // now make an initial allocation
-// int n=p0.size();
- // eventually here one should make a "hard" resize of all the structures
-// log->printf("Reallocating a size of %d atoms to kearsley structure\n",n);
-
-// try{
-// diff.resize(n);
-// }catch(bad_alloc&) {
-// cerr<<"Cannot allocate the vector in Kearsley"<<endl;
-// }
-// try{
-// diff.resize(n);
-// }catch(bad_alloc&) {
-// cerr<<"Cannot allocate the vector in Kearsley"<<endl;
-// }
-
-}
-// do the alignment
-
-double Kearsley::calculate(bool rmsd) {
- // just an ad-hoc scaling factor
- double myscale=1.;
- // basic sanity check
- if(p0.size()!=p1.size() || p1.size()!=align.size()){
- cerr<<"Kearsley: looks like you have not properly allocated the vectors: the two frames have different size"<<endl;
- cerr<<"size of p0 is :"<<p0.size()<<endl;
- cerr<<"size of p1 is :"<<p1.size()<<endl;
- cerr<<"size of align is :"<<align.size()<<endl;
- exit(0);
- }
- if(p0.empty() || p1.empty() ){
- cerr<<"Kearsley: looks like you have not properly allocated the vectors: they do not contain anything"<<endl;
- exit(0);
- }
- Vector rr1,rr0;
- Vector4d q;
- double dddq[4][4][4],gamma[3][3][3],rrsq;
- Matrix<double> m=Matrix<double>(4,4);
-// double dm_r1[4][4][3],dm_r0[4][4][3];
- Vector dm_r1[4][4];
- Vector dm_r0[4][4];
- Tensor pi1,pi0;
- Tensor d;
- Tensor dinv;
- Vector xx;
- double totalign=0., s, tmp1, err;
- unsigned i,j,k,l,ii,ll,jj,mm,n,nn,iii;
- bool verbose=false;
-
- vector<int> alignmap; // on the fly map done for optimization
-
- unsigned natoms=p0.size();
-
-
- // calculate coms
-
- for(i=0;i<natoms;i++){
- if (align[i]>0.){
- alignmap.push_back(i);
- totalign+=align[i];
- }
- if (align[i]<0.){cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl;exit(0);};
-
- }
-
-
- // later will be implemented something for optimizing this piece of crap
-
- com0_is_removed=false;
- com1_is_removed=false;
-
- bool do_center=true; // keep it for legacy code compatibility
-
- if(com0_is_removed==false){
-
- xx.zero();
-
- if(do_center) {// if you dont need to center no prob...
- for(i=0;i<alignmap.size();i++){
- xx+=p0[alignmap[i]]*align[alignmap[i]];
- }
- xx/=(totalign);
- }
-
- com0=xx;
-
- if (p0reset.empty()){p0reset.resize(natoms);}
- for(i=0;i<natoms;i++){
- p0reset[i]=p0[i]-xx;
- }
- com0_is_removed=true;
-
- if (verbose){
- log->printf("P0 RESET\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0reset[i][0]/myscale,p0reset[i][1]/myscale,p0reset[i][2]/myscale);
- }
- log->printf("END\n");
- }
-
- }
-
- if(com1_is_removed==false){
-
- xx.zero();
-
- if(do_center) {// if you dont need to center no prob...
- for(i=0;i<alignmap.size();i++){
- xx+=p1[alignmap[i]]*align[alignmap[i]];
- }
- xx/=(totalign);
- }
-
- com1=xx;
-
- if (p1reset.empty()){p1reset.resize(natoms);}
-
- for(i=0;i<natoms;i++){
- p1reset[i]=p1[i]-xx;
- }
- com1_is_removed=true;
-
- if(verbose){
- log->printf("P1 RESET\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1reset[i][0]/myscale,p1reset[i][1]/myscale,p1reset[i][2]/myscale);
- }
- log->printf("END\n");
- }
-
- }
-
- bool fake=false;
- if(fake){
- // case of trivial alignment
- // set rotmat
- rotmat0on1=Tensor::identity();
- if (p0reset.size()==0){p0reset.resize(natoms);}
- if (p1reset.size()==0){p1reset.resize(natoms);}
-
- derrdp0.resize(natoms);
- derrdp1.resize(natoms);
- dmatdp0.resize(3*3*3*natoms);for(i=0;i<dmatdp0.size();i++)dmatdp0[i]=0.;
- dmatdp1.resize(3*3*3*natoms);for(i=0;i<dmatdp1.size();i++)dmatdp1[i]=0.;
-
- err=0.;
- for(i=0;i<natoms;i++){
- if(align[i]>0.)err+=align[i]*modulo2(p0reset[i]-p1reset[i]);
- }
-
- return 0.;
- }
- //
- // CLEAN M MATRIX
-
- m=0.;
-
- // ASSIGN MATRIX ELEMENTS USING ONLY THE ATOMS INVOLVED IN ALIGNMENT
-
- for(i=0;i<alignmap.size();i++){
-
- k=alignmap[i];
- tmp1=sqrt(align[k]/totalign);
-
- // adopt scaled coordinates
-
- rr1=p1reset[k]*tmp1;
- rr0=p0reset[k]*tmp1;
-
- rrsq=modulo2(rr0)+modulo2(rr1);
-
- m[0][0] += rrsq+2.*(-rr0[0]*rr1[0]-rr0[1]*rr1[1]-rr0[2]*rr1[2]);
- m[1][1] += rrsq+2.*(-rr0[0]*rr1[0]+rr0[1]*rr1[1]+rr0[2]*rr1[2]);
- m[2][2] += rrsq+2.*(+rr0[0]*rr1[0]-rr0[1]*rr1[1]+rr0[2]*rr1[2]);
- m[3][3] += rrsq+2.*(+rr0[0]*rr1[0]+rr0[1]*rr1[1]-rr0[2]*rr1[2]);
- m[0][1] += 2.*(-rr0[1]*rr1[2]+rr0[2]*rr1[1]);
- m[0][2] += 2.*( rr0[0]*rr1[2]-rr0[2]*rr1[0]);
- m[0][3] += 2.*(-rr0[0]*rr1[1]+rr0[1]*rr1[0]);
- m[1][2] -= 2.*( rr0[0]*rr1[1]+rr0[1]*rr1[0]);
- m[1][3] -= 2.*( rr0[0]*rr1[2]+rr0[2]*rr1[0]);
- m[2][3] -= 2.*( rr0[1]*rr1[2]+rr0[2]*rr1[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];
-
- // diagonalize the 4x4 matrix
-
- vector<double> eigenvals;
- Matrix<double> eigenvecs;
-
- int diagerror=diagMat(m, eigenvals, eigenvecs );
-
- if (diagerror!=0){cerr<<"DIAGONALIZATION FAILED WITH ERROR CODE "<<diagerror<<endl;exit(0);}
-
- s=1.0;
- if(eigenvecs(0,0)<0.)s=-1.;//correct for negative values (?)
- // eigenvecs are in rows!!
-
- q[0]=s*eigenvecs(0,0);
- q[1]=s*eigenvecs(0,1);
- q[2]=s*eigenvecs(0,2);
- q[3]=s*eigenvecs(0,3);
- err=eigenvals[0];
-
- //log->printf(" ERR: %20.10f \n",err);
-
- if(verbose){
- log->printf(" ERR: %f \n",err);
- for (i=0;i<4;i++){
- log->printf(" EIGENVALS: %f \n",eigenvals[i]);
- }
- }
-
- if(abs(eigenvals[0]-eigenvals[1])<1.e-8){
- cerr<<"DIAGONALIZATION: NON UNIQUE SOLUTION"<<endl;exit(0);
- }
-
- /*
- * the ROTATION matrix
- */
-
- d[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3] ;
- d[1][0]=2.0*(q[1]*q[2]-q[0]*q[3]);
- d[2][0]=2.0*(q[1]*q[3]+q[0]*q[2]);
- d[0][1]=2.0*(q[1]*q[2]+q[0]*q[3]);
- d[1][1]=q[0]*q[0]+q[2]*q[2]-q[1]*q[1]-q[3]*q[3];
- d[2][1]=2.0*(q[2]*q[3]-q[0]*q[1]);
- d[0][2]=2.0*(q[1]*q[3]-q[0]*q[2]);
- d[1][2]=2.0*(q[2]*q[3]+q[0]*q[1]);
- d[2][2]=q[0]*q[0]+q[3]*q[3]-q[1]*q[1]-q[2]*q[2];
-
- /*
- * first derivative in perturbation theory : derivative of the rotation matrix respect to the
- * quternion vectors
- */
-
- dddq[0][0][0]= 2.0*q[0];
- dddq[1][0][0]=-2.0*q[3];
- dddq[2][0][0]= 2.0*q[2];
- dddq[0][1][0]= 2.0*q[3];
- dddq[1][1][0]= 2.0*q[0];
- dddq[2][1][0]=-2.0*q[1];
- dddq[0][2][0]=-2.0*q[2];
- dddq[1][2][0]= 2.0*q[1];
- dddq[2][2][0]= 2.0*q[0];
-
- dddq[0][0][1]= 2.0*q[1];
- dddq[1][0][1]= 2.0*q[2];
- dddq[2][0][1]= 2.0*q[3];
- dddq[0][1][1]= 2.0*q[2];
- dddq[1][1][1]=-2.0*q[1];
- dddq[2][1][1]=-2.0*q[0];
- dddq[0][2][1]= 2.0*q[3];
- dddq[1][2][1]= 2.0*q[0];
- dddq[2][2][1]=-2.0*q[1];
-
- dddq[0][0][2]=-2.0*q[2];
- dddq[1][0][2]= 2.0*q[1];
- dddq[2][0][2]= 2.0*q[0];
- dddq[0][1][2]= 2.0*q[1];
- dddq[1][1][2]= 2.0*q[2];
- dddq[2][1][2]= 2.0*q[3];
- dddq[0][2][2]=-2.0*q[0];
- dddq[1][2][2]= 2.0*q[3];
- dddq[2][2][2]=-2.0*q[2];
-
- dddq[0][0][3]=-2.0*q[3];
- dddq[1][0][3]=-2.0*q[0];
- dddq[2][0][3]= 2.0*q[1];
- dddq[0][1][3]= 2.0*q[0];
- dddq[1][1][3]=-2.0*q[3];
- dddq[2][1][3]= 2.0*q[2];
- dddq[0][2][3]= 2.0*q[1];
- dddq[1][2][3]= 2.0*q[2];
- dddq[2][2][3]= 2.0*q[3];
-
- /*
- * Build gamma 3x3x3 matrix
- */
- for(i=0;i<3;i++){ //direction
- for(j=0;j<3;j++){ //direction
- for(k=0;k<3;k++){ //eigenvector number
- gamma[i][j][k]=0.0;
- for(l=0;l<4;l++){ //components of each eigenvector in pert. series
- if(abs(eigenvals[0]-eigenvals[k+1])<1.e-8){
- log->printf(" FOUND DEGENERACY IN RMSD_ESS ROUTINE \n");
- log->printf(" I'm DYING....\n");
- log->printf(" COPYING STACK HERE \n");
- log->printf(" P0\n");
- for(ll=0;ll<natoms;ll++)log->printf(" %f %f %f \n",p0reset[ll][0],p0reset[ll][1],p0reset[ll][2]);
- log->printf(" P1\n");
- for(ll=0;ll<natoms;ll++)log->printf(" %f %f %f \n",p1reset[ll][0],p1reset[ll][1],p1reset[ll][2]);
- exit(0);
- }
- else{
- gamma[i][j][k] += dddq[i][j][l]*eigenvecs(k+1,l)/(eigenvals[0]-eigenvals[k+1]);
- }
- }
- //log->printf("GAMMA %2d %2d %2d V %12.6f\n",i,j,k,gamma[i][j][k]);
- }
- }
- }
-
- // allocate various arrays
-
- dmatdp1.resize(3*3*3*natoms);
- for(i=0;i<dmatdp1.size();i++)dmatdp1[i]=0.;
- dmatdp0.resize(3*3*3*natoms);
- for(i=0;i<dmatdp0.size();i++)dmatdp0[i]=0.;
-
- vector<double> dd_dr_temp;dd_dr_temp.resize(natoms);
-
- vector<Vector> derr_dr1;
- derr_dr1.resize(natoms);
- vector<Vector> derr_dr0;
- derr_dr0.resize(natoms);
- vector<Vector> array_3_n;
- array_3_n.resize(natoms);
-
-
- /*
- * Table of Derivative of the quaternion matrix respect to atom position: needed only if simple
- * alignment is required and no correction respect to the rotation matrix is wanted
- */
-
- for(iii=0;iii<alignmap.size();iii++){
-
- i=alignmap[iii];
- tmp1=sqrt(align[i]/totalign);
-
- // once again: derivative respect to scaled distance
-
- rr1=2.*p1reset[i]*tmp1;
- rr0=2.*p0reset[i]*tmp1;
-
-
- dm_r1 [0][0][0]=(rr1[0]-rr0[0]);
- dm_r1 [0][0][1]=(rr1[1]-rr0[1]);
- dm_r1 [0][0][2]=(rr1[2]-rr0[2]);
-
- dm_r1 [0][1][0]=0.;
- dm_r1 [0][1][1]= rr0[2];
- dm_r1 [0][1][2]=-rr0[1];
-
- dm_r1 [0][2][0]=-rr0[2];
- dm_r1 [0][2][1]= 0.;
- dm_r1 [0][2][2]= rr0[0];
-
- dm_r1 [0][3][0]= rr0[1];
- dm_r1 [0][3][1]=-rr0[0];
- dm_r1 [0][3][2]= 0.;
-
- dm_r1 [1][1][0]=(rr1[0]-rr0[0]);
- dm_r1 [1][1][1]=(rr1[1]+rr0[1]);
- dm_r1 [1][1][2]=(rr1[2]+rr0[2]);
-
- dm_r1 [1][2][0]=-rr0[1];
- dm_r1 [1][2][1]=-rr0[0];
- dm_r1 [1][2][2]= 0.;
-
- dm_r1 [1][3][0]=-rr0[2];
- dm_r1 [1][3][1]= 0.;
- dm_r1 [1][3][2]=-rr0[0];
-
- dm_r1 [2][2][0]=(rr1[0]+rr0[0]);
- dm_r1 [2][2][1]=(rr1[1]-rr0[1]);
- dm_r1 [2][2][2]=(rr1[2]+rr0[2]);
-
- dm_r1 [2][3][0]=0.;
- dm_r1 [2][3][1]=-rr0[2];
- dm_r1 [2][3][2]=-rr0[1];
-
- dm_r1 [3][3][0]=(rr1[0]+rr0[0]);
- dm_r1 [3][3][1]=(rr1[1]+rr0[1]);
- dm_r1 [3][3][2]=(rr1[2]-rr0[2]);
- /*
- derivative respec to to the other vector
- */
- dm_r0 [0][0][0]=-(rr1[0]-rr0[0]);
- dm_r0 [0][0][1]=-(rr1[1]-rr0[1]);
- dm_r0 [0][0][2]=-(rr1[2]-rr0[2]);
-
- dm_r0 [0][1][0]=0. ;
- dm_r0 [0][1][1]=-rr1[2];
- dm_r0 [0][1][2]=rr1[1];
-
- dm_r0 [0][2][0]= rr1[2];
- dm_r0 [0][2][1]= 0.;
- dm_r0 [0][2][2]=-rr1[0];
-
- dm_r0 [0][3][0]=-rr1[1] ;
- dm_r0 [0][3][1]= rr1[0];
- dm_r0 [0][3][2]= 0.;
-
- dm_r0 [1][1][0]=-(rr1[0]-rr0[0]);
- dm_r0 [1][1][1]=(rr1[1]+rr0[1]);
- dm_r0 [1][1][2]=(rr1[2]+rr0[2]);
-
- dm_r0 [1][2][0]=-rr1[1];
- dm_r0 [1][2][1]=-rr1[0];
- dm_r0 [1][2][2]= 0.;
-
- dm_r0 [1][3][0]=-rr1[2];
- dm_r0 [1][3][1]= 0.;
- dm_r0 [1][3][2]=-rr1[0];
-
- dm_r0 [2][2][0]=(rr1[0]+rr0[0]);
- dm_r0 [2][2][1]=-(rr1[1]-rr0[1]);
- dm_r0 [2][2][2]=(rr1[2]+rr0[2]);
-
- dm_r0 [2][3][0]=0.;
- dm_r0 [2][3][1]=-rr1[2];
- dm_r0 [2][3][2]=-rr1[1];
-
- dm_r0 [3][3][0]=(rr1[0]+rr0[0]);
- dm_r0 [3][3][1]=(rr1[1]+rr0[1]);
- dm_r0 [3][3][2]=-(rr1[2]-rr0[2]);
- /*
- * write the diagonal
- */
-
- dm_r1[1][0]=dm_r1[0][1];
- dm_r1[2][0]=dm_r1[0][2];
- dm_r1[3][0]=dm_r1[0][3];
- dm_r1[2][1]=dm_r1[1][2];
- dm_r1[3][1]=dm_r1[1][3];
- dm_r1[3][2]=dm_r1[2][3];
-
- dm_r0[1][0]=dm_r0[0][1];
- dm_r0[2][0]=dm_r0[0][2];
- dm_r0[3][0]=dm_r0[0][3];
- dm_r0[2][1]=dm_r0[1][2];
- dm_r0[3][1]=dm_r0[1][3];
- dm_r0[3][2]=dm_r0[2][3];
-
-
- //log->printf("DMDR0 ALIGN %f AT %d VAL %f\n",align[i],i,dm_r0[0][0][j]);
-
-
- /*
- * pi matrix : coefficents in per theory
- */
-
- pi0.zero();
- pi1.zero();
- derr_dr1[i].zero();
- derr_dr0[i].zero();
-
- for(k=0;k<4;k++){
- for(l=0;l<4;l++){
- derr_dr1[i]+=(q[k]*q[l])*dm_r1[l][k];
- derr_dr0[i]+=(q[k]*q[l])*dm_r0[l][k];
- for(mm=0;mm<3;mm++)for(j=0;j<3;j++){
- pi0[mm][j]+=eigenvecs(mm+1,k)*dm_r0[l][k][j]*q[l];
- pi1[mm][j]+=eigenvecs(mm+1,k)*dm_r1[l][k][j]*q[l];
- };
- };
- };
-/*
- derr_dr1[i]/=totalign;
- derr_dr0[i]/=totalign;
-
-
-*/
-
- for(j=0;j<3;j++){
- for (k=0;k<3;k++){
- for(l=0;l<3;l++){
- int ind=j*3*3*natoms+k*3*natoms+l*natoms+i;
- dmatdp0[ind]=0.;
- dmatdp1[ind]=0.;
- for(ii=0;ii<3;ii++){
- dmatdp1[ind]+=gamma[j][k][ii]*pi1[ii][l];
- dmatdp0[ind]+=gamma[j][k][ii]*pi0[ii][l];
- }
-
- }
-
- }
-
- }
- }
-
-
-
-
- // end of the calculation of the derivative of the rotation matrix
-
- /*
- * Now correct for center of mass: only if needed
- *
- */
-
- bool comcorr_r1=true;
-
- if(comcorr_r1){
- for(k=0;k<alignmap.size();k++){
- i=alignmap[k];
- tmp1=sqrt(align[i]/totalign);
- array_3_n[i]=tmp1*derr_dr1[i];
- if(do_center){
- for(jj=0;jj<alignmap.size();jj++){
- j=alignmap[jj];
- array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr1[j];
- }
- }
- }
- for(k=0;k<alignmap.size();k++){
- i=alignmap[k];
- derr_dr1[i]=array_3_n[i];
- }
- }
-
-
- bool do_comcorr_r0=true;
- //
- // correction for r0 frame
- //
- if(do_comcorr_r0){
- for(k=0;k<alignmap.size();k++){
- i=alignmap[k];
- tmp1=sqrt(align[i]/totalign);
- array_3_n[i]=tmp1*derr_dr0[i];
- if(do_center){
- for(jj=0;jj<alignmap.size();jj++){
- j=alignmap[jj];
- array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr0[j];
- }
- }
- }
- for(k=0;k<alignmap.size();k++){
- i=alignmap[k];
- derr_dr0[i]=array_3_n[i];
- }
- }
-
-
- bool do_der_r1=true;
- bool do_der_r0=true;
- bool do_der_rotmat=true;
-
- if(do_der_r1 && do_der_rotmat){
- for(i=0;i<3;i++){
- for(j=0;j<3;j++){
- for(k=0;k<3;k++){
- for(ll=0;ll<alignmap.size();ll++){
- l=alignmap[ll];
- int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
- tmp1=sqrt(align[l]/totalign);
- dd_dr_temp[l]=tmp1*dmatdp1[ind];
- if(do_center){
- for(nn=0;nn<alignmap.size();nn++){
- n=alignmap[nn];
- dd_dr_temp[l]-=dmatdp1[ind-l+n]*tmp1*align[n]/totalign;
- }
- }
-
- }
- for(ll=0;ll<alignmap.size();ll++){
- l=alignmap[ll];
- int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
- dmatdp1[ind]=dd_dr_temp[l];
- }
-
- }
- }
- }
-
- }
-
- if(do_der_r0 && do_der_rotmat){
- for(i=0;i<3;i++){
- for(j=0;j<3;j++){
- for(k=0;k<3;k++){
- for(ll=0;ll<alignmap.size();ll++){
- l=alignmap[ll];
- int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
- tmp1=sqrt(align[l]/totalign);
- dd_dr_temp[l]=tmp1*dmatdp0[ind];
- if(do_center){
- for(nn=0;nn<alignmap.size();nn++){
- n=alignmap[nn];
- dd_dr_temp[l]-=dmatdp0[ind-l+n]*tmp1*align[n]/totalign;
- }
- }
- }
- for(ll=0;ll<alignmap.size();ll++){
- l=alignmap[ll];
- int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
- dmatdp0[ind]=dd_dr_temp[l];
- }
- }
- }
- }
- }
-
-
- bool do_p1rotated=true;
- if (do_p1rotated){
- // resize if not allocated
-
- if(p1.size()!=p1rotated.size())p1rotated.resize(p1.size());
-
-// exit(0);
-
- for(i=0;i<natoms;i++) p1rotated[i]=matmul(d,p1reset[i]);
-
- // reallocate difference vectors
- if(p1.size()!=diff1on0.size())diff1on0.resize(p1.size());
- for(i=0;i<natoms;i++){
- diff1on0[i]=p1rotated[i]-p0reset[i];
- }
-
- if(verbose){
- log->printf("P1-RESET-AND-ROTATED\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1rotated[i][0]/myscale,p1rotated[i][1]/myscale,p1rotated[i][2]/myscale);
- }
- log->printf("END\n");
- log->printf("P0-RESET\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0reset[i][0]/myscale,p0reset[i][1]/myscale,p0reset[i][2]/myscale);
- }
- log->printf("END\n");
- }
-
- }
-
-
-
- dinv=inverse(d);
-
- bool do_p0rotated=true;
- if (do_p0rotated){
- if(p0.size()!=p0rotated.size())p0rotated.resize(p0.size());
- for(i=0;i<natoms;i++) p0rotated[i]=matmul(dinv,p0reset[i]);
- if(p1.size()!=diff0on1.size())diff0on1.resize(p1.size());
- for(i=0;i<natoms;i++) diff0on1[i]=p0rotated[i]-p1reset[i];
- if(verbose){
- log->printf("P0-RESET AND INVERSE ROTATED\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0rotated[i][0]/myscale,p0rotated[i][1]/myscale,p0rotated[i][2]/myscale);
- }
- log->printf("END\n");
- log->printf("P1-RESET\n");
- for(i=0;i<natoms;i++){
- log->printf("ATOM %6u C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1reset[i][0]/myscale,p1reset[i][1]/myscale,p1reset[i][2]/myscale);
- }
- log->printf("END\n");
- }
- }
- // copy on the official vectors:
- rotmat0on1=d;
- rotmat1on0=dinv;
- derrdp0.resize(natoms);
- derrdp0=derr_dr0;
- derrdp1.resize(natoms);
- derrdp1=derr_dr1;
-
- // now rescale accordingly for rmsd instead of msd
- if(rmsd){
- err=sqrt(err);
- double tmp=0.5/err;
- for(ii=0;ii<alignmap.size();ii++){
- i=alignmap[ii];
- derrdp0[i]*=tmp;
- derrdp1[i]*=tmp;
- }
-
- }
-
- return err;
-
-}
-void Kearsley::assignP1(const std::vector<Vector> & p1) {
- this->p1=p1;
- com1_is_removed=false;
-}
-void Kearsley::assignP0(const std::vector<Vector> & p0) {
- this->p0=p0;
- com0_is_removed=false;
-}
-
-void Kearsley::assignAlign(const std::vector<double> & align) {
- this->align=align;
-}
-
-void Kearsley::finiteDifferenceInterface(bool rmsd){
-log->printf("Entering rmsd finite difference test system for kearsley\n");
-log->printf("-------------------------------------------\n");
-log->printf("TEST1: derivative of the value (derr_dr0/derr_dr1)\n");
-//// test 1
-unsigned i,j,l,m;
-double step=1.e-6,olderr,delta;
-// messing up a bit with align weights
-double delta1;
-vector<double> align1;
-align1.resize(p0.size());
-Random rnd;
-for (i=0;i<p0.size();i++){
- // draw a random number
- delta=rnd.RandU01();
- delta1=rnd.RandU01();
- if(delta>delta1){
- //if(delta>0.3){
- align1[i]=delta;
- }else{align1[i]=0.;};
- // log->printf("ALIGN %d IS %8.3f\n",i,align1[i]);
-}
-assignAlign(align1);
-//// get initial value of the error and derivative of it
-olderr=calculate(rmsd);
-log->printf("INITIAL ERROR VALUE: %e\n",olderr);
-// store the matrix
-Tensor old_rotmat0on1=rotmat0on1;
-
-//// get initial value of the error and derivative of it
-
-log->printf("TESTING: derrdp1 \n");
-for(unsigned j=0;j<3;j++){
- for(unsigned i=0;i<derrdp1.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p1[i][j]+=delta;
- com1_is_removed=false; // this is required whenever the assignment is not done with the methods
- com0_is_removed=false; // this is required whenever the assignment is not done with the methods
- err=calculate(rmsd);
- //log->printf("INITIAL ERROR VALUE: %e NEW ERROR %e DELTA %e ELEM %d %d \n",olderr,err,delta,i,j );
- p1[i][j]-=delta;
- switch(j){
- case 0:
- log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
- case 1:
- log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
- case 2:
- log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
-
- }
- }
-}
-//exit(0);
-log->printf("TESTING: derrdp0 \n");
-for(unsigned j=0;j<3;j++){
- for(unsigned i=0;i<derrdp0.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p0[i][j]+=delta;
- com0_is_removed=false; // this is required whenever the assignment is not done with the methods
- com1_is_removed=false; // this is required whenever the assignment is not done with the methods
-
- err=calculate(rmsd);
- p0[i][j]-=delta;
- switch(j){
- case 0:
- log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
- case 1:
- log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
- case 2:
- log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
- }
- }
-}
-
-log->printf("TESTING: dmatdp0 \n");
-for(l=0;l<3;l++){
- for(m=0;m<3;m++){
- for(j=0;j<3;j++){
- for(i=0;i<p0.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p0[i][j]+=delta;
- com0_is_removed=false;
- com1_is_removed=false;
- calculate(rmsd);
- p0[i][j]-=delta;
- int ind=l*3*3*p0.size()+m*3*p0.size()+j*p0.size()+i;
- switch(j){
- case 0:
- log->printf("TESTING: DMATDP0 [ %u ][ %u ]: X %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- case 1:
- log->printf("TESTING: DMATDP0 [ %u ][ %u ]: Y %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- case 2:
- log->printf("TESTING: DMATDP0 [ %u ][ %u ]: Z %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- }
- }
- }
- }
-}
-log->printf("TESTING: dmatdp1 \n");
-for(l=0;l<3;l++){
- for(m=0;m<3;m++){
- for(j=0;j<3;j++){
- for(i=0;i<p1.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p1[i][j]+=delta;
- com0_is_removed=false;
- com1_is_removed=false;
- calculate(rmsd);
- p1[i][j]-=delta;
- int ind=l*3*3*p1.size()+m*3*p1.size()+j*p1.size()+i;
- switch(j){
-
- case 0:
- log->printf("TESTING: DMATDP1 [ %u ][ %u ]: X %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- case 1:
- log->printf("TESTING: DMATDP1 [ %u ][ %u ]: Y %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- case 2:
- log->printf("TESTING: DMATDP1 [ %u ][ %u ]: Z %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
-
- }
- }
- }
- }
-}
-
- exit(0);
-}
-}
diff --git a/src/tools/Kearsley.h b/src/tools/Kearsley.h
deleted file mode 100644
index 1a3de1f1d..000000000
--- a/src/tools/Kearsley.h
+++ /dev/null
@@ -1,239 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- Copyright (c) 2011-2016 The plumed team
- (see the PEOPLE file at the root of the distribution for a list of names)
-
- See http://www.plumed.org for more information.
-
- This file is part of plumed, version 2.
-
- plumed is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- plumed is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with plumed. If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifndef __PLUMED_tools_Kearsley_h
-#define __PLUMED_tools_Kearsley_h
-
-#include "Vector.h"
-#include "Tensor.h"
-#include <vector>
-
-namespace PLMD{
-
-class Log;
-
-/*! A class that implements Kearsley's calculation
- which is optimal alignment via quaternion and
- analytical derivatives via perturbation theory
-
-In Kearsley algorithm (see, S. K. Kearsley, Acta Crystallogr., Sect. A: Found. Crystallogr. 45, 208 1989 ),
-here adapted to included the COM reset,
-one first calculates the COM reset position respect to frame \f$ \{x_0,y_0,z_0\} \f$ (running frame),
-being \f$ w_{tot}^{al}=\sum_i w_{al}^i \f$
-
-\f{eqnarray}
-\tilde x_0^i=x_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_0^i\\
-\tilde y_0^i=y_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_0^i\\
-\tilde z_0^i=z_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_0^i
-\f}
-
-and the same is done with the reference one \f$ \{x_1,y_1,z_1\} \f$
-\f{eqnarray}
-\tilde x_1^i=x_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_1^i\\
-\tilde y_1^i=y_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_1^i\\
-\tilde z_1^i=z_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_1^i
-\f}
-
-Then one can build the \f$ \{x_p,y_p,z_p\} \f$ and \f$ \{x_m,y_m,z_m\} \f$
-vectors of weighted summation and differences:
-\f{eqnarray}
-x_m^i=w_{al}^i(\tilde x_0^i-\tilde x_1^i)\\
-y_m^i=w_{al}^i(\tilde y_0^i-\tilde y_1^i)\\
-z_m^i=w_{al}^i(\tilde z_0^i-\tilde z_1^i)
-\f}
-
-\f{eqnarray}
-x_p^i=w_{al}^i(x_0^i+x_1^i)\\
-y_p^i=w_{al}^i(y_0^i+y_1^i)\\
-z_p^i=w_{al}^i(z_0^i+z_1^i)
-\f}
-
-
-Then one build the COM-resetted matrix
-
-\f{displaymath}
-\mathbf{M}=\left[
-\begin{array}{cccc}
-\sum ( {x}_{m}^{2}+{y}_{m}^{2}+{z}_{m}^{2}) &
- \sum (y_{p}z_{m} -y_{m}z_{p}) &
- \sum ( x_{m}z_{p} -x_{p}z_{m}) &
- \sum (x_{p}y_{m}-x_{m}y_{p} ) \\
- \sum ( y_{p}z_{m} -y_{m}z_{p}) &
-\sum ( {x}_{m}^{2}+{y}_{p}^{2}+{z}_{p}^{2}) &
-\sum ( x_{m}y_{m} -x_{p}y_{p} ) &
- \sum (x_{m}z_{m}-x_{p}z_{p} ) \\
- \sum (x_m z_p - x_p z_m ) &
-\sum ( x_m y_m -x_p y_p) &
-\sum ( {x}_{p}^{2}+{y}_{m}^{2}+{z}_{p}^{2}) &
-\sum ( y_m z_m -y_p z_p) \\
- \sum (x_p y_m -x_m y_p ) &
-\sum (x_m z_m - x_p z_p ) &
-\sum (y_m z_m- y_p z_p ) &
-\sum ( {x}_{p}^{2}+{y}_{p}^{2}+{z}_{m}^{2}) \\
-\end{array}
-\right]
-\f}
-
-by diagonalizing one obtains the mean square deviation by using the lowest eigenvalue \f$ \lambda_0 \f$
-\f{equation}
-MSD= \frac{\lambda_{0}}{w_{tot}^{al}}
-\f}
-
-The rotation matrix is obtained from the eigenvector corresponding to \f$ \lambda_0 \f$ eigenvalue
-having components \f$ q_1, q_2, q_3, q_4 \f$
-
-\f{displaymath}
-\mathbf{R}=\left[
-\begin{array}{ccc}
-q_1 ^2 + q_2 ^2 - q_3 ^2 - q_4^2 &
-2(q_2 q_3 + q_1 q_4) &
-2(q_1 q_4 -q_1 q_3 )\\
-2(q_2 q_3 - q_1 q_4) &
-q_1 ^2 +q_3 ^2 -q_2 ^2 -q_4^2 &
-2(q_3 q_4 - q_1 q_2)\\
-2( q_2 q_4 + q_1 q_3) &
-2( q_3 q_4 - q_1 q_2) &
-q_1^2 +q_4 ^2 - q_2^2 - q_3 ^2 \\
-\end{array}
-\right]
-\f}
-
-by using the perturbation theory one can retrieve the various derivatives:
-
-In derivative calculation we exploited the classical Perturbation Theory up to the first order.
-In extensive manner, we introduce a perturbation over \f$\lambda_{0}\f$ correlated with
-a pertubation of the states \f$\vert q_{0}\rangle \f$ (in bra-ket notation):
-\f{displaymath}
-[\mathbf{M}+d\mathbf{M}][\vert q_{0}\rangle + \vert dq_{0}\rangle ]=
-[\lambda_{0}+d\lambda_{0}][\vert q_{0}\rangle +\vert dq_{0}\rangle ]
-\f}
-Grouping the zero order we recollect the unperturbed equation(see before).
-To the first order:
-\f{displaymath}
-d\mathbf{M}q_{0}+\mathbf{M}\vert dq_{0}\rangle =d\lambda_{0}\vert q_{0}\rangle +\lambda_{0} \vert dq_{0}\rangle
-\f}
-Now we express \f$dq_{0}\f$ as linear combination of the other ortogonal eigenvectors:
-\f{displaymath}
-\vert dq_{0}\rangle =\sum_{j\neq0}c_{j}\vert q_{j}\rangle
-\f}
-thus we have
-\f{displaymath}
-d\mathbf{M}\vert q_{0}\rangle +\sum_{j\neq0}c_{j}\mathbf{M}\vert q_{j}\rangle=
-d\lambda_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\vert q_{j}\rangle
-\f}
-projecting onto the \f$q_{0}\f$ state and deleting the projection onto \f$\vert dq_{0}\rangle\f$ beacuse
-of ortogonality:
-\f{displaymath}
-\langle q_{0}\vert d\mathbf{M}\vert q_{0}\rangle +\sum_{j\neq0}c_{j}\lambda_{j}\langle q_{0} \vert q_{j}\rangle=
-d\lambda_{0}\langle q_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\langle q_{0}\vert q_{j}\rangle
-\f}
-we get
-\f{displaymath}
-\langle q_{0}\vert d\mathbf{M}\vert q_{0}\rangle=d\lambda_{0}
-\f}
-So, using simple chain rules:
-\f{displaymath}
-\langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle
-dr_{k}^{\gamma}=d\lambda_{0}
-\f}
-where here we used the notation \f$r_{k}^{\gamma}\f$ to denote an arbitrary position which can be
-\f$\tilde x_0 ,\tilde y_0,\tilde z_0\f$ or \f$\tilde x_1 ,\tilde y_1,\tilde z_1\f$
-we get
-\f{displaymath}
-\langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle
-=\frac{d\lambda_{0}}{dr_{k}^{\gamma}}
-\f}
-
-The derivatives of the matrix \f$\frac{d\mathbf{M}}{dr_{k}^{\gamma}} \f$ can be readily obtained via the
-chain rule
-\f{displaymath}
-\frac{d\mathbf{M}}{dr_{k}^{\gamma}}=\sum_{\l}^{nat}\sum_{\alpha}^{x,y,z} \frac{d\mathbf{M}}{dP_{l}^{\alpha}}\frac{dP_{l}^{\alpha}}{dr_{k}^{\gamma}} +\\
-\frac{d\mathbf{M}}{dM_{l}^{\alpha}}\frac{dM_{l}^{\alpha}}{dr_{k}^{\gamma}}
-\f}
-
-where \f$ M_{l}^{\alpha} \f$ corresponds to \f$ x_m^{l},y_m^{l},z_m^{l} \f$ and
-\f$ P_{l}^{\alpha} \f$ corresponds to \f$ x_p^{l},y_p^{l},z_p^{l} \f$ according to the \f$ \alpha \f$ component.
-*/
-
-
-class Kearsley
-{
- /// general log reference that needs to be initialized when constructed
- Log* log;
- /// position of atoms (first frame. In md is the running frame)
- std::vector<Vector> p0;
- /// position of atoms (second frame. In md is the reference frame)
- std::vector<Vector> p1;
- /// alignment weight: the rmsd/msd that it provides is only based on this scalar
- std::vector<double> align;
-
- bool com0_is_removed;
- bool com1_is_removed;
-
-public:
- /// error: the distance between two frames (might be rmsd/msd. See below)
- double err;
- /// displacement: the vector that goes from the p0 onto p1
- std::vector<Vector> diff0on1;
- /// displacement: the vector that goes from the p1 onto p0 (via inverse rotation)
- std::vector<Vector> diff1on0;
-
- /// center of mass of p0
- Vector com0;
- /// center of mass of p1
- Vector com1;
- /// position resetted wrt coms p0
- std::vector<Vector> p0reset;
- /// position resetted wrt coms p1
- std::vector<Vector> p1reset;
- /// position rotated: p0
- std::vector<Vector> p0rotated;
- /// position rotated: p1
- std::vector<Vector> p1rotated;
- /// rotation matrices p0 on p1 and reverse (p1 over p0)
- Tensor rotmat0on1,rotmat1on0;
- /// derivatives: derivative of the error respect p0
- std::vector<Vector> derrdp0;
- /// derivatives: derivative of the error respect p1
- std::vector<Vector> derrdp1;
- /// derivative of the rotation matrix
- /// note the dimension 3x3 x 3 x N
- std::vector<double> dmatdp0;
- std::vector<double> dmatdp1;
-
- /// constructor: need the two structure, the alignment vector and the log reference
- Kearsley( const std::vector<Vector> &p0, const std::vector<Vector> &p1, const std::vector<double> &align , Log* &log);
- /// switch the assignment of the structure p0 (e.g. at each md step)
- void assignP0(const std::vector<Vector> & p0);
- /// derivatives: derivative of the error respect p1
- void assignP1(const std::vector<Vector> & p1);
- /// transfer the alignment vector
- void assignAlign(const std::vector<double> & align);
- /// finite differences of all the relevant quantities: takes a bool which decides if giving back rmsd or not (msd in this case)
- void finiteDifferenceInterface(bool rmsd);
- // this makes the real calculation: the rmsd bool decides wether doing rmsd or msd
- double calculate( bool rmsd );
-};
-
-}
-
-#endif
-
diff --git a/src/tools/OptimalAlignment.cpp b/src/tools/OptimalAlignment.cpp
deleted file mode 100644
index eae4aabf6..000000000
--- a/src/tools/OptimalAlignment.cpp
+++ /dev/null
@@ -1,399 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- Copyright (c) 2011-2016 The plumed team
- (see the PEOPLE file at the root of the distribution for a list of names)
-
- See http://www.plumed.org for more information.
-
- This file is part of plumed, version 2.
-
- plumed is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- plumed is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with plumed. If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "OptimalAlignment.h"
-#include "Kearsley.h"
-#include "Log.h"
-#include <cmath>
-#include <iostream>
-#include <cstdlib>
-#include "Random.h"
-
-using namespace std;
-namespace PLMD{
-
-OptimalAlignment::OptimalAlignment( const std::vector<double> & align, const std::vector<double> & displace, const std::vector<Vector> & p0, const std::vector<Vector> & p1 , Log* &log )
-:log(log){
-
- // kearsley init to null
- mykearsley=NULL;
- if (mykearsley==NULL) {
- mykearsley=new Kearsley(p0,p1,align,log);
- }
- // copy the structure into place
- assignP0(p0);
- assignP1(p1);
- assignAlign(align);
- assignDisplace(displace);
-
- // basic check
- if(p0.size() != p1.size()){
- (this->log)->printf("THE SIZE OF THE TWO FRAMES TO BE ALIGNED ARE DIFFERENT\n");
- }
- // fast behaviour: if all the alignment and displacement are 1.0 then go for fast
- fast=true;
- for (unsigned i=0;i<align.size();i++ ){
- if(align[i]!=displace[i] || align[i]!=1.0)fast=false;
- }
-
-}
-
-OptimalAlignment::~OptimalAlignment(){
- if (mykearsley!=NULL)delete mykearsley;
-}
-
-void OptimalAlignment::assignP0( const std::vector<Vector> & p0 ){
- this->p0=p0;
- if(mykearsley!=NULL){mykearsley->assignP0(p0);}else{cerr<<"kearsley is not initialized"<<endl; exit(0);}
-}
-
-void OptimalAlignment::assignP1( const std::vector<Vector> & p1 ){
- this->p1=p1;
- if(mykearsley!=NULL){mykearsley->assignP1(p1);}else{cerr<<"kearsley is not initialized"<<endl; exit(0);}
-}
-
-void OptimalAlignment::assignAlign( const std::vector<double> & align ){
- this->align=align;
- if(mykearsley!=NULL){mykearsley->assignAlign(align);}else{cerr<<"kearsley is not initialized"<<endl; exit(0);}
-}
-
-void OptimalAlignment::assignDisplace( const std::vector<double> & displace ){
- this->displace=displace;
-}
-
-
-double OptimalAlignment::calculate(bool squared, std::vector<Vector> & derivatives){
-
- bool rmsd=!squared ;
-
- double err;
-
- // at this point everything should be already in place for calculating the alignment (p1,p2,align)
- // here everything is done with kearsley algorithm. Extension to other optimal alignment algos is
- // possible here below with a switch
-
- err=mykearsley->calculate(rmsd); // this calculates the MSD: transform into RMSD
-
- // check findiff alignment
- //mykearsley->finiteDifferenceInterface(rmsd);
-
- if(fast){
- //log->printf("Doing fast: ERR %12.6f \n",err);
- derrdp0=mykearsley->derrdp0;
- derrdp1=mykearsley->derrdp1;
- derivatives=derrdp0;
- }else{
- /// TODO this interface really sucks since is strongly asymmetric should be re-engineered.
- err=weightedAlignment(rmsd);
- //log->printf("Doing slow: ERR %12.6f \n",err);
- derivatives=derrdp0;
- }
- // destroy the kearsley object?
-
- return err;
-}
-
-#ifdef __INTEL_COMPILER
-#pragma intel optimization_level 2
-#endif
-/// this does the weighed alignment if the vector of alignment is different from displacement
-double OptimalAlignment::weightedAlignment( bool rmsd){
- double tmp0,tmp1,walign,wdisplace,const1,ret;
- unsigned i,k,l,m,n,o,oo,mm;
-
- unsigned natoms=p0.size();
-
- Kearsley *myk=mykearsley; /// just a shortcut
-
- /// TODO : these two blocks can be calculated once forever after the initialization (exception for certain methods?)
-
- /// clear derivatives
- if (derrdp0.size()!=natoms)derrdp0.resize(natoms);
- if (derrdp1.size()!=natoms)derrdp1.resize(natoms);
-
- // clear the container
- for(i=0;i<natoms;i++){
- derrdp0[i][0]=derrdp0[i][1]=derrdp0[i][2]=0.;
- derrdp1[i][0]=derrdp1[i][1]=derrdp1[i][2]=0.;
- }
-
- walign=0.;
- vector<int> alignmap;
- for(i=0;i<natoms;i++){
- if (align[i]>0.){
- alignmap.push_back(i);
- walign+=align[i];
- }
- if (align[i]<0.){cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl;exit(0);};
- }
-
- wdisplace=0.;
- vector<int> displacemap;
- for(i=0;i<natoms;i++){
- if (displace[i]>0.){
- displacemap.push_back(i);
- wdisplace+=displace[i];
- }
- if (displace[i]<0.){cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl;exit(0);};
- }
-
-
- tmp0=0.;
-
- vector<Vector> array_n_3;
- array_n_3.resize(natoms);
- for(i=0;i<array_n_3.size();i++)array_n_3[i][0]=array_n_3[i][1]=array_n_3[i][2]=0.;
-
- // err= (1/totdisplace) sum_k_l displace_k*((p0reset_k_l- sum_k_m rot_l_m*p1reset_k_m )**2)
-
- //for(kk=0;kk<displacemap.size();kk++){
- // k=displacemap[kk];
- for(k=0;k<natoms;k++){
-
- for(l=0;l<3;l++){
-
- tmp1=0.;
- // contribution from rotated reference frame //
- for(m=0;m<3;m++){
- tmp1-=myk->rotmat0on1[l][m]*myk->p1reset[k][m];
- }
-
- // contribution from running centered frame //
- tmp1+= myk->p0reset[k][l];
-
- array_n_3[k][l]=tmp1; // store coefficents for derivative usage//
- tmp0+=tmp1*tmp1*displace[k]; //squared distance added//
- }
-
- }
-
- tmp0=tmp0/wdisplace;
-
- // log->printf(" ERRR NEW %f \n",tmp0);
-
- ret=tmp0;
-
- /* DERIVATIVE CALCULATION:respect to running frame */
- for(k=0;k<natoms;k++){
- for(l=0;l<3;l++){
-
- tmp1 =2.*array_n_3[k][l]*displace[k]/wdisplace ; //ok
-
- const1=2.*align[k]/(walign*wdisplace);
-
- if(const1>0.){
- for(oo=0;oo<displacemap.size();oo++){
- o=displacemap[oo];
- tmp1 -=const1*array_n_3[o][l]*displace[o]; //ok
- }
- }
-
- for(mm=0;mm<displacemap.size();mm++){
- m=displacemap[mm];
- const1=2.* displace[m]/wdisplace ;
- for(n=0;n<3;n++){
- tmp0=0.;
- for(o=0;o<3;o++){
- int ind=n*3*3*natoms+o*3*natoms+l*natoms+k; //ok
- tmp0+=myk->dmatdp0[ind]*myk->p1reset[m][o];
- }
- tmp0*=-const1*array_n_3[m][n];
-
- tmp1+=tmp0;
- }
- }
-
- derrdp0[k][l]=tmp1;
-
- }
- }
- //exit(0);
-
- //return ret;
- bool do_frameref_der=true;
-
- // derivatives of
- // err= (1/totdisplace) sum_k_l displace_k*((p0reset_k_l- sum_m rot_l_m*p1reset_k_m )**2)
- // respect p1:
- // derr_dp1=(1/totdisplace) sum_k_l 2*displace_k*(p0reset_k_l- sum_m rot_l_m*p1reset_k_m )
- // *d/dp1 ( p0reset_k_l- sum_m rot_l_m*p1reset_k_m)
- // =
- // (1/totdisplace) sum_k_l 2*displace_k*(p0reset_k_l- sum_m rot_l_m*p1reset_k_m )*
- // *(d/dp1 p0reset_k_l
- // - sum_m (d/dp1 rot_l_m)*p1reset_k_m
- // - sum_m rot_l_m*(d/dp1 p1reset_k_m ) )
- // =
- // sum_k_l 2*displace_k/totdisplace* array_n_3_k_l
- // *(- sum_m (d/dp1 rot_l_m)*p1reset_k_m
- // - sum_m rot_l_m*(d/dp1 p1reset_k_m ) )
-
- if(do_frameref_der){
- for(k=0;k<natoms;k++){
-// for(kk=0;kk<displacemap.size();kk++){
-// k=displacemap[kk];
-
- for(l=0;l<3;l++){
-
- tmp1=0.;
- for(mm=0;mm<displacemap.size();mm++){
- m=displacemap[mm];
- const1=2.* displace[m]/wdisplace ;
- for(n=0;n<3;n++){
- tmp0=0.;
- for(o=0;o<3;o++){
- int ind=n*3*3*natoms+o*3*natoms+l*natoms+k;
- tmp0+=myk->dmatdp1[ind]*myk->p1reset[m][o];
- }
- tmp0*=-const1*array_n_3[m][n];
- tmp1+= tmp0;
- }
-
- }
-
- tmp0=0.;
- for(o=0;o<3;o++){
- tmp0+=array_n_3[k][o]*myk->rotmat0on1[o][l];
- }
- tmp1+=-tmp0*2.*displace[k]/wdisplace;
-
- tmp0=0.;
-
- for(mm=0;mm<displacemap.size();mm++){
- m=displacemap[mm];
- for(o=0;o<3;o++){
- tmp0+=array_n_3[m][o]*myk->rotmat0on1[o][l]*displace[m];
- }
- }
- tmp1 += tmp0*2.*align[k]/(walign*wdisplace);
-
- derrdp1[k][l]=tmp1;
-
- }
- }
- }
-
- /// probably not really the way it should be
- if (rmsd){
- ret=sqrt(ret);
- double tmp=0.5/ret;
- for(unsigned i=0;i<natoms;i++){
- derrdp0[i][0]=derrdp0[i][0]*tmp;
- derrdp0[i][1]=derrdp0[i][1]*tmp;
- derrdp0[i][2]=derrdp0[i][2]*tmp;
- derrdp1[i][0]=derrdp1[i][0]*tmp;
- derrdp1[i][1]=derrdp1[i][1]*tmp;
- derrdp1[i][2]=derrdp1[i][2]*tmp;
-
- }
- }
-
- return ret;
-}
-
-double OptimalAlignment::weightedFindiffTest( bool rmsd){
-
- Random rnd;
-
- log->printf("Entering rmsd finite difference test system\n ");
- log->printf("RMSD OR MSD: %s\n",(rmsd)?"rmsd":"msd");
- log->printf("-------------------------------------------\n");
- log->printf("TEST1: derivative of the value (derr_dr0/derr_dr1)\n");
- //// test 1
- double step=1.e-8,olderr,delta,err;
- vector<Vector> fakederivatives;
- fakederivatives.resize(p0.size());
- fast=false;
- // randomizing alignments and displacements
-/* for (i=0;i<p0.size();i++){
- // draw a random number
- delta=drand48();
- delta1=drand48();
- if(delta>delta1){
- align[i]=delta;
- }else{align[i]=0.;};
- delta=drand48();
- delta1=drand48();
- if(delta>delta1){
- displace[i]=delta;
- }else{displace[i]=0.;}
- }*/
- //// get initial value of the error and derivative of it
- assignAlign(align);
- assignDisplace(displace);
- olderr=calculate(rmsd, fakederivatives);
-
- log->printf("INITIAL ERROR VALUE: %e\n",olderr);
-
- // randomizing alignments and displacements
- log->printf("TESTING: derrdp0 \n");
-
- for(unsigned j=0;j<3;j++){
- for(unsigned i=0;i<derrdp0.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p0[i][j]+=delta;
- assignP0( p0 );
- err=calculate(rmsd, fakederivatives);
- p0[i][j]-=delta;
- assignP0( p0 );
- switch(j){
- case 0:
- log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
- case 1:
- log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
- case 2:
- log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
-
- }
- }
- }
-
- log->printf("TESTING: derrdp1 \n");
- for(unsigned j=0;j<3;j++){
- for(unsigned i=0;i<derrdp1.size();i++){
- // random displacement
- delta=(rnd.RandU01()-0.5)*2*step;
- p1[i][j]+=delta;
- assignP1( p1 );
- err=calculate(rmsd, fakederivatives);
- p1[i][j]-=delta;
- assignP1( p1 );
- switch(j){
- case 0:
- log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
- case 1:
- log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
- case 2:
- log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]);break;
-
- }
- }
- }
- exit(0);
-
-// This is to avoid warnings:
- return 0.0;
-
-}
-
-
-
-}
diff --git a/src/tools/OptimalAlignment.h b/src/tools/OptimalAlignment.h
deleted file mode 100644
index 2718d2e8a..000000000
--- a/src/tools/OptimalAlignment.h
+++ /dev/null
@@ -1,82 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- Copyright (c) 2011-2016 The plumed team
- (see the PEOPLE file at the root of the distribution for a list of names)
-
- See http://www.plumed.org for more information.
-
- This file is part of plumed, version 2.
-
- plumed is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- plumed is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with plumed. If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifndef __PLUMED_tools_OptimalAlignment_h
-#define __PLUMED_tools_OptimalAlignment_h
-
-#include "Vector.h"
-#include "Tensor.h"
-#include <vector>
-
-namespace PLMD{
-
-class Log;
-class Kearsley;
-
-/// A class that is intended to include or combine various optimal alignment algorithms
-
-class OptimalAlignment
-{
-private:
- /// a pointer to the object that performs the optimal alignment via quaternions
- Kearsley *mykearsley;
- /// displacement vector : a double that says if the coordinate should be used in calculating the RMSD/MSD
- std::vector<double> displace;
- /// alignment vector: a double that says if the atom has to be used in reset COM and makeing the alignment
- std::vector<double> align;
- /// position of one frame (generally the MD)
- std::vector<Vector> p0;
- /// position of the reference frames
- std::vector<Vector> p1;
- /// derivatives of the error respect to the p0 (MD running frame)
- std::vector<Vector> derrdp0;
- /// derivatives of the error respect to the p1 (static frame, do not remove: useful for SM)
- std::vector<Vector> derrdp1;
- /// the pointer to the logfile
- Log* log;
- /// a bool that decides to make the fast version (alignment vec= displacement vec) or the slower case
- bool fast;
-
-public:
- /// the contructor
- OptimalAlignment( const std::vector<double> & align, const std::vector<double> & displace, const std::vector<Vector> & p0, const std::vector<Vector> & p1 , Log* &log );
- /// the destructor: delete kearsley
- ~OptimalAlignment();
- /// assignment of the running frame p0
- void assignP0( const std::vector<Vector> & p0 );
- /// assignment to the reference frame p1
- void assignP1( const std::vector<Vector> & p1 );
- // this updates align runtime
- void assignAlign( const std::vector<double> & align );
- // this updates displace runtime
- void assignDisplace( const std::vector<double> & displace );
- /// this does the real calculation
- double calculate( bool rmsd, std::vector<Vector> & derivatives);
- /// this should perform the weighted alignment
- double weightedAlignment( bool rmsd);
- // a finite difference test
- double weightedFindiffTest( bool rmsd);
-};
-
-}
-
-#endif
-
--
GitLab