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 COLVAR GYRATION_FAST
27 : /*
28 : Calculate the radius of gyration, or other properties related to it.
29 :
30 : \par Examples
31 :
32 : */
33 : //+ENDPLUMEDOC
34 :
35 : //+PLUMEDOC MCOLVAR GYRATION_TENSOR
36 : /*
37 : Calculate the gyration tensor using a user specified vector of weights
38 :
39 : The elements of the $3 \times 3$ gyration tensor are defined as:
40 :
41 : $$
42 : G_{\alpha\beta} = \frac{\sum_i^{n} w_i (\alpha_i-\alpha_{\rm COM})(\beta_i - \beta_{\rm COM})}{\sum_i^{n} w_i}
43 : $$
44 :
45 : where $\alpha_i$ and $\beta_i$ can be the $x$, $y$ or $z$ coordinates of atom $i$ and $\alpha_{\rm COM}$ and
46 : $\beta_{\rm COM}$ can be the $x$, $y$ or $z$ components of the center, which is calculated using:
47 :
48 : $$
49 : {r}_{\rm COM}=\frac{\sum_i^{n} {r}_i\ w_i }{\sum_i^{n} w_i}
50 : $$
51 :
52 : The following example shows how you can calculate and print the gyration tensor from the positions of the fist 10 atoms using PLUMED
53 :
54 : ```plumed
55 : g: GYRATION_TENSOR ATOMS=1-10
56 : PRINT ARG=g FILE=colvar
57 : ```
58 :
59 : The 9 elements of the gyration matrix will be output to the file `colvar` here.
60 :
61 : Similar functionality to the functionality in the example above is used in the [GYRATION](GYRATION.md) shortcut. There is, however,
62 : 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
63 : used when the weights are all one or when the masses are used as the weights.
64 :
65 : */
66 : //+ENDPLUMEDOC
67 :
68 : class GyrationShortcut : public ActionShortcut {
69 : public:
70 : static void registerKeywords( Keywords& keys );
71 : explicit GyrationShortcut(const ActionOptions&);
72 : };
73 :
74 : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION")
75 : PLUMED_REGISTER_ACTION(GyrationShortcut,"GYRATION_TENSOR")
76 :
77 122 : void GyrationShortcut::registerKeywords( Keywords& keys ) {
78 122 : ActionShortcut::registerKeywords( keys );
79 122 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
80 122 : keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
81 122 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
82 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. "
83 : "If WEIGHTS=@Masses is used the center of mass is computed. If WEIGHTS=@charges the center of charge is computed. If "
84 : "the label of an action is provided PLUMED assumes that that action calculates a list of symmetry functions that can be used "
85 : "as weights. Lastly, an explicit list of numbers to use as weights can be provided");
86 122 : keys.addFlag("PHASES",false,"use trigonometric phases when computing position of center of mass");
87 122 : keys.addFlag("MASS",false,"calculate the center of mass");
88 122 : keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
89 122 : keys.addFlag("UNORMALIZED",false,"do not divide by the sum of the weights");
90 244 : if( keys.getDisplayName()=="GYRATION" ) {
91 232 : keys.setValueDescription("scalar","the radius that was computed from the weights");
92 232 : keys.addActionNameSuffix("_FAST");
93 12 : } else if( keys.getDisplayName()=="GYRATION_TENSOR" ) {
94 12 : keys.setValueDescription("matrix","the gyration tensor that was computed from the weights");
95 : }
96 122 : keys.needsAction("CENTER");
97 122 : keys.needsAction("CONSTANT");
98 122 : keys.needsAction("ONES");
99 122 : keys.needsAction("MASS");
100 122 : keys.needsAction("DISTANCE");
101 122 : keys.needsAction("COVARIANCE_MATRIX");
102 122 : keys.needsAction("SELECT_COMPONENTS");
103 122 : keys.needsAction("SUM");
104 122 : keys.needsAction("CUSTOM");
105 122 : keys.needsAction("DIAGONALIZE");
106 122 : keys.addDOI("10.1021/jp2065612");
107 122 : }
108 :
109 116 : GyrationShortcut::GyrationShortcut(const ActionOptions& ao):
110 : Action(ao),
111 116 : ActionShortcut(ao) {
112 : bool usemass, phases;
113 116 : parseFlag("MASS",usemass);
114 232 : parseFlag("PHASES",phases);
115 : std::vector<std::string> str_weights;
116 232 : parseVector("WEIGHTS",str_weights);
117 : std::string wflab;
118 230 : if( !phases && getName()=="GYRATION" ) {
119 112 : if( usemass || str_weights.size()==0 || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
120 : std::string wt_str;
121 112 : if( str_weights.size()>0 ) {
122 0 : wt_str="WEIGHTS=" + str_weights[0];
123 0 : for(unsigned i=1; i<str_weights.size(); ++i) {
124 0 : wt_str += "," + str_weights[i];
125 : }
126 : }
127 112 : if( usemass || (str_weights.size()==1 && str_weights[0]=="@Masses") ) {
128 : wt_str = "MASS";
129 : }
130 232 : readInputLine( getShortcutLabel() + ": GYRATION_FAST " + wt_str + " " + convertInputLineToString() );
131 : return;
132 : }
133 : }
134 4 : if( usemass ) {
135 0 : str_weights.resize(1);
136 : str_weights[0]="@Masses";
137 : }
138 10 : log<<" Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)")<<"\n";
139 : // Read in the atoms involved
140 : std::vector<std::string> atoms;
141 4 : parseVector("ATOMS",atoms);
142 4 : Tools::interpretRanges(atoms);
143 4 : std::string gtype, atlist=atoms[0];
144 12 : for(unsigned i=1; i<atoms.size(); ++i) {
145 16 : atlist += "," + atoms[i];
146 : }
147 : bool nopbc;
148 8 : parseFlag("NOPBC",nopbc);
149 : std::string pbcstr;
150 4 : if(nopbc) {
151 : pbcstr = " NOPBC";
152 : }
153 : std::string phasestr;
154 4 : if(phases) {
155 : phasestr = " PHASES";
156 : }
157 : // Create the geometric center of the molecule
158 4 : std::string weights_str="";
159 4 : if( str_weights.size()>0 ) {
160 4 : weights_str=" WEIGHTS=" + str_weights[0];
161 8 : for(unsigned i=1; i<str_weights.size(); ++i) {
162 8 : weights_str += "," + str_weights[i];
163 : }
164 : }
165 8 : readInputLine( getShortcutLabel() + "_cent: CENTER ATOMS=" + atlist + pbcstr + phasestr + weights_str );
166 4 : if( str_weights.size()==0 ) {
167 0 : wflab = getShortcutLabel() + "_w";
168 : std::string str_natoms;
169 0 : Tools::convert( atoms.size(), str_natoms );
170 0 : readInputLine( getShortcutLabel() + "_w: ONES SIZE=" + str_natoms );
171 6 : } else if( str_weights.size()==1 && str_weights[0]=="@Masses" ) {
172 0 : wflab = getShortcutLabel() + "_m";
173 0 : readInputLine( getShortcutLabel() + "_m: MASS ATOMS=" + atlist );
174 4 : } else if( str_weights.size()>1 ) {
175 2 : std::string vals=str_weights[0];
176 6 : for(unsigned i=1; i<str_weights.size(); ++i) {
177 8 : vals += "," + str_weights[i];
178 : }
179 4 : readInputLine( getShortcutLabel() + "_w: CONSTANT VALUES=" + vals );
180 4 : wflab=getShortcutLabel() + "_w";
181 : } else {
182 2 : plumed_assert( str_weights.size()==1 );
183 2 : wflab = getShortcutLabel() + "_cent_w";
184 2 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( wflab );
185 2 : if( !av ) {
186 : wflab = str_weights[0];
187 : }
188 : }
189 : // Check for normalisation
190 : bool unorm;
191 8 : parseFlag("UNORMALIZED",unorm);
192 : // Find out the type
193 4 : if( getName()!="GYRATION_TENSOR" ) {
194 0 : parse("TYPE",gtype);
195 0 : if( gtype!="RADIUS" && gtype!="TRACE" && gtype!="GTPC_1" && gtype!="GTPC_2" && gtype!="GTPC_3" && gtype!="ASPHERICITY" && gtype!="ACYLINDRICITY"
196 0 : && gtype!= "KAPPA2" && gtype!="RGYR_1" && gtype!="RGYR_2" && gtype!="RGYR_3" ) {
197 0 : error("type " + gtype + " is invalid");
198 : }
199 : // Check if we need to calculate the unormlised radius
200 0 : if( gtype=="TRACE" || gtype=="KAPPA2" ) {
201 0 : unorm=true;
202 : }
203 : }
204 : // Compute all the vectors separating all the positions from the center
205 4 : std::string distance_act = getShortcutLabel() + "_dists: DISTANCE COMPONENTS" + pbcstr;
206 16 : for(unsigned i=0; i<atoms.size(); ++i) {
207 : std::string num;
208 12 : Tools::convert( i+1, num );
209 24 : distance_act += " ATOMS" + num + "=" + getShortcutLabel() + "_cent," + atoms[i];
210 : }
211 4 : readInputLine( distance_act );
212 : // And calculate the covariance
213 : std::string norm_str;
214 4 : if( unorm ) {
215 : norm_str = " UNORMALIZED";
216 : }
217 4 : if( getName()=="GYRATION_TENSOR" ) {
218 8 : readInputLine( getShortcutLabel() + ": COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
219 : return;
220 : }
221 0 : readInputLine( getShortcutLabel() + "_tensor: COVARIANCE_MATRIX ARG=" + getShortcutLabel() + "_dists.x," + getShortcutLabel() + "_dists.y," + getShortcutLabel() + "_dists.z WEIGHTS=" + wflab + norm_str );
222 : // Pick out the diagonal elements
223 0 : readInputLine( getShortcutLabel() + "_diag_elements: SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_tensor COMPONENTS=1.1,2.2,3.3");
224 0 : if( gtype=="RADIUS") {
225 : // And now we need the average trace for the gyration radius
226 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
227 : // Square root the radius
228 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=sqrt(x) PERIODIC=NO");
229 0 : } else if( gtype=="TRACE" ) {
230 : // Compte the trace of the gyration tensor
231 0 : readInputLine( getShortcutLabel() + "_trace: SUM ARG=" + getShortcutLabel() + "_diag_elements PERIODIC=NO");
232 : // And double it
233 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_trace FUNC=2*x PERIODIC=NO");
234 : } else {
235 : // Diagonalize the gyration tensor
236 0 : readInputLine( getShortcutLabel() + "_diag: DIAGONALIZE ARG=" + getShortcutLabel() + "_tensor VECTORS=all" );
237 0 : if( gtype.find("GTPC")!=std::string::npos ) {
238 0 : std::size_t und=gtype.find_first_of("_");
239 0 : if( und==std::string::npos ) {
240 0 : error( gtype + " is not a valid type for gyration radius");
241 : }
242 0 : std::string num = gtype.substr(und+1);
243 0 : if( num!="1" && num!="2" && num!="3" ) {
244 0 : error( gtype + " is not a valid type for gyration radius");
245 : }
246 : // Now get the appropriate eigenvalue
247 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-" + num + " FUNC=sqrt(x) PERIODIC=NO");
248 0 : } else if( gtype.find("RGYR")!=std::string::npos ) {
249 0 : std::size_t und=gtype.find_first_of("_");
250 0 : if( und==std::string::npos ) {
251 0 : error( gtype + " is not a valid type for gyration radius");
252 : }
253 : unsigned ind;
254 0 : Tools::convert( gtype.substr(und+1), ind );
255 : // Now get the appropriate quantity
256 0 : if( ind==3 ) {
257 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2 FUNC=sqrt(x+y) PERIODIC=NO");
258 0 : } else if( ind==2 ) {
259 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
260 0 : } else if( ind==1 ) {
261 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x+y) PERIODIC=NO");
262 : } else {
263 0 : error( gtype + " is not a valid type for gyration radius");
264 : }
265 0 : } else if( gtype=="ASPHERICITY" ) {
266 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" );
267 0 : } else if( gtype=="ACYLINDRICITY" ) {
268 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=sqrt(x-y) PERIODIC=NO" );
269 0 : } else if( gtype=="KAPPA2" ) {
270 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" );
271 0 : readInputLine( getShortcutLabel() + "_denom: CUSTOM ARG=" + getShortcutLabel() + "_diag.vals-1," + getShortcutLabel() + "_diag.vals-2," + getShortcutLabel() + "_diag.vals-3 FUNC=x+y+z PERIODIC=NO" );
272 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=1-3*(x/(y*y)) PERIODIC=NO");
273 : } else {
274 0 : error( gtype + " is not a valid type for gyration radius");
275 : }
276 : }
277 124 : }
278 :
279 : }
280 : }
|
