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       10483 : PLUMED_REGISTER_ACTION(Gyration,"GYRATION")
      93             : 
      94          34 : void Gyration::registerKeywords(Keywords& keys) {
      95          34 :   Colvar::registerKeywords(keys);
      96          68 :   keys.add("atoms","ATOMS","the group of atoms that you are calculating the Gyration Tensor for");
      97          68 :   keys.add("compulsory","TYPE","RADIUS","The type of calculation relative to the Gyration Tensor you want to perform");
      98          68 :   keys.addFlag("MASS_WEIGHTED",false,"set the masses of all the atoms equal to one");
      99          34 : }
     100             : 
     101          33 : Gyration::Gyration(const ActionOptions&ao):
     102             :   PLUMED_COLVAR_INIT(ao),
     103          33 :   use_masses(false),
     104          33 :   nopbc(false)
     105             : {
     106             :   std::vector<AtomNumber> atoms;
     107          65 :   parseAtomList("ATOMS",atoms);
     108          32 :   if(atoms.size()==0) error("no atoms specified");
     109          66 :   parseFlag("MASS_WEIGHTED",use_masses);
     110             :   std::string Type;
     111          32 :   parse("TYPE",Type);
     112          32 :   parseFlag("NOPBC",nopbc);
     113          32 :   checkRead();
     114             : 
     115          32 :   if(Type=="RADIUS") rg_type=RADIUS;
     116          21 :   else if(Type=="TRACE") rg_type=TRACE;
     117          19 :   else if(Type=="GTPC_1") rg_type=GTPC_1;
     118          17 :   else if(Type=="GTPC_2") rg_type=GTPC_2;
     119          15 :   else if(Type=="GTPC_3") rg_type=GTPC_3;
     120          13 :   else if(Type=="ASPHERICITY") rg_type=ASPHERICITY;
     121          11 :   else if(Type=="ACYLINDRICITY") rg_type=ACYLINDRICITY;
     122           9 :   else if(Type=="KAPPA2") rg_type=KAPPA2;
     123           7 :   else if(Type=="RGYR_3") rg_type=GYRATION_3;
     124           5 :   else if(Type=="RGYR_2") rg_type=GYRATION_2;
     125           3 :   else if(Type=="RGYR_1") rg_type=GYRATION_1;
     126           1 :   else error("Unknown GYRATION type");
     127             : 
     128          31 :   switch(rg_type)
     129             :   {
     130          11 :   case RADIUS:   log.printf("  GYRATION RADIUS (Rg);"); break;
     131           2 :   case TRACE:  log.printf("  TRACE OF THE GYRATION TENSOR;"); break;
     132           2 :   case GTPC_1: log.printf("  THE LARGEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_1);"); break;
     133           2 :   case GTPC_2: log.printf("  THE MIDDLE PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_2);");  break;
     134           2 :   case GTPC_3: log.printf("  THE SMALLEST PRINCIPAL MOMENT OF THE GYRATION TENSOR (S'_3);"); break;
     135           2 :   case ASPHERICITY: log.printf("  THE ASPHERICITY (b');"); break;
     136           2 :   case ACYLINDRICITY: log.printf("  THE ACYLINDRICITY (c');"); break;
     137           2 :   case KAPPA2: log.printf("  THE RELATIVE SHAPE ANISOTROPY (kappa^2);"); break;
     138           2 :   case GYRATION_3: log.printf("  THE SMALLEST PRINCIPAL RADIUS OF GYRATION (r_g3);"); break;
     139           2 :   case GYRATION_2: log.printf("  THE MIDDLE PRINCIPAL RADIUS OF GYRATION (r_g2);"); break;
     140           2 :   case GYRATION_1: log.printf("  THE LARGEST PRINCIPAL RADIUS OF GYRATION (r_g1);"); break;
     141             :   }
     142          68 :   if(rg_type>TRACE) log<<"  Bibliography "<<plumed.cite("Jirí Vymetal and Jirí Vondrasek, J. Phys. Chem. A 115, 11455 (2011)");
     143          31 :   log<<"\n";
     144             : 
     145          31 :   log.printf("  atoms involved : ");
     146         233 :   for(unsigned i=0; i<atoms.size(); ++i) log.printf("%d ",atoms[i].serial());
     147          31 :   log.printf("\n");
     148             : 
     149          31 :   if(nopbc) {
     150           4 :     log<<"  PBC will be ignored\n";
     151             :   } else {
     152          27 :     log<<"  broken molecules will be rebuilt assuming atoms are in the proper order\n";
     153             :   }
     154             : 
     155          31 :   addValueWithDerivatives(); setNotPeriodic();
     156          31 :   requestAtoms(atoms);
     157          35 : }
     158             : 
     159        1195 : void Gyration::calculate() {
     160             : 
     161        1195 :   if(!nopbc) makeWhole();
     162             : 
     163        1195 :   Vector com;
     164             :   double totmass = 0.;
     165        1195 :   if( use_masses ) {
     166           0 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     167           0 :       totmass+=getMass(i);
     168           0 :       com+=getMass(i)*getPosition(i);
     169             :     }
     170             :   } else {
     171        1195 :     totmass = static_cast<double>(getNumberOfAtoms());
     172       10373 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     173        9178 :       com+=getPosition(i);
     174             :     }
     175             :   }
     176        1195 :   com /= totmass;
     177             : 
     178        1195 :   double rgyr=0.;
     179        1195 :   std::vector<Vector> derivatives( getNumberOfAtoms() );
     180        1195 :   Tensor virial;
     181             : 
     182        1195 :   if(rg_type==RADIUS||rg_type==TRACE) {
     183         795 :     if( use_masses ) {
     184           0 :       for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     185           0 :         const Vector diff = delta( com, getPosition(i) );
     186           0 :         rgyr          += getMass(i)*diff.modulo2();
     187           0 :         derivatives[i] = diff*getMass(i);
     188           0 :         virial        -= Tensor(getPosition(i),derivatives[i]);
     189             :       }
     190             :     } else {
     191        7973 :       for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     192        7178 :         const Vector diff = delta( com, getPosition(i) );
     193        7178 :         rgyr          += diff.modulo2();
     194        7178 :         derivatives[i] = diff;
     195        7178 :         virial        -= Tensor(getPosition(i),derivatives[i]);
     196             :       }
     197             :     }
     198             :     double fact;
     199         795 :     if(rg_type==RADIUS) {
     200         665 :       rgyr = std::sqrt(rgyr/totmass);
     201         665 :       fact = 1./(rgyr*totmass);
     202             :     } else {
     203         130 :       rgyr = 2.*rgyr;
     204             :       fact = 4;
     205             :     }
     206         795 :     setValue(rgyr);
     207        7973 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives(i,fact*derivatives[i]);
     208         795 :     setBoxDerivatives(fact*virial);
     209             :     return;
     210             :   }
     211             : 
     212             : 
     213         400 :   Tensor3d gyr_tens;
     214             :   //calculate gyration tensor
     215         400 :   if( use_masses ) {
     216           0 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     217           0 :       const Vector diff=delta( com, getPosition(i) );
     218           0 :       gyr_tens[0][0]+=getMass(i)*diff[0]*diff[0];
     219           0 :       gyr_tens[1][1]+=getMass(i)*diff[1]*diff[1];
     220           0 :       gyr_tens[2][2]+=getMass(i)*diff[2]*diff[2];
     221           0 :       gyr_tens[0][1]+=getMass(i)*diff[0]*diff[1];
     222           0 :       gyr_tens[0][2]+=getMass(i)*diff[0]*diff[2];
     223           0 :       gyr_tens[1][2]+=getMass(i)*diff[1]*diff[2];
     224             :     }
     225             :   } else {
     226        2400 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     227        2000 :       const Vector diff=delta( com, getPosition(i) );
     228        2000 :       gyr_tens[0][0]+=diff[0]*diff[0];
     229        2000 :       gyr_tens[1][1]+=diff[1]*diff[1];
     230        2000 :       gyr_tens[2][2]+=diff[2]*diff[2];
     231        2000 :       gyr_tens[0][1]+=diff[0]*diff[1];
     232        2000 :       gyr_tens[0][2]+=diff[0]*diff[2];
     233        2000 :       gyr_tens[1][2]+=diff[1]*diff[2];
     234             :     }
     235             :   }
     236             : 
     237             :   // first make the matrix symmetric
     238         400 :   gyr_tens[1][0] = gyr_tens[0][1];
     239         400 :   gyr_tens[2][0] = gyr_tens[0][2];
     240         400 :   gyr_tens[2][1] = gyr_tens[1][2];
     241         400 :   Tensor3d ttransf,transf;
     242         400 :   Vector princ_comp,prefactor;
     243             :   //diagonalize gyration tensor
     244         400 :   diagMatSym(gyr_tens, princ_comp, ttransf);
     245         400 :   transf=transpose(ttransf);
     246             :   //sort eigenvalues and eigenvectors
     247         400 :   if (princ_comp[0]<princ_comp[1]) {
     248         400 :     double tmp=princ_comp[0]; princ_comp[0]=princ_comp[1]; princ_comp[1]=tmp;
     249        1600 :     for (unsigned i=0; i<3; i++) {tmp=transf[i][0]; transf[i][0]=transf[i][1]; transf[i][1]=tmp;}
     250             :   }
     251         400 :   if (princ_comp[1]<princ_comp[2]) {
     252         400 :     double tmp=princ_comp[1]; princ_comp[1]=princ_comp[2]; princ_comp[2]=tmp;
     253        1600 :     for (unsigned i=0; i<3; i++) {tmp=transf[i][1]; transf[i][1]=transf[i][2]; transf[i][2]=tmp;}
     254             :   }
     255         400 :   if (princ_comp[0]<princ_comp[1]) {
     256         400 :     double tmp=princ_comp[0]; princ_comp[0]=princ_comp[1]; princ_comp[1]=tmp;
     257        1600 :     for (unsigned i=0; i<3; i++) {tmp=transf[i][0]; transf[i][0]=transf[i][1]; transf[i][1]=tmp;}
     258             :   }
     259             :   //calculate determinant of transformation matrix
     260             :   double det = determinant(transf);
     261             :   // transformation matrix for rotation must have positive determinant, otherwise multiply one column by (-1)
     262         400 :   if(det<0) {
     263        1600 :     for(unsigned j=0; j<3; j++) transf[j][2]=-transf[j][2];
     264         400 :     det = -det;
     265             :   }
     266         400 :   if(std::abs(det-1.)>0.0001) error("Plumed Error: Cannot diagonalize gyration tensor\n");
     267         400 :   switch(rg_type) {
     268         135 :   case GTPC_1:
     269             :   case GTPC_2:
     270             :   case GTPC_3:
     271             :   {
     272         135 :     int pc_index = rg_type-2; //index of principal component
     273         135 :     rgyr=std::sqrt(princ_comp[pc_index]/totmass);
     274         135 :     double rm = rgyr*totmass;
     275         135 :     if(rm>1e-6) prefactor[pc_index]=1.0/rm; //some parts of derivate
     276             :     break;
     277             :   }
     278           0 :   case GYRATION_3:        //the smallest principal radius of gyration
     279             :   {
     280           0 :     rgyr=std::sqrt((princ_comp[1]+princ_comp[2])/totmass);
     281           0 :     double rm = rgyr*totmass;
     282           0 :     if (rm>1e-6) {
     283           0 :       prefactor[1]=1.0/rm;
     284           0 :       prefactor[2]=1.0/rm;
     285             :     }
     286             :     break;
     287             :   }
     288         130 :   case GYRATION_2:       //the midle principal radius of gyration
     289             :   {
     290         130 :     rgyr=std::sqrt((princ_comp[0]+princ_comp[2])/totmass);
     291         130 :     double rm = rgyr*totmass;
     292         130 :     if (rm>1e-6) {
     293         130 :       prefactor[0]=1.0/rm;
     294         130 :       prefactor[2]=1.0/rm;
     295             :     }
     296             :     break;
     297             :   }
     298           0 :   case GYRATION_1:      //the largest principal radius of gyration
     299             :   {
     300           0 :     rgyr=std::sqrt((princ_comp[0]+princ_comp[1])/totmass);
     301           0 :     double rm = rgyr*totmass;
     302           0 :     if (rm>1e-6) {
     303           0 :       prefactor[0]=1.0/rm;
     304           0 :       prefactor[1]=1.0/rm;
     305             :     }
     306             :     break;
     307             :   }
     308           5 :   case ASPHERICITY:
     309             :   {
     310           5 :     rgyr=std::sqrt((princ_comp[0]-0.5*(princ_comp[1]+princ_comp[2]))/totmass);
     311           5 :     double rm = rgyr*totmass;
     312           5 :     if (rm>1e-6) {
     313           5 :       prefactor[0]= 1.0/rm;
     314           5 :       prefactor[1]=-0.5/rm;
     315           5 :       prefactor[2]=-0.5/rm;
     316             :     }
     317             :     break;
     318             :   }
     319           0 :   case ACYLINDRICITY:
     320             :   {
     321           0 :     rgyr=std::sqrt((princ_comp[1]-princ_comp[2])/totmass);
     322           0 :     double rm = rgyr*totmass;
     323           0 :     if (rm>1e-6) {  //avoid division by zero
     324           0 :       prefactor[1]= 1.0/rm;
     325           0 :       prefactor[2]=-1.0/rm;
     326             :     }
     327             :     break;
     328             :   }
     329         130 :   case KAPPA2: // relative shape anisotropy
     330             :   {
     331         130 :     double trace = princ_comp[0]+princ_comp[1]+princ_comp[2];
     332         130 :     double tmp=princ_comp[0]*princ_comp[1]+ princ_comp[1]*princ_comp[2]+ princ_comp[0]*princ_comp[2];
     333         130 :     rgyr=1.0-3*(tmp/(trace*trace));
     334         130 :     if (rgyr>1e-6) {
     335         130 :       prefactor[0]= -3*((princ_comp[1]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
     336         130 :       prefactor[1]= -3*((princ_comp[0]+princ_comp[2])-2*tmp/trace)/(trace*trace) *2;
     337         130 :       prefactor[2]= -3*((princ_comp[0]+princ_comp[1])-2*tmp/trace)/(trace*trace) *2;
     338             :     }
     339             :     break;
     340             :   }
     341             :   }
     342             : 
     343         400 :   if(use_masses) {
     344           0 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     345           0 :       Vector tX;
     346           0 :       const Vector diff=delta( com,getPosition(i) );
     347             :       //project atomic postional vectors to diagonalized frame
     348           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];
     349           0 :       for(unsigned j=0; j<3; j++) derivatives[i][j]=getMass(i)*(prefactor[0]*transf[j][0]*tX[0]+
     350           0 :             prefactor[1]*transf[j][1]*tX[1]+
     351           0 :             prefactor[2]*transf[j][2]*tX[2]);
     352           0 :       setAtomsDerivatives(i,derivatives[i]);
     353             :     }
     354             :   } else {
     355        2400 :     for(unsigned i=0; i<getNumberOfAtoms(); i++) {
     356        2000 :       Vector tX;
     357        2000 :       const Vector diff=delta( com,getPosition(i) );
     358             :       //project atomic postional vectors to diagonalized frame
     359        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];
     360        8000 :       for(unsigned j=0; j<3; j++) derivatives[i][j]=prefactor[0]*transf[j][0]*tX[0]+
     361       12000 :             prefactor[1]*transf[j][1]*tX[1]+
     362        6000 :             prefactor[2]*transf[j][2]*tX[2];
     363        2000 :       setAtomsDerivatives(i,derivatives[i]);
     364             :     }
     365             :   }
     366             : 
     367         400 :   setValue(rgyr);
     368         400 :   setBoxDerivativesNoPbc();
     369             : }
     370             : 
     371             : }
     372             : }