Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2023 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 "ContactMatrix.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC MATRIX CONTACT_MATRIX_PROPER
      26             : /*
      27             : Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
      28             : 
      29             : As discussed [here](module_adjmat.md) a useful tool for developing complex collective variables is the notion of the
      30             : so called adjacency matrix.  An adjacency matrix is an $N \times N$ matrix in which the $i$th, $j$th element tells you whether
      31             : or not the $i$th and $j$th atoms/molecules from a set of $N$ atoms/molecules are adjacent or not.  These matrices can then be further
      32             : analyzed using a number of other algorithms as is detailed in \cite tribello-clustering.
      33             : 
      34             : For this action the elements of the contact matrix are calculated using:
      35             : 
      36             : $$
      37             :  a_{ij} = \sigma( |\mathbf{r}_{ij}| )
      38             : $$
      39             : 
      40             : where $|\mathbf{r}_{ij}|$ is the magnitude of the vector connecting atoms $i$ and $j$ and where $\sigma$ is a switching function.
      41             : 
      42             : ## Examples
      43             : 
      44             : The input shown below calculates a $6 \times 6$ matrix whose elements are equal to one if atom $i$ and atom $j$ are within 0.3 nm
      45             : 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.
      46             : The final quantity output in the colvar file is thus the average coordination number.
      47             : 
      48             : ```plumed
      49             : mat: CONTACT_MATRIX ATOMS=1-6 SWITCH={EXP D_0=0.2 R_0=0.1 D_MAX=0.66}
      50             : ones: ONES SIZE=6
      51             : csums: MATRIX_TIMES_VECTOR ARG=mat,ones
      52             : PRINT ARG=csums FILE=colvar
      53             : ```
      54             : 
      55             : */
      56             : //+ENDPLUMEDOC
      57             : 
      58             : namespace PLMD {
      59             : namespace adjmat {
      60             : 
      61             : typedef AdjacencyMatrixBase<ContactMatrix> cmap;
      62             : PLUMED_REGISTER_ACTION(cmap,"CONTACT_MATRIX_PROPER")
      63             : 
      64         788 : void ContactMatrix::registerKeywords( Keywords& keys ) {
      65        1576 :   keys.setDisplayName("CONTACT_MATRIX");
      66         788 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
      67         788 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
      68         788 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
      69         788 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
      70         788 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
      71             :            "The following provides information on the \\ref switchingfunction that are available. "
      72             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
      73        1576 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
      74         788 : }
      75             : 
      76         219 : void ContactMatrix::parseInput( AdjacencyMatrixBase<ContactMatrix>* action ) {
      77             :   std::string errors;
      78             :   std::string swinput;
      79         438 :   action->parse("SWITCH",swinput);
      80         219 :   if( swinput.length()>0 ) {
      81         192 :     switchingFunction.set( swinput, errors );
      82         192 :     if( errors.length()!=0 ) {
      83           0 :       action->error("problem reading switching function description " + errors);
      84             :     }
      85             :   } else {
      86          27 :     int nn=0;
      87          27 :     int mm=0;
      88          27 :     double r_0=-1.0;
      89          27 :     double d_0= 0.0;
      90          27 :     action->parse("NN",nn);
      91          27 :     action->parse("MM",mm);
      92          27 :     action->parse("R_0",r_0);
      93          27 :     action->parse("D_0",d_0);
      94          27 :     if( r_0<0.0 ) {
      95           0 :       action->error("you must set a value for R_0");
      96             :     }
      97          27 :     switchingFunction.set(nn,mm,r_0,d_0);
      98             :   }
      99             :   // And set the link cell cutoff
     100         219 :   action->log.printf("  switching function cutoff is %s \n",switchingFunction.description().c_str() );
     101         219 :   action->setLinkCellCutoff( true, switchingFunction.get_dmax() );
     102         219 : }
     103             : 
     104    46204029 : void ContactMatrix::calculateWeight( const ContactMatrix& data,
     105             :                                      const AdjacencyMatrixInput& input,
     106             :                                      MatrixOutput output ) {
     107             :   const double mod2 = input.pos.modulo2();
     108    46204029 :   if( mod2<epsilon ) {
     109    42058225 :     return;  // Atoms can't be bonded to themselves
     110             :   }
     111             :   double dfunc;
     112    46202477 :   output.val[0] = data.switchingFunction.calculateSqr( mod2, dfunc );
     113    46202477 :   if( output.val[0]<epsilon ) {
     114    33229841 :     output.val[0] = 0.0;
     115    33229841 :     return;
     116             :   }
     117    12972636 :   if( input.noderiv ) {
     118             :     return;
     119             :   }
     120     4145804 :   const Vector v { (-dfunc)*input.pos[0],
     121             :                    (-dfunc)*input.pos[1],
     122     4145804 :                    (-dfunc)*input.pos[2] };
     123     4145804 :   output.deriv[0] = v[0];
     124     4145804 :   output.deriv[1] = v[1];
     125     4145804 :   output.deriv[2] = v[2];
     126     4145804 :   output.deriv[3] =-v[0];
     127     4145804 :   output.deriv[4] =-v[1];
     128     4145804 :   output.deriv[5] =-v[2];
     129             : 
     130     4145804 :   output.assignOuterProduct( 6, v, input.pos);
     131             : }
     132             : 
     133             : } // namespace adjmat
     134             : } // namespace PLMD