Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "Function.h"
23 : #include "ActionRegister.h"
24 : #include "tools/PDB.h"
25 : #include "reference/MetricRegister.h"
26 : #include "reference/ArgumentOnlyDistance.h"
27 : #include "core/Atoms.h"
28 : #include "core/PlumedMain.h"
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace function {
34 :
35 : //+PLUMEDOC DCOLVAR TARGET
36 : /*
37 : This function measures the pythagorean distance from a particular structure measured in the space defined by some
38 : set of collective variables.
39 :
40 : This collective variable can be used to calculate something akin to:
41 :
42 : \f[
43 : d(X,X') = \vert X - X' \vert
44 : \f]
45 :
46 : where \f$ X \f$ is the instaneous values for a set of collective variables for the system and
47 : \f$ X' \f$ is the values that these self-same set of collective variables take in some reference structure provided as input.
48 : If we call our set of collective variables \f$\{s_i\}\f$ then this CV computes:
49 :
50 : \f[
51 : d = \sqrt{ \sum_{i=1}^N (s_i - s_i^{(ref)})^2 }
52 : \f]
53 :
54 : where \f$s_i^{(ref)}\f$ are the values of the CVs in the reference structure and \f$N\f$ is the number of input CVs.
55 :
56 : We can also calculate normalized euclidean differences using this action and the METRIC=NORM-EUCLIDEAN flag. In other words,
57 : we can compute:
58 :
59 : \f[
60 : d = \sqrt{ \sum_{i=1}^N \sigma_i (s_i - s_i^{(ref)})^2 }
61 : \f]
62 :
63 : where \f$\sigma_i\f$ is a vector of weights. Lastly, by using the METRIC=MAHALONOBIS we can compute mahalonobis distances using:
64 :
65 : \f[
66 : d = \left( \mathbf{s} - \mathbf{s}^{(ref)} \right)^T \mathbf{\Sigma} \left( \mathbf{s} - \mathbf{s}^{(ref)} \right)
67 : \f]
68 :
69 : where \f$\mathbf{s}\f$ is a column vector containing the values of all the CVs and \f$\mathbf{s}^{(ref)}\f$ is a column vector
70 : containg the values of the CVs in the reference configuration. \f$\mathbf{\Sigma}\f$ is then an \f$N \times N\f$ matrix that is
71 : specified in the input.
72 :
73 : \par Examples
74 :
75 : The following input calculates the distance between a reference configuration and the instaneous position of the system in the trajectory.
76 : The position of the reference configuration is specified by providing the values of the distance between atoms 1 and 2 and atoms 3 and 4.
77 :
78 : \plumedfile
79 : d1: DISTANCE ATOMS=1,2
80 : d2: DISTANCE ATOMS=3,4
81 : t1: TARGET REFERENCE=myref.pdb TYPE=EUCLIDEAN
82 : PRINT ARG=t1 FILE=colvar
83 : \endplumedfile
84 :
85 : The contents of the file containing the reference structure (myref.pdb) is shown below. As you can see you must provide information on the
86 : labels of the CVs that are being used to define the position of the reference configuration in this file together with the values that these
87 : quantities take in the reference configuration.
88 :
89 : \verbatim
90 : DESCRIPTION: a reference point.
91 : REMARK WEIGHT=1.0
92 : REMARK ARG=d1,d2
93 : REMARK d1=1.0 d2=1.0
94 : END
95 : \endverbatim
96 :
97 : */
98 : //+ENDPLUMEDOC
99 :
100 : class Target : public Function {
101 : private:
102 : MultiValue myvals;
103 : ReferenceValuePack mypack;
104 : PLMD::ArgumentOnlyDistance* target;
105 : public:
106 : explicit Target(const ActionOptions&);
107 : ~Target();
108 : virtual void calculate();
109 : static void registerKeywords(Keywords& keys );
110 : };
111 :
112 6452 : PLUMED_REGISTER_ACTION(Target,"TARGET")
113 :
114 1 : void Target::registerKeywords(Keywords& keys) {
115 1 : Function::registerKeywords(keys);
116 5 : keys.add("compulsory","TYPE","EUCLIDEAN","the manner in which the distance should be calculated");
117 4 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure. In the PDB file the atomic "
118 : "coordinates and box lengths should be in Angstroms unless you are working with natural units. "
119 : "If you are working with natural units then the coordinates should be in your natural length unit. "
120 : "The charges and masses of the atoms (if required) should be inserted in the beta and occupancy "
121 : "columns respectively. For more details on the PDB file format visit http://www.wwpdb.org/docs.html");
122 1 : }
123 :
124 0 : Target::Target(const ActionOptions&ao):
125 : Action(ao),
126 : Function(ao),
127 : myvals(1,0),
128 0 : mypack(0,0,myvals)
129 : {
130 0 : std::string type; parse("TYPE",type);
131 0 : std::string reference; parse("REFERENCE",reference);
132 0 : checkRead(); PDB pdb;
133 0 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) )
134 0 : error("missing input file " + reference);
135 :
136 : // Use the base ActionWithArguments to expand things like a1.*
137 0 : expandArgKeywordInPDB( pdb );
138 :
139 : // Generate the reference structure
140 0 : target=metricRegister().create<ArgumentOnlyDistance>( type, pdb );
141 :
142 : // Get the argument names
143 0 : std::vector<std::string> args_to_retrieve;
144 0 : target->getArgumentRequests( args_to_retrieve, false );
145 :
146 : // Get the arguments
147 : std::vector<Value*> myargs;
148 0 : interpretArgumentList( args_to_retrieve, myargs );
149 0 : requestArguments( myargs );
150 :
151 : // Now create packs
152 0 : myvals.resize( 1, myargs.size() );
153 0 : mypack.resize( myargs.size(), 0 );
154 :
155 : // Create the value
156 0 : addValueWithDerivatives(); setNotPeriodic();
157 0 : }
158 :
159 0 : Target::~Target() {
160 0 : delete target;
161 0 : }
162 :
163 0 : void Target::calculate() {
164 0 : mypack.clear(); double r=target->calculate( getArguments(), mypack, false ); setValue(r);
165 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) setDerivative( i, mypack.getArgumentDerivative(i) );
166 0 : }
167 :
168 : }
169 4839 : }
|
