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 PLUMED actions
      40             : calculate one scalar quantity or one vector for each of the atoms in the system.  For example
      41             : [COORDINATIONNUMBER](COORDINATIONNUMBER.md) measures the coordination number of each of the atoms in the system and [Q4](Q4.md) 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.  In the paper in the bibliography Lechner and Dellago
      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             : $$
      52             : s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
      53             : $$
      54             : 
      55             : where the $c_i$ and $c_j$ values can be any vector of [symmetry functions](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html)
      56             :  that can be calculated using plumed multicolvars.  The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
      57             : atoms $i$ and $j$.  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 $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
      59             : 
      60             : To see how this works in practice consider the following example input.
      61             : 
      62             : ```plumed
      63             : d1: COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2
      64             : la: LOCAL_AVERAGE SPECIES=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4}
      65             : PRINT ARG=la.* FILE=colvar
      66             : ```
      67             : 
      68             : This input calculates the coordination numbers for all the atoms in the system.  These coordination numbers are then averaged over
      69             : spherical regions.  The number of averaged coordination numbers that are greater than 4 is then output to a file.  Furthermore, if you
      70             : expand the input above you can see how the LOCAL_AVERAGE command is a shortcut action that expands to a longer input that you should be able to
      71             : interpret.
      72             : 
      73             : What Lechner and Dellago did in their paper was a little more complicated than this first example. To reproduce what they did you would use
      74             : an input something like this:
      75             : 
      76             : ```plumed
      77             : Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
      78             : LOCAL_AVERAGE SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
      79             : PRINT ARG=la.* FILE=colvar
      80             : ```
      81             : 
      82             : This example input calculates the [Q4](Q4.md) vectors for each of the atoms in the system.  These vectors are then averaged
      83             : component by component over a spherical region.  The average value for this quantity is then outputeed to a file.  If you want
      84             : to understand more about the shortcut that is used here you can read [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html).
      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          12 :   keys.needsAction("ONES");
     102          12 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
     103          12 :   keys.needsAction("VSTACK");
     104          12 :   keys.needsAction("CUSTOM");
     105          12 :   keys.needsAction("OUTER_PRODUCT");
     106          24 :   keys.setValueDescription("vector","the values of the local averages");
     107          12 :   keys.addDOI("10.1063/1.2977970");
     108          12 : }
     109             : 
     110          10 : LocalAverage::LocalAverage(const ActionOptions&ao):
     111             :   Action(ao),
     112          10 :   ActionShortcut(ao) {
     113             :   std::string sp_str, specA, specB;
     114          10 :   parse("SPECIES",sp_str);
     115          10 :   parse("SPECIESA",specA);
     116          10 :   parse("SPECIESB",specB);
     117          10 :   CoordinationNumbers::expandMatrix( false, getShortcutLabel(), sp_str, specA, specB, this );
     118             :   std::map<std::string,std::string> keymap;
     119          10 :   multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
     120          10 :   if( sp_str.length()>0 ) {
     121             :     specA=specB=sp_str;
     122             :   }
     123             :   // Calculate the coordination numbers
     124          10 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
     125          10 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     126             :   std::string size;
     127          10 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     128          20 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     129          20 :   readInputLine( getShortcutLabel() + "_coord: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones" );
     130             : 
     131          10 :   int l=-1;
     132          10 :   std::vector<ActionShortcut*> shortcuts=plumed.getActionSet().select<ActionShortcut*>();
     133          73 :   for(unsigned i=0; i<shortcuts.size(); ++i) {
     134          72 :     if( specA==shortcuts[i]->getShortcutLabel() ) {
     135          19 :       std::string sname = shortcuts[i]->getName();
     136          70 :       if( sname=="Q1" || sname=="Q3" || sname=="Q4" || sname=="Q6" ) {
     137          18 :         Tools::convert( sname.substr(1), l );
     138             :         break;
     139             :       }
     140             :     }
     141             :   }
     142             : 
     143          10 :   if( l>0 ) {
     144           9 :     std::string vargs, svargs, sargs = "ARG=" + getShortcutLabel() + "_mat";
     145         106 :     for(int i=-l; i<=l; ++i) {
     146          97 :       std::string num = getMomentumSymbol(i);
     147         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_rmn-" + num) ) {
     148         194 :         readInputLine( specB + "_rmn-" + num + ": CUSTOM ARG=" + specB + "_sp.rm-" + num + "," + specB + "_denom FUNC=x/y PERIODIC=NO");
     149             :       }
     150         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_imn-" + num) ) {
     151         194 :         readInputLine( specB  + "_imn-" + num + ": CUSTOM ARG=" + specB + "_sp.im-" + num + "," + specB  + "_denom FUNC=x/y PERIODIC=NO");
     152             :       }
     153          97 :       if( i==-l ) {
     154          18 :         vargs = "ARG=" + specB + "_rmn-" + num + "," + specB + "_imn-" + num;
     155          18 :         svargs = "ARG=" + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num;
     156             :       } else {
     157         176 :         vargs += "," +  specB + "_rmn-" + num + "," + specB + "_imn-" + num;
     158         176 :         svargs += "," + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num;
     159             :       }
     160         194 :       sargs += "," + specB + "_rmn-" + num + "," + specB  + "_imn-" + num;
     161             :     }
     162          18 :     readInputLine( getShortcutLabel() + "_vstack: VSTACK " + vargs );
     163          18 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT " + sargs );
     164          18 :     readInputLine( getShortcutLabel() + "_vpstack: VSTACK " + svargs );
     165             :     std::string twolplusone;
     166           9 :     Tools::convert( 2*(2*l+1), twolplusone );
     167          18 :     readInputLine( getShortcutLabel() + "_lones: ONES SIZE=" + twolplusone );
     168          18 :     readInputLine( getShortcutLabel() + "_unorm: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_coord," + getShortcutLabel() + "_lones" );
     169          18 :     readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "_vpstack," + getShortcutLabel() + "_vstack," + getShortcutLabel() + "_unorm FUNC=(x+y)/(1+z) PERIODIC=NO");
     170          18 :     readInputLine( getShortcutLabel() + "_av2: CUSTOM ARG=" + getShortcutLabel() + "_av FUNC=x*x PERIODIC=NO");
     171          18 :     readInputLine( getShortcutLabel() + "_2: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_av2," + getShortcutLabel() + "_lones");
     172          18 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     173             :   } else {
     174           2 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + sp_str + " " + convertInputLineToString() );
     175           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_prod," + sp_str + "," + getShortcutLabel() + "_coord  FUNC=(x+y)/(1+z) PERIODIC=NO");
     176             :   }
     177          20 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
     178          10 : }
     179             : 
     180          97 : std::string LocalAverage::getMomentumSymbol( const int& m ) const {
     181          97 :   if( m<0 ) {
     182             :     std::string num;
     183          44 :     Tools::convert( -1*m, num );
     184          44 :     return "n" + num;
     185          53 :   } else if( m>0 ) {
     186             :     std::string num;
     187          44 :     Tools::convert( m, num );
     188          44 :     return "p" + num;
     189             :   }
     190           9 :   return "0";
     191             : }
     192             : 
     193             : }
     194             : }