Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2017 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 "core/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/ActionWithValue.h"
      25             : #include "multicolvar/MultiColvarShortcuts.h"
      26             : #include "core/PlumedMain.h"
      27             : #include "core/ActionSet.h"
      28             : #include "CoordinationNumbers.h"
      29             : #include <string>
      30             : #include <cmath>
      31             : 
      32             : namespace PLMD {
      33             : namespace symfunc {
      34             : 
      35             : //+PLUMEDOC MCOLVARF LOCAL_AVERAGE
      36             : /*
      37             : Calculate averages over spherical regions centered on atoms
      38             : 
      39             : As is explained in <a href="http://www.youtube.com/watch?v=iDvZmbWE5ps"> this video </a> certain multicolvars
      40             : calculate one scalar quantity or one vector for each of the atoms in the system.  For example
      41             : \ref COORDINATIONNUMBER measures the coordination number of each of the atoms in the system and \ref Q4 measures
      42             : the 4th order Steinhardt parameter for each of the atoms in the system.  These quantities provide tell us something about
      43             : the disposition of the atoms in the first coordination sphere of each of the atoms of interest.  Lechner and Dellago \cite dellago-q6
      44             : have suggested that one can probe local order in a system by taking the average value of such symmetry functions over
      45             : the atoms within a spherical cutoff of each of these atoms in the systems.  When this is done with Steinhardt parameters
      46             : they claim this gives a coordinate that is better able to distinguish solid and liquid configurations of Lennard-Jones atoms.
      47             : 
      48             : You can calculate such locally averaged quantities within plumed by using the LOCAL_AVERAGE command.  This command calculates
      49             : the following atom-centered quantities:
      50             : 
      51             : \f[
      52             : s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
      53             : \f]
      54             : 
      55             : where the \f$c_i\f$ and \f$c_j\f$ values can be for any one of the symmetry functions that can be calculated using plumed
      56             : multicolvars.  The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
      57             : atoms \f$i\f$ and \f$j\f$.  Lechner and Dellago suggest that the parameters of this function should be set so that it the function is equal to one
      58             : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
      59             : 
      60             : The \f$s_i\f$ quantities calculated using the above command can be again thought of as atom-centred symmetry functions.  They
      61             : thus operate much like multicolvars.  You can thus calculate properties of the distribution of \f$s_i\f$ values using MEAN, LESS_THAN, HISTOGRAM
      62             : and so on.  You can also probe the value of these averaged variables in regions of the box by using the command in tandem with the
      63             : \ref AROUND command.
      64             : 
      65             : \par Examples
      66             : 
      67             : This example input calculates the coordination numbers for all the atoms in the system.  These coordination numbers are then averaged over
      68             : spherical regions.  The number of averaged coordination numbers that are greater than 4 is then output to a file.
      69             : 
      70             : \plumedfile
      71             : COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=d1
      72             : LOCAL_AVERAGE ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4} LABEL=la
      73             : PRINT ARG=la.* FILE=colvar
      74             : \endplumedfile
      75             : 
      76             : This example input calculates the \f$q_4\f$ (see \ref Q4) vectors for each of the atoms in the system.  These vectors are then averaged
      77             : component by component over a spherical region.  The average value for this quantity is then outputeed to a file.  This calculates the
      78             : quantities that were used in the paper by Lechner and Dellago \cite dellago-q6
      79             : 
      80             : \plumedfile
      81             : Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
      82             : LOCAL_AVERAGE ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
      83             : PRINT ARG=la.* FILE=colvar
      84             : \endplumedfile
      85             : 
      86             : */
      87             : //+ENDPLUMEDOC
      88             : 
      89             : class LocalAverage : public ActionShortcut {
      90             : private:
      91             :   std::string getMomentumSymbol( const int& m ) const ;
      92             : public:
      93             :   static void registerKeywords( Keywords& keys );
      94             :   explicit LocalAverage(const ActionOptions&);
      95             : };
      96             : 
      97             : PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
      98             : 
      99          12 : void LocalAverage::registerKeywords( Keywords& keys ) {
     100          12 :   CoordinationNumbers::shortcutKeywords( keys );
     101          24 :   keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
     102          36 :   keys.needsAction("VSTACK"); keys.needsAction("CUSTOM"); keys.needsAction("OUTER_PRODUCT");
     103          12 : }
     104             : 
     105          10 : LocalAverage::LocalAverage(const ActionOptions&ao):
     106             :   Action(ao),
     107          10 :   ActionShortcut(ao)
     108             : {
     109          30 :   std::string sp_str, specA, specB; parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
     110          10 :   CoordinationNumbers::expandMatrix( false, getShortcutLabel(), sp_str, specA, specB, this );
     111          10 :   std::map<std::string,std::string> keymap; multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
     112          10 :   if( sp_str.length()>0 ) specA=specB=sp_str;
     113             :   // Calculate the coordination numbers
     114          10 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
     115          10 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     116          10 :   std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     117          20 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     118          20 :   readInputLine( getShortcutLabel() + "_coord: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones" );
     119             : 
     120          10 :   int l=-1; std::vector<ActionShortcut*> shortcuts=plumed.getActionSet().select<ActionShortcut*>();
     121          73 :   for(unsigned i=0; i<shortcuts.size(); ++i) {
     122          72 :     if( specA==shortcuts[i]->getShortcutLabel() ) {
     123          19 :       std::string sname = shortcuts[i]->getName();
     124          79 :       if( sname=="Q1" || sname=="Q3" || sname=="Q4" || sname=="Q6" ) { Tools::convert( sname.substr(1), l ); break; }
     125             :     }
     126             :   }
     127             : 
     128          10 :   if( l>0 ) {
     129           9 :     std::string vargs, svargs, sargs = "ARG=" + getShortcutLabel() + "_mat";
     130         106 :     for(int i=-l; i<=l; ++i) {
     131          97 :       std::string num = getMomentumSymbol(i);
     132         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_rmn-" + num) ) {
     133         194 :         readInputLine( specB + "_rmn-" + num + ": CUSTOM ARG=" + specB + "_sp.rm-" + num + "," + specB + "_denom FUNC=x/y PERIODIC=NO");
     134             :       }
     135         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_imn-" + num) ) {
     136         194 :         readInputLine( specB  + "_imn-" + num + ": CUSTOM ARG=" + specB + "_sp.im-" + num + "," + specB  + "_denom FUNC=x/y PERIODIC=NO");
     137             :       }
     138         115 :       if( i==-l ) { vargs = "ARG=" + specB + "_rmn-" + num + "," + specB + "_imn-" + num; svargs = "ARG=" + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num; }
     139         264 :       else { vargs += "," +  specB + "_rmn-" + num + "," + specB + "_imn-" + num; svargs += "," + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num; }
     140         194 :       sargs += "," + specB + "_rmn-" + num + "," + specB  + "_imn-" + num;
     141             :     }
     142          18 :     readInputLine( getShortcutLabel() + "_vstack: VSTACK " + vargs );
     143          18 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT " + sargs );
     144          18 :     readInputLine( getShortcutLabel() + "_vpstack: VSTACK " + svargs );
     145          18 :     std::string twolplusone; Tools::convert( 2*(2*l+1), twolplusone ); readInputLine( getShortcutLabel() + "_lones: ONES SIZE=" + twolplusone );
     146          18 :     readInputLine( getShortcutLabel() + "_unorm: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_coord," + getShortcutLabel() + "_lones" );
     147          18 :     readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "_vpstack," + getShortcutLabel() + "_vstack," + getShortcutLabel() + "_unorm FUNC=(x+y)/(1+z) PERIODIC=NO");
     148          18 :     readInputLine( getShortcutLabel() + "_av2: CUSTOM ARG=" + getShortcutLabel() + "_av FUNC=x*x PERIODIC=NO");
     149          18 :     readInputLine( getShortcutLabel() + "_2: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_av2," + getShortcutLabel() + "_lones");
     150          18 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     151             :   } else {
     152           2 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + sp_str + " " + convertInputLineToString() );
     153           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_prod," + sp_str + "," + getShortcutLabel() + "_coord  FUNC=(x+y)/(1+z) PERIODIC=NO");
     154             :   }
     155          20 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
     156          10 : }
     157             : 
     158          97 : std::string LocalAverage::getMomentumSymbol( const int& m ) const {
     159          97 :   if( m<0 ) {
     160          44 :     std::string num; Tools::convert( -1*m, num );
     161          44 :     return "n" + num;
     162          53 :   } else if( m>0 ) {
     163          44 :     std::string num; Tools::convert( m, num );
     164          44 :     return "p" + num;
     165             :   }
     166           9 :   return "0";
     167             : }
     168             : 
     169             : }
     170             : }