Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-2023 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 "CubicHarmonicBase.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/SwitchingFunction.h"
      25             : 
      26             : #include <string>
      27             : #include <cmath>
      28             : 
      29             : using namespace std;
      30             : 
      31             : namespace PLMD {
      32             : namespace crystallization {
      33             : 
      34             : //+PLUMEDOC MCOLVAR FCCUBIC
      35             : /*
      36             : Measure how similar the environment around atoms is to that found in a FCC structure.
      37             : 
      38             : This CV was introduced in this article \cite fcc-michele-1 and again in this article \cite fcc-michele-2
      39             : This CV essentially determines whether the environment around any given atom is similar to that found in
      40             : the FCC structure or not.  The function that is used to make this determination is as follows:
      41             : 
      42             : \f[
      43             : 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}) }
      44             : \f]
      45             : 
      46             : 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
      47             : 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
      48             : 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
      49             : over all of the other atoms in the system ensures that we are calculating an average
      50             : 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$
      51             : 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:
      52             : 
      53             : \f[
      54             : a = \frac{ 80080}{ 2717 + 16 \alpha} \qquad \textrm{and} \qquad b = \frac{ 16(\alpha - 143) }{2717 + 16\alpha}
      55             : \f]
      56             : 
      57             : This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
      58             : 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
      59             : so on.  Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
      60             : 
      61             : \par Examples
      62             : 
      63             : The following input calculates the FCCUBIC parameter for the 64 atoms in the system
      64             : and then calculates and prints the average value for this quantity.
      65             : 
      66             : \plumedfile
      67             : FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN LABEL=d
      68             : PRINT ARG=d.* FILE=colv
      69             : \endplumedfile
      70             : 
      71             : */
      72             : //+ENDPLUMEDOC
      73             : 
      74             : 
      75             : class Fccubic : public CubicHarmonicBase {
      76             : private:
      77             :   double alpha, a1, b1;
      78             : public:
      79             :   static void registerKeywords( Keywords& keys );
      80             :   explicit Fccubic(const ActionOptions&);
      81             :   double calculateCubicHarmonic( const Vector& distance, const double& d2, Vector& myder ) const override;
      82             : };
      83             : 
      84       10427 : PLUMED_REGISTER_ACTION(Fccubic,"FCCUBIC")
      85             : 
      86           5 : void Fccubic::registerKeywords( Keywords& keys ) {
      87           5 :   CubicHarmonicBase::registerKeywords( keys );
      88          10 :   keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function");
      89           5 : }
      90             : 
      91           4 : Fccubic::Fccubic(const ActionOptions&ao):
      92             :   Action(ao),
      93           4 :   CubicHarmonicBase(ao)
      94             : {
      95             :   // Scaling factors such that '1' corresponds to fcc lattice
      96             :   // and '0' corresponds to isotropic (liquid)
      97           4 :   parse("ALPHA",alpha);
      98           4 :   a1 = 80080. / (2717. + 16*alpha); b1 = 16.*(alpha-143)/(2717+16*alpha);
      99           4 :   log.printf("  setting alpha parameter equal to %f \n",alpha);
     100             :   // And setup the ActionWithVessel
     101           4 :   checkRead();
     102           4 : }
     103             : 
     104      320614 : double Fccubic::calculateCubicHarmonic( const Vector& distance, const double& d2, Vector& myder ) const {
     105      320614 :   double x2 = distance[0]*distance[0];
     106      320614 :   double x4 = x2*x2;
     107             : 
     108      320614 :   double y2 = distance[1]*distance[1];
     109      320614 :   double y4 = y2*y2;
     110             : 
     111      320614 :   double z2 = distance[2]*distance[2];
     112      320614 :   double z4 = z2*z2;
     113             : 
     114      320614 :   double r8 = pow( d2, 4 );
     115             :   double r12 = pow( d2, 6 );
     116             : 
     117      320614 :   double tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-alpha*x4*y4*z4/r12;
     118             : 
     119      320614 :   double t0 = (x2*y4+x2*z4)/r8-alpha*x2*y4*z4/r12;
     120      320614 :   double t1 = (y2*x4+y2*z4)/r8-alpha*y2*x4*z4/r12;
     121      320614 :   double t2 = (z2*x4+z2*y4)/r8-alpha*z2*x4*y4/r12;
     122      320614 :   double t3 = (2*tmp-alpha*x4*y4*z4/r12)/d2;
     123             : 
     124      320614 :   myder[0]=4*a1*distance[0]*(t0-t3);
     125      320614 :   myder[1]=4*a1*distance[1]*(t1-t3);
     126      320614 :   myder[2]=4*a1*distance[2]*(t2-t3);
     127             : 
     128      320614 :   return a1*tmp+b1;
     129             : }
     130             : 
     131             : }
     132             : }
     133             :