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 ROPS
31 : /*
32 : Calculate the ROPS order parameter
33 :
34 : \par Examples
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : class RopsShortcut : public ActionShortcut {
40 : public:
41 : static void registerKeywords( Keywords& keys );
42 : explicit RopsShortcut(const ActionOptions&);
43 : };
44 :
45 : PLUMED_REGISTER_ACTION(RopsShortcut,"ROPS")
46 :
47 3 : void RopsShortcut::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_ROPS","the file containing the list of kernel parameters. We expect the normalization factor (height), concentration parameter (kappa), and 4 quaternion pieces of the mean for a bipolar watson distribution (mu_w,mu_i,mu_j,mu_k)): (h*exp(kappa*dot(q_mean,q))), where dot is the 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 ROPS order parameters");
67 3 : keys.needsAction("DISTANCE_MATRIX");
68 3 : keys.needsAction("QUATERNION_PRODUCT_MATRIX");
69 3 : keys.needsAction("ONES");
70 3 : keys.needsAction("CUSTOM");
71 3 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
72 3 : }
73 :
74 1 : RopsShortcut::RopsShortcut(const ActionOptions&ao):
75 : Action(ao),
76 1 : ActionShortcut(ao) {
77 : // Open a file and read in the kernels
78 : double h_dops,h_rops;
79 : std::string kfunc, kfunc_dops,kfunc_rops,fname_dops,fname_rops;
80 1 : parse("KERNELFILE_DOPS",fname_dops);
81 1 : parse("KERNELFILE_ROPS",fname_rops);
82 1 : IFile ifile_dops, ifile_rops;
83 1 : ifile_dops.open(fname_dops);
84 1 : ifile_rops.open(fname_rops);
85 10 : for(unsigned k=0;; ++k) {
86 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_rops.scanField("height",h_rops) ) {
87 : break; //checks eof
88 : }
89 : std::string ktype_dops, ktype_rops;
90 10 : ifile_dops.scanField("kerneltype",ktype_dops);
91 20 : ifile_rops.scanField("kerneltype",ktype_rops);
92 10 : if( ktype_dops!="gaussian" ) {
93 0 : error("cannot process kernels of type " + ktype_dops ); //straightup error
94 : }
95 10 : if( ktype_rops!="gaussian" ) {
96 0 : error("cannot process kernels of type " + ktype_rops );
97 : }
98 :
99 : double mu_dops,mu_w, mu_i, mu_j, mu_k;
100 : std::string hstr_dops, hstr_rops, smu_dops,smu_w, smu_i, smu_j, smu_k, sigmastr,kappastr;
101 :
102 :
103 10 : Tools::convert( h_dops, hstr_dops );
104 10 : Tools::convert( h_rops, hstr_rops );
105 :
106 10 : ifile_dops.scanField("mu",mu_dops);
107 10 : Tools::convert( mu_dops, smu_dops );
108 10 : ifile_rops.scanField("mu_w",mu_w);
109 10 : Tools::convert( mu_w, smu_w );
110 10 : ifile_rops.scanField("mu_i",mu_i);
111 10 : Tools::convert( mu_i, smu_i );
112 10 : ifile_rops.scanField("mu_j",mu_j);
113 10 : Tools::convert( mu_j, smu_j );
114 10 : ifile_rops.scanField("mu_k",mu_k);
115 10 : Tools::convert( mu_k, smu_k );
116 :
117 :
118 : double sigma,kappa;
119 10 : ifile_dops.scanField("sigma",sigma);
120 10 : Tools::convert( sigma, sigmastr );
121 10 : ifile_rops.scanField("kappa",kappa);
122 10 : Tools::convert( kappa, kappastr );
123 :
124 :
125 :
126 10 : ifile_dops.scanField(); /*if( k==0 )*/ kfunc_dops = hstr_dops; //else kfunc_dops += "+" + hstr;
127 10 : ifile_rops.scanField(); /*if( k==0 )*/ kfunc_rops = hstr_rops; //else kfunc_rops += "+" + hstr;
128 :
129 20 : kfunc_rops += "*exp(" + kappastr + "*(w*" + smu_w + "+i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + ")^2)";
130 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
131 10 : if (k==0) {
132 2 : kfunc = kfunc_dops + "*" + kfunc_rops;
133 : } else {
134 18 : kfunc+= "+" + kfunc_dops + "*" + kfunc_rops;
135 : }
136 10 : }
137 : std::string sp_str, specA, specB, grpinfo;
138 : double cutoff;
139 1 : parse("SPECIES",sp_str);
140 1 : parse("SPECIESA",specA);
141 1 : parse("SPECIESB",specB);
142 2 : parse("CUTOFF",cutoff);
143 1 : if( sp_str.length()>0 ) {
144 2 : grpinfo="GROUP=" + sp_str;
145 : } else {//not sure how to use this
146 0 : if( specA.length()==0 || specB.length()==0 ) {
147 0 : error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
148 : }
149 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
150 : }
151 : std::string cutstr;
152 1 : Tools::convert( cutoff, cutstr );
153 : // Setup the contact matrix
154 : // std::string switchstr; parse("SWITCH",switchstr);
155 2 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr);
156 :
157 1 : if( specA.length()==0 ) {
158 : std::string quatstr;
159 1 : parse("QUATERNIONS",quatstr);
160 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_PRODUCT_MATRIX ARG=" + quatstr + ".*," + quatstr + ".*" );
161 : } else {
162 0 : plumed_error();
163 : }
164 : //
165 2 : readInputLine( getShortcutLabel() + "_kfunc: CUSTOM ARG=" + getShortcutLabel() + "_cmat,"+ getShortcutLabel() + "_quatprod.* " + "VAR=x,w,i,j,k PERIODIC=NO FUNC=" + kfunc );
166 : // Element wise product of cmat and kfunc
167 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
168 : // Find the number of ones we need to multiply by
169 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
170 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
171 : std::string size;
172 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
173 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
174 : //
175 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
176 2 : }
177 :
178 : }
179 : }
180 :
181 :
182 :
|
