From 253e263c6541099696117f0997ece060fe7c8e43 Mon Sep 17 00:00:00 2001
From: Jana Pazurikova <pazurikova@ics.muni.cz>
Date: Thu, 14 Sep 2017 15:36:18 +0200
Subject: [PATCH] Close structure - modify PathMSDBase
---
src/colvar/PathMSDBase.cpp | 113 +++++++++++++++++++++++++++++++++++--
src/colvar/PathMSDBase.h | 13 +++++
2 files changed, 120 insertions(+), 6 deletions(-)
diff --git a/src/colvar/PathMSDBase.cpp b/src/colvar/PathMSDBase.cpp
index 69f1aecf4..f908ae4a4 100644
--- a/src/colvar/PathMSDBase.cpp
+++ b/src/colvar/PathMSDBase.cpp
@@ -41,18 +41,27 @@ void PathMSDBase::registerKeywords(Keywords& keys){
keys.add("compulsory","REFERENCE","the pdb is needed to provide the various milestones");
keys.add("optional","NEIGH_SIZE","size of the neighbor list");
keys.add("optional","NEIGH_STRIDE","how often the neighbor list needs to be calculated in time units");
+ keys.add("optional", "EPSILON", "(default=-1) the maximum distance between the close and the current structure, the positive value turn on the close structure method");
+ keys.add("optional", "LOG-CLOSE", "(default=0) value 1 enables logging regarding the close structure");
+ keys.add("optional", "DEBUG-CLOSE", "(default=0) value 1 enables extensive debugging info regarding the close structure, the simulation will run much slower");
}
PathMSDBase::PathMSDBase(const ActionOptions&ao):
PLUMED_COLVAR_INIT(ao),
neigh_size(-1),
neigh_stride(-1),
-nframes(0)
+nframes(0),
+epsilonClose(-1),
+debugClose(0),
+logClose(0)
{
parse("LAMBDA",lambda);
parse("NEIGH_SIZE",neigh_size);
parse("NEIGH_STRIDE",neigh_stride);
parse("REFERENCE",reference);
+ parse("EPSILON", epsilonClose);
+ parse("LOG-CLOSE", logClose);
+ parse("DEBUG-CLOSE", debugClose);
// open the file
FILE* fp=fopen(reference.c_str(),"r");
@@ -85,6 +94,9 @@ nframes(0)
fclose (fp);
log<<"Found TOTAL "<<nframes<< " PDB in the file "<<reference.c_str()<<" \n";
if(nframes==0) error("at least one frame expected");
+ //set up rmsdRefClose, initialize it to the first structure loaded from reference file
+ rmsdPosClose.set(pdbv[0], "OPTIMAL");
+ firstPosClose = true;
}
if(neigh_stride>0 || neigh_size>0){
if(neigh_size>int(nframes)){
@@ -98,6 +110,26 @@ nframes(0)
log.printf(" Neighbor list NOT enabled \n");
}
+ if (epsilonClose > 0) {
+ log.printf(" Computing with the close structure, epsilon = %lf\n", epsilonClose);
+ log << " Bibliography " << plumed.cite("Pazurikova J, Krenek A, Spiwok V, Simkova M J. Chem. Phys. 146, 115101 (2017)") << "\n";
+ }
+ else {
+ debugClose = 0;
+ logClose = 0;
+ }
+ if (debugClose)
+ log.printf(" Extensive debug info regarding close structure turned on\n");
+
+ rotationRefClose = new Tensor[nframes];
+ drotationPosCloseDrr01 = new Tensor[9];
+ savedIndices = vector<unsigned>(nframes);
+
+}
+
+PathMSDBase::~PathMSDBase(){
+ delete[] rotationRefClose;
+ delete[] drotationPosCloseDrr01;
}
void PathMSDBase::calculate(){
@@ -124,21 +156,90 @@ void PathMSDBase::calculate(){
plumed_assert(nframes>0);
plumed_assert(imgVec.size()>0);
+ Tensor* tmp_rotationRefClose = new Tensor[nframes];
+
+ if (epsilonClose > 0){
+ //compute rmsd between positions and close structure, save rotation matrix, drotation_drr01
+ double posclose = rmsdPosClose.calc_Rot_DRotDRr01(getPositions(), rotationPosClose, drotationPosCloseDrr01, true);
+ //if we compute for the first time or the existing close structure is too far from current structure
+ if (firstPosClose || (posclose > epsilonClose)){
+ //set the current structure as close one for a few next steps
+ if (logClose)
+ log << "PLUMED-CLOSE: new close structure, rmsd pos close " << posclose << "\n";
+ rmsdPosClose.clear();
+ rmsdPosClose.setReference(getPositions());
+ //as this is a new close structure, we need to save the rotation matrices fitted to the reference structures
+ // and we need to accurately recalculate for all reference structures
+ computeRefClose = true;
+ imgVec.resize(nframes);
+#pragma simd
+ for(unsigned i=0;i<nframes;i++){
+ imgVec[i].property=indexvec[i];
+ imgVec[i].index=i;
+ }
+ firstPosClose = false;
+ }
+ else{
+ //the current structure is pretty close to the close structure, so we use saved rotation matrices to decrease the complexity of rmsd comuptation
+ if (debugClose)
+ log << "PLUMED-CLOSE: old close structure, rmsd pos close " << posclose << "\n";
+ computeRefClose = false;
+ }
+ }
+
std::vector<double> tmp_distances(imgVec.size(),0.0);
std::vector<Vector> tmp_derivs;
// this array is a merge of all tmp_derivs, so as to allow a single comm.Sum below
std::vector<Vector> tmp_derivs2(imgVec.size()*nat);
// if imgVec.size() is less than nframes, it means that only some msd will be calculated
- for(unsigned i=rank;i<imgVec.size();i+=stride){
-// store temporary local results
- tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
- plumed_assert(tmp_derivs.size()==nat);
- for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
+ if (epsilonClose > 0){
+ if (computeRefClose){
+ //recompute rotation matrices accurately
+ for(unsigned i=rank;i<imgVec.size();i+=stride){
+ tmp_distances[i] = msdv[imgVec[i].index].calc_Rot(getPositions(), tmp_derivs, tmp_rotationRefClose[imgVec[i].index], true);
+ plumed_assert(tmp_derivs.size()==nat);
+ for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
+ }
+ }
+ else{
+ //approximate distance with saved rotation matrices
+ for(unsigned i=rank;i<imgVec.size();i+=stride){
+ tmp_distances[i] = msdv[imgVec[i].index].calculateWithCloseStructure(getPositions(), tmp_derivs, rotationPosClose, rotationRefClose[imgVec[i].index], drotationPosCloseDrr01, true);
+ plumed_assert(tmp_derivs.size()==nat);
+ for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
+ if (debugClose){
+ double withclose = tmp_distances[i];
+ RMSD opt;
+ opt.setType("OPTIMAL");
+ opt.setReference(msdv[imgVec[i].index].getReference());
+ vector<Vector> ders;
+ double withoutclose = opt.calculate(getPositions(), ders, true);
+ float difference = fabs(withoutclose-withclose);
+ log.printf("PLUMED-CLOSE: difference original %lf - with close %lf = %lf, step %d, i %d imgVec[i].index %d \n", withoutclose, withclose, difference, getStep(), i, imgVec[i].index);
+ }
+ }
+ }
}
+ else{
+ // store temporary local results
+ for(unsigned i=rank;i<imgVec.size();i+=stride){
+ tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
+ plumed_assert(tmp_derivs.size()==nat);
+ for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
+ }
+ }
+
// reduce over all processors
comm.Sum(tmp_distances);
comm.Sum(tmp_derivs2);
+ if (epsilonClose > 0 && computeRefClose) {
+ comm.Sum(tmp_rotationRefClose, nframes);
+ for (unsigned i=0; i<nframes;i++) {
+ rotationRefClose[i] = tmp_rotationRefClose[i];
+ }
+ }
+ delete [] tmp_rotationRefClose;
// assign imgVec[i].distance and imgVec[i].distder
for(unsigned i=0;i<imgVec.size();i++){
imgVec[i].distance=tmp_distances[i];
diff --git a/src/colvar/PathMSDBase.h b/src/colvar/PathMSDBase.h
index 4451bffc7..d01c12854 100644
--- a/src/colvar/PathMSDBase.h
+++ b/src/colvar/PathMSDBase.h
@@ -67,6 +67,18 @@ class PathMSDBase : public Colvar {
std::vector<Vector> derivs_s;
std::vector<Vector> derivs_z;
std::vector <ImagePath> imgVec; // this can be used for doing neighlist
+
+ //variables used for the close structure method, i is the number of reference structures
+ double epsilonClose; //the maximal distance between the close and the current structure before reassignment
+ int debugClose; //turns on debug mode
+ int logClose; //turns on logging
+ RMSD rmsdPosClose; //rmsd between the current and the close structure
+ bool firstPosClose; //flag indicating the first time we need to calculate the distance between the close and the current structure
+ bool computeRefClose; //flag indicating necessity to recompute accurately all distances and rotation matrices between the close structure and reference str
+ Tensor *rotationRefClose; //Tensor[i] saved rotation matrices between the close structure and reference structures
+ Tensor rotationPosClose; //rotation matrix between the close and the current structure
+ Tensor *drotationPosCloseDrr01; //Tensor[3][3]; //derivation of the rotation matrix w.r.t rr01, necessary for calculation of derivations
+ std::vector<unsigned> savedIndices; //saved indices of imgVec from previous steps, used for recalculating after neighbourlist update
protected:
std::vector<PDB> pdbv;
std::vector<std::string> labels;
@@ -74,6 +86,7 @@ protected:
unsigned nframes;
public:
explicit PathMSDBase(const ActionOptions&);
+ ~PathMSDBase();
// active methods:
virtual void calculate();
// virtual void prepare();
--
GitLab