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 6 : keys.setValueDescription("vector","the values of the bops order parameters");
67 9 : keys.needsAction("DISTANCE_MATRIX"); keys.needsAction("QUATERNION_BOND_PRODUCT_MATRIX"); keys.needsAction("CUSTOM");
68 6 : keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
69 3 : }
70 :
71 1 : BopsShortcut::BopsShortcut(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_bops; std::string kfunc, kfunc_dops,kfunc_bops,fname_dops,fname_bops;
77 3 : parse("KERNELFILE_DOPS",fname_dops); parse("KERNELFILE_BOPS",fname_bops); IFile ifile_dops, ifile_bops; ifile_dops.open(fname_dops); ifile_bops.open(fname_bops);
78 10 : for(unsigned k=0;; ++k) {
79 21 : if( !ifile_dops.scanField("height",h_dops) || !ifile_bops.scanField("height",h_bops) ) break;//checks eof
80 30 : std::string ktype_dops, ktype_bops; ifile_dops.scanField("kerneltype",ktype_dops); ifile_bops.scanField("kerneltype",ktype_bops);
81 10 : if( ktype_dops!="gaussian" ) error("cannot process kernels of type " + ktype_dops );//straightup error
82 10 : if( ktype_bops!="gaussian" ) error("cannot process kernels of type " + ktype_bops );
83 :
84 : double mu_dops, mu_i, mu_j, mu_k; std::string hstr_dops, hstr_bops, smu_dops,smu_i, smu_j, smu_k, sigmastr,kappastr;
85 :
86 :
87 10 : Tools::convert( h_dops, hstr_dops );
88 10 : Tools::convert( h_bops, hstr_bops );
89 :
90 20 : ifile_dops.scanField("mu",mu_dops); Tools::convert( mu_dops, smu_dops );
91 : //ifile_bops.scanField("mu_w",mu_w); Tools::convert( mu_w, smu_w );
92 20 : ifile_bops.scanField("mu_i",mu_i); Tools::convert( mu_i, smu_i );
93 20 : ifile_bops.scanField("mu_j",mu_j); Tools::convert( mu_j, smu_j );
94 20 : ifile_bops.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_bops.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_bops.scanField(); /*if( k==0 )*/ kfunc_bops = hstr_bops; //else kfunc_bops += "+" + hstr;
105 :
106 20 : kfunc_bops += "*exp(" + kappastr + "*(i*" + smu_i + "+j*" + smu_j + "+k*" + smu_k + "))";
107 20 : kfunc_dops += "*exp(-(x-" + smu_dops +")^2/" + "(2*" + sigmastr +"*" +sigmastr + "))";
108 20 : if (k==0) kfunc = kfunc_dops + "*" + kfunc_bops; else kfunc+= "+" + kfunc_dops + "*" + kfunc_bops;
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 + " COMPONENTS");
123 :
124 1 : if( specA.length()==0 ) {
125 1 : std::string quatstr; parse("QUATERNIONS",quatstr);
126 2 : readInputLine( getShortcutLabel() + "_quatprod: QUATERNION_BOND_PRODUCT_MATRIX ARG=" + quatstr + ".*," + getShortcutLabel() + "_cmat.*" );
127 : } else {
128 0 : plumed_error();
129 : }
130 : //
131 :
132 : ///////////////////
133 : ///replace/////
134 2 : readInputLine( getShortcutLabel() + "_dist: CUSTOM ARG=" + getShortcutLabel() + "_cmat.x," + getShortcutLabel() + "_cmat.y," + getShortcutLabel() + "_cmat.z " +
135 : "FUNC=sqrt((x^2)+(y^2)+(z^2)) PERIODIC=NO");
136 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");
137 :
138 : //replace ^^^ to remove distance hack
139 : //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");
140 : ///end replace////
141 :
142 : // Element wise product of cmat and kfunc
143 : // readInputLine( getShortcutLabel() + "_kdmat: CUSTOM ARG=" + getShortcutLabel() + "_cmat.w," + getShortcutLabel() + "_kfunc FUNC=x*y PERIODIC=NO");
144 : // Find the number of ones we need to multiply by
145 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
146 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
147 1 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
148 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
149 : //
150 2 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kfunc," + getShortcutLabel() + "_ones");
151 2 : }
152 :
153 : }
154 : }
155 :
156 :
157 :
|
