Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 : See http://www.plumed.org for more information.
5 : This file is part of plumed, version 2.
6 : plumed is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 : plumed is distributed in the hope that it will be useful,
11 : but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : GNU Lesser General Public License for more details.
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include "core/ActionRegister.h"
18 : #include "core/PlumedMain.h"
19 : #include "core/ActionSet.h"
20 : #include "core/ActionWithValue.h"
21 : #include "core/ActionShortcut.h"
22 :
23 : namespace PLMD {
24 : namespace colvar {
25 :
26 : //+PLUMEDOC MCOLVAR GYRATION_TENSOR
27 : /*
28 : Calculate the gyration tensor using a user specified vector of weights
29 :
30 : The elements of the $3 \times 3$ gyration tensor are defined as:
31 :
32 : $$
33 : G_{\alpha\beta} = \frac{\sum_i^{n} w_i (\alpha_i-\alpha_{\rm COM})(\beta_i - \beta_{\rm COM})}{\sum_i^{n} w_i}
34 : $$
35 :
36 : where $\alpha_i$ and $\beta_i$ can be the $x$, $y$ or $z$ coordinates of atom $i$ and $\alpha_{\rm COM}$ and
37 : $\beta_{\rm COM}$ can be the $x$, $y$ or $z$ components of the center, which is calculated using:
38 :
39 : $$
40 : {r}_{\rm COM}=\frac{\sum_i^{n} {r}_i\ w_i }{\sum_i^{n} w_i}
41 : $$
42 :
43 : The following example shows how you can calculate and print the gyration tensor from the positions of the fist 10 atoms using PLUMED
44 :
45 : ```plumed
46 : g: GYRATION_TENSOR ATOMS=1-10
47 : PRINT ARG=g FILE=colvar
48 : ```
49 :
50 : In this example input the weights of all the atoms are set equal to one. The 9 elements of the resulting gyration matrix will
51 : be output to the file `colvar` here. If you want the weights of the atoms to be set equal to the masses you can use either of the following
52 : equivalent commands:
53 :
54 : ```plumed
55 : g: GYRATION_TENSOR ATOMS=1-10 MASS
56 : g2: GYRATION_TENSOR ATOMS=1-10 WEIGHTS=@Masses
57 : g3: GYRATION_TENSOR ATOMS=1-10 MASS_WEIGHTED
58 : PRINT ARG=g,g2,g3 FILE=colvar
59 : ```
60 :
61 : This input above should output three identical $3\times 3$ matrices.
62 :
63 : If you want to use an arbitrary set of weights for the atoms you can use the following syntax.
64 :
65 : ```plumed
66 : c: CONSTANT VALUES=1,3,2,4
67 : g: GYRATION_TENSOR ATOMS=1-4 WEIGHTS=c UNORMALIZED
68 : PRINT ARG=g FILE=colvar
69 : ```
70 :
71 : Although the input weights are [CONSTANT](CONSTANT.md) here that is not a requirement. You can use a vector of weights from any action
72 : that outputs a vector in the input for the WEIGHTS keyword here. Notice, also that the denominator in the expression for $G_{\alpha\beta}$
73 : is set equal to one rather than the sum of the weights as we used the UNORMALIZED flag.
74 :
75 : Similar functionality to the functionality in the examples above is used in the [GYRATION](GYRATION.md) shortcut. There is, however,
76 : no fast version of the GYRATION_TENSOR command in the way that there is a fast version of the [GYRATION](GYRATION.md) command that is
77 : used when the weights are all one or when the masses are used as the weights.
78 :
79 : ## A note on periodic boundary conditions
80 :
81 : Calculating the gyration tensor is normally used to determine the shape of a molecule so all the specified atoms
82 : would normally be part of the same molecule. When computing gyration tensors it is important to ensure that the periodic boundaries
83 : are calculated correctly. There are two ways that you can manage periodic boundary conditions when using this action. The
84 : first and simplest is to reconstruct the molecule similarly to the way that [WHOLEMOLECULES](WHOLEMOLECULES.md) operates.
85 : This reconstruction of molecules has been done automatically since PLUMED 2.2. If for some reason you want to turn it off
86 : you can use the NOPBC flag as shown below:
87 :
88 : ```plumed
89 : g: GYRATION_TENSOR ATOMS=1-5 NOPBC
90 : PRINT ARG=g FILE=colvar
91 : ```
92 :
93 : An alternative approach to handling PBC is to use the PHASES keyword. This keyword instructs PLUMED to use the PHASES option
94 : when computing the position of the center using the [CENTER](CENTER.md) command. Distances of atoms from this center are then
95 : computed using PBC as usual. The example shown below shows you how to use this option
96 :
97 : ```plumed
98 : g: GYRATION_TENSOR ATOMS=1-5 PHASES
99 : PRINT ARG=g FILE=colvar
100 : ```
101 :
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : class GyrationShortcut : public ActionShortcut {
107 : public:
108 : static void registerKeywords( Keywords& keys );
109 : explicit GyrationShortcut(const ActionOptions&);
110 : };
111 :
112 : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION")
113 : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION_TENSOR")
114 :
115 156 : void GyrationShortcut::registerKeywords( Keywords& keys ) {
116 156 : ActionShortcut::registerKeywords( keys );
117 156 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
118 156 : if( keys.getDisplayName()=="GYRATION" ) {
119 148 : keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
120 : }
121 156 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
122 156 : keys.add("optional","WEIGHTS","what weights should be used when calculating the center. If this keyword is not present the geometric center is computed. "
123 : "If WEIGHTS=@Masses is used the center of mass is computed. If WEIGHTS=@charges the center of charge is computed. If "
124 : "the label of an action is provided PLUMED assumes that that action calculates a list of symmetry functions that can be used "
125 : "as weights. Lastly, an explicit list of numbers to use as weights can be provided");
126 156 : keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center of mass");
127 156 : keys.addFlag("MASS",false,"calculate the center of mass");
128 156 : keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
129 156 : keys.addFlag("UNORMALIZED",false,"do not divide by the sum of the weights");
130 156 : if( keys.getDisplayName()=="GYRATION" ) {
131 296 : keys.setValueDescription("scalar","the radius that was computed from the weights");
132 148 : keys.addActionNameSuffix("_SCALAR");
133 296 : keys.addActionNameSuffix("_VECTOR");
134 8 : } else if( keys.getDisplayName()=="GYRATION_TENSOR" ) {
135 16 : keys.setValueDescription("matrix","the gyration tensor that was computed from the weights");
136 : }
137 156 : keys.needsAction("CENTER");
138 156 : keys.needsAction("CONSTANT");
139 156 : keys.needsAction("ONES");
140 156 : keys.needsAction("MASS");
141 156 : keys.needsAction("DISTANCE");
142 156 : keys.needsAction("COVARIANCE_MATRIX");
143 156 : keys.needsAction("SELECT_COMPONENTS");
144 156 : keys.needsAction("SUM");
145 156 : keys.needsAction("CUSTOM");
146 156 : keys.needsAction("DIAGONALIZE");
147 156 : keys.addDOI("10.1021/jp2065612");
148 156 : }
149 :
150 146 : GyrationShortcut::GyrationShortcut(const ActionOptions& ao):
151 : Action(ao),
152 146 : ActionShortcut(ao) {
153 : bool usemass, phases, massw;
154 146 : parseFlag("MASS",usemass);
155 146 : parseFlag("MASS_WEIGHTED",massw);
156 146 : if( massw ) {
157 22 : usemass = true;
158 : }
159 146 : parseFlag("PHASES",phases);
160 146 : std::vector<std::string> str_weights;
161 292 : parseVector("WEIGHTS",str_weights);
162 : std::string wflab;
163 146 : if( !phases && getName()=="GYRATION" ) {
164 142 : if( usemass || str_weights.size()==0 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
165 : std::string wt_str;
166 142 : if( str_weights.size()>0 ) {
167 0 : wt_str="WEIGHTS=" + str_weights[0];
168 0 : for(unsigned i=1; i<str_weights.size(); ++i) {
169 2 : wt_str += "," + str_weights[i];
170 : }
171 : }
172 142 : if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
173 : wt_str = "MASS_WEIGHTED";
174 : }
175 : std::string inpt;
176 284 : parseNumbered( "ATOMS", 1, inpt );
177 142 : if( inpt.length()>0 ) {
178 4 : readInputLine( getShortcutLabel() + ": GYRATION_VECTOR " + wt_str + " ATOMS1=" + inpt + " " + convertInputLineToString() );
179 : } else {
180 425 : readInputLine( getShortcutLabel() + ": GYRATION_SCALAR " + wt_str + " " + convertInputLineToString() );
181 : }
182 : return;
183 : }
184 : }
185 4 : if( usemass ) {
186 0 : str_weights.resize(1);
187 : str_weights[0]="@Masses";
188 : }
189 6 : log<<" Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)")<<"\n";
190 : // Read in the atoms involved
191 4 : std::vector<std::string> atoms;
192 4 : parseVector("ATOMS",atoms);
193 4 : Tools::interpretRanges(atoms);
194 4 : if( atoms.empty() ) {
195 0 : plumed_merror("no atoms were specified, use the ATOMS keyword");
196 : }
197 4 : std::string gtype, atlist=atoms[0];
198 12 : for(unsigned i=1; i<atoms.size(); ++i) {
199 16 : atlist += "," + atoms[i];
200 : }
201 : bool nopbc;
202 8 : parseFlag("NOPBC",nopbc);
203 : std::string pbcstr;
204 4 : if(nopbc) {
205 : pbcstr = " NOPBC";
206 : }
207 : std::string phasestr;
208 4 : if(phases) {
209 : phasestr = " PHASES";
210 : }
211 : // Create the geometric center of the molecule
212 4 : std::string weights_str="";
213 4 : if( str_weights.size()>0 ) {
214 4 : weights_str=" WEIGHTS=" + str_weights[0];
215 8 : for(unsigned i=1; i<str_weights.size(); ++i) {
216 8 : weights_str += "," + str_weights[i];
217 : }
218 : }
219 8 : readInputLine( getShortcutLabel() + "_cent: CENTER ATOMS=" + atlist + pbcstr + phasestr + weights_str );
220 4 : if( str_weights.size()==0 ) {
221 0 : wflab = getShortcutLabel() + "_w";
222 : std::string str_natoms;
223 0 : Tools::convert( atoms.size(), str_natoms );
224 0 : readInputLine( getShortcutLabel() + "_w: ONES SIZE=" + str_natoms );
225 4 : } else if( str_weights.size()==1 && str_weights[0]=="@Masses" ) {
226 0 : wflab = getShortcutLabel() + "_m";
227 0 : readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
228 4 : } else if( str_weights.size()>1 ) {
229 2 : std::string vals=str_weights[0];
230 6 : for(unsigned i=1; i<str_weights.size(); ++i) {
231 8 : vals += "," + str_weights[i];
232 : }
233 2 : readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
234 4 : wflab=getShortcutLabel() + "_w";
235 : } else {
236 2 : plumed_assert( str_weights.size()==1 );
237 2 : wflab = getShortcutLabel() + "_cent_w";
238 2 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( wflab );
239 2 : if( !av ) {
240 : wflab = str_weights[0];
241 : }
242 : }
243 : // Check for normalisation
244 : bool unorm;
245 8 : parseFlag("UNORMALIZED",unorm);
246 : // Find out the type
247 4 : if( getName()!="GYRATION_TENSOR" ) {
248 0 : parse("TYPE",gtype);
249 0 : if( gtype!="RADIUS" && gtype!="TRACE" && gtype!="GTPC_1" && gtype!="GTPC_2" && gtype!="GTPC_3" && gtype!="ASPHERICITY" && gtype!="ACYLINDRICITY"
250 0 : && gtype!= "KAPPA2" && gtype!="RGYR_1" && gtype!="RGYR_2" && gtype!="RGYR_3" ) {
251 0 : error("type " + gtype + " is invalid");
252 : }
253 : // Check if we need to calculate the unormlised radius
254 0 : if( gtype=="TRACE" || gtype=="KAPPA2" ) {
255 0 : unorm=true;
256 : }
257 : }
258 : // Compute all the vectors separating all the positions from the center
259 4 : std::string distance_act = getShortcutLabel() + "_dists: DISTANCE COMPONENTS" + pbcstr;
260 16 : for(unsigned i=0; i<atoms.size(); ++i) {
261 : std::string num;
262 12 : Tools::convert( i+1, num );
263 48 : distance_act += " ATOMS" + num + "=" + getShortcutLabel() + "_cent," + atoms[i];
264 : }
265 4 : readInputLine( distance_act );
266 : // And calculate the covariance
267 : std::string norm_str;
268 4 : if( unorm ) {
269 : norm_str = " UNORMALIZED";
270 : }
271 4 : if( getName()=="GYRATION_TENSOR" ) {
272 20 : readInputLine( getShortcutLabel() + ": COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
273 : return;
274 : }
275 0 : readInputLine( getShortcutLabel() + "_tensor: COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
276 : // Pick out the diagonal elements
277 0 : readInputLine( getShortcutLabel() + "_diag_elements: SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_tensor COMPONENTS=1.1,2.2,3.3");
278 0 : if( gtype=="RADIUS") {
279 : // And now we need the average trace for the gyration radius
280 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
281 : // Square root the radius
282 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=sqrt(x) PERIODIC=NO");
283 0 : } else if( gtype=="TRACE" ) {
284 : // Compte the trace of the gyration tensor
285 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
286 : // And double it
287 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=2*x PERIODIC=NO");
288 : } else {
289 : // Diagonalize the gyration tensor
290 0 : readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + getShortcutLabel() + "_tensor VECTORS=all" );
291 0 : if( gtype.find("GTPC")!=std::string::npos ) {
292 0 : std::size_t und=gtype.find_first_of("_");
293 0 : if( und==std::string::npos ) {
294 0 : error( gtype + " is not a valid type for gyration radius");
295 : }
296 0 : std::string num = gtype.substr(und+1);
297 0 : if( num!="1" && num!="2" && num!="3" ) {
298 0 : error( gtype + " is not a valid type for gyration radius");
299 : }
300 : // Now get the appropriate eigenvalue
301 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-" + num + " FUNC=sqrt(x) PERIODIC=NO");
302 0 : } else if( gtype.find("RGYR")!=std::string::npos ) {
303 0 : std::size_t und=gtype.find_first_of("_");
304 0 : if( und==std::string::npos ) {
305 0 : error( gtype + " is not a valid type for gyration radius");
306 : }
307 : unsigned ind;
308 0 : Tools::convert( gtype.substr(und+1), ind );
309 : // Now get the appropriate quantity
310 0 : if( ind==3 ) {
311 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2 FUNC=sqrt(x+y) PERIODIC=NO");
312 0 : } else if( ind==2 ) {
313 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
314 0 : } else if( ind==1 ) {
315 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
316 : } else {
317 0 : error( gtype + " is not a valid type for gyration radius");
318 : }
319 0 : } else if( gtype=="ASPHERICITY" ) {
320 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-0.5*(y+z)) PERIODIC=NO" );
321 0 : } else if( gtype=="ACYLINDRICITY" ) {
322 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-y) PERIODIC=NO" );
323 0 : } else if( gtype=="KAPPA2" ) {
324 0 : readInputLine( getShortcutLabel() + "_numer: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x*y+x*z+y*z PERIODIC=NO" );
325 0 : readInputLine( getShortcutLabel() + "_denom: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x+y+z PERIODIC=NO" );
326 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=1-3*(x/(y*y)) PERIODIC=NO");
327 : } else {
328 0 : error( gtype + " is not a valid type for gyration radius");
329 : }
330 : }
331 154 : }
332 :
333 : }
334 : }
|
