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             : 
      25             : namespace PLMD {
      26             : namespace adjmat {
      27             : 
      28             : //+PLUMEDOC MCOLVAR BRIDGE_MATRIX
      29             : /*
      30             : Calculate the number of atoms that bridge two parts of a structure
      31             : 
      32             : This adjacency matrix is used to implement the [BRIDGE](BRIDGE.md) shortcut. The action outputs a adjacency matrix
      33             : in the same way as [CONTACT_MATRIX](CONTACT_MATRIX.md).  However, the  $j,k$ element of the adjacency matrix is calculated
      34             : using:
      35             : 
      36             : $$
      37             : M_{jk} = \sum_i s_A(r_{ij})s_B(r_{ik})
      38             : $$
      39             : 
      40             : In this expression, the sum runs over all the atoms that were specified using the `BRIDGING_ATOMS` keyword, $s_A$ and
      41             : $s_B$ are switching functions, and $r_{ij}$ and $r_{ik}$ are the distances between atom $i$ and $j$ and between atoms
      42             : $i$ and $k$.  Less formally, this formula ensures that $j,k$ element of the output matrix is one if there is a bridging
      43             : atom between atom $j$ and $k$.
      44             : 
      45             : Notice that the BRIDGE_MATRIX action is also a shortcut. If we have a single group of non-bridging atoms then we can compute
      46             : the contact matrix $\mathbf{C}$ between the non-bridging and bridging atoms. The BRIDGE_MATRIX, $\mathbf{M}$, can by calculated
      47             : from $\mathbf{C}$ by multiplying $\mathbf{C}$ by its transpose.  Similarly, if we have two groups of non-bridging atoms
      48             : we can calculate the BRIDGE_MATRIX by calculating a [CONTACT_MATRIX](CONTACT_MATRIX.md) between the first set of non-bridging atoms and
      49             : a second [CONTACT_MATRIX](CONTACT_MATRIX.md) between the bridging atoms and the bridging atoms and the second set of non-bridging atoms.
      50             : The bridge matrix is the product of these two matrices.
      51             : 
      52             : In the following example input atoms 100-200 can serve as bridging atoms between the atoms in GROUPA and GROUPB and the
      53             : two switching functions $s_A$ and $s_B$ in the formula above are identical.
      54             : 
      55             : ```plumed
      56             : w1: BRIDGE_MATRIX ...
      57             :    BRIDGING_ATOMS=100-200
      58             :    GROUPA=1-10 GROUPB=11-20
      59             :    SWITCH={RATIONAL R_0=0.2}
      60             : ...
      61             : ```
      62             : 
      63             : If you use a single GROUP keyword as in the input below as a single SWITCH keyword the output matrix is symmetric.
      64             : 
      65             : ```plumed
      66             : w2: BRIDGE_MATRIX ...
      67             :    BRIDGING_ATOMS=100-200 GROUP=1-10
      68             :    SWITCH={RATIONAL R_0=0.2}
      69             : ...
      70             : ```
      71             : 
      72             : However, if the two switching functions are not identical, as in the following example, then the output matrix is __not__ symmetric
      73             : even if GROUP is used rather than GROUPA/GROUPB.
      74             : 
      75             : ```plumed
      76             : w2: BRIDGE_MATRIX ...
      77             :    BRIDGING_ATOMS=100-200 GROUP=1-10
      78             :    SWITCHA={RATIONAL R_0=0.2}
      79             :    SWITCHB={RATIONAL R_0=0.4}
      80             : ...
      81             : ```
      82             : 
      83             : Notice that in all the inputs above the $r_{ij}$ and $r_{ik}$ values that enter the formula above are calculated in a way that takes the
      84             : periodic boundary conditions into account.  If you want to ignore the periodic boundary conditions you can use the NOPBC flag as shown below.
      85             : 
      86             : ```plumed
      87             : w2: BRIDGE_MATRIX ...
      88             :    BRIDGING_ATOMS=100-200 GROUP=1-10
      89             :    SWITCH={RATIONAL R_0=0.2}
      90             :    NOPBC
      91             : ...
      92             : ```
      93             : 
      94             : */
      95             : //+ENDPLUMEDOC
      96             : 
      97             : class BridgeMatrix : public ActionShortcut {
      98             : public:
      99             :   static void registerKeywords( Keywords& keys );
     100             :   BridgeMatrix(const ActionOptions&);
     101             : };
     102             : 
     103             : PLUMED_REGISTER_ACTION(BridgeMatrix,"BRIDGE_MATRIX")
     104             : 
     105          24 : void BridgeMatrix::registerKeywords( Keywords& keys ) {
     106          24 :   ActionShortcut::registerKeywords( keys );
     107          24 :   keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
     108          24 :   keys.add("atoms","GROUPA","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB");
     109          24 :   keys.add("atoms","GROUPB","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA");
     110          24 :   keys.add("atoms","BRIDGING_ATOMS","The list of atoms that can form the bridge between the two interesting parts "
     111             :            "of the structure.");
     112          24 :   keys.add("optional","SWITCH","The parameters of the two switching functions in the above formula");
     113          48 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
     114          24 :   keys.add("optional","SWITCHA","The switching function on the distance between bridging atoms and the atoms in "
     115             :            "group A");
     116          48 :   keys.linkActionInDocs("SWITCHA","LESS_THAN");
     117          48 :   keys.setValueDescription("matrix","a matrix containing the weights for the bonds between each pair of atoms");
     118          24 :   keys.add("optional","SWITCHB","The switching function on the distance between the bridging atoms and the atoms in "
     119             :            "group B");
     120          24 :   keys.addFlag("NOPBC", false, "don't use periodic boundary conditions");
     121          48 :   keys.linkActionInDocs("SWITCHB","LESS_THAN");
     122          24 :   keys.needsAction("CONTACT_MATRIX");
     123          24 :   keys.needsAction("TRANSPOSE");
     124          24 :   keys.needsAction("MATRIX_PRODUCT");
     125          24 : }
     126             : 
     127          12 : BridgeMatrix::BridgeMatrix(const ActionOptions&ao):
     128             :   Action(ao),
     129          12 :   ActionShortcut(ao) {
     130             :   std::string b_atoms;
     131          24 :   parse("BRIDGING_ATOMS",b_atoms);
     132             :   std::string swstr, swA_str, swB_str;
     133          24 :   parse("SWITCH",swstr);
     134          12 :   if( swstr.length()==0 ) {
     135           4 :     parse("SWITCHA",swA_str);
     136           2 :     if( swA_str.length()==0 ) {
     137           0 :       error("must set either SWITCH or SWITCHA+SWITCHB");
     138             :     }
     139           4 :     parse("SWITCHB",swB_str);
     140           2 :     if( swB_str.length()==0 ) {
     141           0 :       error("found SWITCHA but no SWITCHB");
     142             :     }
     143             :   }
     144             :   bool nopbc;
     145          12 :   std::string pbcstr="";
     146          12 :   parseFlag("NOPBC",nopbc);
     147          12 :   if( nopbc ) {
     148             :     pbcstr = " NOPBC";
     149             :   }
     150             : 
     151             :   std::string grp_str;
     152          24 :   parse("GROUP",grp_str);
     153          12 :   if( grp_str.length()>0 ) {
     154           2 :     if( swstr.length()>0 ) {
     155           4 :       readInputLine( getShortcutLabel() + "_cmat: CONTACT_MATRIX GROUPA=" + grp_str + " GROUPB=" + b_atoms + " SWITCH={" + swstr + "}" + pbcstr );
     156           3 :       readInputLine( getShortcutLabel() + "_cmatT: TRANSPOSE ARG=" + getShortcutLabel() + "_cmat" );
     157             :     } else {
     158           4 :       readInputLine( getShortcutLabel() + "_cmat: CONTACT_MATRIX GROUPA=" + grp_str + " GROUPB=" + b_atoms + " SWITCH={" + swA_str + "}" + pbcstr );
     159           5 :       readInputLine( getShortcutLabel() + "_cmatT: CONTACT_MATRIX GROUPA=" + b_atoms + " GROUPB=" + grp_str + " SWITCH={" + swB_str + "}" + pbcstr );
     160             :     }
     161           8 :     readInputLine( getShortcutLabel() + ": MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_cmat," + getShortcutLabel() + "_cmatT");
     162             :   } else {
     163          10 :     if( swA_str.length()==0 ) {
     164             :       swA_str = swB_str = swstr;
     165             :     }
     166             :     std::string grpA_str, grpB_str;
     167          10 :     parse("GROUPA",grpA_str);
     168          10 :     parse("GROUPB",grpB_str);
     169          40 :     readInputLine( getShortcutLabel() + "_cmat: CONTACT_MATRIX GROUPA=" + grpA_str + " GROUPB=" + b_atoms + " SWITCH={" + swA_str + "}" + pbcstr );
     170          40 :     readInputLine( getShortcutLabel() + "_cmatT: CONTACT_MATRIX GROUPA=" + b_atoms + " GROUPB=" + grpB_str + " SWITCH={" + swB_str + "}" + pbcstr );
     171          40 :     readInputLine( getShortcutLabel() + ": MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_cmat," + getShortcutLabel() + "_cmatT");
     172             :   }
     173          12 : }
     174             : 
     175             : }
     176             : }