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             : The following is a dummy example which should compute the value 0 because the self interaction
      70             : of atom 1 is skipped. Notice that in plumed 2.0 "self interactions" were not skipped, and the
      71             : same calculation should return 1.
      72             : 
      73             : ```plumed
      74             : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
      75             : PRINT ARG=c STRIDE=10
      76             : ````
      77             : 
      78             : Here's an example that shows what happens when COORDINATION is provided with
      79             : a single group:
      80             : 
      81             : ```plumed
      82             : # define some huge group:
      83             : group: GROUP ATOMS=1-1000
      84             : # Here's coordination of a group against itself:
      85             : c1: COORDINATION GROUPA=group GROUPB=group R_0=0.3
      86             : # Here's coordination within a single group:
      87             : x: COORDINATION GROUPA=group R_0=0.3
      88             : # This is just multiplying times 2 the variable x:
      89             : c2: COMBINE ARG=x COEFFICIENTS=2 PERIODIC=NO
      90             : 
      91             : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
      92             : # since it runs on half of the pairs.
      93             : PRINT ARG=c1,c2 STRIDE=10
      94             : ```
      95             : 
      96             : */
      97             : //+ENDPLUMEDOC
      98             : 
      99             : class Coordination : public CoordinationBase {
     100             :   SwitchingFunction switchingFunction;
     101             : 
     102             : public:
     103             :   explicit Coordination(const ActionOptions&);
     104             : // active methods:
     105             :   static void registerKeywords( Keywords& keys );
     106             :   double pairing(double distance,double&dfunc,unsigned i,unsigned j)const override;
     107             : };
     108             : 
     109             : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
     110             : 
     111         221 : void Coordination::registerKeywords( Keywords& keys ) {
     112         221 :   CoordinationBase::registerKeywords(keys);
     113         221 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
     114         221 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
     115         221 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
     116         221 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
     117         221 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
     118             :            "The following provides information on the \\ref switchingfunction that are available. "
     119             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
     120         442 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
     121         442 :   keys.setValueDescription("scalar","the value of the coordination");
     122         221 : }
     123             : 
     124         219 : Coordination::Coordination(const ActionOptions&ao):
     125             :   Action(ao),
     126         219 :   CoordinationBase(ao) {
     127             : 
     128             :   std::string sw,errors;
     129         438 :   parse("SWITCH",sw);
     130         219 :   if(sw.length()>0) {
     131         170 :     switchingFunction.set(sw,errors);
     132         168 :     if( errors.length()!=0 ) {
     133           1 :       error("problem reading SWITCH keyword : " + errors );
     134             :     }
     135             :   } else {
     136          49 :     int nn=6;
     137          49 :     int mm=0;
     138          49 :     double d0=0.0;
     139          49 :     double r0=0.0;
     140          49 :     parse("R_0",r0);
     141          49 :     if(r0<=0.0) {
     142           0 :       error("R_0 should be explicitly specified and positive");
     143             :     }
     144          49 :     parse("D_0",d0);
     145          49 :     parse("NN",nn);
     146          47 :     parse("MM",mm);
     147          47 :     switchingFunction.set(nn,mm,r0,d0);
     148             :   }
     149             : 
     150         214 :   checkRead();
     151             : 
     152         433 :   log<<"  contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
     153         224 : }
     154             : 
     155    15264886 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
     156             :   (void) i; // avoid warnings
     157             :   (void) j; // avoid warnings
     158    15264886 :   return switchingFunction.calculateSqr(distance,dfunc);
     159             : }
     160             : 
     161             : }
     162             : 
     163             : }