Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "AdjacencyMatrixBase.h"
      23             : #include "tools/SwitchingFunction.h"
      24             : #include "tools/Angle.h"
      25             : #include "core/ActionRegister.h"
      26             : 
      27             : #include <string>
      28             : #include <cmath>
      29             : 
      30             : //+PLUMEDOC MATRIX HBOND_MATRIX
      31             : /*
      32             : Adjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them.
      33             : 
      34             : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
      35             : 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
      36             : 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
      37             : analyzed using a number of other algorithms as is detailed in \cite tribello-clustering.
      38             : 
      39             : For this action the elements of the adjacency matrix are calculated using:
      40             : 
      41             : \f[
      42             : a_{ij} = \sigma_{oo}( |\mathbf{r}_{ij}| ) \sum_{k=1}^N \sigma_{oh}( |\mathbf{r}_{ik}| ) \sigma_{\theta}( \theta_{kij} )
      43             : \f]
      44             : 
      45             : This expression was derived by thinking about how to detect if there is a hydrogen bond between atoms \f$i\f$ and \f$j\f$.  The notion is that
      46             : if the hydrogen bond is present atoms \f$i\f$ and \f$j\f$ should be within a certain cutoff distance.  In addition, there should be a hydrogen
      47             : within a certain cutoff distance of atom \f$i\f$ and this hydrogen should lie on or close to the vector connecting atoms \f$i\f$ and \f$j\f$.
      48             : As such \f$\sigma_{oo}( |\mathbf{r}_{ij}| )\f$ is a \ref switchingfunction that acts on the modulus of the vector connecting atom \f$i\f$ to atom
      49             : \f$j\f$.  The sum over \f$k\f$ then runs over all the hydrogen atoms that are specified using using HYDROGEN keyword.  \f$\sigma_{oh}(|\mathbf{r}_{ik}|)\f$
      50             : is a \ref switchingfunction that acts on the modulus of the vector connecting atom \f$i\f$ to atom \f$k\f$ and \f$\sigma_{\theta}(\theta_{kij})\f$
      51             : is a \ref switchingfunction that acts on the angle between the vector connecting atoms \f$i\f$ and \f$j\f$ and the vector connecting atoms \f$i\f$ and
      52             : \f$k\f$.
      53             : 
      54             : It is important to note that hydrogen bonds, unlike regular bonds, are asymmetric. In other words, the hydrogen atom does not sit at the
      55             : mid point between the two other atoms in this three-center bond.  As a result of this adjacency matrices calculated using \ref HBOND_MATRIX are not
      56             : symmetric like those calculated by \ref CONTACT_MATRIX.  One consequence of this fact is that the quantities found by performing \ref ROWSUMS and
      57             : \ref COLUMNSUMS on a square \ref HBOND_MATRIX are not the same as they would be if you performed \ref ROWSUMS and
      58             : \ref COLUMNSUMS on a square \ref CONTACT_MATRIX.
      59             : 
      60             : \par Examples
      61             : 
      62             : The following input can be used to analyze the number of hydrogen bonds each of the oxygen atoms in a box of water participates in.  Each
      63             : water molecule can participate in a hydrogen bond in one of two ways.  It can either donate one of its hydrogen atom to the neighboring oxygen or
      64             : it can accept a bond between the hydrogen of a neighboring water molecule and its own oxygen.  The input below allows you to output information
      65             : on the number of hydrogen bonds each of the water molecules donates and accepts.  This information is output in two xyz files which each contain
      66             : five columns of data.  The first four of these columns are a label for the atom and the x, y and z position of the oxygen.  The last column is then
      67             : the number of accepted/donated hydrogen bonds.
      68             : 
      69             : \plumedfile
      70             : mat: HBOND_MATRIX ATOMS=1-192:3 HYDROGENS=2-192:3,3-192:3 SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30} ASWITCH={RATIONAL R_0=0.167pi} SUM
      71             : rsums: ROWSUMS MATRIX=mat MEAN
      72             : csums: COLUMNSUMS MATRIX=mat MEAN
      73             : DUMPMULTICOLVAR DATA=rsums FILE=donors.xyz
      74             : DUMPMULTICOLVAR DATA=csums FILE=acceptors.xyz
      75             : \endplumedfile
      76             : 
      77             : */
      78             : //+ENDPLUMEDOC
      79             : 
      80             : namespace PLMD {
      81             : namespace adjmat {
      82             : 
      83             : class HbondMatrix : public AdjacencyMatrixBase {
      84             : private:
      85             :   SwitchingFunction distanceOOSwitch;
      86             :   SwitchingFunction distanceOHSwitch;
      87             :   SwitchingFunction angleSwitch;
      88             : public:
      89             :   static void registerKeywords( Keywords& keys );
      90             :   explicit HbondMatrix(const ActionOptions&);
      91             : // active methods:
      92             :   double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override ;
      93             : };
      94             : 
      95             : PLUMED_REGISTER_ACTION(HbondMatrix,"HBOND_MATRIX")
      96             : 
      97           4 : void HbondMatrix::registerKeywords( Keywords& keys ) {
      98           4 :   AdjacencyMatrixBase::registerKeywords( keys );
      99           8 :   keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond");
     100           8 :   keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond");
     101           8 :   keys.add("atoms","HYDROGENS","The list of atoms that can form the bridge between the two interesting parts "
     102             :            "of the structure.");
     103           8 :   keys.add("numbered","SWITCH","The switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them");
     104           8 :   keys.add("numbered","HSWITCH","The switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be "
     105             :            "considered a hydrogen bond");
     106           8 :   keys.add("numbered","ASWITCH","A switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and "
     107             :            "the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond");
     108           4 : }
     109             : 
     110           2 : HbondMatrix::HbondMatrix(const ActionOptions&ao):
     111             :   Action(ao),
     112           2 :   AdjacencyMatrixBase(ao)
     113             : {
     114           4 :   std::string sfinput,errors; parse("SWITCH",sfinput);
     115           2 :   if( sfinput.length()==0 ) error("could not find SWITCH keyword");
     116           2 :   distanceOOSwitch.set(sfinput,errors);
     117           2 :   if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
     118             : 
     119           4 :   std::string hsfinput; parse("HSWITCH",hsfinput);
     120           2 :   if( hsfinput.length()==0 ) error("could not find HSWITCH keyword");
     121           2 :   distanceOHSwitch.set(hsfinput,errors);
     122           2 :   if( errors.length()!=0 ) error("problem reading HSWITCH keyword : " + errors );
     123             : 
     124           4 :   std::string asfinput; parse("ASWITCH",asfinput);
     125           2 :   if( asfinput.length()==0 ) error("could not find SWITCH keyword");
     126           2 :   angleSwitch.set(asfinput,errors);
     127           2 :   if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
     128             : 
     129             :   // Setup link cells
     130           2 :   setLinkCellCutoff( false, distanceOOSwitch.get_dmax() );
     131             : 
     132             :   // And check everything has been read in correctly
     133           2 :   checkRead();
     134           2 : }
     135             : 
     136        4038 : double HbondMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
     137        4038 :   Vector ood = pos2; double ood_l = ood.modulo2(); // acceptor - donor
     138        4038 :   if( ood_l<epsilon) return 0;
     139        4038 :   double ood_df, ood_sw=distanceOOSwitch.calculateSqr( ood_l, ood_df );
     140             : 
     141             :   double value=0;
     142      520170 :   for(unsigned i=0; i<natoms; ++i) {
     143      516132 :     Vector ohd=getPosition(i,myvals); double ohd_l=ohd.modulo2();
     144      516132 :     double ohd_df, ohd_sw=distanceOHSwitch.calculateSqr( ohd_l, ohd_df );
     145             : 
     146      516132 :     Angle a; Vector ood_adf, ohd_adf; double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
     147      516132 :     double angle_df, angle_sw=angleSwitch.calculate( angle, angle_df );
     148      516132 :     value += ood_sw*ohd_sw*angle_sw;
     149             : 
     150      516132 :     if( !doNotCalculateDerivatives() ) {
     151          36 :       addAtomDerivatives( 0, angle_sw*ohd_sw*(-ood_df)*ood + angle_sw*ood_sw*(-ohd_df)*ohd + ood_sw*ohd_sw*angle_df*angle*(-ood_adf-ohd_adf), myvals );
     152          36 :       addAtomDerivatives( 1, angle_sw*ohd_sw*(+ood_df)*ood + ood_sw*ohd_sw*angle_df*angle*ood_adf, myvals );
     153          36 :       addThirdAtomDerivatives( i, angle_sw*ood_sw*(+ohd_df)*ohd + ood_sw*ohd_sw*angle_df*angle*ohd_adf, myvals );
     154          72 :       addBoxDerivatives( angle_sw*ohd_sw*(-ood_df)*Tensor(ood,ood) + angle_sw*ood_sw*(-ohd_df)*Tensor(ohd,ohd) -
     155         108 :                          ood_sw*ohd_sw*angle_df*angle*(Tensor(ood,ood_adf)+Tensor(ohd,ohd_adf)), myvals );
     156             :     }
     157             :   }
     158        4038 :   return value;
     159             : }
     160             : 
     161             : }
     162             : }