Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-2020 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 "SecondaryStructureRMSD.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionSet.h"
      25             : #include "core/GenericMolInfo.h"
      26             : #include "core/ActionShortcut.h"
      27             : #include "core/ActionRegister.h"
      28             : 
      29             : //+PLUMEDOC MCOLVAR SECONDARY_STRUCTURE_RMSD
      30             : /*
      31             : Calclulate the distance between segments of a protein and a reference structure of interest
      32             : 
      33             : \par Examples
      34             : 
      35             : */
      36             : //+ENDPLUMEDOC
      37             : 
      38             : namespace PLMD {
      39             : namespace secondarystructure {
      40             : 
      41             : PLUMED_REGISTER_ACTION(SecondaryStructureRMSD,"SECONDARY_STRUCTURE_RMSD");
      42             : 
      43          32 : bool SecondaryStructureRMSD::readShortcutWords( std::string& ltmap, ActionShortcut* action ) {
      44          64 :   action->parse("LESS_THAN",ltmap);
      45          32 :   if( ltmap.length()==0 ) {
      46           6 :     std::string nn, mm, d_0, r_0; action->parse("R_0",r_0);
      47           3 :     if( r_0.length()==0 ) r_0="0.08";
      48           9 :     action->parse("NN",nn); action->parse("D_0",d_0); action->parse("MM",mm);
      49           6 :     ltmap = "RATIONAL R_0=" + r_0 + " D_0=" + d_0 + " NN=" + nn + " MM=" + mm;
      50             :     return false;
      51             :   }
      52             :   return true;
      53             : }
      54             : 
      55          32 : void SecondaryStructureRMSD::expandShortcut( const bool& uselessthan, const std::string& labout, const std::string& labin, const std::string& ltmap, ActionShortcut* action ) {
      56          64 :   action->readInputLine( labout + "_lt: LESS_THAN ARG=" + labin + " SWITCH={" + ltmap  +"}");
      57          61 :   if( uselessthan ) action->readInputLine( labout + "_lessthan: SUM ARG=" + labout + "_lt PERIODIC=NO");
      58           6 :   else action->readInputLine( labout + ": SUM ARG=" + labout + "_lt PERIODIC=NO");
      59          32 : }
      60             : 
      61         242 : void SecondaryStructureRMSD::registerKeywords( Keywords& keys ) {
      62         242 :   ActionWithVector::registerKeywords( keys );
      63         484 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
      64         484 :   keys.add("residues","RESIDUES","this command is used to specify the set of residues that could conceivably form part of the secondary structure. "
      65             :            "It is possible to use residues numbers as the various chains and residues should have been identified else using an instance of the "
      66             :            "\\ref MOLINFO action. If you wish to use all the residues from all the chains in your system you can do so by "
      67             :            "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
      68             :            "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain "
      69             :            "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. "
      70             :            "As such if you define portions of the chain with fewer than N residues the code will crash.");
      71         484 :   keys.add("atoms","ATOMS","this is the full list of atoms that we are investigating");
      72         484 :   keys.add("numbered","SEGMENT","this is the lists of atoms in the segment that are being considered");
      73         484 :   keys.add("compulsory","BONDLENGTH","the length to use for bonds");
      74         484 :   keys.add("numbered","STRUCTURE","the reference structure");
      75         484 :   keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. "
      76             :            "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD. For more details on the "
      77             :            "DRMSD method see \\ref DRMSD.");
      78         484 :   keys.add("optional","STRANDS_CUTOFF","If in a segment of protein the two strands are further apart then the calculation "
      79             :            "of the actual RMSD is skipped as the structure is very far from being beta-sheet like. "
      80             :            "This keyword speeds up the calculation enormously when you are using the LESS_THAN option. "
      81             :            "However, if you are using some other option, then this cannot be used");
      82         484 :   keys.add("optional","CUTOFF_ATOMS","the pair of atoms that are used to calculate the strand cutoff");
      83         484 :   keys.addFlag("VERBOSE",false,"write a more detailed output");
      84         484 :   keys.add("optional","LESS_THAN","calculate the number of a residue segments that are within a certain target distance of this secondary structure type. "
      85             :            "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ is a \\ref switchingfunction.");
      86         484 :   keys.add("optional","R_0","The r_0 parameter of the switching function.");
      87         484 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
      88         484 :   keys.add("compulsory","NN","8","The n parameter of the switching function");
      89         484 :   keys.add("compulsory","MM","12","The m parameter of the switching function");
      90         484 :   keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      91         484 :   keys.addOutputComponent("struct","default","the vectors containing the rmsd distances between the residues and each of the reference structures");
      92         484 :   keys.addOutputComponent("lessthan","default","the number blocks of residues that have an RMSD from the secondary structure that is less than the threshold");
      93         726 :   keys.needsAction("SECONDARY_STRUCTURE_RMSD"); keys.needsAction("LESS_THAN"); keys.needsAction("SUM");
      94         242 : }
      95             : 
      96          32 : void SecondaryStructureRMSD::readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& moltype, std::vector<unsigned>& chain_lengths, std::string& all_atoms ) {
      97          32 :   auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
      98          32 :   if( ! moldat ) action->error("Unable to find MOLINFO in input");
      99             : 
     100          64 :   std::vector<std::string> resstrings; action->parseVector( "RESIDUES", resstrings );
     101          32 :   if(resstrings.size()==0) action->error("residues are not defined, check the keyword RESIDUES");
     102          64 :   else if( Tools::caseInSensStringCompare(resstrings[0], "all") ) {
     103             :     resstrings[0]="all";
     104          30 :     action->log.printf("  examining all possible secondary structure combinations\n");
     105             :   } else {
     106           2 :     action->log.printf("  examining secondary structure in residue positions : %s ",resstrings[0].c_str() );
     107           3 :     for(unsigned i=1; i<resstrings.size(); ++i) action->log.printf(", %s",resstrings[i].c_str() );
     108           2 :     action->log.printf("\n");
     109             :   }
     110             :   std::vector< std::vector<AtomNumber> > backatoms;
     111          32 :   moldat->getBackbone( resstrings, moltype, backatoms );
     112             : 
     113          32 :   chain_lengths.resize( backatoms.size() );
     114         473 :   for(unsigned i=0; i<backatoms.size(); ++i) {
     115         441 :     chain_lengths[i]=backatoms[i].size();
     116       18141 :     for(unsigned j=0; j<backatoms[i].size(); ++j) {
     117       17700 :       std::string bat_str; Tools::convert( backatoms[i][j].serial(), bat_str );
     118       17732 :       if( i==0 && j==0 ) all_atoms = "ATOMS=" + bat_str;
     119       35336 :       else all_atoms += "," + bat_str;
     120             :     }
     121             :   }
     122          32 : }
     123             : 
     124             : 
     125          32 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
     126             :   Action(ao),
     127             :   ActionWithVector(ao),
     128          32 :   nopbc(false),
     129          32 :   align_strands(false),
     130          32 :   s_cutoff2(0),
     131          32 :   align_atom_1(0),
     132          32 :   align_atom_2(0)
     133             : {
     134          32 :   if( plumed.usingNaturalUnits() ) error("cannot use this collective variable when using natural units");
     135             : 
     136          64 :   parse("TYPE",alignType); parseFlag("NOPBC",nopbc);
     137          32 :   log.printf("  distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
     138          64 :   log<<"  Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)"); log<<"\n";
     139             : 
     140          32 :   parseFlag("VERBOSE",verbose_output);
     141             : 
     142          64 :   if( keywords.exists("STRANDS_CUTOFF") ) {
     143          32 :     double s_cutoff = 0;
     144          32 :     parse("STRANDS_CUTOFF",s_cutoff); align_strands=true;
     145          32 :     if( s_cutoff>0) {
     146          26 :       log.printf("  ignoring contributions from strands that are more than %f apart\n",s_cutoff);
     147          52 :       std::vector<unsigned> cutatoms; parseVector("CUTOFF_ATOMS",cutatoms);
     148          26 :       if( cutatoms.size()==2 ) {
     149          26 :         align_atom_1=cutatoms[0]; align_atom_2=cutatoms[1];
     150           0 :       } else error("did not find CUTOFF_ATOMS in input");
     151             :     }
     152          32 :     s_cutoff2=s_cutoff*s_cutoff;
     153             :   }
     154             : 
     155             :   // Read in the atoms
     156          64 :   std::vector<AtomNumber> all_atoms; parseAtomList("ATOMS",all_atoms); requestAtoms( all_atoms );
     157             : 
     158      132469 :   for(unsigned i=1;; ++i) {
     159             :     std::vector<unsigned> newatoms;
     160      265002 :     if( !parseNumberedVector("SEGMENT",i,newatoms) ) break;
     161      132469 :     if( verbose_output ) {
     162           0 :       log.printf("  Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
     163           0 :       for(unsigned i=0; i<newatoms.size(); ++i) log.printf("%d ",all_atoms[newatoms[i]].serial() );
     164           0 :       log.printf("\n");
     165             :     }
     166      132469 :     colvar_atoms.push_back( newatoms );
     167      132469 :   }
     168          32 :   if( colvar_atoms.size()==0 ) error("missing SEGMENT keywords -- no atoms have been specified for comparison");
     169             : 
     170          64 :   double bondlength; parse("BONDLENGTH",bondlength); bondlength=bondlength/getUnits().getLength();
     171             : 
     172             :   // Read in the reference structure
     173          47 :   for(unsigned ii=1;; ++ii) {
     174             :     std::vector<double> cstruct;
     175         158 :     if( !parseNumberedVector("STRUCTURE",ii,cstruct) ) break ;
     176          47 :     plumed_assert( cstruct.size()%3==0 && cstruct.size()/3==colvar_atoms[0].size() );
     177          47 :     std::vector<Vector> structure( cstruct.size()/3 );
     178        1457 :     for(unsigned i=0; i<structure.size(); ++i) {
     179        5640 :       for(unsigned j=0; j<3; ++j) structure[i][j] = 0.1*cstruct[3*i+j]/getUnits().getLength();
     180             :     }
     181          47 :     if( alignType=="DRMSD" ) {
     182             :       std::map<std::pair<unsigned,unsigned>, double> targets;
     183         660 :       for(unsigned i=0; i<structure.size()-1; ++i) {
     184       10208 :         for(unsigned j=i+1; j<structure.size(); ++j) {
     185        9570 :           double distance = delta( structure[i], structure[j] ).modulo();
     186        9570 :           if(distance > bondlength) targets[std::make_pair(i,j)] = distance;
     187             :         }
     188             :       }
     189          22 :       drmsd_targets.push_back( targets );
     190             :     } else {
     191          25 :       Vector center; std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
     192         775 :       for(unsigned i=0; i<structure.size(); ++i) center+=structure[i]*align[i];
     193         775 :       for(unsigned i=0; i<structure.size(); ++i) structure[i] -= center;
     194          25 :       RMSD newrmsd; newrmsd.clear();
     195          25 :       newrmsd.set(align,displace,structure,alignType,true,true);
     196          25 :       myrmsd.push_back( newrmsd );
     197          25 :     }
     198          47 :   }
     199             : 
     200             :   // And create values to hold everything
     201          32 :   unsigned nref = myrmsd.size(); if( alignType=="DRMSD" ) nref=drmsd_targets.size();
     202          32 :   plumed_assert( nref>0 );
     203          32 :   std::vector<unsigned> shape(1); shape[0]=colvar_atoms.size();
     204          32 :   if( nref==1 ) { addValue( shape ); setNotPeriodic(); }
     205             :   else {
     206             :     std::string num;
     207          45 :     for(unsigned i=0; i<nref; ++i) {
     208          30 :       Tools::convert( i+1, num ); addComponent( "struct-" + num, shape );
     209          60 :       componentIsNotPeriodic( "struct-" + num );
     210             :     }
     211             :   }
     212          32 : }
     213             : 
     214          79 : void SecondaryStructureRMSD::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
     215          79 :   if( s_cutoff2>0 ) task_reducing_actions.push_back(this);
     216          79 : }
     217             : 
     218      596136 : int SecondaryStructureRMSD::checkTaskStatus( const unsigned& taskno, int& flag ) const {
     219      596136 :   if( s_cutoff2>0 ) {
     220      596136 :     Vector distance=pbcDistance( ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_1) ),
     221             :                                  ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_2) ) );
     222      596136 :     if( distance.modulo2()<s_cutoff2 ) return 1;
     223      541322 :     return 0;
     224           0 :   } return flag;
     225             : }
     226             : 
     227         240 : void SecondaryStructureRMSD::calculate() {
     228         240 :   runAllTasks();
     229         240 : }
     230             : 
     231       55246 : void SecondaryStructureRMSD::performTask( const unsigned& current, MultiValue& myvals ) const {
     232             :   // Resize the derivatives if need be
     233       55246 :   unsigned nderi = 3*getNumberOfAtoms()+9;
     234       55246 :   if( myvals.getNumberOfDerivatives()!=nderi ) myvals.resize( myvals.getNumberOfValues(), nderi, 0, 0 );
     235             :   // Retrieve the positions
     236       55246 :   const unsigned natoms = colvar_atoms[current].size();
     237       55246 :   std::vector<Vector> pos( natoms ), deriv( natoms );
     238     1712626 :   for(unsigned i=0; i<natoms; ++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
     239             : 
     240             :   // This aligns the two strands if this is required
     241       55246 :   Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
     242       55246 :   if( align_strands ) {
     243       55246 :     Vector origin_old, origin_new; origin_old=pos[align_atom_2];
     244       55246 :     origin_new=pos[align_atom_1]+distance;
     245      883936 :     for(unsigned i=15; i<30; ++i) {
     246      828690 :       pos[i]+=( origin_new - origin_old );
     247             :     }
     248           0 :   } else if( !nopbc ) {
     249           0 :     for(unsigned i=0; i<natoms-1; ++i) {
     250           0 :       const Vector & first (pos[i]);
     251           0 :       Vector & second (pos[i+1]);
     252           0 :       second=first+pbcDistance(first,second);
     253             :     }
     254             :   }
     255             :   // Create a holder for the derivatives
     256       55246 :   if( alignType=="DRMSD" ) {
     257             :     // And now calculate the DRMSD
     258       18826 :     const unsigned rs = drmsd_targets.size();
     259       47061 :     for(unsigned i=0; i<rs; ++i) {
     260       28235 :       double drmsd=0; Vector distance; Tensor vir; vir.zero();
     261      875285 :       for(unsigned j=0; j<natoms; ++j) deriv[j].zero();
     262    11519700 :       for(const auto & it : drmsd_targets[i] ) {
     263    11491465 :         const unsigned k=it.first.first;
     264    11491465 :         const unsigned j=it.first.second;
     265             : 
     266    11491465 :         distance=delta( pos[k], pos[j] );
     267    11491465 :         const double len = distance.modulo();
     268    11491465 :         const double diff = len - it.second;
     269    11491465 :         const double der = diff / len;
     270    11491465 :         drmsd += diff*diff;
     271             : 
     272    11491465 :         if( !doNotCalculateDerivatives() ) {
     273    11262001 :           deriv[k] += -der*distance; deriv[j] += der*distance;
     274    11262001 :           vir += -der*Tensor(distance,distance);
     275             :         }
     276             :       }
     277             : 
     278       28235 :       const double inpairs = 1./static_cast<double>(drmsd_targets[i].size());
     279       28235 :       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
     280       28235 :       drmsd = sqrt(inpairs*drmsd); myvals.setValue( ostrn, drmsd );
     281             : 
     282       28235 :       if( !doNotCalculateDerivatives() ) {
     283       27671 :         double scalef = inpairs / drmsd;
     284      857801 :         for(unsigned j=0; j<natoms; ++j) {
     285             :           const unsigned ja = getAtomIndex( current, j );
     286      830130 :           myvals.addDerivative( ostrn, 3*ja + 0, scalef*deriv[j][0] ); myvals.updateIndex( ostrn, 3*ja+0 );
     287      830130 :           myvals.addDerivative( ostrn, 3*ja + 1, scalef*deriv[j][1] ); myvals.updateIndex( ostrn, 3*ja+1 );
     288      830130 :           myvals.addDerivative( ostrn, 3*ja + 2, scalef*deriv[j][2] ); myvals.updateIndex( ostrn, 3*ja+2 );
     289             :         }
     290       27671 :         unsigned nbase = myvals.getNumberOfDerivatives() - 9;
     291      110684 :         for(unsigned k=0; k<3; ++k) {
     292      332052 :           for(unsigned j=0; j<3; ++j) {
     293      249039 :             myvals.addDerivative( ostrn, nbase + 3*k + j, scalef*vir(k,j) );
     294      249039 :             myvals.updateIndex( ostrn, nbase + 3*k + j );
     295             :           }
     296             :         }
     297             :       }
     298             :     }
     299             :   } else {
     300       36420 :     const unsigned rs = myrmsd.size();
     301       91020 :     for(unsigned i=0; i<rs; ++i) {
     302       54600 :       double nr = myrmsd[i].calculate( pos, deriv, false );
     303       54600 :       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
     304       54600 :       myvals.setValue( ostrn, nr );
     305             : 
     306       54600 :       if( !doNotCalculateDerivatives() ) {
     307       54600 :         Tensor vir; vir.zero();
     308     1692600 :         for(unsigned j=0; j<natoms; ++j) {
     309             :           const unsigned ja = getAtomIndex( current, j );
     310     1638000 :           myvals.addDerivative( ostrn, 3*ja + 0, deriv[j][0] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+0 );
     311     1638000 :           myvals.addDerivative( ostrn, 3*ja + 1, deriv[j][1] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+1 );
     312     1638000 :           myvals.addDerivative( ostrn, 3*ja + 2, deriv[j][2] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+2 );
     313     1638000 :           vir+=(-1.0*Tensor( pos[j], deriv[j] ));
     314             :         }
     315       54600 :         unsigned nbase = myvals.getNumberOfDerivatives() - 9;
     316      218400 :         for(unsigned k=0; k<3; ++k) {
     317      655200 :           for(unsigned j=0; j<3; ++j) {
     318      491400 :             myvals.addDerivative( ostrn, nbase + 3*k + j, vir(k,j) );
     319      491400 :             myvals.updateIndex( ostrn, nbase + 3*k + j );
     320             :           }
     321             :         }
     322             :       }
     323             :     }
     324             :   }
     325       55246 :   return;
     326             : }
     327             : 
     328             : }
     329             : }