Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "tools/HistogramBead.h"
24 : #include "MultiColvarFilter.h"
25 :
26 : //+PLUMEDOC MTRANSFORMS MTRANSFORM_BETWEEN
27 : /*
28 : This action can be useed to transform the colvar values calculated by a multicolvar using a \ref histogrambead
29 :
30 : In this action each colvar, \f$s_i\f$, calculated by multicolvar is transformed by a \ref histogrambead function that
31 : is equal to one if the colvar is within a certain range and which is equal to zero otherwise. In other words, we
32 : compute:
33 :
34 : \f[
35 : f_i = \int_a^b K\left( \frac{s-s_i}{w} \right)
36 : \f]
37 :
38 : where \f$a, b\f$ and \f$w\f$ are parameters.
39 :
40 : It is important to understand the distinction between what is done here and what is done by \ref MFILTER_BETWEEN.
41 : In \ref MFILTER_BETWEEN a weight, \f$w_i\f$ for the colvar is calculated using the \ref histogrambead. If one calculates the
42 : MEAN for \ref MFILTER_BETWEEN one is thus calculating:
43 :
44 : \f[
45 : \mu = \frac{ \sum_i f_i s_i }{ \sum_i f_i}
46 : \f]
47 :
48 : In this action by contrast the colvar is being transformed by the \ref histogrambead. If one thus calculates a MEAN for
49 : thia action one computes:
50 :
51 : \f[
52 : \mu = \frac{ \sum_{i=1}^N f_i }{ N }
53 : \f]
54 :
55 : In other words, you are calculating the mean for the transformed colvar.
56 :
57 : \par Examples
58 :
59 : The following input gives an example of how a MTRANSFORM_BETWEEN action can be used to duplicate
60 : functionality that is elsehwere in PLUMED.
61 :
62 : \plumedfile
63 : DISTANCES ...
64 : GROUPA=1-10 GROUPB=11-20
65 : LABEL=d1
66 : ... DISTANCES
67 : MTRANSFORM_BETWEEN DATA=d1 LOWER=1.0 UPPER=2.0 SMEAR=0.5
68 : \endplumedfile
69 :
70 : In this case you can achieve the same result by using:
71 :
72 : \plumedfile
73 : DISTANCES ...
74 : GROUPA=1-10 GROUPB=11-20
75 : BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0}
76 : ... DISTANCES
77 : \endplumedfile
78 : (see \ref DISTANCES)
79 :
80 : The advantage of MTRANSFORM_BETWEEN comes, however, if you want to use transformed colvars as input
81 : for \ref MULTICOLVARDENS
82 :
83 : */
84 : //+ENDPLUMEDOC
85 :
86 : //+PLUMEDOC MFILTERS MFILTER_BETWEEN
87 : /*
88 : This action can be used to filter the colvar values calculated by a multicolvar
89 : so that one can compute the mean and so on for only those multicolvars within a certain range.
90 :
91 : This action can be used to create a dynamic group of atom based on the value of a multicolvar.
92 : In this action a multicolvar is within the dynamic group if its value lies in a particular range.
93 : In practise a weight, \f$w_i\f$ is ascribed to each colvar, \f$s_i\f$ calculated by a multicolvar
94 : and this weight measures the degree to which a colvar is a member of the group. This weight is
95 : calculated using a \ref histogrambead so it is given by:
96 :
97 : \f[
98 : w_i = \int_a^b K\left( \frac{s - s_i}{w} \right)
99 : \f]
100 :
101 : where \f$a, b\f$ and \f$w\f$ are parameters. If one calculates a function of the set of multicolvars
102 : these weights are included in the calculation. As such if one calculates the MEAN, \f$\mu\f$ of a filtered
103 : multicolvar what is computed is the following:
104 :
105 : \f[
106 : \mu = \frac{ \sum_i w_i s_i }{ \sum_i w_i}
107 : \f]
108 :
109 : One is thus calculating the mean for those colvars that are within the range of interest.
110 :
111 : \par Examples
112 :
113 : The example shown below calculates the mean for those distances that are between 0 and 3 nm in length
114 :
115 : \plumedfile
116 : DISTANCES GROUPA=1 GROUPB=2-50 MEAN LABEL=d1
117 : MFILTER_BETWEEN DATA=d1 LOWER=0 UPPER=3.0 SMEAR=0.0001 MEAN LABEL=d4
118 : \endplumedfile
119 :
120 : More complicated things can be done by using the label of a filter as input to a new multicolvar as shown
121 : in the example below. Here the coordination numbers of all atoms are computed. The atoms with a coordination
122 : number between 4 and 6 are then identified using the filter. This reduced list of atoms is then used as input
123 : to a second coordination number calculation. This second coordination number thus measures the number of atoms
124 : 4-6 coordinated atoms each of the 4-6 coordination atoms is bound to.
125 :
126 : \plumedfile
127 : c1: COORDINATIONNUMBER SPECIES=1-150 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
128 : cf: MFILTER_BETWEEN DATA=c1 LOWER=4 UPPER=6 SMEAR=0.5 LOWMEM
129 : c2: COORDINATIONNUMBER SPECIES=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0} MORE_THAN={RATIONAL D_0=2.0 R_0=0.1}
130 : \endplumedfile
131 :
132 : */
133 : //+ENDPLUMEDOC
134 :
135 : namespace PLMD {
136 : namespace multicolvar {
137 :
138 6 : class FilterBetween : public MultiColvarFilter {
139 : private:
140 : HistogramBead hb;
141 : public:
142 : static void registerKeywords( Keywords& keys );
143 : explicit FilterBetween(const ActionOptions& ao);
144 : double applyFilter( const double& val, double& df ) const ;
145 : };
146 :
147 6455 : PLUMED_REGISTER_ACTION(FilterBetween,"MFILTER_BETWEEN")
148 6452 : PLUMED_REGISTER_ACTION(FilterBetween,"MTRANSFORM_BETWEEN")
149 :
150 5 : void FilterBetween::registerKeywords( Keywords& keys ) {
151 5 : MultiColvarFilter::registerKeywords( keys );
152 20 : keys.add("compulsory","LOWER","the lower boundary for the range of interest");
153 20 : keys.add("compulsory","UPPER","the upper boundary for the range of interest");
154 25 : keys.add("compulsory","SMEAR","0.5","the ammount by which to smear the value for kernel density estimation");
155 20 : keys.add("optional","BEAD","This keywords is used if you want to employ an alternative to the function defeind above. "
156 : "The following provides information on the \\ref histogrambead that are available. "
157 : "When this keyword is present you no longer need the LOWER, UPPER and SMEAR keywords.");
158 5 : }
159 :
160 3 : FilterBetween::FilterBetween(const ActionOptions& ao):
161 : Action(ao),
162 3 : MultiColvarFilter(ao)
163 : {
164 : // Read in the switching function
165 6 : std::string sw, errors; parse("BEAD",sw);
166 3 : if( getPntrToMultiColvar()->isPeriodic() ) {
167 0 : std::string min, max; getPntrToMultiColvar()->retrieveDomain( min, max );
168 0 : double mlow, mhigh; Tools::convert( min,mlow ); Tools::convert( max, mhigh);
169 0 : hb.isPeriodic( mlow, mhigh );
170 : } else {
171 : hb.isNotPeriodic();
172 : }
173 :
174 3 : if(sw.length()>0) {
175 0 : hb.set(sw,errors);
176 0 : if( errors.length()!=0 ) error("problem reading BEAD keyword : " + errors );
177 : } else {
178 : double l, u, s; std::string ll, uu, ss;
179 12 : parse("LOWER",l); parse("UPPER",u); parse("SMEAR",s);
180 3 : Tools::convert(l,ll); Tools::convert(u,uu); Tools::convert(s,ss);
181 18 : sw="GAUSSIAN LOWER=" + ll + " UPPER=" + uu + " SMEAR=" + ss;
182 6 : hb.set(sw,errors); plumed_massert(errors.length()==0,"problems with bead" + errors);
183 : }
184 9 : log.printf(" filtering colvar values and focussing only on those values in range %s\n",( hb.description() ).c_str() );
185 :
186 3 : checkRead();
187 3 : }
188 :
189 735 : double FilterBetween::applyFilter( const double& val, double& df ) const {
190 735 : double f = hb.calculate( val, df );
191 735 : return f;
192 : }
193 :
194 : }
195 4839 : }
|
