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 122 : void GyrationShortcut::registerKeywords( Keywords& keys ) {
116 122 : ActionShortcut::registerKeywords( keys );
117 122 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
118 244 : if( keys.getDisplayName()=="GYRATION" ) {
119 116 : keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
120 : }
121 122 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
122 122 : 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 122 : keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center of mass");
127 122 : keys.addFlag("MASS",false,"calculate the center of mass");
128 122 : keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
129 122 : keys.addFlag("UNORMALIZED",false,"do not divide by the sum of the weights");
130 244 : if( keys.getDisplayName()=="GYRATION" ) {
131 232 : keys.setValueDescription("scalar","the radius that was computed from the weights");
132 232 : keys.addActionNameSuffix("_FAST");
133 12 : } else if( keys.getDisplayName()=="GYRATION_TENSOR" ) {
134 12 : keys.setValueDescription("matrix","the gyration tensor that was computed from the weights");
135 : }
136 122 : keys.needsAction("CENTER");
137 122 : keys.needsAction("CONSTANT");
138 122 : keys.needsAction("ONES");
139 122 : keys.needsAction("MASS");
140 122 : keys.needsAction("DISTANCE");
141 122 : keys.needsAction("COVARIANCE_MATRIX");
142 122 : keys.needsAction("SELECT_COMPONENTS");
143 122 : keys.needsAction("SUM");
144 122 : keys.needsAction("CUSTOM");
145 122 : keys.needsAction("DIAGONALIZE");
146 122 : keys.addDOI("10.1021/jp2065612");
147 122 : }
148 :
149 116 : GyrationShortcut::GyrationShortcut(const ActionOptions& ao):
150 : Action(ao),
151 116 : ActionShortcut(ao) {
152 : bool usemass, phases, massw;
153 116 : parseFlag("MASS",usemass);
154 116 : parseFlag("MASS_WEIGHTED",massw);
155 116 : if( massw ) {
156 0 : usemass = true;
157 : }
158 232 : parseFlag("PHASES",phases);
159 : std::vector<std::string> str_weights;
160 232 : parseVector("WEIGHTS",str_weights);
161 : std::string wflab;
162 230 : if( !phases && getName()=="GYRATION" ) {
163 112 : if( usemass || str_weights.size()==0 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
164 : std::string wt_str;
165 112 : if( str_weights.size()>0 ) {
166 0 : wt_str="WEIGHTS=" + str_weights[0];
167 0 : for(unsigned i=1; i<str_weights.size(); ++i) {
168 0 : wt_str += "," + str_weights[i];
169 : }
170 : }
171 112 : if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
172 : wt_str = "MASS_WEIGHTED";
173 : }
174 232 : readInputLine( getShortcutLabel() + ": GYRATION_FAST " + wt_str + " " + convertInputLineToString() );
175 : return;
176 : }
177 : }
178 4 : if( usemass ) {
179 0 : str_weights.resize(1);
180 : str_weights[0]="@Masses";
181 : }
182 10 : log<<" Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)")<<"\n";
183 : // Read in the atoms involved
184 : std::vector<std::string> atoms;
185 4 : parseVector("ATOMS",atoms);
186 4 : Tools::interpretRanges(atoms);
187 4 : std::string gtype, atlist=atoms[0];
188 12 : for(unsigned i=1; i<atoms.size(); ++i) {
189 16 : atlist += "," + atoms[i];
190 : }
191 : bool nopbc;
192 8 : parseFlag("NOPBC",nopbc);
193 : std::string pbcstr;
194 4 : if(nopbc) {
195 : pbcstr = " NOPBC";
196 : }
197 : std::string phasestr;
198 4 : if(phases) {
199 : phasestr = " PHASES";
200 : }
201 : // Create the geometric center of the molecule
202 4 : std::string weights_str="";
203 4 : if( str_weights.size()>0 ) {
204 4 : weights_str=" WEIGHTS=" + str_weights[0];
205 8 : for(unsigned i=1; i<str_weights.size(); ++i) {
206 8 : weights_str += "," + str_weights[i];
207 : }
208 : }
209 8 : readInputLine( getShortcutLabel() + "_cent: CENTER ATOMS=" + atlist + pbcstr + phasestr + weights_str );
210 4 : if( str_weights.size()==0 ) {
211 0 : wflab = getShortcutLabel() + "_w";
212 : std::string str_natoms;
213 0 : Tools::convert( atoms.size(), str_natoms );
214 0 : readInputLine( getShortcutLabel() + "_w: ONES SIZE=" + str_natoms );
215 6 : } else if( str_weights.size()==1 && str_weights[0]=="@Masses" ) {
216 0 : wflab = getShortcutLabel() + "_m";
217 0 : readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
218 4 : } else if( str_weights.size()>1 ) {
219 2 : std::string vals=str_weights[0];
220 6 : for(unsigned i=1; i<str_weights.size(); ++i) {
221 8 : vals += "," + str_weights[i];
222 : }
223 4 : readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
224 4 : wflab=getShortcutLabel() + "_w";
225 : } else {
226 2 : plumed_assert( str_weights.size()==1 );
227 2 : wflab = getShortcutLabel() + "_cent_w";
228 2 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( wflab );
229 2 : if( !av ) {
230 : wflab = str_weights[0];
231 : }
232 : }
233 : // Check for normalisation
234 : bool unorm;
235 8 : parseFlag("UNORMALIZED",unorm);
236 : // Find out the type
237 4 : if( getName()!="GYRATION_TENSOR" ) {
238 0 : parse("TYPE",gtype);
239 0 : if( gtype!="RADIUS" && gtype!="TRACE" && gtype!="GTPC_1" && gtype!="GTPC_2" && gtype!="GTPC_3" && gtype!="ASPHERICITY" && gtype!="ACYLINDRICITY"
240 0 : && gtype!= "KAPPA2" && gtype!="RGYR_1" && gtype!="RGYR_2" && gtype!="RGYR_3" ) {
241 0 : error("type " + gtype + " is invalid");
242 : }
243 : // Check if we need to calculate the unormlised radius
244 0 : if( gtype=="TRACE" || gtype=="KAPPA2" ) {
245 0 : unorm=true;
246 : }
247 : }
248 : // Compute all the vectors separating all the positions from the center
249 4 : std::string distance_act = getShortcutLabel() + "_dists: DISTANCE COMPONENTS" + pbcstr;
250 16 : for(unsigned i=0; i<atoms.size(); ++i) {
251 : std::string num;
252 12 : Tools::convert( i+1, num );
253 24 : distance_act += " ATOMS" + num + "=" + getShortcutLabel() + "_cent," + atoms[i];
254 : }
255 4 : readInputLine( distance_act );
256 : // And calculate the covariance
257 : std::string norm_str;
258 4 : if( unorm ) {
259 : norm_str = " UNORMALIZED";
260 : }
261 4 : if( getName()=="GYRATION_TENSOR" ) {
262 8 : readInputLine( getShortcutLabel() + ": COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
263 : return;
264 : }
265 0 : readInputLine( getShortcutLabel() + "_tensor: COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
266 : // Pick out the diagonal elements
267 0 : readInputLine( getShortcutLabel() + "_diag_elements: SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_tensor COMPONENTS=1.1,2.2,3.3");
268 0 : if( gtype=="RADIUS") {
269 : // And now we need the average trace for the gyration radius
270 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
271 : // Square root the radius
272 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=sqrt(x) PERIODIC=NO");
273 0 : } else if( gtype=="TRACE" ) {
274 : // Compte the trace of the gyration tensor
275 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
276 : // And double it
277 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=2*x PERIODIC=NO");
278 : } else {
279 : // Diagonalize the gyration tensor
280 0 : readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + getShortcutLabel() + "_tensor VECTORS=all" );
281 0 : if( gtype.find("GTPC")!=std::string::npos ) {
282 0 : std::size_t und=gtype.find_first_of("_");
283 0 : if( und==std::string::npos ) {
284 0 : error( gtype + " is not a valid type for gyration radius");
285 : }
286 0 : std::string num = gtype.substr(und+1);
287 0 : if( num!="1" && num!="2" && num!="3" ) {
288 0 : error( gtype + " is not a valid type for gyration radius");
289 : }
290 : // Now get the appropriate eigenvalue
291 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-" + num + " FUNC=sqrt(x) PERIODIC=NO");
292 0 : } else if( gtype.find("RGYR")!=std::string::npos ) {
293 0 : std::size_t und=gtype.find_first_of("_");
294 0 : if( und==std::string::npos ) {
295 0 : error( gtype + " is not a valid type for gyration radius");
296 : }
297 : unsigned ind;
298 0 : Tools::convert( gtype.substr(und+1), ind );
299 : // Now get the appropriate quantity
300 0 : if( ind==3 ) {
301 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2 FUNC=sqrt(x+y) PERIODIC=NO");
302 0 : } else if( ind==2 ) {
303 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
304 0 : } else if( ind==1 ) {
305 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
306 : } else {
307 0 : error( gtype + " is not a valid type for gyration radius");
308 : }
309 0 : } else if( gtype=="ASPHERICITY" ) {
310 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" );
311 0 : } else if( gtype=="ACYLINDRICITY" ) {
312 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-y) PERIODIC=NO" );
313 0 : } else if( gtype=="KAPPA2" ) {
314 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" );
315 0 : readInputLine( getShortcutLabel() + "_denom: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x+y+z PERIODIC=NO" );
316 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=1-3*(x/(y*y)) PERIODIC=NO");
317 : } else {
318 0 : error( gtype + " is not a valid type for gyration radius");
319 : }
320 : }
321 124 : }
322 :
323 : }
324 : }
|
