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 6 : keys.setValueDescription("vector","the values of the ROPS order parameters");
67 6 : keys.needsAction("DISTANCE_MATRIX"); keys.needsAction("QUATERNION_PRODUCT_MATRIX");
68 9 : keys.needsAction("ONES"); keys.needsAction("CUSTOM"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
69 3 : }
70 :
71 1 : RopsShortcut::RopsShortcut(const ActionOptions&ao):
72 : Action(ao),
73 1 : ActionShortcut(ao)
74 : {
75 : // Open a file and read in the kernels
76 : double h_dops,h_rops; std::string kfunc, kfunc_dops,kfunc_rops,fname_dops,fname_rops;
77 3 : parse("KERNELFILE_DOPS",fname_dops); parse("KERNELFILE_ROPS",fname_rops); IFile ifile_dops, ifile_rops; ifile_dops.open(fname_dops); ifile_rops.open(fname_rops);
78 10 : for(unsigned k=0;; ++k) {
79 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_rops.scanField("height",h_rops) ) break;//checks eof
80 30 : std::string ktype_dops, ktype_rops; ifile_dops.scanField("kerneltype",ktype_dops); ifile_rops.scanField("kerneltype",ktype_rops);
81 10 : if( ktype_dops!="gaussian" ) error("cannot process kernels of type " + ktype_dops );//straightup error
82 10 : if( ktype_rops!="gaussian" ) error("cannot process kernels of type " + ktype_rops );
83 :
84 : double mu_dops,mu_w, mu_i, mu_j, mu_k; std::string hstr_dops, hstr_rops, smu_dops,smu_w, smu_i, smu_j, smu_k, sigmastr,kappastr;
85 :
86 :
87 10 : Tools::convert( h_dops, hstr_dops );
88 10 : Tools::convert( h_rops, hstr_rops );
89 :
90 20 : ifile_dops.scanField("mu",mu_dops); Tools::convert( mu_dops, smu_dops );
91 20 : ifile_rops.scanField("mu_w",mu_w); Tools::convert( mu_w, smu_w );
92 20 : ifile_rops.scanField("mu_i",mu_i); Tools::convert( mu_i, smu_i );
93 20 : ifile_rops.scanField("mu_j",mu_j); Tools::convert( mu_j, smu_j );
94 20 : ifile_rops.scanField("mu_k",mu_k); Tools::convert( mu_k, smu_k );
95 :
96 :
97 : double sigma,kappa;
98 20 : ifile_dops.scanField("sigma",sigma); Tools::convert( sigma, sigmastr );
99 20 : ifile_rops.scanField("kappa",kappa); Tools::convert( kappa, kappastr );
100 :
101 :
102 :
103 10 : ifile_dops.scanField(); /*if( k==0 )*/ kfunc_dops = hstr_dops; //else kfunc_dops += "+" + hstr;
104 10 : ifile_rops.scanField(); /*if( k==0 )*/ kfunc_rops = hstr_rops; //else kfunc_rops += "+" + hstr;
105 :
106 20 : kfunc_rops += "*exp(" + kappastr + "*(w*" + smu_w + "+i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + ")^2)";
107 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
108 20 : if (k==0) kfunc = kfunc_dops + "*" + kfunc_rops; else kfunc+= "+" + kfunc_dops + "*" + kfunc_rops;
109 10 : }
110 : std::string sp_str, specA, specB, grpinfo;
111 : double cutoff;
112 5 : parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB); parse("CUTOFF",cutoff);
113 1 : if( sp_str.length()>0 ) {
114 2 : grpinfo="GROUP=" + sp_str;
115 : } else {//not sure how to use this
116 0 : if( specA.length()==0 || specB.length()==0 ) error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
117 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
118 : }
119 1 : std::string cutstr; Tools::convert( cutoff, cutstr );
120 : // Setup the contact matrix
121 : // std::string switchstr; parse("SWITCH",switchstr);
122 2 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr);
123 :
124 1 : if( specA.length()==0 ) {
125 1 : std::string quatstr; parse("QUATERNIONS",quatstr);
126 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_PRODUCT_MATRIX ARG=" + quatstr + ".*," + quatstr + ".*" );
127 : } else {
128 0 : plumed_error();
129 : }
130 : //
131 2 : readInputLine( getShortcutLabel() + "_kfunc: CUSTOM ARG=" + getShortcutLabel() + "_cmat,"+ getShortcutLabel() + "_quatprod.* " + "VAR=x,w,i,j,k PERIODIC=NO FUNC=" + kfunc );
132 : // Element wise product of cmat and kfunc
133 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
134 : // Find the number of ones we need to multiply by
135 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
136 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
137 1 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
138 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
139 : //
140 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
141 2 : }
142 :
143 : }
144 : }
145 :
146 :
147 :
|
