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 {
     133         111 :     return false;
     134             :   }
     135             : };
     136             : 
     137       13887 : PLUMED_REGISTER_ACTION(Distances,"DISTANCES")
     138             : 
     139          55 : void Distances::registerKeywords( Keywords& keys ) {
     140          55 :   MultiColvarBase::registerKeywords( keys );
     141          55 :   keys.use("ALT_MIN");
     142          55 :   keys.use("LOWEST");
     143          55 :   keys.use("HIGHEST");
     144          55 :   keys.use("MEAN");
     145          55 :   keys.use("MIN");
     146          55 :   keys.use("MAX");
     147          55 :   keys.use("LESS_THAN"); // keys.use("DHENERGY");
     148          55 :   keys.use("MORE_THAN");
     149          55 :   keys.use("BETWEEN");
     150          55 :   keys.use("HISTOGRAM");
     151          55 :   keys.use("MOMENTS");
     152         110 :   keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
     153             :            "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
     154             :            "calculated for each ATOM keyword you specify (all ATOM keywords should "
     155             :            "specify the indices of two atoms).  The eventual number of quantities calculated by this "
     156             :            "action will depend on what functions of the distribution you choose to calculate.");
     157         110 :   keys.reset_style("ATOMS","atoms");
     158         110 :   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
     159         110 :   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
     160             :            "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
     161         110 :   keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
     162             :            "in GROUPB. This must be used in conjunction with GROUPA.");
     163          55 : }
     164             : 
     165          51 : Distances::Distances(const ActionOptions&ao):
     166             :   Action(ao),
     167          51 :   MultiColvarBase(ao) {
     168             :   // Read in the atoms
     169             :   std::vector<AtomNumber> all_atoms;
     170         102 :   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
     171          51 :   if( atom_lab.size()==0 ) {
     172          70 :     readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
     173             :   }
     174          51 :   setupMultiColvarBase( all_atoms );
     175             :   // And check everything has been read in correctly
     176          51 :   checkRead();
     177             : 
     178             :   // Now check if we can use link cells
     179          51 :   if( getNumberOfVessels()>0 ) {
     180             :     bool use_link=false;
     181             :     double rcut;
     182          49 :     vesselbase::LessThan* lt=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(0) );
     183          49 :     if( lt ) {
     184             :       use_link=true;
     185          16 :       rcut=lt->getCutoff();
     186             :     } else {
     187          33 :       vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(0) );
     188          33 :       if( bt ) {
     189             :         use_link=true;
     190          30 :         rcut=bt->getCutoff();
     191             :       }
     192             :     }
     193             :     if( use_link ) {
     194          76 :       for(unsigned i=1; i<getNumberOfVessels(); ++i) {
     195          30 :         vesselbase::LessThan* lt2=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(i) );
     196          30 :         vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(i) );
     197          30 :         if( lt2 ) {
     198           0 :           double tcut=lt2->getCutoff();
     199           0 :           if( tcut>rcut ) {
     200           0 :             rcut=tcut;
     201             :           }
     202          30 :         } else if( bt ) {
     203          30 :           double tcut=bt->getCutoff();
     204          30 :           if( tcut>rcut ) {
     205           0 :             rcut=tcut;
     206             :           }
     207             :         } else {
     208             :           use_link=false;
     209             :         }
     210             :       }
     211             :     }
     212          49 :     if( use_link ) {
     213          46 :       setLinkCellCutoff( rcut );
     214             :     }
     215             :   }
     216          51 : }
     217             : 
     218       15275 : double Distances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
     219       15275 :   Vector distance;
     220       15275 :   distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
     221       15275 :   const double value=distance.modulo();
     222       15275 :   const double invvalue=1.0/value;
     223             : 
     224             :   // And finish the calculation
     225       15275 :   addAtomDerivatives( 1, 0,-invvalue*distance, myatoms );
     226       15275 :   addAtomDerivatives( 1, 1, invvalue*distance, myatoms );
     227       15275 :   myatoms.addBoxDerivatives( 1, -invvalue*Tensor(distance,distance) );
     228       15275 :   return value;
     229             : }
     230             : 
     231             : }
     232             : }
     233             :