Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-2020 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/PlumedMain.h"
      25             : #include "core/ActionSet.h"
      26             : #include "core/Group.h"
      27             : #include "AdjacencyMatrixBase.h"
      28             : 
      29             : //+PLUMEDOC MATRIX CONTACT_MATRIX
      30             : /*
      31             : Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
      32             : 
      33             : A useful tool for developing complex collective variables is the notion of the
      34             : so called adjacency matrix. An adjacency matrix can be an $N \times N$ matrix in which the $i$th, $j$th element tells you whether
      35             : or not the $i$th and $j$th atoms/molecules from a set of $N$ atoms/molecules are adjacent or not. Alternatively, you can calculate
      36             : an adjacency matrix between a set of $N$ atoms and a second set of $M$ atoms.  For this type of matrix the $i$th, $j$th element tells you
      37             : whether the whether the $i$th atom in the first group and the $j$th atom in the second group are adjacent or not.  The adjacency matrix in this
      38             : case is thus $N \times M$.
      39             : 
      40             : The simplest type of adjacency matrix is a contact matrix.  The elements of a contact matrix are calculated using:
      41             : 
      42             : $$
      43             : a_{ij} = \sigma( r_{ij} )
      44             : $$
      45             : 
      46             : where $r_{ij}$ is the magnitude of the vector connecting atoms $i$ and $j$ and where $\sigma$ is a switching function. If you want to calculate a
      47             : contact matrix for one group of atoms the input would look something like this:
      48             : 
      49             : ```plumed
      50             : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
      51             : ```
      52             : 
      53             : Alternatively, if you want to calculate the contact matrix between two groups of atoms you would use an input like following:
      54             : 
      55             : ```plumed
      56             : c2: CONTACT_MATRIX GROUPA=1-7 GROUPB=8-20 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
      57             : ```
      58             : 
      59             : Once you have calculated a contact matrix you can perform various matrix operations by using the tools in the matrixtools or clusters modules.
      60             : 
      61             : ## Coordination numbers
      62             : 
      63             : If a contact matrix, $\mathbf{C}$, is multiplied from the back by a vector of ones the $i$th element of the resulting matrix tells you the number of atoms that are
      64             : within $r_c$ of atom $i$.  In other words, the coordination numbers of the atoms can be calculated from the contact matrix by doing matrix vector multiplication.
      65             : 
      66             : This realisation about the relationship between the contact map and the coordination number is heavily used in PLUMED.  For example, to calculate
      67             : and print the coordination numbers of the first 7 atoms in the system with themselves you would use an input something like this:
      68             : 
      69             : ```plumed
      70             : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
      71             : ones: ONES SIZE=7
      72             : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
      73             : PRINT ARG=cc FILE=colvar
      74             : ```
      75             : 
      76             : Implmenting the coordination number this way is useful as there are many different ways to define whether two atoms/molecules and to construct a "contact" matrix based on
      77             : the result.  For example:
      78             : 
      79             : * You could say that two molecules are connected if they are within a certain distance of each other and if they have the same orientation (see [TORSIONS_MATRIX](TORSIONS_MATRIX.md)).
      80             : * You could say that two water molecules are connected if they are hydrogen bonded to each other (see [HBOND_MATRIX](HBOND_MATRIX.md)).
      81             : * You could say that two atoms are connected if they are within a certain distance of each other and if they have similar values for a CV (see [OUTER_PRODUCT](OUTER_PRODUCT.md)).
      82             : 
      83             : When the coordination numbers is implemented in the way described above (by doing the matrix-vector multiplication) you have the flexibility to define the contact matrix that
      84             : is used in the multiplication in whatever way you choose.  In other words, this implementation of the coordination number is much more flexible. For example, suppose you want
      85             : to calculate the number of atoms that have a coordination is greater than 3.0.  You can do this with PLUMED using the following input:
      86             : 
      87             : ```plumed
      88             : # Calculate the contact matrix for the first seven atoms in the system
      89             : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
      90             : # Calculate the coordination numbers for the first seven atoms in the system
      91             : ones: ONES SIZE=7
      92             : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
      93             : # Set the ith element of the vector mtc equal to one if the coordination number of atom i is greater than 3.
      94             : mtc: MORE_THAN ARG=cc SWITCH={RATIONAL D_0=3 R_0=1}
      95             : # Calculate the number of atoms with a coordination number greater than 3.
      96             : s: SUM ARG=mtc PERIODIC=NO
      97             : PRINT ARG=s FILE=colvar
      98             : ```
      99             : 
     100             : Alternatively, consider the CV that was used to study perovskite nucleation in [this paper](https://pubs.acs.org/doi/abs/10.1021/acs.chemmater.9b04259).  The CV here measures the number of
     101             : methylamonium molecules that are attached to at least 5 other methylamoniusm molecules, at least 7 lead atoms and at least 11 ionide ions.  We can calculate something akin to this
     102             : CV and apply a restraint on the resulting quantity by using the following input file:
     103             : 
     104             : ```plumed
     105             : # Lead ions
     106             : Pb: GROUP ATOMS=1-64
     107             : # Iodide atoms
     108             : I: GROUP ATOMS=65-256
     109             : # Methylamonium "atoms" -- in the real CV these are centers of mass rather than single atoms
     110             : cn: GROUP ATOMS=257-320
     111             : 
     112             : ones64: ONES SIZE=64
     113             : # Contact matrix that determines if methylamonium molecules are within 8 A of each other
     114             : cm_cncn: CONTACT_MATRIX GROUP=cn SWITCH={RATIONAL R_0=0.8}
     115             : # Coordination number of methylamounium with methylamonium
     116             : cc_cncn: MATRIX_VECTOR_PRODUCT ARG=cm_cncn,ones64
     117             : # Vector with elements that are one if coordiantion of methylamonium with methylamonium >5
     118             : mt_cncn: MORE_THAN ARG=cc_cncn SWITCH={RATIONAL R_0=5 NN=12 MM=24}
     119             : 
     120             : # Contact matrix that determines if methylamoinium moleulcule and Pb atom are within 7.5 A of each other
     121             : cm_cnpb: CONTACT_MATRIX GROUPA=cn GROUPB=Pb SWITCH={RATIONAL R_0=0.75}
     122             : # Coordination number of methylamonium with Pb
     123             : cc_cnpb: MATRIX_VECTOR_PRODUCT ARG=cm_cnpb,ones64
     124             : # Vector with elements that are one if coordination of methylamounium with lead is >7
     125             : mt_cnpb: MORE_THAN ARG=cc_cnpb SWITCH={RATIONAL R_0=7 NN=12 MM=24}
     126             : 
     127             : ones192: ONES SIZE=192
     128             : # Contact matrix that determines if methylamoinium moleulcule and I atom are within 6.5 A of each other
     129             : cm_cnI: CONTACT_MATRIX GROUPA=cn GROUPB=I SWITCH={RATIONAL R_0=0.65}
     130             : # Coordination number of methylamonium with I
     131             : cc_cnI: MATRIX_VECTOR_PRODUCT ARG=cm_cnI,ones192
     132             : # Vector with elements that are one if coordination of methylamounium with lead is >11
     133             : mt_cnI: MORE_THAN ARG=cc_cnI SWITCH={RATIONAL R_0=11 NN=12 MM=24}
     134             : 
     135             : # Element wise product of these three input vectors.
     136             : # mm[i]==1 if coordination number of corrsponding methylamounium with methylamonium is >5
     137             : # and if coordination of methylamounium with Pb is >7 and if coordination of methylamounium with I > 11
     138             : mm: CUSTOM ARG=mt_cncn,mt_cnpb,mt_cnI FUNC=x*y*z PERIODIC=NO
     139             : 
     140             : # Sum of coordination numbers and thus equal to number of methylamoniums with desired coordination numbers
     141             : ff: SUM ARG=mm PERIODIC=NO
     142             : 
     143             : rr: RESTRAINT ARG=ff AT=62 KAPPA=10
     144             : ```
     145             : 
     146             : ## COMPONENTS flag
     147             : 
     148             : If you add the flag COMPONENTS to the input as shown below:
     149             : 
     150             : ```plumed
     151             : c4: CONTACT_MATRIX GROUP=1-7 COMPONENTS SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
     152             : ```
     153             : 
     154             : then four matrices with the labels `c4.w`, `c4.x`, `c4.y` and `c4.z` are output by the action. The matrix with the label `c4.w` is the adjacency matrix
     155             : that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices `c4.x`, `c4.y` and `c4.z` contain the $x$, $y$ and $z$
     156             : components of the vector connecting atoms $i$ and $j$. Importantly, however, the components of these vectors are only stored in `c4.x`, `c4.y` and `c4.z`
     157             : if the elements of `c4.w` are non-zero.
     158             : 
     159             : ## Optimisation details
     160             : 
     161             : Adjacency matrices are sparse.  Each atom is only be connected to a small number of neighbours and the vast majority of the elements of the contact matrix are thus zero.  To reduce
     162             : the amount of memory that PLUMED requires PLUMED uses sparse matrix storage.  Consequently, whenever you calculate and store a contact matrix only the elements of the matrix that are
     163             : non-zero are stored.  The same thing holds for the additional values that are created when you use the COMPONENTS flag. The components of the vectors connecting atoms are only stored
     164             : when the elements of `c4.w` are non-zero.
     165             : 
     166             : We can also use the sparsity of the adjacency matrix to make the time required to compute a contact matrix scale linearly rather than quadratically with the number of atoms. Element
     167             : $i,j$ of the contact matrix is only non-zero if two atoms are within a cutoff, $r_c$. We can determine that many pairs of atoms are further appart than $r_c$ without computing the
     168             : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists.  __To turn on these features you need to set the `D_MAX` parameter in the
     169             : switching functions.__  The value you pass to the `D_MAX` keyword is used as the cutoff in the link cell algorithm.
     170             : 
     171             : */
     172             : //+ENDPLUMEDOC
     173             : 
     174             : namespace PLMD {
     175             : namespace adjmat {
     176             : 
     177             : class ContactMatrixShortcut : public ActionShortcut {
     178             : public:
     179             :   static void registerKeywords(Keywords& keys);
     180             :   explicit ContactMatrixShortcut(const ActionOptions&);
     181             : };
     182             : 
     183             : PLUMED_REGISTER_ACTION(ContactMatrixShortcut,"CONTACT_MATRIX")
     184             : 
     185         334 : void ContactMatrixShortcut::registerKeywords(Keywords& keys) {
     186         334 :   AdjacencyMatrixBase::registerKeywords( keys );
     187         668 :   keys.remove("GROUP");
     188         334 :   keys.add("numbered","GROUP","specifies the list of atoms that should be assumed indistinguishable");
     189         334 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
     190         334 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
     191         334 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
     192         334 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
     193         334 :   keys.add("numbered","SWITCH","the input for the switching function that acts upon the distance between each pair of atoms");
     194         668 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
     195         334 :   keys.addActionNameSuffix("_PROPER");
     196         334 :   keys.needsAction("TRANSPOSE");
     197         334 :   keys.needsAction("CONCATENATE");
     198         334 : }
     199             : 
     200         208 : ContactMatrixShortcut::ContactMatrixShortcut(const ActionOptions& ao):
     201             :   Action(ao),
     202         208 :   ActionShortcut(ao) {
     203             :   std::vector<std::string> grp_str;
     204         208 :   std::string atomsstr="";
     205             :   std::vector<std::string> atomsvec;
     206         416 :   parseVector("ATOMS",atomsvec);
     207         208 :   if( atomsvec.size()>0 )  {
     208           0 :     for(unsigned i=0; i<atomsvec.size(); ++i) {
     209           0 :       Group* gg = plumed.getActionSet().selectWithLabel<Group*>( atomsvec[i] );
     210           0 :       if( gg ) {
     211           0 :         grp_str.push_back( atomsvec[i] );
     212             :       }
     213             :     }
     214           0 :     if( grp_str.size()!=atomsvec.size() ) {
     215           0 :       grp_str.resize(0);
     216           0 :       atomsstr = " ATOMS=" + atomsvec[0];
     217           0 :       for(unsigned i=1; i<atomsvec.size(); ++i) {
     218           0 :         atomsstr += "," + atomsvec[i];
     219             :       }
     220             :     }
     221             :   } else {
     222             :     std::string grp_inpt;
     223           2 :     for(unsigned i=1;; ++i) {
     224         420 :       if( !parseNumbered("GROUP",i,grp_inpt) ) {
     225             :         break;
     226             :       }
     227           2 :       grp_str.push_back( grp_inpt );
     228             :     }
     229             :   }
     230         208 :   if( grp_str.size()>9 ) {
     231           0 :     error("cannot handle more than 9 groups");
     232             :   }
     233         208 :   if( grp_str.size()==0 )  {
     234         414 :     readInputLine( getShortcutLabel() + ": CONTACT_MATRIX_PROPER " + atomsstr + " " + convertInputLineToString() );
     235             :     return;
     236             :   }
     237             : 
     238           3 :   for(unsigned i=0; i<grp_str.size(); ++i) {
     239             :     std::string sw_str, num;
     240           2 :     Tools::convert( i+1, num );
     241           4 :     parseNumbered("SWITCH", (i+1)*10 + 1 + i,  sw_str );
     242           2 :     if( sw_str.length()==0 ) {
     243           0 :       error("missing SWITCH" + num + num + " keyword");
     244             :     }
     245           4 :     readInputLine( getShortcutLabel() + num +  num + ": CONTACT_MATRIX_PROPER GROUP=" + grp_str[i] + " SWITCH={" + sw_str + "}" );
     246           3 :     for(unsigned j=0; j<i; ++j) {
     247             :       std::string sw_str2, jnum;
     248           1 :       Tools::convert( j+1, jnum );
     249           2 :       parseNumbered("SWITCH", (j+1)*10 + 1 + i, sw_str2);
     250           1 :       if( sw_str2.length()==0 ) {
     251           0 :         error("missing SWITCH" + jnum + num + " keyword");
     252             :       }
     253           2 :       readInputLine( getShortcutLabel() + jnum + num + ": CONTACT_MATRIX_PROPER GROUPA=" + grp_str[j] + " GROUPB=" + grp_str[i] + " SWITCH={" + sw_str2 +"}");
     254           2 :       readInputLine( getShortcutLabel() + num +  jnum + ": TRANSPOSE ARG=" + getShortcutLabel() + jnum + num );
     255             :     }
     256             :   }
     257           1 :   std::string join_matrices = getShortcutLabel() + ": CONCATENATE";
     258           3 :   for(unsigned i=0; i<grp_str.size(); ++i) {
     259             :     std::string inum;
     260           2 :     Tools::convert(i+1,inum);
     261           6 :     for(unsigned j=0; j<grp_str.size(); ++j) {
     262             :       std::string jnum;
     263           4 :       Tools::convert(j+1,jnum);
     264           4 :       if( i>j ) {
     265           2 :         join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum +  jnum;
     266             :       } else {
     267           6 :         join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum +  jnum;
     268             :       }
     269             :     }
     270             :   }
     271           1 :   readInputLine( join_matrices );
     272         416 : }
     273             : 
     274             : }
     275             : }