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 "CoordinationNumbers.h"
23 : #include "multicolvar/MultiColvarShortcuts.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/ActionShortcut.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 :
30 : #include <string>
31 : #include <cmath>
32 :
33 : namespace PLMD {
34 : namespace symfunc {
35 :
36 : //+PLUMEDOC MCOLVAR COORDINATION_SHELL_FUNCTION
37 : /*
38 : Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom
39 :
40 : \par Examples
41 :
42 :
43 : */
44 : //+ENDPLUMEDOC
45 :
46 : //+PLUMEDOC MCOLVAR COORDINATION_SHELL_AVERAGE
47 : /*
48 : Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom and take an average
49 :
50 : \par Examples
51 :
52 :
53 : */
54 : //+ENDPLUMEDOC
55 :
56 : //+PLUMEDOC MCOLVAR SIMPLECUBIC
57 : /*
58 : Calculate whether or not the coordination spheres of atoms are arranged as they would be in a simple cubic structure.
59 :
60 : We can measure how similar the environment around atom \f$i\f$ is to a simple cubic structure is by evaluating
61 : the following quantity:
62 :
63 : \f[
64 : s_i = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \left[ \frac{ x_{ij}^4 + y_{ij}^4 + z_{ij}^4 }{r_{ij}^4} \right] }{ \sum_{i \ne j} \sigma(r_{ij}) }
65 : \f]
66 :
67 : 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
68 : 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 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
69 : over all of the other atoms in the system ensures that we are calculating an average
70 : 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$.
71 : This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
72 : 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
73 : so on. Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
74 :
75 :
76 : \par Examples
77 :
78 : The following input tells plumed to calculate the simple cubic parameter for the atoms 1-100 with themselves.
79 : The mean value is then calculated.
80 : \plumedfile
81 : SIMPLECUBIC SPECIES=1-100 R_0=1.0 MEAN
82 : \endplumedfile
83 :
84 : The following input tells plumed to look at the ways atoms 1-100 are within 3.0 are arranged about atoms
85 : from 101-110. The number of simple cubic parameters that are greater than 0.8 is then output
86 : \plumedfile
87 : SIMPLECUBIC SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=0.8 NN=6 MM=12 D_0=0}
88 : \endplumedfile
89 :
90 : */
91 : //+ENDPLUMEDOC
92 :
93 : //+PLUMEDOC MCOLVAR TETRAHEDRAL
94 : /*
95 : Calculate the degree to which the environment about ions has a tetrahedral order.
96 :
97 : We can measure the degree to which the atoms in the first coordination shell around any atom, \f$i\f$ is
98 : is arranged like a tetrahedron using the following function.
99 :
100 : \f[
101 : s(i) = \frac{1}{\sum_j \sigma( r_{ij} )} \sum_j \sigma( r_{ij} )\left[ \frac{(x_{ij} + y_{ij} + z_{ij})^3}{r_{ij}^3} +
102 : \frac{(x_{ij} - y_{ij} - z_{ij})^3}{r_{ij}^3} +
103 : \frac{(-x_{ij} + y_{ij} - z_{ij})^3}{r_{ij}^3} +
104 : \frac{(-x_{ij} - y_{ij} + z_{ij})^3}{r_{ij}^3} \right]
105 : \f]
106 :
107 : Here \f$r_{ij}\f$ is the magnitude of the vector connecting atom \f$i\f$ to atom \f$j\f$ and \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$
108 : are its three components. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
109 : atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that the function is equal to one
110 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
111 :
112 : \par Examples
113 :
114 : The following command calculates the average value of the TETRAHEDRAL parameter for a set of 64 atoms all of the same type
115 : and outputs this quantity to a file called colvar.
116 :
117 : \plumedfile
118 : tt: TETRAHEDRAL SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
119 : PRINT ARG=tt.mean FILE=colvar
120 : \endplumedfile
121 :
122 : The following command calculates the number of TETRAHEDRAL parameters that are greater than 0.8 in a set of 10 atoms.
123 : In this calculation it is assumed that there are two atom types A and B and that the first coordination sphere of the
124 : 10 atoms of type A contains atoms of type B. The formula above is thus calculated for ten different A atoms and within
125 : it the sum over \f$j\f$ runs over 40 atoms of type B that could be in the first coordination sphere.
126 :
127 : \plumedfile
128 : tt: TETRAHEDRAL SPECIESA=1-10 SPECIESB=11-40 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=0.8}
129 : PRINT ARG=tt.* FILE=colvar
130 : \endplumedfile
131 :
132 : */
133 : //+ENDPLUMEDOC
134 :
135 : class CoordShellVectorFunction : public ActionShortcut {
136 : public:
137 : static void registerKeywords(Keywords& keys);
138 : explicit CoordShellVectorFunction(const ActionOptions&);
139 : };
140 :
141 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"FCCUBIC")
142 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"TETRAHEDRAL")
143 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"SIMPLECUBIC")
144 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_FUNCTION")
145 : PLUMED_REGISTER_ACTION(CoordShellVectorFunction,"COORDINATION_SHELL_AVERAGE")
146 :
147 23 : void CoordShellVectorFunction::registerKeywords( Keywords& keys ) {
148 23 : CoordinationNumbers::shortcutKeywords( keys );
149 46 : keys.add("compulsory","FUNCTION","the function of the bond vectors that you would like to evaluate");
150 46 : keys.add("compulsory","PHI","0.0","The Euler rotational angle phi");
151 46 : keys.add("compulsory","THETA","0.0","The Euler rotational angle theta");
152 46 : keys.add("compulsory","PSI","0.0","The Euler rotational angle psi");
153 46 : keys.add("compulsory","ALPHA","3.0","The alpha parameter of the angular function that is used for FCCUBIC");
154 46 : keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
155 69 : keys.needsAction("CONTACT_MATRIX"); keys.needsAction("FCCUBIC_FUNC"); keys.needsAction("CUSTOM");
156 46 : keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
157 23 : }
158 :
159 9 : CoordShellVectorFunction::CoordShellVectorFunction(const ActionOptions& ao):
160 : Action(ao),
161 9 : ActionShortcut(ao)
162 : {
163 9 : std::string matlab, sp_str, specA, specB; bool lowmem; parseFlag("LOWMEM",lowmem);
164 9 : if( lowmem ) warning("LOWMEM flag is deprecated and is no longer required for this action");
165 36 : parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
166 9 : if( sp_str.length()>0 || specA.length()>0 ) {
167 9 : matlab = getShortcutLabel() + "_mat";
168 9 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this );
169 0 : } else error("found no input atoms use SPECIES/SPECIESA");
170 27 : double phi, theta, psi; parse("PHI",phi); parse("THETA",theta); parse("PSI",psi);
171 9 : std::vector<std::string> rotelements(9); std::string xvec = matlab + ".x", yvec = matlab + ".y", zvec = matlab + ".z";
172 9 : if( phi!=0 || theta!=0 || psi!=0 ) {
173 0 : Tools::convert( std::cos(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::sin(psi), rotelements[0] );
174 0 : Tools::convert( std::cos(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::sin(psi), rotelements[1] );
175 0 : Tools::convert( std::sin(psi)*std::sin(theta), rotelements[2] );
176 :
177 0 : Tools::convert( -std::sin(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::cos(psi), rotelements[3] );
178 0 : Tools::convert( -std::sin(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::cos(psi), rotelements[4] );
179 0 : Tools::convert( std::cos(psi)*std::sin(theta), rotelements[5] );
180 :
181 0 : Tools::convert( std::sin(theta)*std::sin(phi), rotelements[6] );
182 0 : Tools::convert( -std::sin(theta)*std::cos(phi), rotelements[7] );
183 0 : Tools::convert( std::cos(theta), rotelements[8] );
184 0 : readInputLine( getShortcutLabel() + "_xrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[0] + "*x+" + rotelements[1] + "*y+" + rotelements[2] + "*z PERIODIC=NO");
185 0 : readInputLine( getShortcutLabel() + "_yrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[3] + "*x+" + rotelements[4] + "*y+" + rotelements[5] + "*z PERIODIC=NO");
186 0 : readInputLine( getShortcutLabel() + "_zrot: CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z FUNC=" + rotelements[6] + "*x+" + rotelements[7] + "*y+" + rotelements[8] + "*z PERIODIC=NO");
187 : }
188 : // Calculate FCC cubic function from bond vectors
189 9 : if( getName()=="FCCUBIC" ) {
190 5 : std::string alpha; parse("ALPHA",alpha);
191 10 : readInputLine( getShortcutLabel() + "_vfunc: FCCUBIC_FUNC ARG=" + xvec + "," + yvec + "," + zvec+ " ALPHA=" + alpha);
192 4 : } else if( getName()=="TETRAHEDRAL" ) {
193 2 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
194 3 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r"
195 2 : + " VAR=x,y,z,r PERIODIC=NO FUNC=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3" );
196 3 : } else if( getName()=="SIMPLECUBIC" ) {
197 2 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
198 3 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r"
199 2 : + " VAR=x,y,z,r PERIODIC=NO FUNC=(x^4+y^4+z^4)/(r^4)" );
200 : } else {
201 2 : std::string myfunc; parse("FUNCTION",myfunc);
202 2 : if( myfunc.find("r")!=std::string::npos ) {
203 4 : readInputLine( getShortcutLabel() + "_r: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=sqrt(x*x+y*y+z*z)");
204 4 : readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + "," + getShortcutLabel() + "_r VAR=x,y,z,r PERIODIC=NO FUNC=" + myfunc );
205 0 : } else readInputLine( getShortcutLabel() + "_vfunc: CUSTOM ARG=" + xvec + "," + yvec + "," + zvec + " PERIODIC=NO FUNC=" + myfunc );
206 : }
207 : // Hadamard product of function above and weights
208 18 : readInputLine( getShortcutLabel() + "_wvfunc: CUSTOM ARG=" + getShortcutLabel() + "_vfunc," + matlab + ".w FUNC=x*y PERIODIC=NO");
209 : // And coordination numbers
210 9 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
211 9 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
212 9 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
213 18 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
214 18 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_wvfunc," + getShortcutLabel() + "_ones");
215 9 : std::string olab=getShortcutLabel();
216 9 : if( getName()!="COORDINATION_SHELL_FUNCTION" ) {
217 7 : olab = getShortcutLabel() + "_n";
218 : // Calculate coordination numbers for denominator
219 14 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + matlab + ".w," + getShortcutLabel() + "_ones");
220 : // And normalise
221 14 : readInputLine( getShortcutLabel() + "_n: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
222 : }
223 : // And expand the functions
224 9 : std::map<std::string,std::string> keymap; multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
225 18 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), olab, "", keymap, this );
226 18 : }
227 :
228 : }
229 : }
230 :
|
