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.setValueDescription("vector","a vector containing the rmsd distance between each of the residue segments and the reference structure");
      92         484 :   keys.addOutputComponent("struct","default","vector","the vectors containing the rmsd distances between the residues and each of the reference structures");
      93         484 :   keys.addOutputComponent("lessthan","default","scalar","the number blocks of residues that have an RMSD from the secondary structure that is less than the threshold");
      94         726 :   keys.needsAction("SECONDARY_STRUCTURE_RMSD"); keys.needsAction("LESS_THAN"); keys.needsAction("SUM");
      95         242 : }
      96             : 
      97          32 : void SecondaryStructureRMSD::readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& moltype, std::vector<unsigned>& chain_lengths, std::string& all_atoms ) {
      98          32 :   auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
      99          32 :   if( ! moldat ) action->error("Unable to find MOLINFO in input");
     100             : 
     101          64 :   std::vector<std::string> resstrings; action->parseVector( "RESIDUES", resstrings );
     102          32 :   if(resstrings.size()==0) action->error("residues are not defined, check the keyword RESIDUES");
     103          64 :   else if( Tools::caseInSensStringCompare(resstrings[0], "all") ) {
     104             :     resstrings[0]="all";
     105          30 :     action->log.printf("  examining all possible secondary structure combinations\n");
     106             :   } else {
     107           2 :     action->log.printf("  examining secondary structure in residue positions : %s ",resstrings[0].c_str() );
     108           3 :     for(unsigned i=1; i<resstrings.size(); ++i) action->log.printf(", %s",resstrings[i].c_str() );
     109           2 :     action->log.printf("\n");
     110             :   }
     111             :   std::vector< std::vector<AtomNumber> > backatoms;
     112          32 :   moldat->getBackbone( resstrings, moltype, backatoms );
     113             : 
     114          32 :   chain_lengths.resize( backatoms.size() );
     115         473 :   for(unsigned i=0; i<backatoms.size(); ++i) {
     116         441 :     chain_lengths[i]=backatoms[i].size();
     117       18141 :     for(unsigned j=0; j<backatoms[i].size(); ++j) {
     118       17700 :       std::string bat_str; Tools::convert( backatoms[i][j].serial(), bat_str );
     119       17732 :       if( i==0 && j==0 ) all_atoms = "ATOMS=" + bat_str;
     120       35336 :       else all_atoms += "," + bat_str;
     121             :     }
     122             :   }
     123          32 : }
     124             : 
     125             : 
     126          32 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
     127             :   Action(ao),
     128             :   ActionWithVector(ao),
     129          32 :   nopbc(false),
     130          32 :   align_strands(false),
     131          32 :   s_cutoff2(0),
     132          32 :   align_atom_1(0),
     133          32 :   align_atom_2(0)
     134             : {
     135          32 :   if( plumed.usingNaturalUnits() ) error("cannot use this collective variable when using natural units");
     136             : 
     137          64 :   parse("TYPE",alignType); parseFlag("NOPBC",nopbc);
     138          32 :   log.printf("  distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
     139          64 :   log<<"  Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)"); log<<"\n";
     140             : 
     141          32 :   parseFlag("VERBOSE",verbose_output);
     142             : 
     143          64 :   if( keywords.exists("STRANDS_CUTOFF") ) {
     144          32 :     double s_cutoff = 0;
     145          32 :     parse("STRANDS_CUTOFF",s_cutoff); align_strands=true;
     146          32 :     if( s_cutoff>0) {
     147          26 :       log.printf("  ignoring contributions from strands that are more than %f apart\n",s_cutoff);
     148          52 :       std::vector<unsigned> cutatoms; parseVector("CUTOFF_ATOMS",cutatoms);
     149          26 :       if( cutatoms.size()==2 ) {
     150          26 :         align_atom_1=cutatoms[0]; align_atom_2=cutatoms[1];
     151           0 :       } else error("did not find CUTOFF_ATOMS in input");
     152             :     }
     153          32 :     s_cutoff2=s_cutoff*s_cutoff;
     154             :   }
     155             : 
     156             :   // Read in the atoms
     157          64 :   std::vector<AtomNumber> all_atoms; parseAtomList("ATOMS",all_atoms); requestAtoms( all_atoms );
     158             : 
     159      132469 :   for(unsigned i=1;; ++i) {
     160             :     std::vector<unsigned> newatoms;
     161      265002 :     if( !parseNumberedVector("SEGMENT",i,newatoms) ) break;
     162      132469 :     if( verbose_output ) {
     163           0 :       log.printf("  Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
     164           0 :       for(unsigned i=0; i<newatoms.size(); ++i) log.printf("%d ",all_atoms[newatoms[i]].serial() );
     165           0 :       log.printf("\n");
     166             :     }
     167      132469 :     colvar_atoms.push_back( newatoms );
     168      132469 :   }
     169          32 :   if( colvar_atoms.size()==0 ) error("missing SEGMENT keywords -- no atoms have been specified for comparison");
     170             : 
     171          64 :   double bondlength; parse("BONDLENGTH",bondlength); bondlength=bondlength/getUnits().getLength();
     172             : 
     173             :   // Read in the reference structure
     174          47 :   for(unsigned ii=1;; ++ii) {
     175             :     std::vector<double> cstruct;
     176         158 :     if( !parseNumberedVector("STRUCTURE",ii,cstruct) ) break ;
     177          47 :     plumed_assert( cstruct.size()%3==0 && cstruct.size()/3==colvar_atoms[0].size() );
     178          47 :     std::vector<Vector> structure( cstruct.size()/3 );
     179        1457 :     for(unsigned i=0; i<structure.size(); ++i) {
     180        5640 :       for(unsigned j=0; j<3; ++j) structure[i][j] = 0.1*cstruct[3*i+j]/getUnits().getLength();
     181             :     }
     182          47 :     if( alignType=="DRMSD" ) {
     183             :       std::map<std::pair<unsigned,unsigned>, double> targets;
     184         660 :       for(unsigned i=0; i<structure.size()-1; ++i) {
     185       10208 :         for(unsigned j=i+1; j<structure.size(); ++j) {
     186        9570 :           double distance = delta( structure[i], structure[j] ).modulo();
     187        9570 :           if(distance > bondlength) targets[std::make_pair(i,j)] = distance;
     188             :         }
     189             :       }
     190          22 :       drmsd_targets.push_back( targets );
     191             :     } else {
     192          25 :       Vector center; std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
     193         775 :       for(unsigned i=0; i<structure.size(); ++i) center+=structure[i]*align[i];
     194         775 :       for(unsigned i=0; i<structure.size(); ++i) structure[i] -= center;
     195          25 :       RMSD newrmsd; newrmsd.clear();
     196          25 :       newrmsd.set(align,displace,structure,alignType,true,true);
     197          25 :       myrmsd.push_back( newrmsd );
     198          25 :     }
     199          47 :   }
     200             : 
     201             :   // And create values to hold everything
     202          32 :   unsigned nref = myrmsd.size(); if( alignType=="DRMSD" ) nref=drmsd_targets.size();
     203          32 :   plumed_assert( nref>0 );
     204          32 :   std::vector<unsigned> shape(1); shape[0]=colvar_atoms.size();
     205          32 :   if( nref==1 ) { addValue( shape ); setNotPeriodic(); }
     206             :   else {
     207             :     std::string num;
     208          45 :     for(unsigned i=0; i<nref; ++i) {
     209          30 :       Tools::convert( i+1, num ); addComponent( "struct-" + num, shape );
     210          60 :       componentIsNotPeriodic( "struct-" + num );
     211             :     }
     212             :   }
     213          32 : }
     214             : 
     215          79 : void SecondaryStructureRMSD::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
     216          79 :   if( s_cutoff2>0 ) task_reducing_actions.push_back(this);
     217          79 : }
     218             : 
     219      596136 : int SecondaryStructureRMSD::checkTaskStatus( const unsigned& taskno, int& flag ) const {
     220      596136 :   if( s_cutoff2>0 ) {
     221      596136 :     Vector distance=pbcDistance( ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_1) ),
     222             :                                  ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_2) ) );
     223      596136 :     if( distance.modulo2()<s_cutoff2 ) return 1;
     224      541322 :     return 0;
     225           0 :   } return flag;
     226             : }
     227             : 
     228         240 : void SecondaryStructureRMSD::calculate() {
     229         240 :   runAllTasks();
     230         240 : }
     231             : 
     232       55246 : void SecondaryStructureRMSD::performTask( const unsigned& current, MultiValue& myvals ) const {
     233             :   // Resize the derivatives if need be
     234       55246 :   unsigned nderi = 3*getNumberOfAtoms()+9;
     235       55246 :   if( myvals.getNumberOfDerivatives()!=nderi ) myvals.resize( myvals.getNumberOfValues(), nderi, 0, 0 );
     236             :   // Retrieve the positions
     237       55246 :   const unsigned natoms = colvar_atoms[current].size();
     238       55246 :   std::vector<Vector> pos( natoms ), deriv( natoms );
     239     1712626 :   for(unsigned i=0; i<natoms; ++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
     240             : 
     241             :   // This aligns the two strands if this is required
     242       55246 :   Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
     243       55246 :   if( align_strands ) {
     244       55246 :     Vector origin_old, origin_new; origin_old=pos[align_atom_2];
     245       55246 :     origin_new=pos[align_atom_1]+distance;
     246      883936 :     for(unsigned i=15; i<30; ++i) {
     247      828690 :       pos[i]+=( origin_new - origin_old );
     248             :     }
     249           0 :   } else if( !nopbc ) {
     250           0 :     for(unsigned i=0; i<natoms-1; ++i) {
     251           0 :       const Vector & first (pos[i]);
     252           0 :       Vector & second (pos[i+1]);
     253           0 :       second=first+pbcDistance(first,second);
     254             :     }
     255             :   }
     256             :   // Create a holder for the derivatives
     257       55246 :   if( alignType=="DRMSD" ) {
     258             :     // And now calculate the DRMSD
     259       18826 :     const unsigned rs = drmsd_targets.size();
     260       47061 :     for(unsigned i=0; i<rs; ++i) {
     261       28235 :       double drmsd=0; Vector distance; Tensor vir; vir.zero();
     262      875285 :       for(unsigned j=0; j<natoms; ++j) deriv[j].zero();
     263    11519700 :       for(const auto & it : drmsd_targets[i] ) {
     264    11491465 :         const unsigned k=it.first.first;
     265    11491465 :         const unsigned j=it.first.second;
     266             : 
     267    11491465 :         distance=delta( pos[k], pos[j] );
     268    11491465 :         const double len = distance.modulo();
     269    11491465 :         const double diff = len - it.second;
     270    11491465 :         const double der = diff / len;
     271    11491465 :         drmsd += diff*diff;
     272             : 
     273    11491465 :         if( !doNotCalculateDerivatives() ) {
     274    11262001 :           deriv[k] += -der*distance; deriv[j] += der*distance;
     275    11262001 :           vir += -der*Tensor(distance,distance);
     276             :         }
     277             :       }
     278             : 
     279       28235 :       const double inpairs = 1./static_cast<double>(drmsd_targets[i].size());
     280       28235 :       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
     281       28235 :       drmsd = sqrt(inpairs*drmsd); myvals.setValue( ostrn, drmsd );
     282             : 
     283       28235 :       if( !doNotCalculateDerivatives() ) {
     284       27671 :         double scalef = inpairs / drmsd;
     285      857801 :         for(unsigned j=0; j<natoms; ++j) {
     286             :           const unsigned ja = getAtomIndex( current, j );
     287      830130 :           myvals.addDerivative( ostrn, 3*ja + 0, scalef*deriv[j][0] ); myvals.updateIndex( ostrn, 3*ja+0 );
     288      830130 :           myvals.addDerivative( ostrn, 3*ja + 1, scalef*deriv[j][1] ); myvals.updateIndex( ostrn, 3*ja+1 );
     289      830130 :           myvals.addDerivative( ostrn, 3*ja + 2, scalef*deriv[j][2] ); myvals.updateIndex( ostrn, 3*ja+2 );
     290             :         }
     291       27671 :         unsigned nbase = myvals.getNumberOfDerivatives() - 9;
     292      110684 :         for(unsigned k=0; k<3; ++k) {
     293      332052 :           for(unsigned j=0; j<3; ++j) {
     294      249039 :             myvals.addDerivative( ostrn, nbase + 3*k + j, scalef*vir(k,j) );
     295      249039 :             myvals.updateIndex( ostrn, nbase + 3*k + j );
     296             :           }
     297             :         }
     298             :       }
     299             :     }
     300             :   } else {
     301       36420 :     const unsigned rs = myrmsd.size();
     302       91020 :     for(unsigned i=0; i<rs; ++i) {
     303       54600 :       double nr = myrmsd[i].calculate( pos, deriv, false );
     304       54600 :       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
     305       54600 :       myvals.setValue( ostrn, nr );
     306             : 
     307       54600 :       if( !doNotCalculateDerivatives() ) {
     308       54600 :         Tensor vir; vir.zero();
     309     1692600 :         for(unsigned j=0; j<natoms; ++j) {
     310             :           const unsigned ja = getAtomIndex( current, j );
     311     1638000 :           myvals.addDerivative( ostrn, 3*ja + 0, deriv[j][0] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+0 );
     312     1638000 :           myvals.addDerivative( ostrn, 3*ja + 1, deriv[j][1] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+1 );
     313     1638000 :           myvals.addDerivative( ostrn, 3*ja + 2, deriv[j][2] ); myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+2 );
     314     1638000 :           vir+=(-1.0*Tensor( pos[j], deriv[j] ));
     315             :         }
     316       54600 :         unsigned nbase = myvals.getNumberOfDerivatives() - 9;
     317      218400 :         for(unsigned k=0; k<3; ++k) {
     318      655200 :           for(unsigned j=0; j<3; ++j) {
     319      491400 :             myvals.addDerivative( ostrn, nbase + 3*k + j, vir(k,j) );
     320      491400 :             myvals.updateIndex( ostrn, nbase + 3*k + j );
     321             :           }
     322             :         }
     323             :       }
     324             :     }
     325             :   }
     326       55246 :   return;
     327             : }
     328             : 
     329             : }
     330             : }