Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2017 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 "function/FunctionShortcut.h"
23 : #include "function/FunctionOfScalar.h"
24 : #include "function/FunctionOfVector.h"
25 : #include "core/ActionRegister.h"
26 : #include "function/FunctionSetup.h"
27 :
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace refdist {
32 :
33 : //+PLUMEDOC FUNCTION DIFFERENCE
34 : /*
35 : Calculate the differences between two scalars or vectors
36 :
37 : This action can be used to calculate the difference between two values. For example, if we want to calculate the difference
38 : between the distances between atoms 1 and 2 and 3 and 4 we can use an input like the one shown below:
39 :
40 : ```plumed
41 : d1: DISTANCE ATOMS=1,2
42 : d2: DISTANCE ATOMS=3,4
43 : diff: DIFFERENCE ARG=d1,d2
44 : PRINT ARG=diff FILE=colvar
45 : ```
46 :
47 : At first sight this action appears pointless as you can achieve the same result with the following input:
48 :
49 : ```plumed
50 : d1: DISTANCE ATOMS=1,2
51 : d2: DISTANCE ATOMS=3,4
52 : diff: CUSTOM ARG=d1,d2 FUNC=x-y PERIODIC=NO
53 : PRINT ARG=diff FILE=colvar
54 : ```
55 :
56 : However, this second example will not give an equivalent result to the following input:
57 :
58 : ```plumed
59 : t1: TORSION ATOMS=1,2,3,4
60 : t2: TORSION ATOMS=5,6,7,8
61 : diff: DIFFERENCE ARG=t1,t2
62 : PRINT ARG=diff FILE=colvar
63 : ```
64 :
65 : as the [CUSTOM](CUSTOM.md) command cannot caclulate differences between periodic input variables. In this case you have to
66 : use DIFFERENCE to ensure that the periodic boundary conditions are applied when the final difference is calculated.
67 :
68 : ## Difference from reference value/s
69 :
70 : The DIFFERENCE command is frequently used to determine the difference between the instantaneous value of a quantity and
71 : a reference value as illustrated in the example below:
72 :
73 : ```plumed
74 : ref_psi: CONSTANT VALUES=2.25029540
75 : t: TORSION ATOMS=1,2,3,4
76 : diff: DIFFERENCE ARG=t,ref_psi
77 : PRINT ARG=diff FILE=colvar
78 : ```
79 :
80 : You can also use this action to calculate the difference between the instantaneous values of a vector of quantities and a
81 : reference vector as illustrated in the example below:
82 :
83 : ```plumed
84 : ref: CONSTANT VALUES=1.387980,1.126120,1.269380,1.321120,1.212420
85 : d: DISTANCES ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8 ATOMS5=9,10
86 : diff: DIFFERENCE ARG=d,ref
87 : PRINT ARG=diff FILE=colvar
88 : ```
89 :
90 : The output in this case is a five dimensional vector. Each element of this vector contains the difference between the
91 : corresponding elements of the two input vectors.
92 :
93 : Notice that you __cannot__ use a mixture of scalars and vectors in the input for this action. Only two values can be input
94 : to this action and these two input values __must__ have the same rank and the same number of elements.
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 : class Difference {
100 : public:
101 : bool periodic;
102 : double max_minus_min;
103 : double inv_max_minus_min;
104 : static void registerKeywords(Keywords& keys);
105 : static void read( Difference& func, ActionWithArguments* action, function::FunctionOptions& options );
106 : static void calc( const Difference& func, bool noderiv, const View<const double,helpers::dynamic_extent>& args, function::FunctionOutput& funcout );
107 : };
108 :
109 :
110 : typedef function::FunctionShortcut<Difference> DifferenceShortcut;
111 : PLUMED_REGISTER_ACTION(DifferenceShortcut,"DIFFERENCE")
112 : typedef function::FunctionOfScalar<Difference> ScalarDifference;
113 : PLUMED_REGISTER_ACTION(ScalarDifference,"DIFFERENCE_SCALAR")
114 : typedef function::FunctionOfVector<Difference> VectorDifference;
115 : PLUMED_REGISTER_ACTION(VectorDifference,"DIFFERENCE_VECTOR")
116 :
117 743 : void Difference::registerKeywords(Keywords& keys) {
118 1486 : keys.setValueDescription("scalar/vector","a function that measures the difference");
119 743 : }
120 :
121 184 : void Difference::read( Difference& func, ActionWithArguments* action, function::FunctionOptions& options ) {
122 184 : if( action->getNumberOfArguments()!=2 ) {
123 0 : action->error("should be two arguments to this action");
124 : }
125 184 : if( action->getPntrToArgument(0)->getRank()==action->getPntrToArgument(1)->getRank() ) {
126 51 : std::vector<std::size_t> shape( action->getPntrToArgument(0)->getShape() );
127 61 : for(unsigned i=0; i<shape.size(); ++i) {
128 10 : if( shape[i]!=action->getPntrToArgument(1)->getShape()[i] ) {
129 0 : action->error("shapes of input actions do not match");
130 : }
131 : }
132 51 : }
133 :
134 184 : func.periodic=false;
135 : std::string min0, max0;
136 184 : if( action->getPntrToArgument(0)->isPeriodic() ) {
137 38 : func.periodic=true;
138 38 : action->getPntrToArgument(0)->getDomain( min0, max0 );
139 38 : if( !action->getPntrToArgument(1)->isConstant() && !action->getPntrToArgument(1)->isPeriodic() ) {
140 0 : action->error("period for input variables " + action->getPntrToArgument(0)->getName() + " and " + action->getPntrToArgument(1)->getName() + " should be the same 0");
141 : }
142 38 : if( !action->getPntrToArgument(1)->isConstant() ) {
143 : std::string min1, max1;
144 0 : action->getPntrToArgument(1)->getDomain( min1, max1 );
145 0 : if( min0!=min1 || max0!=max1 ) {
146 0 : action->error("domain for input variables should be the same");
147 : }
148 : } else {
149 38 : action->getPntrToArgument(1)->setDomain( min0, max0 );
150 : }
151 146 : } else if( action->getPntrToArgument(1)->isPeriodic() ) {
152 4 : func.periodic=true;
153 4 : action->getPntrToArgument(1)->getDomain( min0, max0 );
154 4 : if( !action->getPntrToArgument(0)->isConstant() ) {
155 0 : action->error("period for input variables " + action->getPntrToArgument(0)->getName() + " and " + action->getPntrToArgument(1)->getName() + " should be the same 1");
156 : } else {
157 4 : action->getPntrToArgument(0)->setDomain( min0, max0 );
158 : }
159 : }
160 184 : if( func.periodic ) {
161 : double min, max;
162 42 : Tools::convert( min0, min );
163 42 : Tools::convert( max0, max );
164 42 : func.max_minus_min=max-min;
165 42 : func.inv_max_minus_min=1.0/func.max_minus_min;
166 : }
167 184 : }
168 :
169 3916895 : void Difference::calc( const Difference& func, bool noderiv, const View<const double,helpers::dynamic_extent>& args, function::FunctionOutput& funcout ) {
170 : plumed_dbg_assert( args.size()==2 );
171 3916895 : if( func.periodic ) {
172 262908 : funcout.values[0] = func.max_minus_min*Tools::pbc( func.inv_max_minus_min*(args[0] - args[1]) );
173 : } else {
174 3653987 : funcout.values[0] = args[0] - args[1];
175 : }
176 3916895 : if( !noderiv ) {
177 1921829 : funcout.derivs[0][0] = 1.0;
178 1921829 : funcout.derivs[0][1]=-1;
179 : }
180 3916895 : }
181 :
182 : }
183 : }
184 :
185 :
|
