Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 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 "Colvar.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 :
26 : namespace PLMD {
27 : namespace colvar {
28 :
29 : //+PLUMEDOC COLVAR GYRATION
30 : /*
31 : Calculate the radius of gyration, or other properties related to it.
32 :
33 : The radius of gyration is calculated using:
34 :
35 : $$
36 : s_{\rm Gyr}=\Big ( \frac{\sum_i^{n} w_i \vert {r}_i -{r}_{\rm COM} \vert ^2 }{\sum_i^{n} w_i} \Big)^{1/2}
37 : $$
38 :
39 : with the position of the center ${r}_{\rm COM}$ given by:
40 :
41 : $$
42 : {r}_{\rm COM}=\frac{\sum_i^{n} {r}_i\ w_i }{\sum_i^{n} w_i}
43 : $$
44 :
45 : In these expressions $r_i$ indicates the position of atom $i$ and $w_i$ is the weight for atom $i$. The following input
46 : shows how you can calculate the expressions for a set of atoms by using PLUMED:
47 :
48 : ```plumed
49 : w: CONSTANT VALUES=1,2,2,3,4
50 : g: GYRATION ATOMS=1-5 TYPE=RADIUS WEIGHTS=w
51 : PRINT ARG=g FILE=colvar
52 : ```
53 :
54 : In the above input the weights in the expressions above are set equal to the eleemnts of a constant vector, `w`. A more common
55 : approach is to use the masses of the atoms as in the input below:
56 :
57 : ```plumed
58 : g: GYRATION ATOMS=1-5 TYPE=RADIUS MASS_WEIGHTED
59 : # This input is equivalent
60 : # g: GYRATION ATOMS=1-5 TYPE=RADIUS WEIGHTS=@Masses
61 : PRINT ARG=g FILE=colvar
62 : ```
63 :
64 : or to set all the input weights equal to one as is done in the input below;
65 :
66 : ```plumed
67 : g: GYRATION ATOMS=1-5 TYPE=RADIUS
68 : PRINT ARG=g FILE=colvar
69 : ```
70 :
71 : Notice that in the first input above GYRATION is a shortcut and that in the second two inputs it is not. This is because there
72 : is a faster (non-shortcutted) version that can be used for these two more-commonplace inputs. The shortcut input is still useful
73 : as it allows one to do less comonplace analyses such as this one where a clustering is performed and the radius of gyration for the
74 : largest cluster is determined.
75 :
76 : ```plumed
77 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1 D_MAX=0.5}
78 : dfs: DFSCLUSTERING ARG=c1
79 : w: CLUSTER_WEIGHTS CLUSTERS=dfs CLUSTER=1
80 : g: GYRATION ATOMS=1-100 WEIGHTS=w TYPE=RADIUS
81 : ```
82 :
83 : In calculating the gyration radius you first calculate the [GYRATION_TENSOR](GYRATION_TENSOR.md). Instad of computing the radius from this
84 : tensor you can instead compute the trace of the tensor by using the following input:
85 :
86 : ```plumed
87 : w: CONSTANT VALUES=1,2,2,3,4
88 : g: GYRATION ATOMS=1-5 TYPE=TRACE WEIGHTS=w
89 : PRINT ARG=g FILE=colvar
90 : ```
91 :
92 : You can also calculate the largest, middle and smallest principal moments of the tensor:
93 :
94 : ```plumed
95 : w: CONSTANT VALUES=1,2,2,3,4
96 : # This computes the largest principal moment
97 : l: GYRATION ATOMS=1-5 TYPE=GTPC_1 WEIGHTS=w
98 : # This computes the middle principal moment
99 : m: GYRATION ATOMS=1-5 TYPE=GTPC_2 WEIGHTS=w
100 : # This computes the smallest principal moment
101 : s: GYRATION ATOMS=1-5 TYPE=GTPC_3 WEIGHTS=w
102 : PRINT ARG=l,m,s FILE=colvar
103 : ```
104 :
105 : From these principal moments you can calculate the Asphericiry, Acylindricity, or the Relative Shape Anisotropy by using the following input:
106 :
107 : ```plumed
108 : w: CONSTANT VALUES=1,2,2,3,4
109 : # This computes the Asphericiry
110 : l: GYRATION ATOMS=1-5 TYPE=ASPHERICITY WEIGHTS=w
111 : # This computes the Acylindricity
112 : m: GYRATION ATOMS=1-5 TYPE=ACYLINDRICITY WEIGHTS=w
113 : # This computes the Relative Shape Anisotropy
114 : s: GYRATION ATOMS=1-5 TYPE=KAPPA2 WEIGHTS=w
115 : PRINT ARG=l,m,s FILE=colvar
116 : ```
117 :
118 : Lastly, you can calculate the largest, middle and smallest principal radii of gyration by using the following input:
119 :
120 : ```plumed
121 : w: CONSTANT VALUES=1,2,2,3,4
122 : # This computes the largest principal radius of gyration
123 : l: GYRATION ATOMS=1-5 TYPE=RGYR_1 WEIGHTS=w
124 : # This computes the middle principal radius of gyration
125 : m: GYRATION ATOMS=1-5 TYPE=RGYR_2 WEIGHTS=w
126 : # This computes the smallest principal radius of gyration
127 : s: GYRATION ATOMS=1-5 TYPE=RGYR_3 WEIGHTS=w
128 : PRINT ARG=l,m,s FILE=colvar
129 : ```
130 :
131 : By expanding the shortcuts in the inputs above or by reading the papers cited below you can see how these quantities are computed from the [GYRATION_TENSOR](GYRATION_TENSOR.md).
132 : Notice, however, that as with the radius if you use the masses or a vector of ones as weights a faster implemention of this colvar that
133 : calculates these quantities directly will be used.
134 :
135 : ## A note on periodic boundary conditions
136 :
137 : Calculating the radius of gyration is normally used to determine the shape of a molecule so all the specified atoms
138 : would normally be part of the same molecule. When computing this CV it is important to ensure that the periodic boundaries
139 : are calculated correctly. There are two ways that you can manage periodic boundary conditions when using this action. The
140 : first and simplest is to reconstruct the molecule similarly to the way that [WHOLEMOLECULES](WHOLEMOLECULES.md) operates.
141 : This reconstruction of molecules has been done automatically since PLUMED 2.2. If for some reason you want to turn it off
142 : you can use the NOPBC flag.
143 :
144 : An alternative approach to handling PBC is to use the PHASES keyword. This keyword instructs PLUMED to use the PHASES option
145 : when computing the position of the center using the [CENTER](CENTER.md) command. Distances of atoms from this center are then
146 : computed using PBC as usual.
147 :
148 : */
149 : //+ENDPLUMEDOC
150 :
151 : class Gyration : public Colvar {
152 : private:
153 : enum CV_TYPE {RADIUS, TRACE, GTPC_1, GTPC_2, GTPC_3, ASPHERICITY, ACYLINDRICITY, KAPPA2, GYRATION_3, GYRATION_2, GYRATION_1, TOT};
154 : int rg_type;
155 : bool use_masses;
156 : bool nopbc;
157 : std::vector<Vector> derivatives;
158 : public:
159 : static void registerKeywords(Keywords& keys);
160 : explicit Gyration(const ActionOptions&);
161 : void calculate() override;
162 : };
163 :
164 : PLUMED_REGISTER_ACTION(Gyration,"GYRATION_FAST")
165 :
166 224 : void Gyration::registerKeywords(Keywords& keys) {
167 224 : Colvar::registerKeywords(keys);
168 448 : keys.setDisplayName("GYRATION");
169 224 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
170 224 : keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
171 224 : keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
172 448 : keys.setValueDescription("scalar","the radius of gyration");
173 224 : }
174 :
175 112 : Gyration::Gyration(const ActionOptions&ao):
176 : PLUMED_COLVAR_INIT(ao),
177 112 : use_masses(false),
178 112 : nopbc(false) {
179 : std::vector<AtomNumber> atoms;
180 223 : parseAtomList("ATOMS",atoms);
181 111 : if(atoms.size()==0) {
182 0 : error("no atoms specified");
183 : }
184 224 : parseFlag("MASS_WEIGHTED",use_masses);
185 : std::string Type;
186 111 : parse("TYPE",Type);
187 111 : parseFlag("NOPBC",nopbc);
188 111 : checkRead();
189 :
190 111 : if(Type=="RADIUS") {
191 90 : rg_type=RADIUS;
192 21 : } else if(Type=="TRACE") {
193 2 : rg_type=TRACE;
194 19 : } else if(Type=="GTPC_1") {
195 2 : rg_type=GTPC_1;
196 17 : } else if(Type=="GTPC_2") {
197 2 : rg_type=GTPC_2;
198 15 : } else if(Type=="GTPC_3") {
199 2 : rg_type=GTPC_3;
200 13 : } else if(Type=="ASPHERICITY") {
201 2 : rg_type=ASPHERICITY;
202 11 : } else if(Type=="ACYLINDRICITY") {
203 2 : rg_type=ACYLINDRICITY;
204 9 : } else if(Type=="KAPPA2") {
205 2 : rg_type=KAPPA2;
206 7 : } else if(Type=="RGYR_3") {
207 2 : rg_type=GYRATION_3;
208 5 : } else if(Type=="RGYR_2") {
209 2 : rg_type=GYRATION_2;
210 3 : } else if(Type=="RGYR_1") {
211 2 : rg_type=GYRATION_1;
212 : } else {
213 1 : error("Unknown GYRATION type");
214 : }
215 :
216 110 : switch(rg_type) {
217 90 : case RADIUS:
218 90 : log.printf(" GYRATION RADIUS (Rg);");
219 : break;
220 2 : case TRACE:
221 2 : log.printf(" TRACE OF THE GYRATION TENSOR;");
222 : break;
223 2 : case GTPC_1:
224 2 : log.printf(" THE LARGEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_1);");
225 : break;
226 2 : case GTPC_2:
227 2 : log.printf(" THE MIDDLE PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_2);");
228 : break;
229 2 : case GTPC_3:
230 2 : log.printf(" THE SMALLEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_3);");
231 : break;
232 2 : case ASPHERICITY:
233 2 : log.printf(" THE ASPHERICITY (b');");
234 : break;
235 2 : case ACYLINDRICITY:
236 2 : log.printf(" THE ACYLINDRICITY (c');");
237 : break;
238 2 : case KAPPA2:
239 2 : log.printf(" THE RELATIVE SHAPE ANISOTROPY (kappa^2);");
240 : break;
241 2 : case GYRATION_3:
242 2 : log.printf(" THE SMALLEST PRINCIPAL RADIUS OF GYRATION (r_g3);");
243 : break;
244 2 : case GYRATION_2:
245 2 : log.printf(" THE MIDDLE PRINCIPAL RADIUS OF GYRATION (r_g2);");
246 : break;
247 2 : case GYRATION_1:
248 2 : log.printf(" THE LARGEST PRINCIPAL RADIUS OF GYRATION (r_g1);");
249 : break;
250 : }
251 110 : if(rg_type>TRACE) {
252 36 : log<<" Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)");
253 : }
254 110 : log<<"\n";
255 :
256 110 : log.printf(" atoms involved : ");
257 1023 : for(unsigned i=0; i<atoms.size(); ++i) {
258 913 : log.printf("%d ",atoms[i].serial());
259 : }
260 110 : log.printf("\n");
261 :
262 110 : if(nopbc) {
263 4 : log<<" PBC will be ignored\n";
264 : } else {
265 106 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
266 : }
267 :
268 111 : addValueWithDerivatives();
269 110 : setNotPeriodic();
270 110 : requestAtoms(atoms);
271 114 : }
272 :
273 1238 : void Gyration::calculate() {
274 :
275 1238 : if(!nopbc) {
276 580 : makeWhole();
277 : }
278 :
279 1238 : Vector com;
280 : double totmass = 0.;
281 1238 : if( use_masses ) {
282 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
283 0 : totmass+=getMass(i);
284 0 : com+=getMass(i)*getPosition(i);
285 : }
286 : } else {
287 1238 : totmass = static_cast<double>(getNumberOfAtoms());
288 10803 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
289 9565 : com+=getPosition(i);
290 : }
291 : }
292 1238 : com /= totmass;
293 :
294 1238 : double rgyr=0.;
295 1238 : derivatives.resize(getNumberOfAtoms());
296 :
297 1238 : if(rg_type==RADIUS||rg_type==TRACE) {
298 838 : if( use_masses ) {
299 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
300 0 : const Vector diff = delta( com, getPosition(i) );
301 0 : rgyr += getMass(i)*diff.modulo2();
302 0 : derivatives[i] = diff*getMass(i);
303 : }
304 : } else {
305 8403 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
306 7565 : const Vector diff = delta( com, getPosition(i) );
307 7565 : rgyr += diff.modulo2();
308 7565 : derivatives[i] = diff;
309 : }
310 : }
311 : double fact;
312 838 : if(rg_type==RADIUS) {
313 708 : rgyr = std::sqrt(rgyr/totmass);
314 708 : fact = 1./(rgyr*totmass);
315 : } else {
316 130 : rgyr = 2.*rgyr;
317 : fact = 4;
318 : }
319 838 : setValue(rgyr);
320 8403 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
321 7565 : setAtomsDerivatives(i,fact*derivatives[i]);
322 : }
323 838 : setBoxDerivativesNoPbc();
324 838 : return;
325 : }
326 :
327 :
328 400 : Tensor3d gyr_tens;
329 : //calculate gyration tensor
330 400 : if( use_masses ) {
331 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
332 0 : const Vector diff=delta( com, getPosition(i) );
333 0 : gyr_tens[0][0]+=getMass(i)*diff[0]*diff[0];
334 0 : gyr_tens[1][1]+=getMass(i)*diff[1]*diff[1];
335 0 : gyr_tens[2][2]+=getMass(i)*diff[2]*diff[2];
336 0 : gyr_tens[0][1]+=getMass(i)*diff[0]*diff[1];
337 0 : gyr_tens[0][2]+=getMass(i)*diff[0]*diff[2];
338 0 : gyr_tens[1][2]+=getMass(i)*diff[1]*diff[2];
339 : }
340 : } else {
341 2400 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
342 2000 : const Vector diff=delta( com, getPosition(i) );
343 2000 : gyr_tens[0][0]+=diff[0]*diff[0];
344 2000 : gyr_tens[1][1]+=diff[1]*diff[1];
345 2000 : gyr_tens[2][2]+=diff[2]*diff[2];
346 2000 : gyr_tens[0][1]+=diff[0]*diff[1];
347 2000 : gyr_tens[0][2]+=diff[0]*diff[2];
348 2000 : gyr_tens[1][2]+=diff[1]*diff[2];
349 : }
350 : }
351 :
352 : // first make the matrix symmetric
353 400 : gyr_tens[1][0] = gyr_tens[0][1];
354 400 : gyr_tens[2][0] = gyr_tens[0][2];
355 400 : gyr_tens[2][1] = gyr_tens[1][2];
356 400 : Tensor3d ttransf,transf;
357 400 : Vector princ_comp,prefactor;
358 : //diagonalize gyration tensor
359 400 : diagMatSym(gyr_tens, princ_comp, ttransf);
360 400 : transf=transpose(ttransf);
361 : //sort eigenvalues and eigenvectors
362 400 : if (princ_comp[0]<princ_comp[1]) {
363 : double tmp=princ_comp[0];
364 400 : princ_comp[0]=princ_comp[1];
365 400 : princ_comp[1]=tmp;
366 1600 : for (unsigned i=0; i<3; i++) {
367 1200 : tmp=transf[i][0];
368 1200 : transf[i][0]=transf[i][1];
369 1200 : transf[i][1]=tmp;
370 : }
371 : }
372 400 : if (princ_comp[1]<princ_comp[2]) {
373 : double tmp=princ_comp[1];
374 400 : princ_comp[1]=princ_comp[2];
375 400 : princ_comp[2]=tmp;
376 1600 : for (unsigned i=0; i<3; i++) {
377 1200 : tmp=transf[i][1];
378 1200 : transf[i][1]=transf[i][2];
379 1200 : transf[i][2]=tmp;
380 : }
381 : }
382 400 : if (princ_comp[0]<princ_comp[1]) {
383 : double tmp=princ_comp[0];
384 400 : princ_comp[0]=princ_comp[1];
385 400 : princ_comp[1]=tmp;
386 1600 : for (unsigned i=0; i<3; i++) {
387 1200 : tmp=transf[i][0];
388 1200 : transf[i][0]=transf[i][1];
389 1200 : transf[i][1]=tmp;
390 : }
391 : }
392 : //calculate determinant of transformation matrix
393 : double det = determinant(transf);
394 : // transformation matrix for rotation must have positive determinant, otherwise multiply one column by (-1)
395 400 : if(det<0) {
396 1600 : for(unsigned j=0; j<3; j++) {
397 1200 : transf[j][2]=-transf[j][2];
398 : }
399 400 : det = -det;
400 : }
401 400 : if(std::abs(det-1.)>0.0001) {
402 0 : error("Plumed Error: Cannot diagonalize gyration tensor\n");
403 : }
404 400 : switch(rg_type) {
405 135 : case GTPC_1:
406 : case GTPC_2:
407 : case GTPC_3: {
408 135 : int pc_index = rg_type-2; //index of principal component
409 135 : rgyr=std::sqrt(princ_comp[pc_index]/totmass);
410 135 : double rm = rgyr*totmass;
411 135 : if(rm>1e-6) {
412 135 : prefactor[pc_index]=1.0/rm; //some parts of derivate
413 : }
414 : break;
415 : }
416 0 : case GYRATION_3: { //the smallest principal radius of gyration
417 0 : rgyr=std::sqrt((princ_comp[1]+princ_comp[2])/totmass);
418 0 : double rm = rgyr*totmass;
419 0 : if (rm>1e-6) {
420 0 : prefactor[1]=1.0/rm;
421 0 : prefactor[2]=1.0/rm;
422 : }
423 : break;
424 : }
425 130 : case GYRATION_2: { //the midle principal radius of gyration
426 130 : rgyr=std::sqrt((princ_comp[0]+princ_comp[2])/totmass);
427 130 : double rm = rgyr*totmass;
428 130 : if (rm>1e-6) {
429 130 : prefactor[0]=1.0/rm;
430 130 : prefactor[2]=1.0/rm;
431 : }
432 : break;
433 : }
434 0 : case GYRATION_1: { //the largest principal radius of gyration
435 0 : rgyr=std::sqrt((princ_comp[0]+princ_comp[1])/totmass);
436 0 : double rm = rgyr*totmass;
437 0 : if (rm>1e-6) {
438 0 : prefactor[0]=1.0/rm;
439 0 : prefactor[1]=1.0/rm;
440 : }
441 : break;
442 : }
443 5 : case ASPHERICITY: {
444 5 : rgyr=std::sqrt((princ_comp[0]-0.5*(princ_comp[1]+princ_comp[2]))/totmass);
445 5 : double rm = rgyr*totmass;
446 5 : if (rm>1e-6) {
447 5 : prefactor[0]= 1.0/rm;
448 5 : prefactor[1]=-0.5/rm;
449 5 : prefactor[2]=-0.5/rm;
450 : }
451 : break;
452 : }
453 0 : case ACYLINDRICITY: {
454 0 : rgyr=std::sqrt((princ_comp[1]-princ_comp[2])/totmass);
455 0 : double rm = rgyr*totmass;
456 0 : if (rm>1e-6) { //avoid division by zero
457 0 : prefactor[1]= 1.0/rm;
458 0 : prefactor[2]=-1.0/rm;
459 : }
460 : break;
461 : }
462 130 : case KAPPA2: { // relative shape anisotropy
463 130 : double trace = princ_comp[0]+princ_comp[1]+princ_comp[2];
464 130 : double tmp=princ_comp[0]*princ_comp[1]+ princ_comp[1]*princ_comp[2]+ princ_comp[0]*princ_comp[2];
465 130 : rgyr=1.0-3*(tmp/(trace*trace));
466 130 : if (rgyr>1e-6) {
467 130 : prefactor[0]= -3*((princ_comp[1]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
468 130 : prefactor[1]= -3*((princ_comp[0]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
469 130 : prefactor[2]= -3*((princ_comp[0]+princ_comp[1])-2*tmp/trace)/(trace*trace) *2;
470 : }
471 : break;
472 : }
473 : }
474 :
475 400 : if(use_masses) {
476 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
477 0 : Vector tX;
478 0 : const Vector diff=delta( com,getPosition(i) );
479 : //project atomic postional vectors to diagonalized frame
480 0 : for(unsigned j=0; j<3; j++) {
481 0 : tX[j]=transf[0][j]*diff[0]+transf[1][j]*diff[1]+transf[2][j]*diff[2];
482 : }
483 0 : for(unsigned j=0; j<3; j++)
484 0 : derivatives[i][j]=getMass(i)*(prefactor[0]*transf[j][0]*tX[0]+
485 0 : prefactor[1]*transf[j][1]*tX[1]+
486 0 : prefactor[2]*transf[j][2]*tX[2]);
487 0 : setAtomsDerivatives(i,derivatives[i]);
488 : }
489 : } else {
490 2400 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
491 2000 : Vector tX;
492 2000 : const Vector diff=delta( com,getPosition(i) );
493 : //project atomic postional vectors to diagonalized frame
494 8000 : for(unsigned j=0; j<3; j++) {
495 6000 : tX[j]=transf[0][j]*diff[0]+transf[1][j]*diff[1]+transf[2][j]*diff[2];
496 : }
497 8000 : for(unsigned j=0; j<3; j++)
498 6000 : derivatives[i][j]=prefactor[0]*transf[j][0]*tX[0]+
499 6000 : prefactor[1]*transf[j][1]*tX[1]+
500 6000 : prefactor[2]*transf[j][2]*tX[2];
501 2000 : setAtomsDerivatives(i,derivatives[i]);
502 : }
503 : }
504 :
505 400 : setValue(rgyr);
506 400 : setBoxDerivativesNoPbc();
507 : }
508 :
509 : }
510 : }
|
