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