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 CV was introduced in this article \cite fcc-michele-1 and again in this article \cite fcc-michele-2
      56             : This CV essentially determines whether the environment around any given atom is similar to that found in
      57             : the FCC structure or not.  The function that is used to make this determination is as follows:
      58             : 
      59             : \f[
      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             : \f]
      62             : 
      63             : In this expression \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ are the \f$x\f$, \f$y\f$ and \f$z\f$ components of the vector connecting atom \f$i\f$ to
      64             : atom \f$j\f$ and \f$r_{ij}\f$ is the magnitude of this vector.  \f$\sigma(r_{ij})\f$ is a \ref switchingfunction that acts on the distance between
      65             : atom \f$i\f$ and atom \f$j\f$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
      66             : over all of the other atoms in the system ensures that we are calculating an average
      67             : of the function of \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ for the atoms in the first coordination sphere around atom \f$i\f$.  Lastly, \f$\alpha\f$
      68             : is a parameter that can be set by the user, which by default is equal to three.  The values of \f$a\f$ and \f$b\f$ are calculated from \f$\alpha\f$ using:
      69             : 
      70             : \f[
      71             : a = \frac{ 80080}{ 2717 + 16 \alpha} \qquad \textrm{and} \qquad b = \frac{ 16(\alpha - 143) }{2717 + 16\alpha}
      72             : \f]
      73             : 
      74             : This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
      75             : the average value for the atoms in your system, the number of atoms that have an \f$s_i\f$ value that is more that some target and
      76             : so on.  Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
      77             : 
      78             : \par Examples
      79             : 
      80             : The following input calculates the FCCUBIC parameter for the 64 atoms in the system
      81             : and then calculates and prints the average value for this quantity.
      82             : 
      83             : \plumedfile
      84             : FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN LABEL=d
      85             : PRINT ARG=d.* FILE=colv
      86             : \endplumedfile
      87             : 
      88             : */
      89             : //+ENDPLUMEDOC
      90             : 
      91             : 
      92          40 : class Fccubic : public function::FunctionTemplateBase {
      93             : private:
      94             :   double alpha, a1, b1;
      95             : public:
      96             :   void registerKeywords( Keywords& keys ) override;
      97             :   void read( ActionWithArguments* action ) override;
      98             :   void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
      99             : };
     100             : 
     101             : typedef function::FunctionShortcut<Fccubic> FccubicShortcut;
     102             : PLUMED_REGISTER_ACTION(FccubicShortcut,"FCCUBIC_FUNC")
     103             : typedef function::FunctionOfMatrix<Fccubic> MatrixFccubic;
     104             : PLUMED_REGISTER_ACTION(MatrixFccubic,"FCCUBIC_FUNC_MATRIX")
     105             : 
     106          28 : void Fccubic::registerKeywords( Keywords& keys ) {
     107          56 :   keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function");
     108          56 :   keys.setValueDescription("matrix","a function that measures the similarity with an fcc environment");
     109          28 : }
     110             : 
     111           6 : void Fccubic::read( ActionWithArguments* action ) {
     112             :   // Scaling factors such that '1' corresponds to fcc lattice
     113             :   // and '0' corresponds to isotropic (liquid)
     114           6 :   parse(action,"ALPHA",alpha);
     115           6 :   a1 = 80080. / (2717. + 16*alpha); b1 = 16.*(alpha-143)/(2717+16*alpha);
     116           6 :   action->log.printf("  setting alpha paramter equal to %f \n",alpha);
     117           6 : }
     118             : 
     119     2486314 : void Fccubic::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
     120     2486314 :   double x2 = args[0]*args[0];
     121     2486314 :   double x4 = x2*x2;
     122             : 
     123     2486314 :   double y2 = args[1]*args[1];
     124     2486314 :   double y4 = y2*y2;
     125             : 
     126     2486314 :   double z2 = args[2]*args[2];
     127     2486314 :   double z4 = z2*z2;
     128             : 
     129     2486314 :   double d2 = x2 + y2 + z2;
     130     2486314 :   double r8 = pow( d2, 4 );
     131     2486314 :   double r12 = pow( d2, 6 );
     132             : 
     133     2486314 :   double tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-alpha*x4*y4*z4/r12;
     134             : 
     135     2486314 :   double t0 = (x2*y4+x2*z4)/r8-alpha*x2*y4*z4/r12;
     136     2486314 :   double t1 = (y2*x4+y2*z4)/r8-alpha*y2*x4*z4/r12;
     137     2486314 :   double t2 = (z2*x4+z2*y4)/r8-alpha*z2*x4*y4/r12;
     138     2486314 :   double t3 = (2*tmp-alpha*x4*y4*z4/r12)/d2;
     139             : 
     140     2486314 :   derivatives(0,0)=4*a1*args[0]*(t0-t3);
     141     2486314 :   derivatives(0,1)=4*a1*args[1]*(t1-t3);
     142     2486314 :   derivatives(0,2)=4*a1*args[2]*(t2-t3);
     143             : 
     144             :   // Set the value and the derivatives
     145     2486314 :   vals[0] = (a1*tmp+b1);
     146     2486314 : }
     147             : 
     148             : }
     149             : }
     150             :