Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-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 "function/FunctionTemplateBase.h"
      23             : #include "function/FunctionShortcut.h"
      24             : #include "function/FunctionOfMatrix.h"
      25             : #include "core/ActionRegister.h"
      26             : 
      27             : #include <string>
      28             : #include <cmath>
      29             : 
      30             : namespace PLMD {
      31             : namespace symfunc {
      32             : 
      33             : //+PLUMEDOC MCOLVAR FCCUBIC_FUNC
      34             : /*
      35             : Measure how similar the environment around atoms is to that found in a FCC structure.
      36             : 
      37             : \par Examples
      38             : 
      39             : */
      40             : //+ENDPLUMEDOC
      41             : 
      42             : //+PLUMEDOC MCOLVAR FCCUBIC_FUNC_MATRIX
      43             : /*
      44             : Measure how similar the environment around atoms is to that found in a FCC structure.
      45             : 
      46             : \par Examples
      47             : 
      48             : */
      49             : //+ENDPLUMEDOC
      50             : 
      51             : //+PLUMEDOC MCOLVAR FCCUBIC
      52             : /*
      53             : Measure how similar the environment around atoms is to that found in a FCC structure.
      54             : 
      55             : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
      56             : which we can use to measure how similar the environment around atom $i$ is to an fcc structure.
      57             : The function that is used to make this determination is as follows:
      58             : 
      59             : $$
      60             : s_i = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \left\{ a\left[ \frac{(x_{ij}y_{ij})^4 + (x_{ij}z_{ij})^4 + (y_{ij}z_{ij})^4}{r_{ij}^8} - \frac{\alpha (x_{ij}y_{ij}z_{ij})^4}{r_{ij}^{12}} \right] + b \right\} }{ \sum_{i \ne j} \sigma(r_{ij}) }
      61             : $$
      62             : 
      63             : 
      64             : In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
      65             : atom $j$ and $r_{ij}$ is the magnitude of this vector.  $\sigma(r_{ij})$ is a switching function that acts on the distance between
      66             : atom $i$ and atom $j$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
      67             : over all of the other atoms in the system ensures that we are calculating an average
      68             : of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$.  Lastly, $\alpha$
      69             : is a parameter that can be set by the user, which by default is equal to three.  The values of $a$ and $b$ are calculated from $\alpha$ using:
      70             : 
      71             : $$
      72             : a = \frac{ 80080}{ 2717 + 16 \alpha} \qquad \textrm{and} \qquad b = \frac{ 16(\alpha - 143) }{2717 + 16\alpha}
      73             : $$
      74             : 
      75             : This action was been used in all the articles in the bibliography. We thus wrote an explict
      76             : action to calculate it in PLUMED instead of using a shortcut as we did for [SIMPLECUBIC](SIMPLECUBIC.md) so that we could get
      77             : good computational performance.
      78             : 
      79             : The following input calculates the FCCUBIC parameter for the 64 atoms in the system
      80             : and then calculates and prints the average value for this quantity.
      81             : 
      82             : ```plumed
      83             : d: FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN
      84             : PRINT ARG=d.* FILE=colv
      85             : ```
      86             : 
      87             : Notice that you can you can rotate the bond vectors before computing the
      88             : function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
      89             : 
      90             : */
      91             : //+ENDPLUMEDOC
      92             : 
      93             : 
      94          40 : class Fccubic : public function::FunctionTemplateBase {
      95             : private:
      96             :   double alpha, a1, b1;
      97             : public:
      98             :   void registerKeywords( Keywords& keys ) override;
      99             :   void read( ActionWithArguments* action ) override;
     100             :   void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
     101             : };
     102             : 
     103             : typedef function::FunctionShortcut<Fccubic> FccubicShortcut;
     104             : PLUMED_REGISTER_ACTION(FccubicShortcut,"FCCUBIC_FUNC")
     105             : typedef function::FunctionOfMatrix<Fccubic> MatrixFccubic;
     106             : PLUMED_REGISTER_ACTION(MatrixFccubic,"FCCUBIC_FUNC_MATRIX")
     107             : 
     108          28 : void Fccubic::registerKeywords( Keywords& keys ) {
     109          28 :   keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function");
     110          56 :   keys.setValueDescription("matrix","a function that measures the similarity with an fcc environment");
     111          28 : }
     112             : 
     113           6 : void Fccubic::read( ActionWithArguments* action ) {
     114             :   // Scaling factors such that '1' corresponds to fcc lattice
     115             :   // and '0' corresponds to isotropic (liquid)
     116           6 :   parse(action,"ALPHA",alpha);
     117           6 :   a1 = 80080. / (2717. + 16*alpha);
     118           6 :   b1 = 16.*(alpha-143)/(2717+16*alpha);
     119           6 :   action->log.printf("  setting alpha paramter equal to %f \n",alpha);
     120           6 : }
     121             : 
     122     2486314 : void Fccubic::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
     123     2486314 :   double x2 = args[0]*args[0];
     124     2486314 :   double x4 = x2*x2;
     125             : 
     126     2486314 :   double y2 = args[1]*args[1];
     127     2486314 :   double y4 = y2*y2;
     128             : 
     129     2486314 :   double z2 = args[2]*args[2];
     130     2486314 :   double z4 = z2*z2;
     131             : 
     132     2486314 :   double d2 = x2 + y2 + z2;
     133     2486314 :   double r8 = pow( d2, 4 );
     134     2486314 :   double r12 = pow( d2, 6 );
     135             : 
     136     2486314 :   double tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-alpha*x4*y4*z4/r12;
     137             : 
     138     2486314 :   double t0 = (x2*y4+x2*z4)/r8-alpha*x2*y4*z4/r12;
     139     2486314 :   double t1 = (y2*x4+y2*z4)/r8-alpha*y2*x4*z4/r12;
     140     2486314 :   double t2 = (z2*x4+z2*y4)/r8-alpha*z2*x4*y4/r12;
     141     2486314 :   double t3 = (2*tmp-alpha*x4*y4*z4/r12)/d2;
     142             : 
     143     2486314 :   derivatives(0,0)=4*a1*args[0]*(t0-t3);
     144     2486314 :   derivatives(0,1)=4*a1*args[1]*(t1-t3);
     145     2486314 :   derivatives(0,2)=4*a1*args[2]*(t2-t3);
     146             : 
     147             :   // Set the value and the derivatives
     148     2486314 :   vals[0] = (a1*tmp+b1);
     149     2486314 : }
     150             : 
     151             : }
     152             : }
     153             :