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 "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 different properties can be calculated and selected by the TYPE keyword:
34 : the Radius of Gyration (RADIUS); the Trace of the Gyration Tensor (TRACE);
35 : the Largest Principal Moment of the Gyration Tensor (GTPC_1); the middle Principal Moment of the Gyration Tensor (GTPC_2);
36 : the Smallest Principal Moment of the Gyration Tensor (GTPC_3); the Asphericiry (ASPHERICITY); the Acylindricity (ACYLINDRICITY);
37 : the Relative Shape Anisotropy (KAPPA2); the Smallest Principal Radius Of Gyration (GYRATION_3);
38 : the Middle Principal Radius of Gyration (GYRATION_2); the Largest Principal Radius of Gyration (GYRATION_1).
39 : A derivation of all these different variants can be found in \cite Vymetal:2011gv
40 :
41 : The radius of gyration is calculated using:
42 :
43 : \f[
44 : s_{\rm Gyr}=\Big ( \frac{\sum_i^{n}
45 : m_i \vert {r}_i -{r}_{\rm COM} \vert ^2 }{\sum_i^{n} m_i} \Big)^{1/2}
46 : \f]
47 :
48 : with the position of the center of mass \f${r}_{\rm COM}\f$ given by:
49 :
50 : \f[
51 : {r}_{\rm COM}=\frac{\sum_i^{n} {r}_i\ m_i }{\sum_i^{n} m_i}
52 : \f]
53 :
54 : The radius of gyration usually makes sense when atoms used for the calculation
55 : are all part of the same molecule.
56 : When running with periodic boundary conditions, the atoms should be
57 : in the proper periodic image. This is done automatically since PLUMED 2.2,
58 : by considering the ordered list of atoms and rebuilding the broken entities using a procedure
59 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
60 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
61 : which actually modifies the coordinates stored in PLUMED.
62 :
63 : In case you want to recover the old behavior you should use the NOPBC flag.
64 : In that case you need to take care that atoms are in the correct
65 : periodic image.
66 :
67 :
68 : \par Examples
69 :
70 : The following input tells plumed to print the radius of gyration of the
71 : chain containing atoms 10 to 20.
72 : \plumedfile
73 : GYRATION TYPE=RADIUS ATOMS=10-20 LABEL=rg
74 : PRINT ARG=rg STRIDE=1 FILE=colvar
75 : \endplumedfile
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 : class Gyration : public Colvar {
81 : private:
82 : enum CV_TYPE {RADIUS, TRACE, GTPC_1, GTPC_2, GTPC_3, ASPHERICITY, ACYLINDRICITY, KAPPA2, GYRATION_3, GYRATION_2, GYRATION_1, TOT};
83 : int rg_type;
84 : bool use_masses;
85 : bool nopbc;
86 : public:
87 : static void registerKeywords(Keywords& keys);
88 : explicit Gyration(const ActionOptions&);
89 : void calculate() override;
90 : };
91 :
92 13849 : PLUMED_REGISTER_ACTION(Gyration,"GYRATION")
93 :
94 37 : void Gyration::registerKeywords(Keywords& keys) {
95 37 : Colvar::registerKeywords(keys);
96 74 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
97 74 : keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
98 74 : keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
99 37 : }
100 :
101 33 : Gyration::Gyration(const ActionOptions&ao):
102 : PLUMED_COLVAR_INIT(ao),
103 33 : use_masses(false),
104 33 : nopbc(false) {
105 : std::vector<AtomNumber> atoms;
106 65 : parseAtomList("ATOMS",atoms);
107 32 : if(atoms.size()==0) {
108 0 : error("no atoms specified");
109 : }
110 66 : parseFlag("MASS_WEIGHTED",use_masses);
111 : std::string Type;
112 32 : parse("TYPE",Type);
113 32 : parseFlag("NOPBC",nopbc);
114 32 : checkRead();
115 :
116 32 : if(Type=="RADIUS") {
117 11 : rg_type=RADIUS;
118 21 : } else if(Type=="TRACE") {
119 2 : rg_type=TRACE;
120 19 : } else if(Type=="GTPC_1") {
121 2 : rg_type=GTPC_1;
122 17 : } else if(Type=="GTPC_2") {
123 2 : rg_type=GTPC_2;
124 15 : } else if(Type=="GTPC_3") {
125 2 : rg_type=GTPC_3;
126 13 : } else if(Type=="ASPHERICITY") {
127 2 : rg_type=ASPHERICITY;
128 11 : } else if(Type=="ACYLINDRICITY") {
129 2 : rg_type=ACYLINDRICITY;
130 9 : } else if(Type=="KAPPA2") {
131 2 : rg_type=KAPPA2;
132 7 : } else if(Type=="RGYR_3") {
133 2 : rg_type=GYRATION_3;
134 5 : } else if(Type=="RGYR_2") {
135 2 : rg_type=GYRATION_2;
136 3 : } else if(Type=="RGYR_1") {
137 2 : rg_type=GYRATION_1;
138 : } else {
139 1 : error("Unknown GYRATION type");
140 : }
141 :
142 31 : switch(rg_type) {
143 11 : case RADIUS:
144 11 : log.printf(" GYRATION RADIUS (Rg);");
145 : break;
146 2 : case TRACE:
147 2 : log.printf(" TRACE OF THE GYRATION TENSOR;");
148 : break;
149 2 : case GTPC_1:
150 2 : log.printf(" THE LARGEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_1);");
151 : break;
152 2 : case GTPC_2:
153 2 : log.printf(" THE MIDDLE PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_2);");
154 : break;
155 2 : case GTPC_3:
156 2 : log.printf(" THE SMALLEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_3);");
157 : break;
158 2 : case ASPHERICITY:
159 2 : log.printf(" THE ASPHERICITY (b');");
160 : break;
161 2 : case ACYLINDRICITY:
162 2 : log.printf(" THE ACYLINDRICITY (c');");
163 : break;
164 2 : case KAPPA2:
165 2 : log.printf(" THE RELATIVE SHAPE ANISOTROPY (kappa^2);");
166 : break;
167 2 : case GYRATION_3:
168 2 : log.printf(" THE SMALLEST PRINCIPAL RADIUS OF GYRATION (r_g3);");
169 : break;
170 2 : case GYRATION_2:
171 2 : log.printf(" THE MIDDLE PRINCIPAL RADIUS OF GYRATION (r_g2);");
172 : break;
173 2 : case GYRATION_1:
174 2 : log.printf(" THE LARGEST PRINCIPAL RADIUS OF GYRATION (r_g1);");
175 : break;
176 : }
177 31 : if(rg_type>TRACE) {
178 37 : log<<" Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)");
179 : }
180 31 : log<<"\n";
181 :
182 31 : log.printf(" atoms involved : ");
183 233 : for(unsigned i=0; i<atoms.size(); ++i) {
184 202 : log.printf("%d ",atoms[i].serial());
185 : }
186 31 : log.printf("\n");
187 :
188 31 : if(nopbc) {
189 4 : log<<" PBC will be ignored\n";
190 : } else {
191 27 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
192 : }
193 :
194 31 : addValueWithDerivatives();
195 31 : setNotPeriodic();
196 31 : requestAtoms(atoms);
197 35 : }
198 :
199 1195 : void Gyration::calculate() {
200 :
201 1195 : if(!nopbc) {
202 537 : makeWhole();
203 : }
204 :
205 1195 : Vector com;
206 : double totmass = 0.;
207 1195 : if( use_masses ) {
208 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
209 0 : totmass+=getMass(i);
210 0 : com+=getMass(i)*getPosition(i);
211 : }
212 : } else {
213 1195 : totmass = static_cast<double>(getNumberOfAtoms());
214 10373 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
215 9178 : com+=getPosition(i);
216 : }
217 : }
218 1195 : com /= totmass;
219 :
220 1195 : double rgyr=0.;
221 1195 : std::vector<Vector> derivatives( getNumberOfAtoms() );
222 1195 : Tensor virial;
223 :
224 1195 : if(rg_type==RADIUS||rg_type==TRACE) {
225 795 : if( use_masses ) {
226 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
227 0 : const Vector diff = delta( com, getPosition(i) );
228 0 : rgyr += getMass(i)*diff.modulo2();
229 0 : derivatives[i] = diff*getMass(i);
230 0 : virial -= Tensor(getPosition(i),derivatives[i]);
231 : }
232 : } else {
233 7973 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
234 7178 : const Vector diff = delta( com, getPosition(i) );
235 7178 : rgyr += diff.modulo2();
236 7178 : derivatives[i] = diff;
237 7178 : virial -= Tensor(getPosition(i),derivatives[i]);
238 : }
239 : }
240 : double fact;
241 795 : if(rg_type==RADIUS) {
242 665 : rgyr = std::sqrt(rgyr/totmass);
243 665 : fact = 1./(rgyr*totmass);
244 : } else {
245 130 : rgyr = 2.*rgyr;
246 : fact = 4;
247 : }
248 795 : setValue(rgyr);
249 7973 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
250 7178 : setAtomsDerivatives(i,fact*derivatives[i]);
251 : }
252 795 : setBoxDerivatives(fact*virial);
253 : return;
254 : }
255 :
256 :
257 400 : Tensor3d gyr_tens;
258 : //calculate gyration tensor
259 400 : if( use_masses ) {
260 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
261 0 : const Vector diff=delta( com, getPosition(i) );
262 0 : gyr_tens[0][0]+=getMass(i)*diff[0]*diff[0];
263 0 : gyr_tens[1][1]+=getMass(i)*diff[1]*diff[1];
264 0 : gyr_tens[2][2]+=getMass(i)*diff[2]*diff[2];
265 0 : gyr_tens[0][1]+=getMass(i)*diff[0]*diff[1];
266 0 : gyr_tens[0][2]+=getMass(i)*diff[0]*diff[2];
267 0 : gyr_tens[1][2]+=getMass(i)*diff[1]*diff[2];
268 : }
269 : } else {
270 2400 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
271 2000 : const Vector diff=delta( com, getPosition(i) );
272 2000 : gyr_tens[0][0]+=diff[0]*diff[0];
273 2000 : gyr_tens[1][1]+=diff[1]*diff[1];
274 2000 : gyr_tens[2][2]+=diff[2]*diff[2];
275 2000 : gyr_tens[0][1]+=diff[0]*diff[1];
276 2000 : gyr_tens[0][2]+=diff[0]*diff[2];
277 2000 : gyr_tens[1][2]+=diff[1]*diff[2];
278 : }
279 : }
280 :
281 : // first make the matrix symmetric
282 400 : gyr_tens[1][0] = gyr_tens[0][1];
283 400 : gyr_tens[2][0] = gyr_tens[0][2];
284 400 : gyr_tens[2][1] = gyr_tens[1][2];
285 400 : Tensor3d ttransf,transf;
286 400 : Vector princ_comp,prefactor;
287 : //diagonalize gyration tensor
288 400 : diagMatSym(gyr_tens, princ_comp, ttransf);
289 400 : transf=transpose(ttransf);
290 : //sort eigenvalues and eigenvectors
291 400 : if (princ_comp[0]<princ_comp[1]) {
292 : double tmp=princ_comp[0];
293 400 : princ_comp[0]=princ_comp[1];
294 400 : princ_comp[1]=tmp;
295 1600 : for (unsigned i=0; i<3; i++) {
296 1200 : tmp=transf[i][0];
297 1200 : transf[i][0]=transf[i][1];
298 1200 : transf[i][1]=tmp;
299 : }
300 : }
301 400 : if (princ_comp[1]<princ_comp[2]) {
302 : double tmp=princ_comp[1];
303 400 : princ_comp[1]=princ_comp[2];
304 400 : princ_comp[2]=tmp;
305 1600 : for (unsigned i=0; i<3; i++) {
306 1200 : tmp=transf[i][1];
307 1200 : transf[i][1]=transf[i][2];
308 1200 : transf[i][2]=tmp;
309 : }
310 : }
311 400 : if (princ_comp[0]<princ_comp[1]) {
312 : double tmp=princ_comp[0];
313 400 : princ_comp[0]=princ_comp[1];
314 400 : princ_comp[1]=tmp;
315 1600 : for (unsigned i=0; i<3; i++) {
316 1200 : tmp=transf[i][0];
317 1200 : transf[i][0]=transf[i][1];
318 1200 : transf[i][1]=tmp;
319 : }
320 : }
321 : //calculate determinant of transformation matrix
322 : double det = determinant(transf);
323 : // transformation matrix for rotation must have positive determinant, otherwise multiply one column by (-1)
324 400 : if(det<0) {
325 1600 : for(unsigned j=0; j<3; j++) {
326 1200 : transf[j][2]=-transf[j][2];
327 : }
328 400 : det = -det;
329 : }
330 400 : if(std::abs(det-1.)>0.0001) {
331 0 : error("Plumed Error: Cannot diagonalize gyration tensor\n");
332 : }
333 400 : switch(rg_type) {
334 135 : case GTPC_1:
335 : case GTPC_2:
336 : case GTPC_3: {
337 135 : int pc_index = rg_type-2; //index of principal component
338 135 : rgyr=std::sqrt(princ_comp[pc_index]/totmass);
339 135 : double rm = rgyr*totmass;
340 135 : if(rm>1e-6) {
341 135 : prefactor[pc_index]=1.0/rm; //some parts of derivate
342 : }
343 : break;
344 : }
345 0 : case GYRATION_3: { //the smallest principal radius of gyration
346 0 : rgyr=std::sqrt((princ_comp[1]+princ_comp[2])/totmass);
347 0 : double rm = rgyr*totmass;
348 0 : if (rm>1e-6) {
349 0 : prefactor[1]=1.0/rm;
350 0 : prefactor[2]=1.0/rm;
351 : }
352 : break;
353 : }
354 130 : case GYRATION_2: { //the midle principal radius of gyration
355 130 : rgyr=std::sqrt((princ_comp[0]+princ_comp[2])/totmass);
356 130 : double rm = rgyr*totmass;
357 130 : if (rm>1e-6) {
358 130 : prefactor[0]=1.0/rm;
359 130 : prefactor[2]=1.0/rm;
360 : }
361 : break;
362 : }
363 0 : case GYRATION_1: { //the largest principal radius of gyration
364 0 : rgyr=std::sqrt((princ_comp[0]+princ_comp[1])/totmass);
365 0 : double rm = rgyr*totmass;
366 0 : if (rm>1e-6) {
367 0 : prefactor[0]=1.0/rm;
368 0 : prefactor[1]=1.0/rm;
369 : }
370 : break;
371 : }
372 5 : case ASPHERICITY: {
373 5 : rgyr=std::sqrt((princ_comp[0]-0.5*(princ_comp[1]+princ_comp[2]))/totmass);
374 5 : double rm = rgyr*totmass;
375 5 : if (rm>1e-6) {
376 5 : prefactor[0]= 1.0/rm;
377 5 : prefactor[1]=-0.5/rm;
378 5 : prefactor[2]=-0.5/rm;
379 : }
380 : break;
381 : }
382 0 : case ACYLINDRICITY: {
383 0 : rgyr=std::sqrt((princ_comp[1]-princ_comp[2])/totmass);
384 0 : double rm = rgyr*totmass;
385 0 : if (rm>1e-6) { //avoid division by zero
386 0 : prefactor[1]= 1.0/rm;
387 0 : prefactor[2]=-1.0/rm;
388 : }
389 : break;
390 : }
391 130 : case KAPPA2: { // relative shape anisotropy
392 130 : double trace = princ_comp[0]+princ_comp[1]+princ_comp[2];
393 130 : double tmp=princ_comp[0]*princ_comp[1]+ princ_comp[1]*princ_comp[2]+ princ_comp[0]*princ_comp[2];
394 130 : rgyr=1.0-3*(tmp/(trace*trace));
395 130 : if (rgyr>1e-6) {
396 130 : prefactor[0]= -3*((princ_comp[1]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
397 130 : prefactor[1]= -3*((princ_comp[0]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
398 130 : prefactor[2]= -3*((princ_comp[0]+princ_comp[1])-2*tmp/trace)/(trace*trace) *2;
399 : }
400 : break;
401 : }
402 : }
403 :
404 400 : if(use_masses) {
405 0 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
406 0 : Vector tX;
407 0 : const Vector diff=delta( com,getPosition(i) );
408 : //project atomic postional vectors to diagonalized frame
409 0 : for(unsigned j=0; j<3; j++) {
410 0 : tX[j]=transf[0][j]*diff[0]+transf[1][j]*diff[1]+transf[2][j]*diff[2];
411 : }
412 0 : for(unsigned j=0; j<3; j++)
413 0 : derivatives[i][j]=getMass(i)*(prefactor[0]*transf[j][0]*tX[0]+
414 0 : prefactor[1]*transf[j][1]*tX[1]+
415 0 : prefactor[2]*transf[j][2]*tX[2]);
416 0 : setAtomsDerivatives(i,derivatives[i]);
417 : }
418 : } else {
419 2400 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
420 2000 : Vector tX;
421 2000 : const Vector diff=delta( com,getPosition(i) );
422 : //project atomic postional vectors to diagonalized frame
423 8000 : for(unsigned j=0; j<3; j++) {
424 6000 : tX[j]=transf[0][j]*diff[0]+transf[1][j]*diff[1]+transf[2][j]*diff[2];
425 : }
426 8000 : for(unsigned j=0; j<3; j++)
427 6000 : derivatives[i][j]=prefactor[0]*transf[j][0]*tX[0]+
428 12000 : prefactor[1]*transf[j][1]*tX[1]+
429 6000 : prefactor[2]*transf[j][2]*tX[2];
430 2000 : setAtomsDerivatives(i,derivatives[i]);
431 : }
432 : }
433 :
434 400 : setValue(rgyr);
435 400 : setBoxDerivativesNoPbc();
436 : }
437 :
438 : }
439 : }
|
