Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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 "Steinhardt.h"
      23             : #include "LocalSteinhardt.h"
      24             : #include "core/ActionRegister.h"
      25             : 
      26             : //+PLUMEDOC MCOLVAR Q6
      27             : /*
      28             : Calculate sixth order Steinhardt parameters.
      29             : 
      30             : The sixth order Steinhardt parameters allow us to measure the degree to which the first coordination shell
      31             : around an atom is ordered.  The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
      32             : calculated using the following formula:
      33             : 
      34             : \f[
      35             : q_{6m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{6m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
      36             : \f]
      37             : 
      38             : where \f$Y_{6m}\f$ is one of the sixth order spherical harmonics so \f$m\f$ is a number that runs from \f$-6\f$ to
      39             : \f$+6\f$.  The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
      40             : atoms \f$i\f$ and \f$j\f$.  The parameters of this function should be set so that it the function is equal to one
      41             : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
      42             : 
      43             : The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways.  The
      44             : simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
      45             : 
      46             : \f[
      47             : Q_6(i) = \sqrt{ \sum_{m=-6}^6 q_{6m}(i)^{*} q_{6m}(i) }
      48             : \f]
      49             : 
      50             : This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
      51             : the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
      52             : that is investigated.
      53             : 
      54             : Other measures of order can be taken by averaging the components of the individual \f$q_6\f$ vectors individually or by taking dot products of
      55             : the \f$q_{6}\f$ vectors on adjacent atoms.  More information on these variables can be found in the documentation for \ref LOCAL_Q6,
      56             : \ref LOCAL_AVERAGE and \ref NLINKS.
      57             : 
      58             : \par Examples
      59             : 
      60             : The following command calculates the average Q6 parameter for the 64 atoms in a box of Lennard Jones and prints this
      61             : quantity to a file called colvar:
      62             : 
      63             : \plumedfile
      64             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q6
      65             : PRINT ARG=q6.mean FILE=colvar
      66             : \endplumedfile
      67             : 
      68             : The following command calculates the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones and prints these
      69             : quantities to a file called colvar:
      70             : 
      71             : \plumedfile
      72             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q6
      73             : PRINT ARG=q6.* FILE=colvar
      74             : \endplumedfile
      75             : 
      76             : The following command could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
      77             : sodium atoms in sodium chloride.  The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
      78             : with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions.  Once again the average Q6 parameter is calculated and output to a
      79             : file called colvar
      80             : 
      81             : \plumedfile
      82             : Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q6
      83             : PRINT ARG=q6.mean FILE=colvar
      84             : \endplumedfile
      85             : 
      86             : If you simply want to examine the values of the Q6 parameters for each of the atoms in your system you can do so by exploiting the
      87             : command \ref DUMPMULTICOLVAR as shown in the example below.  The following output file will output a file in an extended xyz format
      88             : called q6.xyz for each frame of the analyzed MD trajectory.  The first column in this file will contain a dummy name for each of the
      89             : atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q6 parameter, columns
      90             : 6-19 will contain the real parts of the director of the \f$q_{6m}\f$ vector while columns 20-33 will contain the imaginary parts of this director.
      91             : 
      92             : \plumedfile
      93             : q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
      94             : DUMPMULTICOLVAR DATA=q6 FILE=q6.xyz
      95             : \endplumedfile
      96             : 
      97             : */
      98             : //+ENDPLUMEDOC
      99             : 
     100             : //+PLUMEDOC MCOLVARF LOCAL_Q6
     101             : /*
     102             : Calculate the local degree of order around an atoms by taking the average dot product between the \f$q_6\f$ vector on the central atom and the \f$q_6\f$ vector on the atoms in the first coordination sphere.
     103             : 
     104             : The \ref Q6 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
     105             : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
     106             : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
     107             : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
     108             : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
     109             : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
     110             : because the number of atoms is relatively small.
     111             : 
     112             : When the average \ref Q6 parameter is used to bias the dynamics a problems
     113             : can occur. These averaged coordinates cannot distinguish between the correct,
     114             : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
     115             : themselves into their solid-like configuration simultaneously. This second type
     116             : of pathway would be impossible in reality because there is a large entropic
     117             : barrier that prevents concerted processes like this from happening.  However,
     118             : in the finite sized systems that are commonly simulated this barrier is reduced
     119             : substantially. As a result in simulations where average Steinhardt parameters
     120             : are biased there are often quite dramatic system size effects
     121             : 
     122             : If one wants to simulate nucleation using some form on
     123             : biased dynamics what is really required is an order parameter that measures:
     124             : 
     125             : - Whether or not the coordination spheres around atoms are ordered
     126             : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
     127             : 
     128             : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
     129             : LOCAL_Q6 is another variable that can be used in these sorts of calculations. The LOCAL_Q6 parameter for a particular atom is a number that measures the extent to
     130             : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom.  It does this by calculating the following
     131             : quantity for each of the atoms in the system:
     132             : 
     133             : \f[
     134             :  s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
     135             : \f]
     136             : 
     137             : where \f$q_{6m}(i)\f$ and \f$q_{6m}(j)\f$ are the sixth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
     138             : complex conjugation.  The function
     139             : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$.  The parameters of this function should be set
     140             : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.  The sum in the numerator
     141             : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
     142             : adjacent atoms is correlated.
     143             : 
     144             : \par Examples
     145             : 
     146             : The following command calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
     147             : quantity to a file called colvar.
     148             : 
     149             : \plumedfile
     150             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
     151             : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq6
     152             : PRINT ARG=lq6.mean FILE=colvar
     153             : \endplumedfile
     154             : 
     155             : The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
     156             : 
     157             : \plumedfile
     158             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
     159             : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq6
     160             : PRINT ARG=lq6.* FILE=colvar
     161             : \endplumedfile
     162             : 
     163             : The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
     164             : are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
     165             : 
     166             : \plumedfile
     167             : Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6a
     168             : Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6b
     169             : 
     170             : LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w6
     171             : PRINT ARG=w6.* FILE=colvar
     172             : \endplumedfile
     173             : 
     174             : */
     175             : //+ENDPLUMEDOC
     176             : 
     177             : namespace PLMD {
     178             : namespace crystallization {
     179             : 
     180             : class Q6 : public Steinhardt {
     181             : public:
     182             :   static void registerKeywords( Keywords& keys );
     183             :   explicit Q6( const ActionOptions& ao );
     184             : };
     185             : 
     186             : // For some reason, this is not seen correctly by cppcheck
     187             : // cppcheck-suppress unknownMacro
     188       10449 : PLUMED_REGISTER_ACTION(Q6,"Q6")
     189             : typedef LocalSteinhardt<Q6> LOCAL_Q6;
     190       10429 : PLUMED_REGISTER_ACTION(LOCAL_Q6,"LOCAL_Q6")
     191             : 
     192          16 : void Q6::registerKeywords( Keywords& keys ) {
     193          16 :   Steinhardt::registerKeywords( keys );
     194          16 : }
     195             : 
     196          15 : Q6::Q6(const ActionOptions& ao ):
     197             :   Action(ao),
     198          15 :   Steinhardt(ao)
     199             : {
     200          15 :   setAngularMomentum(6);
     201             : 
     202          15 :   normaliz.resize( 7 );
     203          15 :   normaliz[0] = sqrt( ( 13.0*720.0 ) / (4.0*pi*720.0) );
     204          15 :   normaliz[1] = -sqrt( ( 13.0*120.0 ) / (4.0*pi*5040) );
     205          15 :   normaliz[2] = sqrt( ( 13.0*24) / (4.0*pi*40320) );
     206          15 :   normaliz[3] = -sqrt( ( 13.0*6) / (4.0*pi*362880) );
     207          15 :   normaliz[4] = sqrt( (13.0*2) / (4.0*pi*3628800) );
     208          15 :   normaliz[5] = -sqrt( (13.0*1) / (4.0*pi*39916800) );
     209          15 :   normaliz[6] = sqrt( (13.0*1) / (4.0*pi*479001600) );
     210             : 
     211          15 :   coeff_poly.resize( 7 );
     212          15 :   coeff_poly[0]=-0.3125; coeff_poly[1]=0.0;
     213          15 :   coeff_poly[2]=6.5625; coeff_poly[3]=0.0;
     214          15 :   coeff_poly[4]=-19.6875; coeff_poly[5]=0.0;
     215          15 :   coeff_poly[6]=14.4375;
     216          15 : }
     217             : 
     218             : }
     219             : }
     220             :