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 the BOPS order parameter
33 :
34 : \par Examples
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : class BopsShortcut : public ActionShortcut {
40 : public:
41 : static void registerKeywords( Keywords& keys );
42 : explicit BopsShortcut(const ActionOptions&);
43 : };
44 :
45 : PLUMED_REGISTER_ACTION(BopsShortcut,"BOPS")
46 :
47 3 : void BopsShortcut::registerKeywords( Keywords& keys ) {
48 3 : ActionShortcut::registerKeywords( keys );
49 6 : keys.add("atoms","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
50 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
51 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
52 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
53 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
54 : "coordination number more than four for example");
55 6 : keys.add("atoms-2","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
56 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
57 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
58 : "using the label of another multicolvar");
59 6 : keys.add("atoms-2","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
60 : "the documentation for that keyword");
61 6 : keys.add("compulsory","QUATERNIONS","the label of the action that computes the quaternions that should be used");
62 6 : 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)");
63 6 : 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.");
64 6 : keys.add("compulsory", "CUTOFF", "cutoff for the distance matrix");
65 : // keys.add("compulsory","SWITCH","the switching function that acts on the distances between points)");
66 3 : keys.setValueDescription("the values of the bops order parameters");
67 3 : keys.needsAction("DISTANCE_MATRIX");
68 3 : keys.needsAction("QUATERNION_BOND_PRODUCT_MATRIX");
69 3 : keys.needsAction("CUSTOM");
70 3 : keys.needsAction("ONES");
71 3 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
72 3 : }
73 :
74 1 : BopsShortcut::BopsShortcut(const ActionOptions&ao):
75 : Action(ao),
76 1 : ActionShortcut(ao) {
77 : // Open a file and read in the kernels
78 : double h_dops,h_bops;
79 : std::string kfunc, kfunc_dops,kfunc_bops,fname_dops,fname_bops;
80 1 : parse("KERNELFILE_DOPS",fname_dops);
81 1 : parse("KERNELFILE_BOPS",fname_bops);
82 1 : IFile ifile_dops, ifile_bops;
83 1 : ifile_dops.open(fname_dops);
84 1 : ifile_bops.open(fname_bops);
85 10 : for(unsigned k=0;; ++k) {
86 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_bops.scanField("height",h_bops) ) {
87 : break; //checks eof
88 : }
89 : std::string ktype_dops, ktype_bops;
90 10 : ifile_dops.scanField("kerneltype",ktype_dops);
91 20 : ifile_bops.scanField("kerneltype",ktype_bops);
92 10 : if( ktype_dops!="gaussian" ) {
93 0 : error("cannot process kernels of type " + ktype_dops ); //straightup error
94 : }
95 10 : if( ktype_bops!="gaussian" ) {
96 0 : error("cannot process kernels of type " + ktype_bops );
97 : }
98 :
99 : double mu_dops, mu_i, mu_j, mu_k;
100 : std::string hstr_dops, hstr_bops, smu_dops,smu_i, smu_j, smu_k, sigmastr,kappastr;
101 :
102 :
103 10 : Tools::convert( h_dops, hstr_dops );
104 10 : Tools::convert( h_bops, hstr_bops );
105 :
106 10 : ifile_dops.scanField("mu",mu_dops);
107 10 : Tools::convert( mu_dops, smu_dops );
108 : //ifile_bops.scanField("mu_w",mu_w); Tools::convert( mu_w, smu_w );
109 10 : ifile_bops.scanField("mu_i",mu_i);
110 10 : Tools::convert( mu_i, smu_i );
111 10 : ifile_bops.scanField("mu_j",mu_j);
112 10 : Tools::convert( mu_j, smu_j );
113 10 : ifile_bops.scanField("mu_k",mu_k);
114 10 : Tools::convert( mu_k, smu_k );
115 :
116 :
117 : double sigma,kappa;
118 10 : ifile_dops.scanField("sigma",sigma);
119 10 : Tools::convert( sigma, sigmastr );
120 10 : ifile_bops.scanField("kappa",kappa);
121 10 : Tools::convert( kappa, kappastr );
122 :
123 :
124 :
125 10 : ifile_dops.scanField(); /*if( k==0 )*/ kfunc_dops = hstr_dops; //else kfunc_dops += "+" + hstr;
126 10 : ifile_bops.scanField(); /*if( k==0 )*/ kfunc_bops = hstr_bops; //else kfunc_bops += "+" + hstr;
127 :
128 20 : kfunc_bops += "*exp(" + kappastr + "*(i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + "))";
129 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
130 10 : if (k==0) {
131 2 : kfunc = kfunc_dops + "*" + kfunc_bops;
132 : } else {
133 18 : kfunc+= "+" + kfunc_dops + "*" + kfunc_bops;
134 : }
135 10 : }
136 : std::string sp_str, specA, specB, grpinfo;
137 : double cutoff;
138 1 : parse("SPECIES",sp_str);
139 1 : parse("SPECIESA",specA);
140 1 : parse("SPECIESB",specB);
141 2 : parse("CUTOFF",cutoff);
142 1 : if( sp_str.length()>0 ) {
143 2 : grpinfo="GROUP=" + sp_str;
144 : } else {//not sure how to use this
145 0 : if( specA.length()==0 || specB.length()==0 ) {
146 0 : error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
147 : }
148 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
149 : }
150 : std::string cutstr;
151 1 : Tools::convert( cutoff, cutstr );
152 : // Setup the contact matrix
153 : // std::string switchstr; parse("SWITCH",switchstr);
154 2 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr + " COMPONENTS");
155 :
156 1 : if( specA.length()==0 ) {
157 : std::string quatstr;
158 1 : parse("QUATERNIONS",quatstr);
159 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_BOND_PRODUCT_MATRIX ARG=" + quatstr + ".*," + getShortcutLabel() + "_cmat.*" );
160 : } else {
161 0 : plumed_error();
162 : }
163 : //
164 :
165 : ///////////////////
166 : ///replace/////
167 2 : readInputLine( getShortcutLabel() + "_dist: CUSTOM ARG=" + getShortcutLabel() + "_cmat.x," + getShortcutLabel() + "_cmat.y," + getShortcutLabel() + "_cmat.z " +
168 : "FUNC=sqrt((x^2)+(y^2)+(z^2)) PERIODIC=NO");
169 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");
170 :
171 : //replace ^^^ to remove distance hack
172 : //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");
173 : ///end replace////
174 :
175 : // Element wise product of cmat and kfunc
176 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
177 : // Find the number of ones we need to multiply by
178 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
179 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
180 : std::string size;
181 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
182 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
183 : //
184 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
185 2 : }
186 :
187 : }
188 : }
189 :
190 :
191 :
|
