Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2018 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 : \verbatim
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 : \endverbatim
70 :
71 : */
72 : //+ENDPLUMEDOC
73 :
74 :
75 6 : 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 ;
82 : };
83 :
84 2526 : PLUMED_REGISTER_ACTION(Fccubic,"FCCUBIC")
85 :
86 4 : void Fccubic::registerKeywords( Keywords& keys ) {
87 4 : CubicHarmonicBase::registerKeywords( keys );
88 4 : keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function");
89 4 : }
90 :
91 3 : Fccubic::Fccubic(const ActionOptions&ao):
92 : Action(ao),
93 3 : CubicHarmonicBase(ao)
94 : {
95 : // Scaling factors such that '1' corresponds to fcc lattice
96 : // and '0' corresponds to isotropic (liquid)
97 3 : parse("ALPHA",alpha);
98 3 : a1 = 80080. / (2717. + 16*alpha); b1 = 16.*(alpha-143)/(2717+16*alpha);
99 3 : log.printf(" setting alpha paramter equal to %f \n",alpha);
100 : // And setup the ActionWithVessel
101 3 : checkRead();
102 3 : }
103 :
104 320184 : double Fccubic::calculateCubicHarmonic( const Vector& distance, const double& d2, Vector& myder ) const {
105 320184 : double x2 = distance[0]*distance[0];
106 320225 : double x4 = x2*x2;
107 :
108 320225 : double y2 = distance[1]*distance[1];
109 320225 : double y4 = y2*y2;
110 :
111 320225 : double z2 = distance[2]*distance[2];
112 320224 : double z4 = z2*z2;
113 :
114 320224 : double r8 = pow( d2, 4 );
115 320188 : double r12 = pow( d2, 6 );
116 :
117 320196 : double tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-alpha*x4*y4*z4/r12;
118 :
119 320196 : double t0 = (x2*y4+x2*z4)/r8-alpha*x2*y4*z4/r12;
120 320196 : double t1 = (y2*x4+y2*z4)/r8-alpha*y2*x4*z4/r12;
121 320196 : double t2 = (z2*x4+z2*y4)/r8-alpha*z2*x4*y4/r12;
122 320196 : double t3 = (2*tmp-alpha*x4*y4*z4/r12)/d2;
123 :
124 320196 : myder[0]=4*a1*distance[0]*(t0-t3);
125 320242 : myder[1]=4*a1*distance[1]*(t1-t3);
126 320233 : myder[2]=4*a1*distance[2]*(t2-t3);
127 :
128 320234 : return a1*tmp+b1;
129 : }
130 :
131 : }
132 2523 : }
133 :
|
