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 "core/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             : 
      36             : $$
      37             : \sum_{i\in A} \sum_{i\in B} s_{ij}
      38             : $$
      39             : 
      40             : where $s_{ij}$ is 1 if the contact between atoms $i$ and $j$ is formed and zero otherwise.
      41             : In actuality, $s_{ij}$ is replaced with a switching function so as to ensure that the calculated CV has continuous derivatives.
      42             : 
      43             : The default switching function is:
      44             : 
      45             : $$
      46             : s_{ij} = \frac{ 1 - \left(\frac{r_{ij}-d_0}{r_0}\right)^n } { 1 - \left(\frac{r_{ij}-d_0}{r_0}\right)^m }
      47             : $$
      48             : 
      49             : but it can be changed using the optional SWITCH option.  You can find more information about the various switching functions that you can
      50             : use with this action in the documentation for [LESS_THAN](LESS_THAN.md).
      51             : 
      52             : To make your __calculation faster you can use a neighbor list__, which makes it that only a relevant subset of the pairwise distance are calculated at every step.
      53             : 
      54             : If GROUPB is empty, the coordination  number will be calculated based on the $\frac{N(N-1)}{2}$ pairs in GROUPA. This avoids computing
      55             : permuted indexes (e.g. pair (i,j) and (j,i)) twice and ensures that the calculation runs at twice the speed.
      56             : 
      57             : Notice that if there are common atoms between GROUPA and GROUPB the switching function should be
      58             : equal to one. These "self contacts" are discarded by plumed (since version 2.1),
      59             : so that they actually count as "zero".
      60             : 
      61             : ## Examples
      62             : 
      63             : 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.
      64             : 
      65             : ```plumed
      66             : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
      67             : ```
      68             : 
      69             : ```plumed
      70             : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLISTCELLS NL_CUTOFF=0.5 NL_STRIDE=80
      71             : ```
      72             : 
      73             : The following is a dummy example which should compute the value 0 because the self interaction
      74             : of atom 1 is skipped. Notice that in plumed 2.0 "self interactions" were not skipped, and the
      75             : same calculation should return 1.
      76             : 
      77             : ```plumed
      78             : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
      79             : PRINT ARG=c STRIDE=10
      80             : ````
      81             : 
      82             : Here's an example that shows what happens when COORDINATION is provided with
      83             : a single group:
      84             : 
      85             : ```plumed
      86             : # define some huge group:
      87             : group: GROUP ATOMS=1-1000
      88             : # Here's coordination of a group against itself:
      89             : c1: COORDINATION GROUPA=group GROUPB=group R_0=0.3
      90             : # Here's coordination within a single group:
      91             : x: COORDINATION GROUPA=group R_0=0.3
      92             : # This is just multiplying times 2 the variable x:
      93             : c2: COMBINE ARG=x COEFFICIENTS=2 PERIODIC=NO
      94             : 
      95             : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
      96             : # since it runs on half of the pairs.
      97             : PRINT ARG=c1,c2 STRIDE=10
      98             : ```
      99             : 
     100             : */
     101             : //+ENDPLUMEDOC
     102             : 
     103             : class Coordination : public CoordinationBase {
     104             :   SwitchingFunction switchingFunction;
     105             : 
     106             : public:
     107             :   explicit Coordination(const ActionOptions&);
     108             : // active methods:
     109             :   static void registerKeywords( Keywords& keys );
     110             :   double pairing(double distance,double&dfunc,unsigned i,unsigned j)const override;
     111             : };
     112             : 
     113             : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
     114             : 
     115         227 : void Coordination::registerKeywords( Keywords& keys ) {
     116         227 :   CoordinationBase::registerKeywords(keys);
     117         227 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
     118         227 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
     119         227 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
     120         227 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
     121         227 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
     122             :            "The following provides information on the \\ref switchingfunction that are available. "
     123             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
     124         454 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
     125         454 :   keys.setValueDescription("scalar","the value of the coordination");
     126         227 : }
     127             : 
     128         225 : Coordination::Coordination(const ActionOptions&ao):
     129             :   Action(ao),
     130         225 :   CoordinationBase(ao) {
     131             : 
     132             :   std::string sw,errors;
     133         450 :   parse("SWITCH",sw);
     134         225 :   if(sw.length()>0) {
     135         176 :     switchingFunction.set(sw,errors);
     136         174 :     if( errors.length()!=0 ) {
     137           1 :       error("problem reading SWITCH keyword : " + errors );
     138             :     }
     139             :   } else {
     140          49 :     int nn=6;
     141          49 :     int mm=0;
     142          49 :     double d0=0.0;
     143          49 :     double r0=0.0;
     144          49 :     parse("R_0",r0);
     145          49 :     if(r0<=0.0) {
     146           0 :       error("R_0 should be explicitly specified and positive");
     147             :     }
     148          49 :     parse("D_0",d0);
     149          49 :     parse("NN",nn);
     150          47 :     parse("MM",mm);
     151          47 :     switchingFunction.set(nn,mm,r0,d0);
     152             :   }
     153             : 
     154         220 :   checkRead();
     155             : 
     156         445 :   log<<"  contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
     157         235 : }
     158             : 
     159    16358386 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
     160             :   (void) i; // avoid warnings
     161             :   (void) j; // avoid warnings
     162    16358386 :   return switchingFunction.calculateSqr(distance,dfunc);
     163             : }
     164             : 
     165             : }
     166             : 
     167             : }