Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-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 "MultiColvarBase.h"
      23             : #include "AtomValuePack.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "vesselbase/LessThan.h"
      26             : #include "vesselbase/Between.h"
      27             : 
      28             : #include <string>
      29             : #include <cmath>
      30             : 
      31             : namespace PLMD {
      32             : namespace multicolvar {
      33             : 
      34             : //+PLUMEDOC MCOLVAR DISTANCES
      35             : /*
      36             : Calculate the distances between one or many pairs of atoms.  You can then calculate functions of the distribution of
      37             :  distances such as the minimum, the number less than a certain quantity and so on.
      38             : 
      39             : \par Examples
      40             : 
      41             : The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2 and to
      42             : print the minimum for these two distances.
      43             : \plumedfile
      44             : d1: DISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1}
      45             : PRINT ARG=d1.min
      46             : \endplumedfile
      47             : (See also \ref PRINT).
      48             : 
      49             : The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2
      50             : and then to calculate the number of these distances that are less than 0.1 nm.  The number of distances
      51             : less than 0.1nm is then printed to a file.
      52             : \plumedfile
      53             : d1: DISTANCES ATOMS1=3,5 ATOMS2=1,2 LESS_THAN={RATIONAL R_0=0.1}
      54             : PRINT ARG=d1.lessthan
      55             : \endplumedfile
      56             : (See also \ref PRINT \ref switchingfunction).
      57             : 
      58             : The following input tells plumed to calculate all the distances between atoms 1, 2 and 3 (i.e. the distances between atoms
      59             : 1 and 2, atoms 1 and 3 and atoms 2 and 3).  The average of these distances is then calculated.
      60             : \plumedfile
      61             : d1: DISTANCES GROUP=1-3 MEAN
      62             : PRINT ARG=d1.mean
      63             : \endplumedfile
      64             : (See also \ref PRINT)
      65             : 
      66             : The following input tells plumed to calculate all the distances between the atoms in GROUPA and the atoms in GROUPB.
      67             : In other words the distances between atoms 1 and 2 and the distance between atoms 1 and 3.  The number of distances
      68             : more than 0.1 is then printed to a file.
      69             : \plumedfile
      70             : d1: DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
      71             : PRINT ARG=d1.morethan
      72             : \endplumedfile
      73             : (See also \ref PRINT \ref switchingfunction)
      74             : 
      75             : 
      76             : \par Calculating minimum distances
      77             : 
      78             : To calculate and print the minimum distance between two groups of atoms you use the following commands
      79             : 
      80             : \plumedfile
      81             : d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MIN={BETA=500.}
      82             : PRINT ARG=d1.min FILE=colvar STRIDE=10
      83             : \endplumedfile
      84             : (see \ref DISTANCES and \ref PRINT)
      85             : 
      86             : In order to ensure that the minimum value has continuous derivatives we use the following function:
      87             : 
      88             : \f[
      89             : s = \frac{\beta}{ \log \sum_i \exp\left( \frac{\beta}{s_i} \right) }
      90             : \f]
      91             : 
      92             : where \f$\beta\f$ is a user specified parameter.
      93             : 
      94             : This input is used rather than a separate MINDIST colvar so that the same routine and the same input style can be
      95             : used to calculate minimum coordination numbers (see \ref COORDINATIONNUMBER), minimum
      96             : angles (see \ref ANGLES) and many other variables.
      97             : 
      98             : This new way of calculating mindist is part of plumed 2's multicolvar functionality.  These special actions
      99             : allow you to calculate multiple functions of a distribution of simple collective variables.  As an example you
     100             : can calculate the number of distances less than 1.0, the minimum distance, the number of distances more than
     101             : 2.0 and the number of distances between 1.0 and 2.0 by using the following command:
     102             : 
     103             : \plumedfile
     104             : d1: DISTANCES ...
     105             :  GROUPA=1-10 GROUPB=11-20
     106             :  LESS_THAN={RATIONAL R_0=1.0}
     107             :  MORE_THAN={RATIONAL R_0=2.0}
     108             :  BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0}
     109             :  MIN={BETA=500.}
     110             : ...
     111             : PRINT ARG=d1.lessthan,d1.morethan,d1.between,d1.min FILE=colvar STRIDE=10
     112             : \endplumedfile
     113             : (see \ref DISTANCES and \ref PRINT)
     114             : 
     115             : A calculation performed this way is fast because the expensive part of the calculation - the calculation of all the distances - is only
     116             : done once per step.  Furthermore, it can be made faster by using the TOL keyword to discard those distance that make only a small contributions
     117             : to the final values together with the NL_STRIDE keyword, which ensures that the distances that make only a small contribution to the final values aren't
     118             : calculated at every step.
     119             : 
     120             : */
     121             : //+ENDPLUMEDOC
     122             : 
     123             : 
     124             : class Distances : public MultiColvarBase {
     125             : private:
     126             : public:
     127             :   static void registerKeywords( Keywords& keys );
     128             :   explicit Distances(const ActionOptions&);
     129             : // active methods:
     130             :   double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
     131             : /// Returns the number of coordinates of the field
     132         111 :   bool isPeriodic() override { return false; }
     133             : };
     134             : 
     135       10521 : PLUMED_REGISTER_ACTION(Distances,"DISTANCES")
     136             : 
     137          52 : void Distances::registerKeywords( Keywords& keys ) {
     138          52 :   MultiColvarBase::registerKeywords( keys );
     139         156 :   keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
     140         208 :   keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN"); // keys.use("DHENERGY");
     141         208 :   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
     142         104 :   keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
     143             :            "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
     144             :            "calculated for each ATOM keyword you specify (all ATOM keywords should "
     145             :            "specify the indices of two atoms).  The eventual number of quantities calculated by this "
     146             :            "action will depend on what functions of the distribution you choose to calculate.");
     147         104 :   keys.reset_style("ATOMS","atoms");
     148         104 :   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
     149         104 :   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
     150             :            "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
     151         104 :   keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
     152             :            "in GROUPB. This must be used in conjunction with GROUPA.");
     153          52 : }
     154             : 
     155          51 : Distances::Distances(const ActionOptions&ao):
     156             :   Action(ao),
     157          51 :   MultiColvarBase(ao)
     158             : {
     159             :   // Read in the atoms
     160             :   std::vector<AtomNumber> all_atoms;
     161         102 :   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
     162          86 :   if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
     163          51 :   setupMultiColvarBase( all_atoms );
     164             :   // And check everything has been read in correctly
     165          51 :   checkRead();
     166             : 
     167             :   // Now check if we can use link cells
     168          51 :   if( getNumberOfVessels()>0 ) {
     169             :     bool use_link=false; double rcut;
     170          49 :     vesselbase::LessThan* lt=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(0) );
     171          49 :     if( lt ) {
     172          16 :       use_link=true; rcut=lt->getCutoff();
     173             :     } else {
     174          33 :       vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(0) );
     175          33 :       if( bt ) { use_link=true; rcut=bt->getCutoff(); }
     176             :     }
     177             :     if( use_link ) {
     178          76 :       for(unsigned i=1; i<getNumberOfVessels(); ++i) {
     179          30 :         vesselbase::LessThan* lt2=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(i) );
     180          30 :         vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(i) );
     181          30 :         if( lt2 ) {
     182           0 :           double tcut=lt2->getCutoff();
     183           0 :           if( tcut>rcut ) rcut=tcut;
     184          30 :         } else if( bt ) {
     185          30 :           double tcut=bt->getCutoff();
     186          30 :           if( tcut>rcut ) rcut=tcut;
     187             :         } else {
     188             :           use_link=false;
     189             :         }
     190             :       }
     191             :     }
     192          49 :     if( use_link ) setLinkCellCutoff( rcut );
     193             :   }
     194          51 : }
     195             : 
     196       15275 : double Distances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
     197       15275 :   Vector distance;
     198       15275 :   distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
     199       15275 :   const double value=distance.modulo();
     200       15275 :   const double invvalue=1.0/value;
     201             : 
     202             :   // And finish the calculation
     203       15275 :   addAtomDerivatives( 1, 0,-invvalue*distance, myatoms );
     204       15275 :   addAtomDerivatives( 1, 1, invvalue*distance, myatoms );
     205       15275 :   myatoms.addBoxDerivatives( 1, -invvalue*Tensor(distance,distance) );
     206       15275 :   return value;
     207             : }
     208             : 
     209             : }
     210             : }
     211             :