Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "MultiColvarShortcuts.h"
25 :
26 : //+PLUMEDOC MCOLVAR COORD_ANGLES
27 : /*
28 : Calculate functions of the distribution of angles between bonds in the first coordination spheres of a set of atoms
29 :
30 : This action ncan be used to calculate functions, $g$, such as:
31 :
32 : $$
33 : f(x) = \sum_{ijk} s(r_{ij})s(r_{jk}) g(\theta_{ijk})
34 : $$
35 :
36 : where $\theta_{ijk}$ is the angle between the vector connecting atom $i$ and $j$ and the vector connecting atom $j$ and atom $k$ and
37 : where $s(r)$ is a switching function. The switching functions in the expression above ensure that we can calculate all the angles in the first coordination
38 : spheres of an atom using an input like the one shown below:
39 :
40 : ```plumed
41 : c1: COORD_ANGLES ...
42 : CATOMS=1 GROUP=2-100 SWITCH={RATIONAL R_0=1.0} SUM
43 : ...
44 : PRINT ARG=c1.sum FILE=colvar
45 : ```
46 :
47 : The input above will output the sum of all the angles in the first coordination sphere.
48 :
49 : __As you can see if you expand the shortcut in the input above, the calculation of the above sum is obtained by calculating a [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices))
50 : of an [OUTER_PRODUCT](OUTER_PRODUCT.md) matrix and a [MATRIX_PRODUCT](MATRIX_PRODUCT.md). We have thus deprecated this action as we believe that there is value in learning to use the more complicated syntax
51 : directly.__
52 :
53 : */
54 : //+ENDPLUMEDOC
55 :
56 : namespace PLMD {
57 : namespace multicolvar {
58 :
59 : class CoordAngles : public ActionShortcut {
60 : public:
61 : static void registerKeywords(Keywords& keys);
62 : static void pruneShortcuts(Keywords& keys);
63 : explicit CoordAngles(const ActionOptions&);
64 : private:
65 : /// Create the actions that calculate the angles between the bonds that connect catom to the atoms in group
66 : void createAngles( const std::string& lab, const std::string& catom, const std::vector<std::string>& group, const std::string& switch_input );
67 : };
68 :
69 : PLUMED_REGISTER_ACTION(CoordAngles,"COORD_ANGLES")
70 :
71 18 : void CoordAngles::registerKeywords(Keywords& keys) {
72 18 : ActionShortcut::registerKeywords( keys );
73 18 : keys.add("atoms","CATOMS","all the angles between the bonds that radiate out from these central atom are computed");
74 18 : keys.add("atoms","GROUP","a list of angls between pairs of bonds connecting one of the atoms specified using the CATOM command and two of the atoms specified here are computed");
75 18 : keys.add("compulsory","SWITCH","the switching function specifies that only those bonds that have a length that is less than a certain threshold are considered");
76 18 : MultiColvarShortcuts::shortcutKeywords( keys );
77 18 : pruneShortcuts( keys );
78 18 : keys.needsAction("DISTANCE");
79 18 : keys.needsAction("COMBINE");
80 18 : keys.needsAction("CUSTOM");
81 18 : keys.needsAction("VSTACK");
82 18 : keys.needsAction("TRANSPOSE");
83 18 : keys.needsAction("OUTER_PRODUCT");
84 18 : keys.needsAction("MATRIX_PRODUCT");
85 18 : keys.setDeprecated("MATRIX_PRODUCT");
86 18 : }
87 :
88 29 : void CoordAngles::pruneShortcuts(Keywords& keys) {
89 29 : keys.remove("ALT_MIN");
90 29 : keys.remove("MIN");
91 29 : keys.remove("MAX");
92 29 : keys.remove("HIGHEST");
93 29 : keys.remove("LOWEST");
94 29 : }
95 :
96 15 : void CoordAngles::createAngles( const std::string& lab, const std::string& catom, const std::vector<std::string>& group, const std::string& switch_input ) {
97 : // Create the list of atoms. An atom is not bonded to itself so the central atom is not used as one of its own neighbours
98 : std::string atlist;
99 15 : unsigned k=1;
100 295 : for(unsigned j=0; j<group.size(); ++j) {
101 280 : if( group[j]==catom ) {
102 4 : continue;
103 : }
104 : std::string num;
105 276 : Tools::convert( k, num );
106 828 : atlist += " ATOMS" + num + "=" + catom + "," + group[j];
107 276 : k++;
108 : }
109 15 : if( k==1 ) {
110 1 : error("there are no atoms in GROUP that are different from the central atom " + catom );
111 : }
112 : // Calculate the distances
113 14 : readInputLine( lab + "_dd: DISTANCE" + atlist );
114 : // Transform with the switching function
115 42 : readInputLine( lab + "_sw: LESS_THAN ARG=" + lab + "_dd SWITCH={" + switch_input +"}");
116 : // Now get the normalised vectors
117 28 : readInputLine( lab + "_comp: DISTANCE" + atlist + " COMPONENTS");
118 56 : readInputLine( lab + "_norm2: COMBINE ARG=" + lab + "_comp.x" + "," + lab + "_comp.y," + lab + "_comp.z POWERS=2,2,2 PERIODIC=NO");
119 28 : readInputLine( lab + "_norm: CUSTOM ARG=" + lab + "_norm2 FUNC=sqrt(x) PERIODIC=NO");
120 42 : readInputLine( lab + "_norm_x: CUSTOM ARG=" + lab + "_comp.x," + lab + "_norm FUNC=x/y PERIODIC=NO");
121 42 : readInputLine( lab + "_norm_y: CUSTOM ARG=" + lab + "_comp.y," + lab + "_norm FUNC=x/y PERIODIC=NO");
122 42 : readInputLine( lab + "_norm_z: CUSTOM ARG=" + lab + "_comp.z," + lab + "_norm FUNC=x/y PERIODIC=NO");
123 56 : readInputLine( lab + "_stack: VSTACK ARG=" + lab + "_norm_x" + "," + lab + "_norm_y," + lab + "_norm_z");
124 28 : readInputLine( lab + "_stackT: TRANSPOSE ARG=" + lab + "_stack");
125 : // Create the matrix of weights
126 42 : readInputLine( lab + "_swd: OUTER_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=" + lab + "_sw," + lab + "_sw");
127 : // Avoid double counting
128 28 : readInputLine( lab + "_wmat: CUSTOM ARG=" + lab + "_swd FUNC=0.5*x PERIODIC=NO");
129 : // And the matrix of dot products and the angles
130 42 : readInputLine( lab + "_dpmat: MATRIX_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=" + lab + "_stack," + lab + "_stackT");
131 43 : readInputLine( lab + "_angles: CUSTOM ARG=" + lab + "_dpmat FUNC=acos(x) PERIODIC=NO");
132 14 : }
133 :
134 11 : CoordAngles::CoordAngles(const ActionOptions& ao):
135 : Action(ao),
136 11 : ActionShortcut(ao) {
137 : // Parse the central atoms
138 11 : std::vector<std::string> catoms;
139 11 : parseVector("CATOMS",catoms);
140 11 : Tools::interpretRanges(catoms);
141 11 : if( catoms.size()==0 ) {
142 0 : error("no central atoms were found in CATOMS");
143 : }
144 : // Parse the coordination sphere
145 11 : std::vector<std::string> group;
146 11 : parseVector("GROUP",group);
147 11 : Tools::interpretRanges(group);
148 11 : if( group.size()==0 ) {
149 0 : error("no atoms were found in GROUP");
150 : }
151 : std::string switch_input;
152 12 : parse("SWITCH",switch_input);
153 : // Read the input
154 11 : Keywords keys;
155 11 : MultiColvarShortcuts::shortcutKeywords( keys );
156 11 : pruneShortcuts( keys );
157 : bool do_mean;
158 22 : parseFlag("MEAN",do_mean);
159 : std::map<std::string,std::string> keymap;
160 11 : readShortcutKeywords( keys, keymap );
161 11 : if( do_mean ) {
162 11 : keymap.insert(std::pair<std::string,std::string>("SUM",""));
163 : }
164 11 : const std::string lab = getShortcutLabel();
165 11 : if( catoms.size()==1 ) {
166 8 : createAngles( lab, catoms[0], group, switch_input );
167 14 : MultiColvarShortcuts::expandFunctions( lab, lab + "_angles", lab + "_wmat", keymap, this );
168 7 : if( do_mean ) {
169 8 : readInputLine( lab + "_denom: SUM ARG=" + lab + "_wmat PERIODIC=NO");
170 17 : readInputLine( lab + "_mean: CUSTOM ARG=" + lab + "_sum," + lab + "_denom FUNC=x/y PERIODIC=NO");
171 : }
172 : return;
173 : }
174 : // The angles are calculated separately for each of the central atoms, as the angles between bonds
175 : // that start from different atoms are not wanted. Everything that we output is a sum over the angles
176 : // so the final quantities are the sums of the quantities for each central atom.
177 3 : std::vector<std::string> outputs, prefixes;
178 10 : for(unsigned i=0; i<catoms.size(); ++i) {
179 : std::string num;
180 7 : Tools::convert( i+1, num );
181 14 : prefixes.push_back( lab + "_c" + num );
182 7 : createAngles( prefixes[i], catoms[i], group, switch_input );
183 7 : std::vector<std::string> thisoutputs;
184 14 : MultiColvarShortcuts::expandFunctions( prefixes[i], prefixes[i] + "_angles", prefixes[i] + "_wmat", keymap, this, &thisoutputs );
185 7 : if( i==0 ) {
186 3 : outputs = thisoutputs;
187 : }
188 7 : if( do_mean ) {
189 6 : readInputLine( prefixes[i] + "_denom: SUM ARG=" + prefixes[i] + "_wmat PERIODIC=NO");
190 : }
191 7 : }
192 10 : for(const auto& suffix : outputs ) {
193 7 : std::string args = prefixes[0] + suffix;
194 15 : for(unsigned i=1; i<prefixes.size(); ++i) {
195 16 : args += "," + prefixes[i] + suffix;
196 : }
197 28 : readInputLine( lab + suffix + ": COMBINE ARG=" + args + " PERIODIC=NO");
198 : }
199 3 : if( do_mean ) {
200 1 : std::string args = prefixes[0] + "_denom";
201 2 : for(unsigned i=1; i<prefixes.size(); ++i) {
202 2 : args += "," + prefixes[i] + "_denom";
203 : }
204 2 : readInputLine( lab + "_denom: COMBINE ARG=" + args + " PERIODIC=NO");
205 4 : readInputLine( lab + "_mean: CUSTOM ARG=" + lab + "_sum," + lab + "_denom FUNC=x/y PERIODIC=NO");
206 : }
207 28 : }
208 :
209 : }
210 : }
|
