From 360b98c146ce520d591beb40bbc3398feaa05951 Mon Sep 17 00:00:00 2001
From: Sandro <sbottaro@sissa.it>
Date: Tue, 13 Sep 2016 10:47:49 +0200
Subject: [PATCH] Added ERMSD implementation
Notice that as of now it is not fully integrated in the plumed
reference/ module. However, it is possible to compute
ERMSD from a given structure and to use it as a CV
---
src/colvar/ERMSD.cpp | 187 +++++++++++++++++++++++
src/tools/ERMSD.cpp | 305 +++++++++++++++++++++++++++++++++++++
src/tools/ERMSD.h | 72 +++++++++
src/tools/MolDataClass.cpp | 9 ++
src/tools/PDB.cpp | 9 ++
src/tools/PDB.h | 2 +
src/tools/Tensor.h | 31 ++++
user-doc/bibliography.bib | 10 ++
8 files changed, 625 insertions(+)
create mode 100644 src/colvar/ERMSD.cpp
create mode 100644 src/tools/ERMSD.cpp
create mode 100644 src/tools/ERMSD.h
diff --git a/src/colvar/ERMSD.cpp b/src/colvar/ERMSD.cpp
new file mode 100644
index 000000000..4c7677599
--- /dev/null
+++ b/src/colvar/ERMSD.cpp
@@ -0,0 +1,187 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2015 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ 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 "Colvar.h"
+#include "core/PlumedMain.h"
+#include "ActionRegister.h"
+#include "tools/PDB.h"
+#include "tools/ERMSD.h"
+#include "core/Atoms.h"
+#include <iostream>
+
+using namespace std;
+
+namespace PLMD{
+namespace colvar{
+
+
+//+PLUMEDOC COLVAR ERMSD
+/*
+Calculate eRMSD with respect to a reference structure.
+
+eRMSD is a metric developed for measuring distances between three-dimensional RNA structures.
+The standard RMSD measure is highly inaccurate when measuring distances among three-dimensional
+structures of nucleic acids.
+It is not unusual, for example, that two RNA structures with low RMSD (i.e. less than 0.4nm) display a completely different network of base-base interactions.
+
+eRMSD measures the distance between structures by considering only the relative positions and orientations of nucleobases. The eRMSD can be considered as a vectorial version of contact maps and it is calculated as follows:
+
+1. Set up a local reference system in the center of the six-membered ring of each nucleobase in a molecule.
+ The xy plane lies on the plane of the nucleobase, and it is oriented such that the Watson-Crick interaction is always at \f$ \theta \approx 60^{\circ} \f$.
+
+2. Calculate all pairwise distance vectors \f$ \vec{r}_{i,j} \f$ among base centers.
+
+3. Rescale distance vectors as \f$ \tilde{\vec{r}}_{i,j} = (r_x/a,r_y/a,r_z/b) \f$, where \f$ a=b=5 \AA, c= 3 \AA\f$. This rescaling has the effect of weghting more deviations on the z-axis with respect to the x/y directions.
+
+4. Calculate the G vectors
+
+\f[
+\vec{G}(\tilde{\vec{r}}) = (\sin(\gamma \tilde{r}) \tilde{r}_x/\tilde{r},\sin(\gamma \tilde{r}) \tilde{r}_y/\tilde{r},\sin(\gamma \tilde{r}) \tilde{r}_z/\tilde{r}, 1+\cos(\gamma \tilde{r})) \times \Theta(\tilde{r}_{cutoff}-\tilde{r})
+\f]
+
+Here, \f$ \gamma = \pi/\tilde{r}_{cutoff}\f$ and \f$ \Theta \f$ is the Heaviside step function. The default cutoff is set to 2.4.
+
+5. The eRMSD between two structures \f$ \alpha \f$ and \f$ \beta \f$ reads
+
+\f[
+eRMSD = \sqrt{\frac{1}{N} \sum_{j,k} \vert \vec{G}(\tilde{\vec{r}}_{jk}^{\alpha}) - \vec{G}(\tilde{\vec{r}}_{jk}^{\beta}) \vert^2 }
+\f]
+
+Using the default cutoff, two structures with eRMSD of 0.7 or lower can be considered as significantly similar. A full description of the eRMSD can be found in \cite bott14
+
+\par Examples
+
+Calculate the eRMSD from reference structure reference.pdb using the default cutoff (2.4). The list of residues involved in the calculation has to be specified. In this example, the eRMSD is calculated
+considering residues 1,2,3,4,5,6.
+
+\verbatim
+eRMSD1: ERMSD REFERENCE=reference.pdb ATOMS=@lcs-1,@lcs-2,@lcs-3,@lcs-4,@lcs-5,@lcs-6
+\endverbatim
+
+*/
+//+ENDPLUMEDOC
+
+
+class ERMSD : public Colvar {
+
+
+ vector<Vector> derivs;
+ PLMD::ERMSD ermsd;
+ bool pbc;
+
+public:
+ ERMSD(const ActionOptions&);
+ virtual void calculate();
+ static void registerKeywords(Keywords& keys);
+};
+
+PLUMED_REGISTER_ACTION(ERMSD,"ERMSD")
+
+void ERMSD::registerKeywords(Keywords& keys){
+ Colvar::registerKeywords(keys);
+
+ keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
+ keys.add("compulsory","CUTOFF","2.4","only pairs of atoms closer than CUTOFF are considered in the calculation.");
+ keys.add("atoms","ATOMS","the list of atoms (use @lcs-nn)");
+ keys.add("optional","PAIRS","List of pairs considered. All pairs are considered if this value is not specified.");
+
+}
+
+ERMSD::ERMSD(const ActionOptions&ao):
+PLUMED_COLVAR_INIT(ao), pbc(true)
+{
+ string reference;
+ parse("REFERENCE",reference);
+ double cutoff=2.4;
+ parse("CUTOFF",cutoff);
+
+
+ bool nopbc(false);
+ parseFlag("NOPBC",nopbc);
+ pbc=!nopbc;
+
+ vector<AtomNumber> atoms_;
+ parseAtomList("ATOMS",atoms_);
+
+ vector<unsigned> pairs_;
+ parseVector("PAIRS",pairs_);
+ checkRead();
+
+ addValueWithDerivatives(); setNotPeriodic();
+
+ if(atoms_.size()<6) error("at least six atoms should be specified");
+ if(atoms_.size()%3!=0) error("Atoms are not multiple of 3");
+ if(pairs_.size()%2!=0) error("pairs are not multiple of 2");
+
+
+ //checkRead();
+ //log.printf(" of atoms");
+ //for(unsigned i=0;i<atoms.size();++i) log.printf(" %d",atoms[i].serial());
+ //requestAtoms(atoms);
+
+ // read everything in ang and transform to nm if we are not in natural units
+ PDB pdb;
+ if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
+ error("missing input file " + reference );
+ // store target_ distance
+ vector <Vector> reference_positions;
+ unsigned natoms = atoms_.size();
+ log.printf("Read %d atoms\n",natoms);
+
+ reference_positions.resize(natoms);
+ for(unsigned i=0;i<natoms;i++){
+ reference_positions[i] = pdb.getPosition(atoms_[i]);
+ //log.printf("%f %f %f \n",reference_positions[i][0],reference_positions[i][1],reference_positions[i][2]);
+ }
+
+// shift to count from zero
+ for(unsigned i=0;i<pairs_.size();++i) pairs_[i]--;
+
+ ermsd.setReference(reference_positions,pairs_,cutoff);
+
+ requestAtoms(atoms_);
+ derivs.resize(natoms);
+
+ log.printf(" reference from file %s\n",reference.c_str());
+ log.printf(" which contains %d atoms\n",natoms);
+
+}
+
+// calculator
+void ERMSD::calculate(){
+// set derivatives to zero
+ for(unsigned i=0;i<derivs.size();++i) {derivs[i].zero();}
+ double ermsdist;
+ Tensor virial;
+ makeWhole();
+ ermsdist=ermsd.calculate(getPositions(),getPbc(),derivs,virial);
+ setValue(ermsdist);
+
+ for(unsigned i=0;i<derivs.size();++i) {setAtomsDerivatives(i,derivs[i]);}
+
+ setBoxDerivativesNoPbc();
+
+ //setBoxDerivatives(virial);
+
+ }
+
+}
+}
diff --git a/src/tools/ERMSD.cpp b/src/tools/ERMSD.cpp
new file mode 100644
index 000000000..3a5b27876
--- /dev/null
+++ b/src/tools/ERMSD.cpp
@@ -0,0 +1,305 @@
+#include "ERMSD.h"
+#include "PDB.h"
+#include "Matrix.h"
+#include "Tensor.h"
+
+#include "Pbc.h"
+#include <cmath>
+#include <iostream>
+
+
+using namespace std;
+namespace PLMD{
+
+
+ void ERMSD::clear(){
+ reference_mat.clear();
+ }
+
+ //void ERMSD::calcLcs(const vector<Vector> & positions, vector<Vector> &)
+
+ void ERMSD::setReference(const vector<Vector> & reference, const std::vector<unsigned> & pairs_vec, double mycutoff){
+
+
+ natoms = reference.size();
+ nresidues = natoms/3;
+ unsigned npairs = pairs_vec.size()/2;
+ pairs.resize(npairs);
+ //for(unsigned i=0;i<2*npairs;++i) {
+ // std::cout << "CCC " << pairs_vec[i] << " ";
+ //}
+ for(unsigned i=0;i<npairs;++i) {
+
+ pairs[i].first = pairs_vec[2*i];
+ pairs[i].second = pairs_vec[2*i+1];
+ }
+
+ cutoff = mycutoff;
+ std::vector<TensorGeneric<4,3> > deri;
+ deri.resize(natoms*natoms);
+ reference_mat.resize(nresidues*nresidues);
+ Pbc fake_pbc;
+
+
+ calcMat(reference,fake_pbc,reference_mat,deri);
+
+ }
+
+
+ bool ERMSD::inPair(unsigned i, unsigned j){
+
+ //return true;
+ if(pairs.size()==0) return true;
+ for(unsigned idx=0;idx<pairs.size();idx++){
+ //std::cout << "AAA " << pairs[idx][0] << " " << pairs[idx][1] << "\n";
+ if(pairs[idx].first == i && pairs[idx].second == j) return true;
+ if(pairs[idx].second == i && pairs[idx].first == j) return true;
+ }
+ return false;
+ }
+//double ERMSD::calculate(const std::vector<Vector> & positions,
+// std::vector<Vector> &derivatives, Tensor& virial) const {
+
+// Pbc fake_pbc;
+// return ERMSD::calculate(positions,fake_pbc,derivatives,virial,false);
+// }
+
+
+
+ void ERMSD::calcMat(const std::vector<Vector> & positions,const Pbc& pbc, std::vector<Vector4d> &mat, std::vector<TensorGeneric<4,3> > &Gderi) {
+
+ std::vector<Vector3d> pos;
+ pos.resize(3*nresidues);
+
+ std::vector<Tensor3d> deri;
+ deri.resize(nresidues*9);
+
+ std::vector<Vector> centers;
+ centers.resize(nresidues);
+
+ unsigned at_idx = 0;
+ unsigned idx_deri = 0;
+
+ Tensor da_dxa = (2./3.)*Tensor::identity();
+ Tensor da_dxb = -(1./3.)*Tensor::identity();
+ Tensor da_dxc = -(1./3.)*Tensor::identity();
+
+ Tensor db_dxa = -(1./3.)*Tensor::identity();
+ Tensor db_dxb = (2./3.)*Tensor::identity();
+ Tensor db_dxc = -(1./3.)*Tensor::identity();
+
+ // Form factors - should this be somewhere else?
+
+ double w = 1./3.;
+ Vector form_factor = Vector(2.0,2.0,1.0/0.3);
+
+ for(unsigned res_idx=0;res_idx<natoms/3;res_idx++){
+
+
+ at_idx = 3*res_idx;
+ //center
+ for (unsigned j=0;j<3;j++){
+ centers[res_idx] += w*positions[at_idx+j];
+ }
+
+ Vector3d a = delta(centers[res_idx],positions[at_idx]);
+ Vector3d b = delta(centers[res_idx],positions[at_idx+1]);
+ Vector3d d = crossProduct(a,b);
+ double ianorm = 1./a.modulo();
+ double idnorm = 1./d.modulo();
+
+ // X vector: COM-C2
+ pos[at_idx] = a*ianorm;
+ // Z versor: C2 x (COM-C4/C6)
+ pos[at_idx+2] = d*idnorm;
+ // Y versor: Z x Y
+ pos[at_idx+1] = crossProduct(pos[at_idx+2],pos[at_idx]);
+
+ // Derivatives ////////
+ Tensor3d t1 = ianorm*(Tensor::identity()-extProduct(pos[at_idx],pos[at_idx]));
+ // dv1/dxa
+ deri[idx_deri] = (2./3. )*t1;
+ // dv1/dxb
+ deri[idx_deri+3] = -(1./3.)*t1;
+ // dv1/dxc
+ deri[idx_deri+6] = -(1./3.)*t1;
+
+ Tensor dd_dxa = VcrossTensor(a,db_dxa) -VcrossTensor(b,da_dxa);
+ Tensor dd_dxb = VcrossTensor(a,db_dxb)-VcrossTensor(b,da_dxb);
+ Tensor dd_dxc = VcrossTensor(a,db_dxc)-VcrossTensor(b,da_dxc);
+
+ // dv3/dxa
+ deri[idx_deri+2] = deriNorm(d,dd_dxa);
+ // dv3/dxb
+ deri[idx_deri+5] = deriNorm(d,dd_dxb);
+ // dv3/dxc
+ deri[idx_deri+8] = deriNorm(d,dd_dxc);
+
+ // dv2/dxa = dv3/dxa cross v1 + v3 cross dv1/dxa
+ deri[idx_deri+1] = (VcrossTensor(deri[idx_deri+2],pos[at_idx]) + \
+ VcrossTensor(pos[at_idx+2],deri[idx_deri]));
+ // dv2/dxb
+ deri[idx_deri+4] = (VcrossTensor(deri[idx_deri+5],pos[at_idx]) + \
+ VcrossTensor(pos[at_idx+2],deri[idx_deri+3]));
+ // dv2/dxc
+ deri[idx_deri+7] = (VcrossTensor(deri[idx_deri+8],pos[at_idx]) + \
+ VcrossTensor(pos[at_idx+2],deri[idx_deri+6]));
+
+ idx_deri += 9;
+ // End derivatives ///////
+
+ }
+
+
+ // Initialization (unnecessary?)
+ for (unsigned i1=0;i1<nresidues*nresidues;i1++){
+ for (unsigned i2=0;i2<4;i2++){
+ mat[i1][i2] = 0.0;
+ }
+ }
+
+ double maxdist = cutoff/form_factor[0];
+ double gamma = pi/cutoff;
+ unsigned idx;
+ unsigned idx1 = 0;
+ // Calculate mat
+ for (unsigned i=0;i<nresidues;i++){
+ for (unsigned j=0;j<nresidues;j++){
+
+ // skip i==j
+ if(inPair(i,j) and i != j){
+ //if(i!=j){
+
+
+ // Calculate normal distance first
+ Vector diff = delta(centers[i],centers[j]);
+ double d1 = diff.modulo();
+ //std::cout << inPair(i,j) << " " << i << " " << j << " "<< d1 <<"\n";
+ //std::cout << inPair(i,j) << " " << i << " " << j << " "<< d1 <<"\n";
+ if(d1<maxdist){
+
+ // calculate r_tilde_ij
+ Vector3d rtilde;
+ for (unsigned k=0;k<3;k++){
+ for (unsigned l=0;l<3;l++){
+ rtilde[l] += pos[3*i+l][k]*diff[k]*form_factor[l];
+ }
+ }
+ double rtilde_norm = rtilde.modulo();
+
+ double irnorm = 1./rtilde_norm;
+
+ // ellipsoidal cutoff
+ if(rtilde_norm < cutoff){
+ idx = i*nresidues + j;
+ //std::cout << i << " " << j << " " << rtilde_norm << " " << idx <<"\n";
+
+
+ // fill 4d matrix
+ double dummy = sin(gamma*rtilde_norm)/(rtilde_norm*gamma);
+ mat[idx][0] = dummy*rtilde[0];
+ mat[idx][1] = dummy*rtilde[1];
+ mat[idx][2] = dummy*rtilde[2];
+ mat[idx][3] = (1.+ cos(gamma*rtilde_norm))/gamma;
+
+ // Derivative (drtilde_dx)
+ std::vector<Tensor3d> drtilde_dx;
+ drtilde_dx.resize(6);
+ unsigned pos_idx = 3*i;
+ unsigned deri_idx = 9*i;
+ for (unsigned at=0;at<3;at++){
+ for (unsigned l=0;l<3;l++){
+ Vector3d rvec = form_factor[l]*((pos[pos_idx+l])/3.);
+ Vector3d vvec = form_factor[l]*(matmul(deri[deri_idx+3*at+l],diff));
+ drtilde_dx[at].setRow(l,vvec-rvec);
+ drtilde_dx[at+3].setRow(l,rvec);
+ }
+ }
+
+ //std::vector<TensorGeneric<4,3> > dG_dx;
+ //dG_dx.resize(6);
+
+ double dummy1 = (cos(gamma*rtilde_norm) - dummy);
+
+ idx1 = i*nresidues*6 + j*6;
+
+ for (unsigned l=0;l<6;l++){
+ //std::cout << i << " " << j << " " << idx1 << " " << idx1+l << "\n";
+
+ // components 1,2,3
+ // sin(gamma*|rtilde|)/gamma*|rtilde|*d_rtilde +
+ // + ((d_rtilde*r_tilde/r_tilde^2) out r_tilde)*
+ // (cos(gamma*|rtilde| - sin(gamma*|rtilde|)/gamma*|rtilde|))
+ Vector3d rdr = matmul(rtilde,drtilde_dx[l]);
+ Tensor tt = dummy*drtilde_dx[l] + (dummy1*irnorm*irnorm)*Tensor(rtilde,rdr);
+ for (unsigned m=0;m<3;m++){
+ // Transpose here
+ //dG_dx[l].setRow(m,tt.getRow(m));
+ Gderi[idx1+l].setRow(m,tt.getRow(m));
+ }
+ // component 4
+ // - sin(gamma*|rtilde|)/|rtilde|*(r_tilde*d_rtilde)
+ //dG_dx[l].setRow(3,-dummy*gamma*rdr);
+ Gderi[idx1+l].setRow(3,-dummy*gamma*rdr);
+ }
+
+
+
+
+ }
+ }
+ }
+
+ }
+ }
+
+ }
+
+
+ double ERMSD::calculate(const std::vector<Vector> & positions,const Pbc& pbc,\
+ std::vector<Vector> &derivatives, Tensor& virial) {
+
+
+ double ermsd=0.;
+ std::vector<Vector4d> mat;
+ mat.resize(nresidues*nresidues);
+
+ std::vector<TensorGeneric<4,3> > Gderi;
+ Gderi.resize(natoms*natoms);
+
+ calcMat(positions,pbc,mat,Gderi);
+
+ unsigned idx1 = 0;
+ for(unsigned i=0;i<nresidues;i++){
+ for(unsigned j=0;j<nresidues;j++){
+ unsigned ii = i*nresidues + j;
+
+ Vector4d dd = delta(reference_mat[ii],mat[ii]);
+ double val = dd.modulo2();
+ //std::cout << "AAA " << i << " " << j << " " << ii << " "<< val << "\n";
+
+ if(val>0.0 && i!=j){
+
+ for(unsigned k=0;k<3;k++){
+ idx1 = i*nresidues*6 + j*6 + k;
+
+ derivatives[3*i+k] += matmul(dd,Gderi[idx1]);
+ derivatives[3*j+k] += matmul(dd,Gderi[idx1+3]);
+ }
+ ermsd += val;
+ }
+ }
+ }
+
+ //std::cout << ermsd << " ";
+ //if(pairs.size()!=0) nresidues=pairs.size();
+ //std::cout << ermsd << " " << nresidues;
+ ermsd = sqrt(ermsd/nresidues);
+ double iermsd = 1.0/(ermsd*nresidues);
+ for(unsigned i=0;i<natoms;++i){derivatives[i] *= iermsd;}
+
+ return ermsd;
+}
+
+
+}
diff --git a/src/tools/ERMSD.h b/src/tools/ERMSD.h
new file mode 100644
index 000000000..5ba98f844
--- /dev/null
+++ b/src/tools/ERMSD.h
@@ -0,0 +1,72 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ 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_ERMSD_h
+#define __PLUMED_tools_ERMSD_h
+
+#include "Tensor.h"
+#include "Vector.h"
+#include <vector>
+#include <limits>
+#include <string>
+#include <map>
+#include <utility>
+
+namespace PLMD{
+
+class PDB;
+class Pbc;
+
+/// A class that implements ERMSD calculations
+class ERMSD {
+ //std::map< std::pair <unsigned,unsigned> , double> targets;
+ //unsigned natoms;
+ std::vector<Vector4d> reference_mat;
+ unsigned natoms;
+ unsigned nresidues;
+ std::vector<std::pair <unsigned,unsigned> > pairs;
+ double cutoff;
+
+ public:
+/// Constructor
+ ERMSD(): natoms(0),nresidues(0) {}
+
+/// clear the structure
+ void clear();
+
+ bool inPair(unsigned i, unsigned j);
+
+/// set reference coordinates
+ void setReference(const std::vector<Vector> & reference, const std::vector<unsigned> & pairs_vec,double mycutoff=0.24);
+
+ void calcMat(const std::vector<Vector> & positions, const Pbc& pbc,std::vector<Vector4d> &mat,std::vector<TensorGeneric<4,3> > & Gderivatives);
+
+/// Compute ermsd ( no pbc )
+// double calculate(const std::vector<Vector> & positions,
+// std::vector<Vector> &derivatives, Tensor& virial) const ;
+/// Compute ermsd ( with pbc )
+ double calculate(const std::vector<Vector>& positions, const Pbc& pbc,std::vector<Vector> &derivatives, Tensor& virial);
+};
+
+}
+
+#endif
+
diff --git a/src/tools/MolDataClass.cpp b/src/tools/MolDataClass.cpp
index 01adb71fb..49f0e0fac 100644
--- a/src/tools/MolDataClass.cpp
+++ b/src/tools/MolDataClass.cpp
@@ -368,6 +368,15 @@ void MolDataClass::specialSymbol( const std::string& type, const std::string& sy
numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N7",resnum,chainid));
numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C8",resnum,chainid));
} else plumed_error();
+ } else if( name=="lcs" ) {
+ numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
+ if(basetype=="T" || basetype=="U" || basetype=="C"){
+ numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
+ numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
+ } else if(basetype=="G" || basetype=="A"){
+ numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
+ numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
+ } else plumed_error();
} else numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
}
}
diff --git a/src/tools/PDB.cpp b/src/tools/PDB.cpp
index 38f113605..3cd3de176 100644
--- a/src/tools/PDB.cpp
+++ b/src/tools/PDB.cpp
@@ -296,5 +296,14 @@ Log& operator<<(Log& ostr, const PDB& pdb){
return ostr;
}
+Vector PDB::getPosition(AtomNumber a)const{
+ std::map<AtomNumber,unsigned>::const_iterator p;
+ p=number2index.find(a);
+ if(p==number2index.end()) plumed_merror("atom not available");
+ else return positions[p->second];
+}
+
+
+
}
diff --git a/src/tools/PDB.h b/src/tools/PDB.h
index 5d1ee0e85..fdea48ce0 100644
--- a/src/tools/PDB.h
+++ b/src/tools/PDB.h
@@ -108,6 +108,8 @@ public:
unsigned getNumberOfAtomBlocks() const ;
/// Set the position array
void setPositions(const std::vector<Vector> &v);
+/// Access to the position array
+ Vector getPosition(AtomNumber a)const;
};
}
diff --git a/src/tools/Tensor.h b/src/tools/Tensor.h
index ecb30d0bd..a1f15859d 100644
--- a/src/tools/Tensor.h
+++ b/src/tools/Tensor.h
@@ -171,6 +171,10 @@ public:
friend TensorGeneric<n_,m_> extProduct(const VectorGeneric<n>&,const VectorGeneric<m>&);
friend TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&,const VectorGeneric<3>&);
friend TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&,const VectorGeneric<3>&);
+ friend TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
+ friend TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&,const VectorGeneric<3>&);
+/// Derivative of a normalized vector
+ friend TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
/// << operator.
/// Allows printing tensor `t` with `std::cout<<t;`
template<unsigned n_,unsigned m_>
@@ -468,6 +472,33 @@ typedef TensorGeneric<4,4> Tensor4d;
/// \ingroup TOOLBOX
typedef Tensor3d Tensor;
+inline
+TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2){
+
+ TensorGeneric<3,3> t;
+ for(unsigned i=0;i<3;i++){
+ t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
+ }
+ return t;
+}
+
+inline
+TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&v2,const VectorGeneric<3>&v1){
+ TensorGeneric<3,3> t;
+ for(unsigned i=0;i<3;i++){
+ t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
+ }
+ return t;
+}
+
+
+inline
+TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2){
+ // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
+ double over_norm = 1./v1.modulo();
+ return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
+}
+
diff --git a/user-doc/bibliography.bib b/user-doc/bibliography.bib
index b15ee2505..222f48f16 100644
--- a/user-doc/bibliography.bib
+++ b/user-doc/bibliography.bib
@@ -399,6 +399,16 @@ journal = { J. Chem. Phys. }
Volume = {126},
Year = {2007}}
+@article{bott14,
+ Title={The role of nucleobase interactions in RNA structure and dynamics},
+ Author={Bottaro, Sandro and Di Palma, Francesco and Bussi, Giovanni},
+ Journal={Nucleic acids research},
+ Pages={13306--13314},
+ Number={42},
+ Volume={21},
+ Year={2014},
+ Publisher={Oxford Univ Press}}
+
@inproceedings{92avconAPAF,
Address = {La Jolla, CA},
Author = {G. Abdulla and S. Patel and M. Abrams and E. A. Fox},
--
GitLab