Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "CoordinationBase.h"
      23             : #include "tools/SwitchingFunction.h"
      24             : #include "ActionRegister.h"
      25             : 
      26             : #include <string>
      27             : 
      28             : using namespace std;
      29             : 
      30             : namespace PLMD {
      31             : namespace colvar {
      32             : 
      33             : //+PLUMEDOC COLVAR COORDINATION
      34             : /*
      35             : Calculate coordination numbers.
      36             : 
      37             : This keyword can be used to calculate the number of contacts between two groups of atoms
      38             : and is defined as
      39             : \f[
      40             : \sum_{i\in A} \sum_{i\in B} s_{ij}
      41             : \f]
      42             : where \f$s_{ij}\f$ is 1 if the contact between atoms \f$i\f$ and \f$j\f$ is formed,
      43             : zero otherwise.
      44             : In actuality, \f$s_{ij}\f$ is replaced with a switching function so as to ensure that the calculated CV has continuous derivatives.
      45             : The default switching function is:
      46             : \f[
      47             : 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 }
      48             : \f]
      49             : but it can be changed using the optional SWITCH option.
      50             : 
      51             : To make your calculation faster you can use a neighbor list, which makes it that only a
      52             : relevant subset of the pairwise distance are calculated at every step.
      53             : 
      54             : If GROUPB is empty, it will sum the \f$\frac{N(N-1)}{2}\f$ pairs in GROUPA. This avoids computing
      55             : twice permuted indexes (e.g. pair (i,j) and (j,i)) thus running 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             : 
      62             : \par Examples
      63             : 
      64             : 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.
      65             : \plumedfile
      66             : COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
      67             : \endplumedfile
      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             : \plumedfile
      73             : c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
      74             : PRINT ARG=c STRIDE=10
      75             : \endplumedfile
      76             : 
      77             : Here's an example that shows what happens when providing COORDINATION with
      78             : a single group:
      79             : \plumedfile
      80             : # define some huge group:
      81             : group: GROUP ATOMS=1-1000
      82             : # Here's coordination of a group against itself:
      83             : c1: COORDINATION GROUPA=group GROUPB=group R_0=0.3
      84             : # Here's coordination within a single group:
      85             : x: COORDINATION GROUPA=group R_0=0.3
      86             : # This is just multiplying times 2 the variable x:
      87             : c2: COMBINE ARG=x COEFFICIENTS=2 PERIODIC=NO
      88             : 
      89             : # the two variables c1 and c2 should be identical, but the calculation of c2 is twice faster
      90             : # since it runs on half of the pairs.
      91             : PRINT ARG=c1,c2 STRIDE=10
      92             : \endplumedfile
      93             : 
      94             : 
      95             : 
      96             : */
      97             : //+ENDPLUMEDOC
      98             : 
      99         352 : 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             :   virtual double pairing(double distance,double&dfunc,unsigned i,unsigned j)const;
     107             : };
     108             : 
     109        7711 : PLUMED_REGISTER_ACTION(Coordination,"COORDINATION")
     110             : 
     111         180 : void Coordination::registerKeywords( Keywords& keys ) {
     112         180 :   CoordinationBase::registerKeywords(keys);
     113         900 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
     114         900 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
     115         900 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
     116         720 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
     117         720 :   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         180 : }
     121             : 
     122         179 : Coordination::Coordination(const ActionOptions&ao):
     123             :   Action(ao),
     124         182 :   CoordinationBase(ao)
     125             : {
     126             : 
     127             :   string sw,errors;
     128         358 :   parse("SWITCH",sw);
     129         179 :   if(sw.length()>0) {
     130         130 :     switchingFunction.set(sw,errors);
     131         133 :     if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
     132             :   } else {
     133          49 :     int nn=6;
     134          49 :     int mm=0;
     135          49 :     double d0=0.0;
     136          49 :     double r0=0.0;
     137          98 :     parse("R_0",r0);
     138          49 :     if(r0<=0.0) error("R_0 should be explicitly specified and positive");
     139          98 :     parse("D_0",d0);
     140          98 :     parse("NN",nn);
     141          94 :     parse("MM",mm);
     142          47 :     switchingFunction.set(nn,mm,r0,d0);
     143             :   }
     144             : 
     145         176 :   checkRead();
     146             : 
     147         355 :   log<<"  contacts are counted with cutoff "<<switchingFunction.description()<<"\n";
     148         176 : }
     149             : 
     150    15251910 : double Coordination::pairing(double distance,double&dfunc,unsigned i,unsigned j)const {
     151             :   (void) i; // avoid warnings
     152             :   (void) j; // avoid warnings
     153    15251910 :   return switchingFunction.calculateSqr(distance,dfunc);
     154             : }
     155             : 
     156             : }
     157             : 
     158        5517 : }