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 DOPS
31 : /*
32 : Calculate DOPS order parameter for molecules.
33 :
34 : DOPS is a shortcut to calculate the Distance Order Parameters that are described in the paper below.
35 : As arguments, DOPS takes a list of atoms, a cutoff, and a kernel file detailing the means, variances, and normalization factors of probability distributions. DOPS returns a vector of order parameters.
36 :
37 : 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.
38 :
39 : DOPS returns one order parameter per atom given, evaluated over each atom's neighbors within the cutoff given. The kernel file defines a radial distribution function. The order parameter is obtained by evaluating the RDF at each neighbor's position within the cutoff, and summing them up.
40 :
41 : This example file calculates the DOPS for a system of 3 molecules.
42 :
43 : ```plumed
44 : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-dops-forces/kernels.dat
45 : dops: DOPS SPECIES=1,4,7 CUTOFF=7.4 KERNELFILE=regtest/crystdistrib/rt-dops-forces/kernels.dat
46 : ```
47 :
48 : If you want to calculate DOPS from the radial distribution function of the GROUPB atoms around the GROUPA atoms you use an input like the one shown below:
49 :
50 : ```plumed
51 : #SETTINGS INPUTFILES=regtest/crystdistrib/rt-dops-forces/kernels.dat
52 : dops: DOPS ...
53 : SPECIESA=1,4,7 SPECIESB=10,13,16,19 CUTOFF=7.4
54 : KERNELFILE=regtest/crystdistrib/rt-dops-forces/kernels.dat
55 : ...
56 : ```
57 :
58 : */
59 : //+ENDPLUMEDOC
60 :
61 :
62 : class DopsShortcut : public ActionShortcut {
63 : public:
64 : static void registerKeywords( Keywords& keys );
65 : explicit DopsShortcut(const ActionOptions&);
66 : };
67 :
68 : PLUMED_REGISTER_ACTION(DopsShortcut,"DOPS")
69 :
70 4 : void DopsShortcut::registerKeywords( Keywords& keys ) {
71 4 : ActionShortcut::registerKeywords( keys );
72 4 : keys.add("atoms","SPECIES","the list of atoms for which the DOPS are being calculated and the atoms that can be in the environments");
73 4 : keys.add("atoms-2","SPECIESA","the list of atoms for which DOPS are being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment");
74 4 : keys.add("atoms-2","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the DOPS are being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which DOPS are being calculated.");
75 4 : keys.add("compulsory","KERNELFILE","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)");
76 4 : keys.add("compulsory","CUTOFF","6.25","to make the calculation faster we calculate a cutoff value on the distances. The input to this keyword determines x in this expreession max(mu + sqrt(2*x)/sigma)");
77 8 : keys.setValueDescription("vector","the values of the DOPS order parameters");
78 4 : keys.needsAction("DISTANCE_MATRIX");
79 4 : keys.needsAction("CUSTOM");
80 4 : keys.needsAction("ONES");
81 4 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
82 4 : keys.addDOI("10.1063/1.3548889");
83 4 : }
84 :
85 2 : DopsShortcut::DopsShortcut(const ActionOptions&ao):
86 : Action(ao),
87 2 : ActionShortcut(ao) {
88 : // Open a file and read in the kernels
89 2 : double cutoff=0, h;
90 : std::string kfunc,fname;
91 : double dp2cutoff;
92 2 : parse("CUTOFF",dp2cutoff);
93 2 : parse("KERNELFILE",fname);
94 2 : IFile ifile;
95 2 : ifile.open(fname);
96 20 : for(unsigned k=0;; ++k) {
97 44 : if( !ifile.scanField("height",h) ) {
98 : break;
99 : }
100 : std::string ktype;
101 40 : ifile.scanField("kerneltype",ktype);
102 20 : if( ktype!="gaussian" ) {
103 0 : error("cannot process kernels of type " + ktype );
104 : }
105 : double mu, sigma;
106 20 : ifile.scanField("mu",mu);
107 20 : ifile.scanField("sigma",sigma);
108 20 : ifile.scanField();
109 : std::string hstr, mustr, sigmastr;
110 20 : Tools::convert( h, hstr );
111 20 : Tools::convert( 2*sigma*sigma, sigmastr );
112 20 : Tools::convert( mu, mustr );
113 : // Get a sensible value for the cutoff
114 20 : double support = sqrt(2.0*dp2cutoff)*(1.0/sigma);
115 20 : if( mu+support>cutoff ) {
116 6 : cutoff= mu + support;
117 : }
118 : // And make the kernel
119 20 : if( k==0 ) {
120 : kfunc = hstr;
121 : } else {
122 36 : kfunc += "+" + hstr;
123 : }
124 40 : kfunc += "*exp(-(x-" + mustr +")^2/" + sigmastr + ")";
125 20 : }
126 : std::string sp_str, specA, specB, grpinfo;
127 2 : parse("SPECIES",sp_str);
128 2 : parse("SPECIESA",specA);
129 4 : parse("SPECIESB",specB);
130 2 : if( sp_str.length()>0 ) {
131 4 : grpinfo="GROUP=" + sp_str;
132 : } else {
133 0 : if( specA.length()==0 || specB.length()==0 ) {
134 0 : error("no atoms were specified in input use either SPECIES or SPECIESA + SPECIESB");
135 : }
136 0 : grpinfo="GROUPA=" + specA + " GROUPB=" + specB;
137 : }
138 : std::string cutstr;
139 2 : Tools::convert( cutoff, cutstr );
140 : // Setup the contact matrix
141 4 : readInputLine( getShortcutLabel() + "_cmat: DISTANCE_MATRIX " + grpinfo + " CUTOFF=" + cutstr);
142 : // And the kernels
143 4 : readInputLine( getShortcutLabel() + "_kval: CUSTOM ARG=" + getShortcutLabel() + "_cmat PERIODIC=NO FUNC=" + kfunc );
144 : // Find the number of ones we need to multiply by
145 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmat");
146 2 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
147 : std::string size;
148 2 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
149 4 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
150 : // And the final order parameters
151 4 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_kval," + getShortcutLabel() + "_ones");
152 4 : }
153 :
154 : }
155 : }
156 :
157 :
158 :
|
