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          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          12 : }
     107             : 
     108          10 : LocalAverage::LocalAverage(const ActionOptions&ao):
     109             :   Action(ao),
     110          10 :   ActionShortcut(ao) {
     111             :   std::string sp_str, specA, specB;
     112          10 :   parse("SPECIES",sp_str);
     113          10 :   parse("SPECIESA",specA);
     114          10 :   parse("SPECIESB",specB);
     115          10 :   CoordinationNumbers::expandMatrix( false, getShortcutLabel(), sp_str, specA, specB, this );
     116             :   std::map<std::string,std::string> keymap;
     117          10 :   multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
     118          10 :   if( sp_str.length()>0 ) {
     119             :     specA=specB=sp_str;
     120             :   }
     121             :   // Calculate the coordination numbers
     122          10 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
     123          10 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     124             :   std::string size;
     125          10 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     126          20 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     127          20 :   readInputLine( getShortcutLabel() + "_coord: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_ones" );
     128             : 
     129          10 :   int l=-1;
     130          10 :   std::vector<ActionShortcut*> shortcuts=plumed.getActionSet().select<ActionShortcut*>();
     131          73 :   for(unsigned i=0; i<shortcuts.size(); ++i) {
     132          72 :     if( specA==shortcuts[i]->getShortcutLabel() ) {
     133          19 :       std::string sname = shortcuts[i]->getName();
     134          70 :       if( sname=="Q1" || sname=="Q3" || sname=="Q4" || sname=="Q6" ) {
     135          18 :         Tools::convert( sname.substr(1), l );
     136             :         break;
     137             :       }
     138             :     }
     139             :   }
     140             : 
     141          10 :   if( l>0 ) {
     142           9 :     std::string vargs, svargs, sargs = "ARG=" + getShortcutLabel() + "_mat";
     143         106 :     for(int i=-l; i<=l; ++i) {
     144          97 :       std::string num = getMomentumSymbol(i);
     145         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_rmn-" + num) ) {
     146         194 :         readInputLine( specB + "_rmn-" + num + ": CUSTOM ARG=" + specB + "_sp.rm-" + num + "," + specB + "_denom FUNC=x/y PERIODIC=NO");
     147             :       }
     148         194 :       if( !plumed.getActionSet().selectWithLabel<ActionWithValue*>(specB + "_imn-" + num) ) {
     149         194 :         readInputLine( specB  + "_imn-" + num + ": CUSTOM ARG=" + specB + "_sp.im-" + num + "," + specB  + "_denom FUNC=x/y PERIODIC=NO");
     150             :       }
     151          97 :       if( i==-l ) {
     152          18 :         vargs = "ARG=" + specB + "_rmn-" + num + "," + specB + "_imn-" + num;
     153          18 :         svargs = "ARG=" + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num;
     154             :       } else {
     155         176 :         vargs += "," +  specB + "_rmn-" + num + "," + specB + "_imn-" + num;
     156         176 :         svargs += "," + getShortcutLabel() + "_prod." + specB + "_rmn-" + num + "," + getShortcutLabel() + "_prod." + specB + "_imn-" + num;
     157             :       }
     158         194 :       sargs += "," + specB + "_rmn-" + num + "," + specB  + "_imn-" + num;
     159             :     }
     160          18 :     readInputLine( getShortcutLabel() + "_vstack: VSTACK " + vargs );
     161          18 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT " + sargs );
     162          18 :     readInputLine( getShortcutLabel() + "_vpstack: VSTACK " + svargs );
     163             :     std::string twolplusone;
     164           9 :     Tools::convert( 2*(2*l+1), twolplusone );
     165          18 :     readInputLine( getShortcutLabel() + "_lones: ONES SIZE=" + twolplusone );
     166          18 :     readInputLine( getShortcutLabel() + "_unorm: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_coord," + getShortcutLabel() + "_lones" );
     167          18 :     readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "_vpstack," + getShortcutLabel() + "_vstack," + getShortcutLabel() + "_unorm FUNC=(x+y)/(1+z) PERIODIC=NO");
     168          18 :     readInputLine( getShortcutLabel() + "_av2: CUSTOM ARG=" + getShortcutLabel() + "_av FUNC=x*x PERIODIC=NO");
     169          18 :     readInputLine( getShortcutLabel() + "_2: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_av2," + getShortcutLabel() + "_lones");
     170          18 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
     171             :   } else {
     172           2 :     readInputLine( getShortcutLabel() + "_prod: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat," + sp_str + " " + convertInputLineToString() );
     173           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_prod," + sp_str + "," + getShortcutLabel() + "_coord  FUNC=(x+y)/(1+z) PERIODIC=NO");
     174             :   }
     175          20 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
     176          10 : }
     177             : 
     178          97 : std::string LocalAverage::getMomentumSymbol( const int& m ) const {
     179          97 :   if( m<0 ) {
     180             :     std::string num;
     181          44 :     Tools::convert( -1*m, num );
     182          44 :     return "n" + num;
     183          53 :   } else if( m>0 ) {
     184             :     std::string num;
     185          44 :     Tools::convert( m, num );
     186          44 :     return "p" + num;
     187             :   }
     188           9 :   return "0";
     189             : }
     190             : 
     191             : }
     192             : }