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/FunctionSetup.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
      34             : /*
      35             : Measure how similar the environment around atoms is to that found in a FCC structure.
      36             : 
      37             : This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
      38             : which we can use to measure how similar the environment around atom $i$ is to an fcc structure.
      39             : The function that is used to make this determination is as follows:
      40             : 
      41             : $$
      42             : 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}) }
      43             : $$
      44             : 
      45             : 
      46             : In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
      47             : atom $j$ and $r_{ij}$ is the magnitude of this vector.  $\sigma(r_{ij})$ is a switching function that acts on the distance between
      48             : 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
      49             : over all of the other atoms in the system ensures that we are calculating an average
      50             : of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$.  Lastly, $\alpha$
      51             : 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:
      52             : 
      53             : $$
      54             : a = \frac{ 80080}{ 2717 + 16 \alpha} \qquad \textrm{and} \qquad b = \frac{ 16(\alpha - 143) }{2717 + 16\alpha}
      55             : $$
      56             : 
      57             : This action was been used in all the articles in the bibliography. We thus wrote an explict
      58             : action to calculate it in PLUMED instead of using a shortcut as we did for [SIMPLECUBIC](SIMPLECUBIC.md) so that we could get
      59             : good computational performance.  Notice also that you can you can rotate the bond vectors before computing the
      60             : function in the above expression by using the PHI, THETA and PSI keywords as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md).
      61             : 
      62             : The following input calculates the FCCUBIC parameter for the 64 atoms in the system
      63             : and then calculates and prints the average value for this quantity.
      64             : 
      65             : ```plumed
      66             : d: FCCUBIC SPECIES=1-64 D_0=3.0 R_0=1.5 NN=6 MM=12
      67             : dm: MEAN ARG=d PERIODIC=NO
      68             : PRINT ARG=dm FILE=colv
      69             : ```
      70             : 
      71             : In the input above we use a rational [switching function](LESS_THAN.md) with the parameters above. We would recommend using SWITCH syntax
      72             : rather than the syntax above when giving the parameters for the switching function as you can then use any of the switching functions described
      73             : in the documentation for [LESS_THAN](LESS_THAN.md).  More importantly, however, using this syntax allows you to set the D_MAX parameter for the
      74             : switching function as demonstrated below:
      75             : 
      76             : ```plumed
      77             : d: FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
      78             : dm: MEAN ARG=d PERIODIC=NO
      79             : PRINT ARG=dm FILE=colv
      80             : ```
      81             : 
      82             : Setting the `D_MAX` can substantially improve PLUMED performance as it turns on the linked list algorithm that is discussed in the optimisation details part
      83             : of the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
      84             : 
      85             : ## Working with two types of atom
      86             : 
      87             : If you would like to calculate whether the atoms in GROUPB are arranged around the atoms in GROUPA as they in in an FCC structure you use an input like the one
      88             : shown below:
      89             : 
      90             : ```plumed
      91             : d: FCCUBIC SPECIESA=1-64 SPECIESB=65-200 SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
      92             : lt: MORE_THAN ARG=d SWITCH={RATIONAL R_0=0.5}
      93             : s: SUM ARG=lt PERIODIC=NO
      94             : PRINT ARG=s FILE=colv
      95             : ```
      96             : 
      97             : This input calculates how many of the 64 atoms that were input to the SPECIESA keyword have an fcc value that is greater than 0.5.
      98             : 
      99             : ## The MASK keyword
     100             : 
     101             : You can use the MASK keyword with this action in the same way that it is used with [COORDINATIONNUMBER](COORDINATIONNUMBER.md).  This keyword thus expects a vector in
     102             : input, which tells FCCUBIC that it is safe not to calculate the FCCUBIC parameter for some of the atoms.  As illustrated below, this is useful if you are using functionality
     103             : from the [volumes module](module_volumes.md) to calculate the average value of the FCCUBIC parameter for only those atoms that lie in a certain part of the simulation box.
     104             : 
     105             : ```plumed
     106             : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
     107             : center: FIXEDATOM AT=2.5,2.5,2.5
     108             : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
     109             : sphere: INSPHERE ATOMS=1-400 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
     110             : # Calculate the fccubic parameter of the atoms
     111             : cc: FCCUBIC ...
     112             :   SPECIES=1-400 MASK=sphere
     113             :   SWITCH={RATIONAL D_0=3.0 R_0=1.5 D_MAX=6.0}
     114             : ...
     115             : # Multiply fccubic parameters numbers by sphere vector
     116             : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
     117             : # Sum of coordination numbers for atoms that are in the sphere of interest
     118             : numer: SUM ARG=prod PERIODIC=NO
     119             : # Number of atoms that are in sphere of interest
     120             : denom: SUM ARG=sphere PERIODIC=NO
     121             : # Average coordination number for atoms in sphere of interest
     122             : av: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
     123             : # And print out final CV to a file
     124             : PRINT ARG=av FILE=colvar STRIDE=1
     125             : ```
     126             : 
     127             : This input calculate the average value of the FCCUBIC parameter for only those atoms that are within a spherical region that is centered on the point
     128             : $(2.5,2.5,2.5)$.
     129             : 
     130             : ## Deprecated syntax
     131             : 
     132             : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
     133             : 
     134             : */
     135             : //+ENDPLUMEDOC
     136             : 
     137             : //+PLUMEDOC MCOLVAR FCCUBIC_FUNC
     138             : /*
     139             : Measure how similar the environment around atoms is to that found in a FCC structure.
     140             : 
     141             : This is the function that is used in the [FCCUBIC](FCCUBIC.md) shortcut. You can see an example that shows how it is used
     142             : if you expand the FCCUBIC shortcut in the following example input:
     143             : 
     144             : ```plumed
     145             : d: FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN
     146             : PRINT ARG=d.* FILE=colv
     147             : ```
     148             : 
     149             : Notice that the decomposition of the function that is illustrated above allows you to do things like the calculation in the following input:
     150             : 
     151             : ```plumed
     152             : # Calculate the distances between atoms
     153             : d_mat: DISTANCE_MATRIX GROUP=1-64 CUTOFF=4.5 COMPONENTS
     154             : # Find the six nearest atom to each of the coordinates
     155             : nn: NEIGHBORS ARG=d_mat.w NLOWEST=6
     156             : # Evalulate the FCC function
     157             : d_vfunc: FCCUBIC_FUNC ARG=d_mat.x,d_mat.y,d_mat.z MASK=nn ALPHA=3.0
     158             : # Take the product of the weights with the fcc function evaluations
     159             : d_wvfunc: CUSTOM ARG=d_vfunc,nn FUNC=x*y PERIODIC=NO
     160             : # Calculate the sum of fcc cubic function values for each atom
     161             : d_ones: ONES SIZE=64
     162             : d: MATRIX_VECTOR_PRODUCT ARG=d_wvfunc,d_ones
     163             : # Calculate the number of neighbours
     164             : d_denom: MATRIX_VECTOR_PRODUCT ARG=nn,d_ones
     165             : # Calculate the average value of the fcc cubic function per bonds
     166             : d_n: CUSTOM ARG=d,d_denom FUNC=x/y PERIODIC=NO
     167             : d_mean: MEAN ARG=d_n PERIODIC=NO
     168             : PRINT ARG=d_mean FILE=colv
     169             : ```
     170             : 
     171             : In this input we evaluate the [FCCUBIC](FCCUBIC.md) function for the six nearest neighbours to each atom rather than the atoms that are within a certain cutoff.
     172             : By using the MASK keyword in the input to FCCUBIC_FUNC we ensure that the [FCCUBIC](FCCUBIC.md) function is only evaluated for the six nearest neighbors.
     173             : 
     174             : */
     175             : //+ENDPLUMEDOC
     176             : 
     177             : class Fccubic {
     178             : public:
     179             :   double alpha, a1, b1;
     180             :   static void registerKeywords( Keywords& keys );
     181             :   static void read( Fccubic& func, ActionWithArguments* action, function::FunctionOptions& funcout );
     182             :   static void calc( const Fccubic& func, bool noderiv, View<const double> args, function::FunctionOutput& funcout );
     183             :   Fccubic& operator=(const Fccubic& m) {
     184             :     alpha = m.alpha;
     185             :     a1 = m.a1;
     186             :     b1 = m.b1;
     187             :     return *this;
     188             :   }
     189             : };
     190             : 
     191             : typedef function::FunctionShortcut<Fccubic> FccubicShortcut;
     192             : PLUMED_REGISTER_ACTION(FccubicShortcut,"FCCUBIC_FUNC")
     193             : typedef function::FunctionOfMatrix<Fccubic> MatrixFccubic;
     194             : PLUMED_REGISTER_ACTION(MatrixFccubic,"FCCUBIC_FUNC_MATRIX")
     195             : 
     196          28 : void Fccubic::registerKeywords( Keywords& keys ) {
     197          28 :   keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function");
     198          56 :   keys.setValueDescription("matrix","a function that measures the similarity with an fcc environment");
     199          28 : }
     200             : 
     201           6 : void Fccubic::read( Fccubic& func, ActionWithArguments* action, function::FunctionOptions& funcout ) {
     202             :   // Scaling factors such that '1' corresponds to fcc lattice
     203             :   // and '0' corresponds to isotropic (liquid)
     204           6 :   action->parse("ALPHA",func.alpha);
     205           6 :   func.a1 = 80080. / (2717. + 16*func.alpha);
     206           6 :   func.b1 = 16.*(func.alpha-143)/(2717+16*func.alpha);
     207           6 :   action->log.printf("  setting alpha paramter equal to %f \n",func.alpha);
     208           6 : }
     209             : 
     210     2490346 : void Fccubic::calc( const Fccubic& func, bool noderiv, const View<const double> args, function::FunctionOutput& funcout ) {
     211     2490346 :   double x2 = args[0]*args[0];
     212     2490346 :   double x4 = x2*x2;
     213             : 
     214     2490346 :   double y2 = args[1]*args[1];
     215     2490346 :   double y4 = y2*y2;
     216             : 
     217     2490346 :   double z2 = args[2]*args[2];
     218     2490346 :   double z4 = z2*z2;
     219             : 
     220     2490346 :   double d2 = x2 + y2 + z2;
     221     2490346 :   if( d2 < epsilon ) {
     222             :     d2 = 1;
     223             :   }
     224     2490346 :   double r8 = pow( d2, 4 );
     225     2490346 :   double r12 = pow( d2, 6 );
     226             : 
     227     2490346 :   double tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-func.alpha*x4*y4*z4/r12;
     228             : 
     229     2490346 :   double t0 = (x2*y4+x2*z4)/r8-func.alpha*x2*y4*z4/r12;
     230     2490346 :   double t1 = (y2*x4+y2*z4)/r8-func.alpha*y2*x4*z4/r12;
     231     2490346 :   double t2 = (z2*x4+z2*y4)/r8-func.alpha*z2*x4*y4/r12;
     232     2490346 :   double t3 = (2*tmp-func.alpha*x4*y4*z4/r12)/d2;
     233             : 
     234     2490346 :   if( !noderiv ) {
     235        5057 :     funcout.derivs[0][0]=4*func.a1*args[0]*(t0-t3);
     236        5057 :     funcout.derivs[0][1]=4*func.a1*args[1]*(t1-t3);
     237        5057 :     funcout.derivs[0][2]=4*func.a1*args[2]*(t2-t3);
     238             :   }
     239             : 
     240             :   // Set the value and the derivatives
     241     2490346 :   funcout.values[0] = (func.a1*tmp+func.b1);
     242     2490346 : }
     243             : 
     244             : }
     245             : }
     246             :