Enable EMST in all colvars using makeWhole()

CHANGES/v2.10.md

@@ -41,6 +41,7 @@ Before switching to version 2.10, users are invited to carefully read the follow
  - Action CENTER_OF_MULTICOLVAR no longer exists. You now simply use \ref CENTER with the PHASES option and a vector as input for the weights.
  - The rational switching function when used with MM!=2*MM now returns a more correct derivative around `dist~=R0`
  - Dimensionality reduction methods and landmark selection actions have a new syntax. A good introduction that explains how to use these actions can be found in [this tutorial](https://plumed-school.github.io/lessons/21/006/data/DIMENSIONALITY.html)
+ - When using \ref MOLINFO with the `WHOLE` flag, PBCs in the following actions will be reconstructed using a minimum spanning tree based on the coordinates stored in the MOLFILE reference pdb.
 
 - Places where we strongly recommend using the new sytax:
  - If you are using \ref DFSCLUSTERING and the \ref CLUSTER_PROPERTIES or \ref CLUSTER_DISTRIBUTION actions you are strongly encouraged to read: https://plumed-school.github.io/lessons/23/001/data/Clusters.html

regtest/basic/rt-emst/COLVAR.reference

@@ -1,3 +1,3 @@
-#! FIELDS time r1 r2 r3 r4 r5 r6
- 0.000000 12.624019 0.000000 12.624019 12.624019 12.624019 0.000000
- 1.000000 12.624019 8.981462 12.624019 12.624019 12.624019 0.000000
+#! FIELDS time r1 r2 r3 r4 r5 r6 r7
+ 0.000000 12.624019 0.000000 12.624019 12.624019 0.000000 0.000000 12.624019
+ 1.000000 12.624019 8.981462 12.624019 12.624019 0.000000 0.000000 12.624019

regtest/basic/rt-emst/plumed.dat

@@ -13,17 +13,22 @@ r3: RMSD REFERENCE=ref.pdb
 r4: RMSD REFERENCE=ref.pdb NOPBC
 
 MOLINFO STRUCTURE=ref.pdb WHOLE
+# here I leave the EMST flag which is not necessary anymore and redundant
 WHOLEMOLECULES ENTITY0=1-123 EMST
 DUMPATOMS ATOMS=1-123 FILE=right.xyz PRECISION=3
-# this is non zero (incorrectly reconstructed)
+
 r5: RMSD REFERENCE=ref.pdb
 # finally, this is always zero
 r6: RMSD REFERENCE=ref.pdb NOPBC
 
+MOLINFO STRUCTURE=ref.pdb # this is to make sure RMSD does not see this
+# this is non zero (incorrectly reconstructed)
+r7: RMSD REFERENCE=ref.pdb
 
 # check EMST with a subset of the atoms in the PDB
 MOLINFO STRUCTURE=ref_even.pdb WHOLE
-WHOLEMOLECULES ENTITY0=2-123 EMST
+# here I remove the EMST flag which is not necessary anymore and redundant
+WHOLEMOLECULES ENTITY0=2-123 # EMST 
 
 # note that a complete PDB is required now to avoid errors in DUMPATOMS
 MOLINFO STRUCTURE=ref.pdb

src/colvar/Dipole.cpp

@@ -33,7 +33,7 @@ Calculate the dipole moment for a group of atoms.
 
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.5,
-by considering the ordered list of atoms and rebuilding the molecule with a procedure
+by rebuilding the molecule with a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/colvar/Gyration.cpp

@@ -55,7 +55,7 @@ The radius of gyration usually makes sense when atoms used for the calculation
 are all part of the same molecule.
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.2,
-by considering the ordered list of atoms and rebuilding the broken entities using a procedure
+by rebuilding the broken entities using a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/colvar/PathMSD.cpp

@@ -37,7 +37,7 @@ in input ("sss" component) and the distance from them ("zzz" component).
 
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.5,
-by considering the ordered list of atoms and rebuilding molecules with a procedure
+by rebuilding molecules with a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/colvar/PropertyMap.cpp

@@ -46,7 +46,7 @@ where the parameters \f$X_i\f$ and \f$Y_i\f$ are provided in the input pdb (al
 
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.5,
-by considering the ordered list of atoms and rebuilding molecules using a procedure
+by rebuilding molecules using a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/colvar/RMSD.cpp

@@ -129,7 +129,7 @@ that are available in plumed. More information on these various methods can be
 
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.5,
-by considering the ordered list of atoms and rebuilding molecules using a procedure
+by rebuilding molecules using a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/core/ActionAtomistic.cpp

@@ -32,6 +32,7 @@
 #include "tools/Exception.h"
 #include "tools/Pbc.h"
 #include "tools/PDB.h"
+#include "tools/Tree.h"
 
 namespace PLMD {
 
@@ -49,6 +50,7 @@ ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
  chargesWereSet(false)
 {
  ActionWithValue* bv = plumed.getActionSet().selectWithLabel<ActionWithValue*>("Box");
+ moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
  if( bv ) boxValue=bv->copyOutput(0);
  // We now get all the information about atoms that are lying about
  std::vector<ActionWithValue*> vatoms = plumed.getActionSet().select<ActionWithValue*>();
@@ -82,6 +84,8 @@ void ActionAtomistic::registerKeywords( Keywords& keys ) {
 
 void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
  plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
+ // this makes the EMST invalid so that it will be regenerated at next request
+ tree.reset();
  int nat=a.size();
  indexes=a;
  positions.resize(nat);
@@ -275,7 +279,6 @@ void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::
  if(!groupfound) plumed_error()<<"group has not been found in index file";
  ok=true;
  } else {
- auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
  if( moldat ) {
  std::vector<AtomNumber> atom_list; moldat->interpretSymbol( symbol, atom_list );
  if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
@@ -432,10 +435,22 @@ unsigned ActionAtomistic::getTotAtoms()const {
 }
 
 void ActionAtomistic::makeWhole() {
- for(unsigned j=0; j<positions.size()-1; ++j) {
- const Vector & first (positions[j]);
- Vector & second (positions[j+1]);
- second=first+pbcDistance(first,second);
+ if(moldat && moldat->isWhole()) {
+ // make sure the tree has been constructed
+ if(!tree) tree=std::make_unique<Tree>(moldat);
+ const auto & tree_indexes=tree->getTreeIndexes();
+ const auto & root_indexes=tree->getRootIndexes();
+ for(unsigned j=0; j<root_indexes.size(); j++) {
+ const Vector & first (positions[root_indexes[j]]);
+ Vector & second (positions[tree_indexes[j]]);
+ second=first+pbcDistance(first,second);
+ }
+ } else {
+ for(unsigned j=0; j<positions.size()-1; ++j) {
+ const Vector & first (positions[j]);
+ Vector & second (positions[j+1]);
+ second=first+pbcDistance(first,second);
+ }
  }
 }
 

src/core/ActionAtomistic.h

@@ -34,6 +34,8 @@ namespace PLMD {
 
 class Pbc;
 class PDB;
+class GenericMolInfo;
+class Tree;
 
 namespace colvar {
 class SelectMassCharge;
@@ -76,6 +78,10 @@ class ActionAtomistic :
  bool donotretrieve;
  bool donotforce;
 
+ // EMST
+ GenericMolInfo* moldat{nullptr};
+ std::unique_ptr<Tree> tree;
+
 /// Values that hold information about atom positions and charges
  std::vector<Value*> xpos, ypos, zpos, masv, chargev;
  void updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l );

src/generic/FitToTemplate.cpp

@@ -78,7 +78,7 @@ this action is performed at every MD step.
 
 When running with periodic boundary conditions, the atoms should be
 in the proper periodic image. This is done automatically since PLUMED 2.5,
-by considering the ordered list of atoms and rebuilding the molecules using a procedure
+by rebuilding the molecules using a procedure
 that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
 rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
 which actually modifies the coordinates stored in PLUMED.

src/generic/MolInfo.cpp

@@ -48,6 +48,16 @@ generated with `gmx editconf -f topol.tpr -o reference.pdb`.
 
 More information of the PDB parser implemented in PLUMED can be found \ref pdbreader "at this page".
 
+If the flag `WHOLE` is used, the reference structure is assumed to be whole, i.e. not broken by PBC.
+**This will impact PDB calculations in all the actions following the MOLINFO action!**.
+In particular, actions that reconstruct molecules locally will do the reconstruction using a
+minimum spanning tree rather than the order in which atoms are provided.
+Notice that PLUMED 2.8 and 2.9 this would only have affected the behavior of \ref WHOLEMOLECULES
+actions, whereas starting with PLUMED 2.10 all actions that reconstruct PBCs are affected.
+If you want a variable to ignore this, you can just repeate the MOLINFO action without
+the `WHOLE` flag: all the following actions will use the traditional reconstruction algorithm,
+based on the order in which atoms are listed.
+
 Providing `MOLTYPE=protein`, `MOLTYPE=rna`, or `MOLTYPE=dna` will instruct plumed to look
 for known residues from these three types of molecule. In other words, this is available for
 historical reasons and to allow future extensions where alternative lists will be provided.

src/generic/WholeMolecules.cpp

@@ -57,7 +57,9 @@
 know exactly what you are doing, leave the default stride (1), so that
 this action is performed at every MD step.
 
-The way WHOLEMOLECULES modifies each of the listed entities is this:
+The behavior of WHOLEMOLECULES is affected by the last \ref MOLINFO action
+present in the input file before WHOLEMOLECULES. Specifically, if the
+\ref MOLINFO action does not have a `WHOLE` flag, then the behavior is the following:
 - First atom of the list is left in place
 - Each atom of the list is shifted by a lattice vectors so that it becomes as close as possible
  to the previous one, iteratively.
@@ -67,6 +69,11 @@
 This can be usually achieved selecting consecutive atoms (1-100), but it is also possible
 to skip some atoms, provided consecutive chosen atoms are close enough.
 
+If instead the \ref MOLINFO action does have a `WHOLE` flag, then a minimum spanning tree
+is built based on the atoms passed to WHOLEMOLECULES using the coordinates in the PDB
+passed to \ref MOLINFO as a reference, and this tree is used to reconstruct PBCs.
+This approach is more robust when dealing with complexes of multiple molecules.
+
 \par Examples
 
 This command instructs plumed to reconstruct the molecule containing atoms 1-20
@@ -129,7 +136,7 @@
  "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
  "you are interested in as a list of numbers");
  keys.add("optional","MOLTYPE","the type of molecule that is under study. This is used to define the backbone atoms");
- keys.addFlag("EMST", false, "Define atoms sequence in entities using an Euclidean minimum spanning tree");
+ keys.addFlag("EMST", false, "only for backward compatibility, as of PLUMED 2.10 this is the default when using MOLINFO with WHOLE");
  keys.addFlag("ADDREFERENCE", false, "Define the reference position of the first atom of each entity using a PDB file");
 }
 
@@ -142,7 +149,9 @@
  std::vector<std::vector<AtomNumber> > groups;
  std::vector<std::vector<AtomNumber> > roots;
  // parse optional flags
- parseFlag("EMST", doemst);
+ bool doemst_tmp;
+ parseFlag("EMST", doemst_tmp);
+ if(doemst_tmp) log << "EMST option is not needed any more as of PLUMED 2.10\n";
  parseFlag("ADDREFERENCE", addref);
 
  // create groups from ENTITY
@@ -174,8 +183,12 @@
  if(groups.size()==0) error("no atoms found for WHOLEMOLECULES!");
 
  // if using PDBs reorder atoms in groups based on proximity in PDB file
+ auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
+ if(moldat && moldat->isWhole()) doemst=true;
+
+ if(doemst_tmp && ! doemst) error("cannot enable EMST if MOLINFO is not WHOLE");
+
  if(doemst) {
- auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
  if( !moldat ) error("MOLINFO is required to use EMST");
  // initialize tree
  Tree tree = Tree(moldat);
@@ -197,7 +210,7 @@
 
  // adding reference if needed
  if(addref) {
  auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
  if( !moldat ) error("MOLINFO is required to use ADDREFERENCE");
  for(unsigned i=0; i<groups.size(); ++i) {
  // add reference position of first atom in entity

src/sasa/sasa_HASEL.cpp

@@ -87,7 +87,7 @@ A Python script for the computation of free energy of transfer values to describ
 
 If the DELTAGFILE is not provided, the program computes the free energy of transfer values as if they had to take into account the effect of temperature according to approaches 2 or 3 in the paper \cite Arsiccio-T-SASA-2021. Please read and cite this paper if using the transfer model for computing the effect of temperature in implicit solvent simulations. For this purpose, the keyword APPROACH should be added, and set to either 2 or 3.
 
-The SASA usually makes sense when atoms used for the calculation are all part of the same molecule. When running with periodic boundary conditions, the atoms should be in the proper periodic image. This is done automatically since PLUMED 2.2, by considering the ordered list of atoms and rebuilding the broken entities using a procedure that is equivalent to that done in \ref WHOLEMOLECULES. Notice that rebuilding is local to this action. This is different from \ref WHOLEMOLECULES which actually modifies the coordinates stored in PLUMED.
+The SASA usually makes sense when atoms used for the calculation are all part of the same molecule. When running with periodic boundary conditions, the atoms should be in the proper periodic image. This is done automatically since PLUMED 2.2, by rebuilding the broken entities using a procedure that is equivalent to that done in \ref WHOLEMOLECULES. Notice that rebuilding is local to this action. This is different from \ref WHOLEMOLECULES which actually modifies the coordinates stored in PLUMED.
 
 In case you want to recover the old behavior you should use the NOPBC flag.
 In that case you need to take care that atoms are in the correct periodic image.

src/sasa/sasa_LCPO.cpp

@@ -85,7 +85,7 @@ A Python script for the computation of free energy of transfer values to describ
 
 If the DELTAGFILE is not provided, the program computes the free energy of transfer values as if they had to take into account the effect of temperature according to approaches 2 or 3 in the paper \cite Arsiccio-T-SASA-2021. Please read and cite this paper if using the transfer model for computing the effect of temperature in implicit solvent simulations. For this purpose, the keyword APPROACH should be added, and set to either 2 or 3.
 
-The SASA usually makes sense when atoms used for the calculation are all part of the same molecule. When running with periodic boundary conditions, the atoms should be in the proper periodic image. This is done automatically since PLUMED 2.2, by considering the ordered list of atoms and rebuilding the broken entities using a procedure that is equivalent to that done in \ref WHOLEMOLECULES. Notice that rebuilding is local to this action. This is different from \ref WHOLEMOLECULES which actually modifies the coordinates stored in PLUMED.
+The SASA usually makes sense when atoms used for the calculation are all part of the same molecule. When running with periodic boundary conditions, the atoms should be in the proper periodic image. This is done automatically since PLUMED 2.2, by rebuilding the broken entities using a procedure that is equivalent to that done in \ref WHOLEMOLECULES. Notice that rebuilding is local to this action. This is different from \ref WHOLEMOLECULES which actually modifies the coordinates stored in PLUMED.
 
 In case you want to recover the old behavior you should use the NOPBC flag.
 In that case you need to take care that atoms are in the correct periodic image.