Skip to content

Commit

Permalink
Changed SPRINT shortcut so that you can use both the new syntax for t…
Browse files Browse the repository at this point in the history
…his action and the old one
  • Loading branch information
Gareth Aneurin Tribello authored and Gareth Aneurin Tribello committed Mar 24, 2024
1 parent 3baf976 commit 4ceefe1
Show file tree
Hide file tree
Showing 14 changed files with 589 additions and 456 deletions.
4 changes: 2 additions & 2 deletions regtest/adjmat/rt-graph-2/graph.md.reference
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
flowchart BT
11(["label=#64;11
12(["label=#64;12
BIASVALUE
"])
11 -- s --> s
12 -- s --> s
subgraph subc1 [c1]
subgraph subc1_mat [c1]
c1(["label=c1
Expand Down
8 changes: 4 additions & 4 deletions regtest/adjmat/rt-graph-3/graph.md.reference
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,10 @@ ff -- ff --> rr
rr(["label=rr
RESTRAINT
"])
ff -- ff --> 25
rr -- rr.bias --> 25
rr -- rr.force2 --> 25
25("label=#64;25
ff -- ff --> 28
rr -- rr.bias --> 28
rr -- rr.force2 --> 28
28("label=#64;28
PRINT
FILE=colvar
")
16 changes: 8 additions & 8 deletions regtest/adjmat/rt-graph/graph.md.reference
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,9 @@ s(["label=s
SUM
"])
end
c1 -- c1 --> 1
c1 -- c1 --> 2
linkStyle 2 stroke:red,color:red;
1("label=#64;1
2("label=#64;2
PRINT
FILE=colvar
")
Expand All @@ -37,22 +37,22 @@ c1 -- c1 --> cc
linkStyle 3 stroke:red,color:red;
ones -- ones --> cc
linkStyle 4 stroke:blue,color:blue;
cc -- cc --> 5
cc -- cc --> 6
linkStyle 5 stroke:blue,color:blue;
5("label=#64;5
6("label=#64;6
PRINT
FILE=coords
")
cc -- cc --> mtc
linkStyle 6 stroke:blue,color:blue;
mtc -- mtc --> s
linkStyle 7 stroke:blue,color:blue;
s -- s --> 10
10("label=#64;10
s -- s --> 11
11("label=#64;11
PRINT
FILE=fcolvar
")
s -- s --> 11
11(["label=#64;11
s -- s --> 12
12(["label=#64;12
BIASVALUE
"])
2 changes: 1 addition & 1 deletion regtest/clusters/rt-dfg2/dfs.ndx.reference
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
[ @38 step 0 ]
[ @40 step 0 ]
349 351 369 370 371 549 550 551 568 569 570 571 589 769
2 changes: 1 addition & 1 deletion regtest/clusters/rt-dfg2/dfs2.ndx.reference
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
[ @59 step 0 ]
[ @62 step 0 ]
86 87 88 89 90 105 106 107 108 109 110 111 125 126 127
128 129 130 131 147 149 306 307 308 326 328 1505 1683 1685 1687
1703 1704 1705 1706 1707 1723 1724 1725 1726 1727 1865 1883 1884 1885 1886
Expand Down
200 changes: 100 additions & 100 deletions regtest/dimred/rt-mds/analysis.0.embed.reference

Large diffs are not rendered by default.

200 changes: 100 additions & 100 deletions regtest/dimred/rt-mds/embed.reference

Large diffs are not rendered by default.

200 changes: 100 additions & 100 deletions regtest/dimred/rt-mds2/analysis.0.embed.reference

Large diffs are not rendered by default.

200 changes: 100 additions & 100 deletions regtest/dimred/rt-mds2/embed.reference

Large diffs are not rendered by default.

7 changes: 5 additions & 2 deletions src/adjmat/AdjacencyMatrixBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
keys.add("atoms","GROUPA","");
keys.add("atoms","GROUPB","");
keys.add("atoms-2","ATOMS","the atoms for which you would like to calculate the adjacency matrix (basically equivalent to GROUP)");
keys.add("atoms-2","ATOMS","the atoms for which you would like to calculate the adjacency matrix. This is a depracated syntax that is equivalent to GROUP. You are strongly recommened to use GROUP instead of ATOMS.");
keys.reserve("atoms","GROUPC","");
keys.addFlag("COMPONENTS",false,"also calculate the components of the vector connecting the atoms in the contact matrix");
keys.addFlag("NOPBC",false,"don't use pbc");
Expand All @@ -54,7 +54,10 @@ AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
natoms_per_list(0)
{
std::vector<unsigned> shape(2); std::vector<AtomNumber> t; parseAtomList("GROUP", t );
if( t.size()==0 ) parseAtomList("ATOMS", t);
if( t.size()==0 ) {
parseAtomList("ATOMS", t);
if( t.size()>0 ) warning("using depracated syntax for contact matrix. You are strongly recommended to use GROUP instead of ATOMS");
}

if( t.size()==0 ) {
std::vector<AtomNumber> ta; parseAtomList("GROUPA",ta);
Expand Down
6 changes: 4 additions & 2 deletions src/adjmat/ContactMatrix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
#include "tools/SwitchingFunction.h"
#include "tools/Matrix.h"

//+PLUMEDOC MATRIX CONTACT_MATRIX
//+PLUMEDOC MATRIX CONTACT_MATRIX_PROPER
/*
Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
Expand Down Expand Up @@ -71,9 +71,11 @@ class ContactMatrix : public AdjacencyMatrixBase {
explicit ContactMatrix(const ActionOptions&);
/// This does nothing
double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override;
/// Override this so we write the graph properly
std::string writeInGraph() const override { return "CONTACT_MATRIX"; }
};

PLUMED_REGISTER_ACTION(ContactMatrix,"CONTACT_MATRIX")
PLUMED_REGISTER_ACTION(ContactMatrix,"CONTACT_MATRIX_PROPER")

void ContactMatrix::registerKeywords( Keywords& keys ) {
AdjacencyMatrixBase::registerKeywords( keys );
Expand Down
138 changes: 138 additions & 0 deletions src/adjmat/ContactMatrixShortcut.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Copyright (c) 2013-2020 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 "core/ActionShortcut.h"
#include "core/ActionRegister.h"
#include "core/PlumedMain.h"
#include "core/ActionSet.h"
#include "core/Group.h"

//+PLUMEDOC MATRIX CONTACT_MATRIX
/*
Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
so called adjacency matrix. An adjacency matrix is an \f$N \times N\f$ matrix in which the \f$i\f$th, \f$j\f$th element tells you whether
or not the \f$i\f$th and \f$j\f$th atoms/molecules from a set of \f$N\f$ atoms/molecules are adjacent or not. These matrices can then be further
analyzed using a number of other algorithms as is detailed in \cite tribello-clustering.
For this action the elements of the contact matrix are calculated using:
\f[
a_{ij} = \sigma( |\mathbf{r}_{ij}| )
\f]
where \f$|\mathbf{r}_{ij}|\f$ is the magnitude of the vector connecting atoms \f$i\f$ and \f$j\f$ and where \f$\sigma\f$ is a \ref switchingfunction.
\par Examples
The input shown below calculates a \f$6 \times 6\f$ matrix whose elements are equal to one if atom \f$i\f$ and atom \f$j\f$ are within 0.3 nm
of each other and which is zero otherwise. The columns in this matrix are then summed so as to give the coordination number for each atom.
The final quantity output in the colvar file is thus the average coordination number.
\plumedfile
mat: CONTACT_MATRIX ATOMS=1-6 SWITCH={EXP D_0=0.2 R_0=0.1 D_MAX=0.66}
COLUMNSUMS MATRIX=mat MEAN LABEL=csums
PRINT ARG=csums.* FILE=colvar
\endplumedfile
*/
//+ENDPLUMEDOC

namespace PLMD {
namespace adjmat {

class ContactMatrixShortcut : public ActionShortcut {
public:
static void registerKeywords(Keywords& keys);
explicit ContactMatrixShortcut(const ActionOptions&);
};

PLUMED_REGISTER_ACTION(ContactMatrixShortcut,"CONTACT_MATRIX")

void ContactMatrixShortcut::registerKeywords(Keywords& keys) {
ActionShortcut::registerKeywords( keys );
keys.add("atoms","GROUPA","");
keys.add("atoms","GROUPB","");
keys.add("atoms-2","ATOMS","the atoms for which you would like to calculate the adjacency matrix (basically equivalent to GROUP)");
keys.add("numbered","GROUP","specifies the list of atoms that should be assumed indistinguishable");
keys.add("compulsory","NN","6","The n parameter of the switching function ");
keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
keys.add("compulsory","R_0","The r_0 parameter of the switching function");
keys.add("numbered","SWITCH","specify the switching function to use between two sets of indistinguishable atoms");
keys.add("compulsory","NL_CUTOFF","0.0","The cutoff for the neighbor list. A value of 0 means we are not using a neighbor list");
keys.add("compulsory","NL_STRIDE","1","The frequency with which we are updating the atoms in the neighbor list");
keys.addFlag("COMPONENTS",false,"also calculate the components of the vector connecting the atoms in the contact matrix");
keys.addFlag("NOPBC",false,"don't use pbc");
keys.addActionNameSuffix("_PROPER"); keys.needsAction("TRANSPOSE"); keys.needsAction("CONCATENATE");
}

ContactMatrixShortcut::ContactMatrixShortcut(const ActionOptions& ao):
Action(ao),
ActionShortcut(ao)
{
std::vector<std::string> grp_str; std::string atomsstr="";
std::vector<std::string> atomsvec; parseVector("ATOMS",atomsvec);
if( atomsvec.size()>0 ) {
for(unsigned i=0; i<atomsvec.size(); ++i) {
Group* gg = plumed.getActionSet().selectWithLabel<Group*>( atomsvec[i] );
if( gg ) grp_str.push_back( atomsvec[i] );
}
if( grp_str.size()!=atomsvec.size() ) {
grp_str.resize(0);
atomsstr = " ATOMS=" + atomsvec[0]; for(unsigned i=1; i<atomsvec.size(); ++i) atomsstr += "," + atomsvec[i];
}
} else {
std::string grp_inpt;
for(unsigned i=1;; ++i) {
if( !parseNumbered("GROUP",i,grp_inpt) ) break;
grp_str.push_back( grp_inpt );
}
}
if( grp_str.size()>9 ) error("cannot handle more than 9 groups");
if( grp_str.size()==0 ) { readInputLine( getShortcutLabel() + ": CONTACT_MATRIX_PROPER " + atomsstr + " " + convertInputLineToString() ); return; }

for(unsigned i=0; i<grp_str.size(); ++i) {
std::string sw_str, num; Tools::convert( i+1, num ); parseNumbered("SWITCH", (i+1)*10 + 1 + i, sw_str );
if( sw_str.length()==0 ) error("missing SWITCH" + num + num + " keyword");
readInputLine( getShortcutLabel() + num + num + ": CONTACT_MATRIX_PROPER GROUP=" + grp_str[i] + " SWITCH={" + sw_str + "}" );
for(unsigned j=0; j<i; ++j) {
std::string sw_str2, jnum; Tools::convert( j+1, jnum ); parseNumbered("SWITCH", (j+1)*10 + 1 + i, sw_str2);
if( sw_str2.length()==0 ) error("missing SWITCH" + jnum + num + " keyword");
readInputLine( getShortcutLabel() + jnum + num + ": CONTACT_MATRIX_PROPER GROUPA=" + grp_str[j] + " GROUPB=" + grp_str[i] + " SWITCH={" + sw_str2 +"}");
readInputLine( getShortcutLabel() + num + jnum + ": TRANSPOSE ARG=" + getShortcutLabel() + jnum + num );
}
}
std::string join_matrices = getShortcutLabel() + ": CONCATENATE";
for(unsigned i=0; i<grp_str.size(); ++i) {
std::string inum; Tools::convert(i+1,inum);
for(unsigned j=0; j<grp_str.size(); ++j) {
std::string jnum; Tools::convert(j+1,jnum);
if( i>j ) join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum + jnum;
else join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum + jnum;
}
}
readInputLine( join_matrices );
}

}
}
58 changes: 24 additions & 34 deletions src/adjmat/Sprint.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,9 @@
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "core/ActionShortcut.h"
#include "core/ActionRegister.h"
#include "core/ActionWithValue.h"
#include "core/PlumedMain.h"
#include "core/ActionSet.h"

//+PLUMEDOC MATRIXF SPRINT
/*
Expand Down Expand Up @@ -88,56 +91,39 @@ PLUMED_REGISTER_ACTION(Sprint,"SPRINT")

void Sprint::registerKeywords(Keywords& keys) {
ActionShortcut::registerKeywords( keys );
keys.add("optional","MATRIX","the matrix that you would like to perform SPRINT on");
keys.add("numbered","GROUP","specifies the list of atoms that should be assumed indistinguishable");
keys.add("numbered","SWITCH","specify the switching function to use between two sets of indistinguishable atoms");
keys.needsAction("CONTACT_MATRIX"); keys.needsAction("TRANSPOSE"); keys.needsAction("CONCATENATE");
keys.needsAction("DIAGONALIZE"); keys.needsAction("CUSTOM");
keys.needsAction("SELECT_COMPONENTS"); keys.needsAction("SORT");
keys.needsAction("CONTACT_MATRIX"); keys.needsAction("DIAGONALIZE"); keys.needsAction("CUSTOM");
keys.needsAction("SELECT_COMPONENTS"); keys.needsAction("SORT"); keys.needsAction("COMBINE");
keys.addOutputComponent("coord","default","the sprint coordinates");
}

Sprint::Sprint(const ActionOptions& ao):
Action(ao),
ActionShortcut(ao)
{
std::vector<std::string> grp_str; std::string grp_inpt;
for(unsigned i=1;; ++i) {
if( !parseNumbered("GROUP",i,grp_inpt) ) break;
grp_str.push_back( grp_inpt );
}
if( grp_str.size()>9 ) error("cannot handle more than 9 groups");
std::string matinp; parse("MATRIX",matinp);
if( matinp.length()==0 ) {
readInputLine( getShortcutLabel() + "_jmat: CONTACT_MATRIX " + convertInputLineToString() );
matinp = getShortcutLabel() + "_jmat";
}
std::vector<unsigned> nin_group; unsigned ntot_atoms=0;
for(unsigned i=0; i<grp_str.size(); ++i) {
std::string sw_str, num; Tools::convert( i+1, num ); parseNumbered("SWITCH", (i+1)*10 + 1 + i, sw_str );
if( sw_str.length()==0 ) error("missing SWITCH" + num + num + " keyword");
readInputLine( getShortcutLabel() + "_mat" + num + num + ": CONTACT_MATRIX GROUP=" + grp_str[i] + " SWITCH={" + sw_str + "}" );
// Get number of atoms in each group
std::vector<std::string> words=Tools::getWords(grp_str[i],"\t\n ,"); Tools::interpretRanges(words);
nin_group.push_back( words.size() ); ntot_atoms += words.size();
for(unsigned j=0; j<i; ++j) {
std::string sw_str2, jnum; Tools::convert( j+1, jnum ); parseNumbered("SWITCH", (j+1)*10 + 1 + i, sw_str2);
if( sw_str2.length()==0 ) error("missing SWITCH" + jnum + num + " keyword");
readInputLine( getShortcutLabel() + "_mat" + jnum + num + ": CONTACT_MATRIX GROUPA=" + grp_str[j] + " GROUPB=" + grp_str[i] + " SWITCH={" + sw_str2 +"}");
readInputLine( getShortcutLabel() + "_mat" + num + jnum + ": TRANSPOSE ARG=" + getShortcutLabel() + "_mat" + jnum + num );
}
}
std::string join_matrices = getShortcutLabel() + "_jmat: CONCATENATE";
for(unsigned i=0; i<grp_str.size(); ++i) {
std::string inum; Tools::convert(i+1,inum);
for(unsigned j=0; j<grp_str.size(); ++j) {
std::string jnum; Tools::convert(j+1,jnum);
if( i>j ) join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + "_mat" + inum + jnum;
else join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + "_mat" + inum + jnum;
}
for(unsigned i=1;; ++i) {
std::string inum; Tools::convert( i, inum );
ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( matinp + inum + inum );
if( !av ) break ;
unsigned natoms = (av->copyOutput(0))->getShape()[0]; nin_group.push_back( natoms ); ntot_atoms += natoms;
}
readInputLine( join_matrices );

// Diagonalization
readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + getShortcutLabel() + "_jmat VECTORS=1");
readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + matinp + " VECTORS=1");
// Compute sprint coordinates as product of eigenvalue and eigenvector times square root of number of atoms in all groups
std::string str_natoms; Tools::convert( ntot_atoms, str_natoms );
readInputLine( getShortcutLabel() + "_sp: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() +
"_diag.vecs-1 FUNC=sqrt(" + str_natoms + ")*x*y PERIODIC=NO");
// Sort sprint coordinates for each group of atoms
unsigned k=0;
unsigned k=0, kk=0;
for(unsigned j=0; j<nin_group.size(); ++j) {
std::string jnum, knum; Tools::convert( j+1, jnum ); Tools::convert(k+1, knum); k++;
std::string sort_act = getShortcutLabel() + "_selection" + jnum + ": SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_sp COMPONENTS=" + knum;
Expand All @@ -146,6 +132,10 @@ Sprint::Sprint(const ActionOptions& ao):
}
readInputLine( sort_act );
readInputLine( getShortcutLabel() + jnum + ": SORT ARG=" + getShortcutLabel() + "_selection" + jnum );
for(unsigned n=0; n<nin_group[j]; ++n) {
std::string knum, nnum; Tools::convert( kk, knum ); Tools::convert( n+1, nnum ); kk++;
readInputLine( getShortcutLabel() + "_coord-" + knum + ": COMBINE ARG=" + getShortcutLabel() + jnum + "." + nnum + " PERIODIC=NO" );
}
}
}

Expand Down
4 changes: 2 additions & 2 deletions src/valtools/Concatenate.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -108,11 +108,11 @@ Concatenate::Concatenate(const ActionOptions& ao):
std::vector<unsigned> shape(2); shape[0]=0; unsigned k=0;
row_starts.resize( arglist.size() ); col_starts.resize( arglist.size() );
for(unsigned i=0; i<nrows; ++i) {
unsigned cstart = 0, nr = 1; if( arglist[k]->getRank()==2 ) nr=arglist[k]->getShape()[1];
unsigned cstart = 0, nr = 1; if( arglist[k]->getRank()==2 ) nr=arglist[k]->getShape()[0];
for(unsigned j=0; j<ncols; ++j) {
if( arglist[k]->getRank()==0 ) {
if( nr!=1 ) error("mismatched matrix sizes");
} else if( nrows>1 && arglist[k]->getShape()[1]!=nr ) error("mismatched matrix sizes");
} else if( nrows>1 && arglist[k]->getShape()[0]!=nr ) error("mismatched matrix sizes");
row_starts[k] = shape[0]; col_starts[k] = cstart;
if( arglist[k]->getRank()==0 ) cstart += 1;
else cstart += arglist[k]->getShape()[1];
Expand Down

0 comments on commit 4ceefe1

Please sign in to comment.