Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "CoordinationBase.h"
      23             : #include "tools/SwitchingFunction.h"
      24             : #include "ActionRegister.h"
      25             : 
      26             : namespace PLMD {
      27             : namespace colvar {
      28             : 
      29             : //+PLUMEDOC COLVAR COORDINATION
      30             : /*
      31             : Calculate coordination numbers.
      32             : 
      33             : This keyword can be used to calculate the number of contacts between two groups of atoms
      34             : and is defined as
      35             : \f[
      36             : \sum_{i\in A} \sum_{i\in B} s_{ij}
      37             : \f]
      38             : where \f$s_{ij}\f$ is 1 if the contact between atoms \f$i\f$ and \f$j\f$ is formed,
      39             : zero otherwise.
      40             : In actuality, \f$s_{ij}\f$ is replaced with a switching function so as to ensure that the calculated CV has continuous derivatives.
      41             : The default switching function is:
      42             : \f[
      43             : s_{ij} = \frac{ 1 - \left(\frac{{\bf r}_{ij}-d_0}{r_0}\right)^n } { 1 - \left(\frac{{\bf r}_{ij}-d_0}{r_0}\right)^m }
      44             : \f]
      45             : but it can be changed using the optional SWITCH option.
      46             : 
      47             : To make your calculation faster you can use a neighbor list, which makes it that only a
      48             : relevant subset of the pairwise distance are calculated at every step.
      49             : 
      50             : If GROUPB is empty, it will sum the \f$\frac{N(N-1)}{2}\f$ pairs in GROUPA. This avoids computing
      51             : twice permuted indexes (e.g. pair (i,j) and (j,i)) thus running at twice the speed.
      52             : 
      53             : Notice that if there are common atoms between GROUPA and GROUPB the switching function should be
      54             : equal to one. These "self contacts" are discarded by plumed (since version 2.1),
      55             : so that they actually count as "zero".
      56             : 
      57             : 
      58             : \par Examples
      59             : 
      60             : The following example instructs plumed to calculate the total coordination number of the atoms in group 1-10 with the atoms in group 20-100.  For atoms 1-10 coordination numbers are calculated that count the number of atoms from the second group that are within 0.3 nm of the central atom.  A neighbor list is used to make this calculation faster, this neighbor list is updated every 100 steps.
      61             : \plumedfile
      62             : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
      63             : \endplumedfile
      64             : 
      65             : The following is a dummy example which should compute the value 0 because the self interaction
      66             : of atom 1 is skipped. Notice that in plumed 2.0 "self interactions" were not skipped, and the
      67             : same calculation should return 1.
      68             : \plumedfile
      69             : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
      70             : PRINT ARG=c STRIDE=10
      71             : \endplumedfile
      72             : 
      73             : Here's an example that shows what happens when providing COORDINATION with
      74             : a single group:
      75             : \plumedfile
      76             : # define some huge group:
      77             : group: GROUP ATOMS=1-1000
      78             : # Here's coordination of a group against itself:
      79             : c1: COORDINATION GROUPA=group GROUPB=group R_0=0.3
      80             : # Here's coordination within a single group:
      81             : x: COORDINATION GROUPA=group R_0=0.3
      82             : # This is just multiplying times 2 the variable x:
      83             : c2: COMBINE ARG=x COEFFICIENTS=2 PERIODIC=NO
      84             : 
      85             : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
      86             : # since it runs on half of the pairs.
      87             : PRINT ARG=c1,c2 STRIDE=10
      88             : \endplumedfile
      89             : 
      90             : 
      91             : 
      92             : */
      93             : //+ENDPLUMEDOC
      94             : 
      95             : class Coordination : public CoordinationBase {
      96             :   SwitchingFunction switchingFunction;
      97             : 
      98             : public:
      99             :   explicit Coordination(const ActionOptions&);
     100             : // active methods:
     101             :   static void registerKeywords( Keywords& keys );
     102             :   double pairing(double distance,double&dfunc,unsigned i,unsigned j)const override;
     103             : };
     104             : 
     105       10784 : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
     106             : 
     107         186 : void Coordination::registerKeywords( Keywords& keys ) {
     108         186 :   CoordinationBase::registerKeywords(keys);
     109         372 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
     110         372 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
     111         372 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
     112         372 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
     113         372 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
     114             :            "The following provides information on the \\ref switchingfunction that are available. "
     115             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
     116         186 : }
     117             : 
     118         185 : Coordination::Coordination(const ActionOptions&ao):
     119             :   Action(ao),
     120         185 :   CoordinationBase(ao)
     121             : {
     122             : 
     123             :   std::string sw,errors;
     124         370 :   parse("SWITCH",sw);
     125         185 :   if(sw.length()>0) {
     126         136 :     switchingFunction.set(sw,errors);
     127         134 :     if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
     128             :   } else {
     129          49 :     int nn=6;
     130          49 :     int mm=0;
     131          49 :     double d0=0.0;
     132          49 :     double r0=0.0;
     133          49 :     parse("R_0",r0);
     134          49 :     if(r0<=0.0) error("R_0 should be explicitly specified and positive");
     135          49 :     parse("D_0",d0);
     136          49 :     parse("NN",nn);
     137          47 :     parse("MM",mm);
     138          47 :     switchingFunction.set(nn,mm,r0,d0);
     139             :   }
     140             : 
     141         180 :   checkRead();
     142             : 
     143         365 :   log<<"  contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
     144         195 : }
     145             : 
     146    15263796 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
     147             :   (void) i; // avoid warnings
     148             :   (void) j; // avoid warnings
     149    15263796 :   return switchingFunction.calculateSqr(distance,dfunc);
     150             : }
     151             : 
     152             : }
     153             : 
     154             : }