Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) crystdistrib 2023-2023 The code team
3 : (see the PEOPLE-crystdistrib file at the root of this folder for a list of names)
4 :
5 : This file is part of crystdistrib code module.
6 :
7 : The crystdistrib code module is free software: you can redistribute it and/or modify
8 : it under the terms of the GNU Lesser General Public License as published by
9 : the Free Software Foundation, either version 3 of the License, or
10 : (at your option) any later version.
11 :
12 : The crystdistrib code module is distributed in the hope that it will be useful,
13 : but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : GNU Lesser General Public License for more details.
16 :
17 : You should have received a copy of the GNU Lesser General Public License
18 : along with the crystdistrib code module. If not, see <http://www.gnu.org/licenses/>.
19 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
20 : #include "core/ActionShortcut.h"
21 : #include "core/ActionRegister.h"
22 : #include "core/PlumedMain.h"
23 : #include "core/ActionSet.h"
24 : #include "core/ActionWithValue.h"
25 : #include "tools/IFile.h"
26 :
27 : namespace PLMD {
28 : namespace crystdistrib {
29 :
30 : //+PLUMEDOC COLVAR BOPS
31 : /*
32 : Calculate Bond orientational order parameters for molecules.
33 :
34 : BOPS is a shortcut to calculate the Bond-orientational Order Parameters detailed that are described in the paper cited below.
35 : As arguments, BOPS takes a list of atoms (corresponding to molecules), a vector of quaternions, a cutoff distance, and two kernel files
36 : detailing the means, variances, and normalization factors of probability distributions. BOPS returns a vector of order parameters.
37 :
38 : The DOPS kernel file has FIELDS height, mu, and sigma corresponding to the normalization factor, mean, and variance of the gaussian distributions used in the order parameters.
39 : The SET kerneltype is gaussian.
40 :
41 : The BOPS kernel file has FIELDS height, kappa, mu\_i, mu\_j, and mu\_k, which correspond to the normalization factor, reciprocal variance, and imaginary components of the
42 : mean quaternion frame of the fisher distribution used in the order parameters. The SET kerneltype is gaussian.
43 :
44 : BOPS returns one order parameter per atom given, evaluated over each atom's neighbors within the cutoff given. The distribution defined by the kernel files, analogous to a radial distribution function, is defined over all possible unit vectors which could be drawn between two atoms. The order parameter is obtained by evaluating the distribution at each unit vector pointing to all neighbors within the cutoff, and summing them up.
45 :
46 :
47 : This example file calculates the BOPS for a system of 3 molecules.
48 :
49 : ```plumed
50 : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-bops-shortcut/kernels.dat,regtest/crystdistrib/rt-bops-shortcut/kernels2.dat
51 : quat: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
52 : bops: BOPS ...
53 : SPECIES=1,4,7 QUATERNIONS=quat CUTOFF=100.0
54 : KERNELFILE_DOPS=regtest/crystdistrib/rt-bops-shortcut/kernels.dat
55 : KERNELFILE_BOPS=regtest/crystdistrib/rt-bops-shortcut/kernels2.dat
56 : ...
57 : ```
58 :
59 : To calculate the BOPS between the orientation of the molecules in GROUPA and the bonds to the atoms in GROUPB you use an input like the one shown below:
60 :
61 : ```plumed
62 : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-bops-shortcut/kernels.dat,regtest/crystdistrib/rt-bops-shortcut/kernels2.dat
63 : quat: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
64 : bops: BOPS ...
65 : QUATERNIONS=quat CUTOFF=100.0
66 : SPECIESA=1,4,7 SPECIESB=10,11,12,13,14
67 : KERNELFILE_DOPS=regtest/crystdistrib/rt-bops-shortcut/kernels.dat
68 : KERNELFILE_BOPS=regtest/crystdistrib/rt-bops-shortcut/kernels2.dat
69 : ...
70 : ```
71 :
72 : */
73 : //+ENDPLUMEDOC
74 :
75 : class BopsShortcut : public ActionShortcut {
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit BopsShortcut(const ActionOptions&);
79 : };
80 :
81 : PLUMED_REGISTER_ACTION(BopsShortcut,"BOPS")
82 :
83 3 : void BopsShortcut::registerKeywords( Keywords& keys ) {
84 3 : ActionShortcut::registerKeywords( keys );
85 3 : keys.add("atoms","SPECIES","the list of atoms for which the BOPS are being calculated and the atoms that can be in the environments");
86 3 : keys.add("atoms-2","SPECIESA","the list of atoms for which BOPS are being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment");
87 3 : keys.add("atoms-2","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the BOPS are being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which BOPS are being calculated.");
88 3 : keys.add("compulsory","QUATERNIONS","the label of the action that computes the quaternions that should be used");
89 3 : keys.add("compulsory","KERNELFILE_DOPS","the file containing the list of kernel parameters. We expect h, mu and sigma parameters for a 1D Gaussian kernel of the form h*exp(-(x-mu)^2/2sigma^2)");
90 3 : keys.add("compulsory","KERNELFILE_BOPS","the second file containing the list of kernel parameters. Expecting a normalization factor (height), concentration parameter (kappa), and 3 norm vector pieces of the mean (mu_i, mu_j, mu_k )for a fisher distribution. of the form h*exp(kappa*dot(r_mean,r)), where dot is a standard dot product.");
91 3 : keys.add("compulsory", "CUTOFF", "cutoff for the distance matrix");
92 : // keys.add("compulsory","SWITCH","the switching function that acts on the distances between points)");
93 6 : keys.setValueDescription("vector","the values of the bops order parameters");
94 3 : keys.needsAction("DISTANCE_MATRIX");
95 3 : keys.needsAction("QUATERNION_BOND_PRODUCT_MATRIX");
96 3 : keys.needsAction("CUSTOM");
97 3 : keys.needsAction("ONES");
98 3 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
99 3 : keys.addDOI("10.1063/1.3548889");
100 3 : }
101 :
102 1 : BopsShortcut::BopsShortcut(const ActionOptions&ao):
103 : Action(ao),
104 1 : ActionShortcut(ao) {
105 : // Open a file and read in the kernels
106 : double h_dops,h_bops;
107 : std::string kfunc, kfunc_dops,kfunc_bops,fname_dops,fname_bops;
108 1 : parse("KERNELFILE_DOPS",fname_dops);
109 1 : parse("KERNELFILE_BOPS",fname_bops);
110 1 : IFile ifile_dops, ifile_bops;
111 1 : ifile_dops.open(fname_dops);
112 1 : ifile_bops.open(fname_bops);
113 10 : for(unsigned k=0;; ++k) {
114 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_bops.scanField("height",h_bops) ) {
115 : break; //checks eof
116 : }
117 : std::string ktype_dops, ktype_bops;
118 10 : ifile_dops.scanField("kerneltype",ktype_dops);
119 20 : ifile_bops.scanField("kerneltype",ktype_bops);
120 10 : if( ktype_dops!="gaussian" ) {
121 0 : error("cannot process kernels of type " + ktype_dops ); //straightup error
122 : }
123 10 : if( ktype_bops!="gaussian" ) {
124 0 : error("cannot process kernels of type " + ktype_bops );
125 : }
126 :
127 : double mu_dops, mu_i, mu_j, mu_k;
128 : std::string hstr_dops, hstr_bops, smu_dops,smu_i, smu_j, smu_k, sigmastr,kappastr;
129 :
130 :
131 10 : Tools::convert( h_dops, hstr_dops );
132 10 : Tools::convert( h_bops, hstr_bops );
133 :
134 10 : ifile_dops.scanField("mu",mu_dops);
135 10 : Tools::convert( mu_dops, smu_dops );
136 : //ifile_bops.scanField("mu_w",mu_w); Tools::convert( mu_w, smu_w );
137 10 : ifile_bops.scanField("mu_i",mu_i);
138 10 : Tools::convert( mu_i, smu_i );
139 10 : ifile_bops.scanField("mu_j",mu_j);
140 10 : Tools::convert( mu_j, smu_j );
141 10 : ifile_bops.scanField("mu_k",mu_k);
142 10 : Tools::convert( mu_k, smu_k );
143 :
144 :
145 : double sigma,kappa;
146 10 : ifile_dops.scanField("sigma",sigma);
147 10 : Tools::convert( sigma, sigmastr );
148 10 : ifile_bops.scanField("kappa",kappa);
149 10 : Tools::convert( kappa, kappastr );
150 :
151 :
152 :
153 10 : ifile_dops.scanField(); /*if( k==0 )*/ kfunc_dops = hstr_dops; //else kfunc_dops += "+" + hstr;
154 10 : ifile_bops.scanField(); /*if( k==0 )*/ kfunc_bops = hstr_bops; //else kfunc_bops += "+" + hstr;
155 :
156 20 : kfunc_bops += "*exp(" + kappastr + "*(i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + "))";
157 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
158 10 : if (k==0) {
159 2 : kfunc = kfunc_dops + "*" + kfunc_bops;
160 : } else {
161 18 : kfunc+= "+" + kfunc_dops + "*" + kfunc_bops;
162 : }
163 10 : }
164 : std::string sp_str, specA, specB, grpinfo;
165 : double cutoff;
166 1 : parse("SPECIES",sp_str);
167 1 : parse("SPECIESA",specA);
168 1 : parse("SPECIESB",specB);
169 2 : parse("CUTOFF",cutoff);
170 1 : if( sp_str.length()>0 ) {
171 2 : grpinfo="GROUP=" + sp_str;
172 : } else {//not sure how to use this
173 0 : if( specA.length()==0 || specB.length()==0 ) {
174 0 : error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
175 : }
176 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
177 : }
178 : std::string cutstr;
179 1 : Tools::convert( cutoff, cutstr );
180 : // Setup the contact matrix
181 : // std::string switchstr; parse("SWITCH",switchstr);
182 2 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr + " COMPONENTS");
183 :
184 : // if( specA.length()==0 ) {
185 : std::string quatstr;
186 1 : parse("QUATERNIONS",quatstr);
187 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_BOND_PRODUCT_MATRIX ARG=" + quatstr + ".*," + getShortcutLabel() + "_cmat.*" );
188 : // } else {
189 : // plumed_error();
190 : // }
191 : //
192 :
193 : ///////////////////
194 : ///replace/////
195 2 : readInputLine( getShortcutLabel() + "_dist: CUSTOM ARG=" + getShortcutLabel() + "_cmat.x," + getShortcutLabel() + "_cmat.y," + getShortcutLabel() + "_cmat.z " +
196 : "FUNC=sqrt((x^2)+(y^2)+(z^2)) PERIODIC=NO");
197 2 : readInputLine( getShortcutLabel() + "_kfunc: CUSTOM ARG=" + getShortcutLabel() + "_quatprod.i," + getShortcutLabel() + "_quatprod.j," + getShortcutLabel() + "_quatprod.k,"+ getShortcutLabel() + "_dist " + "VAR=i,j,k,x FUNC=" + kfunc + " PERIODIC=NO");
198 :
199 : //replace ^^^ to remove distance hack
200 : //readInputLine( getShortcutLabel() + "_kfunc: CUSTOM ARG=" + getShortcutLabel() + "_quatprod.i," + getShortcutLabel() + "_quatprod.j," + getShortcutLabel() + "_quatprod.k,"+ getShortcutLabel() + "_cmat.w " + "VAR=i,j,k,x FUNC=" + kfunc + " PERIODIC=NO");
201 : ///end replace////
202 :
203 : // Element wise product of cmat and kfunc
204 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
205 : // Find the number of ones we need to multiply by
206 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
207 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
208 : std::string size;
209 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
210 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
211 : //
212 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
213 2 : }
214 :
215 : }
216 : }
217 :
218 :
219 :
|
