Line data    Source code

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