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 : ROPS is a shortcut to calculate the Relative-Orientational Order Parameters that are described in the paper that is referenced below. As arguments, ROPS takes a list of atoms (corresponding to molecules), a vector of quaternions, a cutoff distance, and two kernel files detailing the means, variances, and normalization factors of the probability distributions. ROPS returns a vector of order parameters.
35 :
36 : 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. The SET kerneltype is gaussian.
37 : The ROPS kernel file has FIELDS height, kappa, mu\_w, mu\_i, mu\_j, and mu\_k, which correspond to the normalization factor, reciprocal variance, and components of the mean quaternion frame of the Bipolar Watson distribution used in the order parameters. The SET kerneltype is gaussian.
38 :
39 : ROPS returns one order parameter per atom given, evaluated over each atom's neighbors within the cutoff. The distribution defined by the kernel files, analogous to a radial distribution function, defined over all possible spatial orientations two molecules could take relative to one another. The domain is the set of all unit quaternions. The order parameter is obtained by evaluating the distribution at each neighbor's unit quaternion, and summing them.
40 :
41 : This example file calculates the ROPS for a system of 3 molecules.
42 :
43 : ```plumed
44 : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-rops-shortcut/kernels.dat,regtest/crystdistrib/rt-rops-shortcut/kernels2.dat
45 : quat: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
46 : bops: ROPS ...
47 : SPECIES=1,4,7 QUATERNIONS=quat CUTOFF=100.0
48 : KERNELFILE_DOPS=regtest/crystdistrib/rt-rops-shortcut/kernels.dat
49 : KERNELFILE_ROPS=regtest/crystdistrib/rt-rops-shortcut/kernels2.dat
50 : ...
51 : ```
52 :
53 : */
54 : //+ENDPLUMEDOC
55 :
56 : class RopsShortcut : public ActionShortcut {
57 : public:
58 : static void registerKeywords( Keywords& keys );
59 : explicit RopsShortcut(const ActionOptions&);
60 : };
61 :
62 : PLUMED_REGISTER_ACTION(RopsShortcut,"ROPS")
63 :
64 3 : void RopsShortcut::registerKeywords( Keywords& keys ) {
65 3 : ActionShortcut::registerKeywords( keys );
66 3 : keys.add("atoms","SPECIES","the list of atoms for which the ROPS are being calculated and the atoms that can be in the environments");
67 3 : keys.add("compulsory","QUATERNIONS","the label of the action that computes the quaternions that should be used");
68 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)");
69 3 : 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 ");
70 3 : keys.add("compulsory", "CUTOFF", "cutoff for the distance matrix");
71 : // keys.add("compulsory","SWITCH","the switching function that acts on the distances between points)");
72 6 : keys.setValueDescription("vector","the values of the ROPS order parameters");
73 3 : keys.needsAction("DISTANCE_MATRIX");
74 3 : keys.needsAction("QUATERNION_PRODUCT_MATRIX");
75 3 : keys.needsAction("ONES");
76 3 : keys.needsAction("CUSTOM");
77 3 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
78 3 : keys.addDOI("10.1063/1.3548889");
79 3 : }
80 :
81 1 : RopsShortcut::RopsShortcut(const ActionOptions&ao):
82 : Action(ao),
83 1 : ActionShortcut(ao) {
84 : // Open a file and read in the kernels
85 : double h_dops,h_rops;
86 : std::string kfunc, kfunc_dops,kfunc_rops,fname_dops,fname_rops;
87 1 : parse("KERNELFILE_DOPS",fname_dops);
88 1 : parse("KERNELFILE_ROPS",fname_rops);
89 1 : IFile ifile_dops, ifile_rops;
90 1 : ifile_dops.open(fname_dops);
91 1 : ifile_rops.open(fname_rops);
92 10 : for(unsigned k=0;; ++k) {
93 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_rops.scanField("height",h_rops) ) {
94 : break; //checks eof
95 : }
96 : std::string ktype_dops, ktype_rops;
97 10 : ifile_dops.scanField("kerneltype",ktype_dops);
98 20 : ifile_rops.scanField("kerneltype",ktype_rops);
99 10 : if( ktype_dops!="gaussian" ) {
100 0 : error("cannot process kernels of type " + ktype_dops ); //straightup error
101 : }
102 10 : if( ktype_rops!="gaussian" ) {
103 0 : error("cannot process kernels of type " + ktype_rops );
104 : }
105 :
106 : double mu_dops,mu_w, mu_i, mu_j, mu_k;
107 : std::string hstr_dops, hstr_rops, smu_dops,smu_w, smu_i, smu_j, smu_k, sigmastr,kappastr;
108 :
109 :
110 10 : Tools::convert( h_dops, hstr_dops );
111 10 : Tools::convert( h_rops, hstr_rops );
112 :
113 10 : ifile_dops.scanField("mu",mu_dops);
114 10 : Tools::convert( mu_dops, smu_dops );
115 10 : ifile_rops.scanField("mu_w",mu_w);
116 10 : Tools::convert( mu_w, smu_w );
117 10 : ifile_rops.scanField("mu_i",mu_i);
118 10 : Tools::convert( mu_i, smu_i );
119 10 : ifile_rops.scanField("mu_j",mu_j);
120 10 : Tools::convert( mu_j, smu_j );
121 10 : ifile_rops.scanField("mu_k",mu_k);
122 10 : Tools::convert( mu_k, smu_k );
123 :
124 :
125 : double sigma,kappa;
126 10 : ifile_dops.scanField("sigma",sigma);
127 10 : Tools::convert( sigma, sigmastr );
128 10 : ifile_rops.scanField("kappa",kappa);
129 10 : Tools::convert( kappa, kappastr );
130 :
131 :
132 :
133 10 : ifile_dops.scanField(); /*if( k==0 )*/ kfunc_dops = hstr_dops; //else kfunc_dops += "+" + hstr;
134 10 : ifile_rops.scanField(); /*if( k==0 )*/ kfunc_rops = hstr_rops; //else kfunc_rops += "+" + hstr;
135 :
136 20 : kfunc_rops += "*exp(" + kappastr + "*(w*" + smu_w + "+i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + ")^2)";
137 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
138 10 : if (k==0) {
139 2 : kfunc = kfunc_dops + "*" + kfunc_rops;
140 : } else {
141 18 : kfunc+= "+" + kfunc_dops + "*" + kfunc_rops;
142 : }
143 10 : }
144 : std::string sp_str, specA, specB, grpinfo;
145 : double cutoff;
146 1 : parse("SPECIES",sp_str);
147 : //parse("SPECIESA",specA);
148 : //parse("SPECIESB",specB);
149 2 : parse("CUTOFF",cutoff);
150 1 : if( sp_str.length()>0 ) {
151 2 : grpinfo="GROUP=" + sp_str;
152 : } else {//not sure how to use this
153 0 : if( specA.length()==0 || specB.length()==0 ) {
154 0 : error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
155 : }
156 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
157 : }
158 : std::string cutstr;
159 1 : Tools::convert( cutoff, cutstr );
160 : // Setup the contact matrix
161 : // std::string switchstr; parse("SWITCH",switchstr);
162 2 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr);
163 :
164 1 : if( specA.length()==0 ) {
165 : std::string quatstr;
166 1 : parse("QUATERNIONS",quatstr);
167 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_PRODUCT_MATRIX MASK=" + getShortcutLabel() + "_cmat ARG=" + quatstr + ".*," + quatstr + ".*" );
168 : } else {
169 0 : plumed_error();
170 : }
171 : //
172 2 : readInputLine( getShortcutLabel() + "_kfunc: CUSTOM MASK=" + getShortcutLabel() + "_cmat ARG=" + getShortcutLabel() + "_cmat,"+ getShortcutLabel() + "_quatprod.* " + "VAR=x,w,i,j,k PERIODIC=NO FUNC=" + kfunc );
173 : // Element wise product of cmat and kfunc
174 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
175 : // Find the number of ones we need to multiply by
176 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
177 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
178 : std::string size;
179 1 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
180 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
181 : //
182 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
183 2 : }
184 :
185 : }
186 : }
187 :
188 :
189 :
|
