Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "MetainferenceBase.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/NeighborList.h"
      25             : #include "tools/Pbc.h"
      26             : #include <memory>
      27             : 
      28             : namespace PLMD {
      29             : namespace isdb {
      30             : 
      31             : //+PLUMEDOC ISDB_COLVAR PRE
      32             : /*
      33             : Calculates the Paramagnetic Resonance Enhancement intensity ratio between a spin label atom and a list of atoms .
      34             : 
      35             : The reference atom for the spin label is added with SPINLABEL, the affected atom(s)
      36             : are give as numbered GROUPA1, GROUPA2, ...
      37             : The additional parameters needed for the calculation are given as INEPT, the inept
      38             : time, TAUC the correlation time, OMEGA, the Larmor frequency and RTWO for the relaxation
      39             : time.
      40             : 
      41             : \ref METAINFERENCE can be activated using DOSCORE and the other relevant keywords.
      42             : 
      43             : \par Examples
      44             : 
      45             : In the following example five PRE intensities are calculated using the distance between the
      46             : oxygen of the spin label and the backbone hydrogen atoms. Omega is the NMR frequency, RTWO the
      47             : R2 for the hydrogen atoms, INEPT of 8 ms for the experiment and a TAUC of 1.21 ns
      48             : 
      49             : \plumedfile
      50             : PRE ...
      51             : LABEL=HN_pre
      52             : INEPT=8
      53             : TAUC=1.21
      54             : OMEGA=900
      55             : SPINLABEL=1818
      56             : GROUPA1=86  RTWO1=0.0120272827
      57             : GROUPA2=177 RTWO2=0.0263953158
      58             : GROUPA3=285 RTWO3=0.0058899829
      59             : GROUPA4=335 RTWO4=0.0102072646
      60             : GROUPA5=451 RTWO5=0.0086341843
      61             : ... PRE
      62             : 
      63             : PRINT ARG=HN_pre.* FILE=PRE.dat STRIDE=1
      64             : 
      65             : \endplumedfile
      66             : 
      67             : */
      68             : //+ENDPLUMEDOC
      69             : 
      70             : class PRE :
      71             :   public MetainferenceBase {
      72             : private:
      73             :   bool             pbc;
      74             :   bool             doratio;
      75             :   double           constant;
      76             :   double           inept;
      77             :   std::vector<double>   rtwo;
      78             :   std::vector<unsigned> nga;
      79             :   std::unique_ptr<NeighborList> nl;
      80             :   unsigned         tot_size;
      81             : public:
      82             :   static void registerKeywords( Keywords& keys );
      83             :   explicit PRE(const ActionOptions&);
      84             :   void calculate() override;
      85             :   void update() override;
      86             : };
      87             : 
      88       13793 : PLUMED_REGISTER_ACTION(PRE,"PRE")
      89             : 
      90           8 : void PRE::registerKeywords( Keywords& keys ) {
      91           8 :   componentsAreNotOptional(keys);
      92           8 :   MetainferenceBase::registerKeywords(keys);
      93          16 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
      94          16 :   keys.addFlag("NORATIO",false,"Set to TRUE if you want to compute PRE without Intensity Ratio");
      95          16 :   keys.add("compulsory","INEPT","is the INEPT time (in ms).");
      96          16 :   keys.add("compulsory","TAUC","is the correlation time (in ns) for this electron-nuclear interaction.");
      97          16 :   keys.add("compulsory","OMEGA","is the Larmor frequency of the nuclear spin (in MHz).");
      98          16 :   keys.add("atoms","SPINLABEL","The atom to be used as the paramagnetic center.");
      99          16 :   keys.add("numbered","GROUPA","the atoms involved in each of the contacts you wish to calculate. "
     100             :            "Keywords like GROUPA1, GROUPA2, GROUPA3,... should be listed and one contact will be "
     101             :            "calculated for each ATOM keyword you specify.");
     102          16 :   keys.reset_style("GROUPA","atoms");
     103          16 :   keys.add("numbered","RTWO","The relaxation of the atom/atoms in the corresponding GROUPA of atoms. "
     104             :            "Keywords like RTWO1, RTWO2, RTWO3,... should be listed.");
     105          16 :   keys.add("numbered","PREINT","Add an experimental value for each PRE.");
     106          16 :   keys.addOutputComponent("pre","default","the # PRE");
     107          16 :   keys.addOutputComponent("exp","PREINT","the # PRE experimental intensity");
     108           8 : }
     109             : 
     110           4 : PRE::PRE(const ActionOptions&ao):
     111             :   PLUMED_METAINF_INIT(ao),
     112           4 :   pbc(true),
     113           4 :   doratio(true) {
     114           4 :   bool nopbc=!pbc;
     115           4 :   parseFlag("NOPBC",nopbc);
     116           4 :   pbc=!nopbc;
     117             : 
     118           4 :   bool noratio=!doratio;
     119           4 :   parseFlag("NORATIO",noratio);
     120           4 :   doratio=!noratio;
     121             : 
     122             :   std::vector<AtomNumber> atom;
     123           8 :   parseAtomList("SPINLABEL",atom);
     124           4 :   if(atom.size()!=1) {
     125           0 :     error("Number of specified atom should be 1");
     126             :   }
     127             : 
     128             :   // Read in the atoms
     129             :   std::vector<AtomNumber> t, ga_lista, gb_lista;
     130          12 :   for(int i=1;; ++i ) {
     131          32 :     parseAtomList("GROUPA", i, t );
     132          16 :     if( t.empty() ) {
     133             :       break;
     134             :     }
     135          28 :     for(unsigned j=0; j<t.size(); j++) {
     136          16 :       ga_lista.push_back(t[j]);
     137          16 :       gb_lista.push_back(atom[0]);
     138             :     }
     139          12 :     nga.push_back(t.size());
     140          12 :     t.resize(0);
     141          12 :   }
     142             : 
     143             :   // Read in reference values
     144           4 :   rtwo.resize( nga.size() );
     145           4 :   if(doratio) {
     146             :     unsigned ntarget=0;
     147           4 :     for(unsigned i=0; i<nga.size(); ++i) {
     148           8 :       if( !parseNumbered( "RTWO", i+1, rtwo[i] ) ) {
     149             :         break;
     150             :       }
     151           0 :       ntarget++;
     152             :     }
     153           4 :     if( ntarget==0 ) {
     154           4 :       parse("RTWO",rtwo[0]);
     155          12 :       for(unsigned i=1; i<nga.size(); ++i) {
     156           8 :         rtwo[i]=rtwo[0];
     157             :       }
     158           0 :     } else if( ntarget!=nga.size() ) {
     159           0 :       error("found wrong number of RTWO values");
     160             :     }
     161             :   }
     162             : 
     163           4 :   double tauc=0.;
     164           4 :   parse("TAUC",tauc);
     165           4 :   if(tauc==0.) {
     166           0 :     error("TAUC must be set");
     167             :   }
     168             : 
     169           4 :   double omega=0.;
     170           4 :   parse("OMEGA",omega);
     171           4 :   if(omega==0.) {
     172           0 :     error("OMEGA must be set");
     173             :   }
     174             : 
     175           4 :   inept=0.;
     176           4 :   if(doratio) {
     177           4 :     parse("INEPT",inept);
     178           4 :     if(inept==0.) {
     179           0 :       error("INEPT must be set");
     180             :     }
     181           4 :     inept *= 0.001; // ms2s
     182             :   }
     183             : 
     184             :   const double ns2s   = 0.000000001;
     185             :   const double MHz2Hz = 1000000.;
     186             :   const double Kappa  = 12300000000.00; // this is 1/15*S*(S+1)*gamma^2*g^2*beta^2
     187             :   // where gamma is the nuclear gyromagnetic ratio,
     188             :   // g is the electronic g factor, and beta is the Bohr magneton
     189             :   // in nm^6/s^2
     190           4 :   constant = (4.*tauc*ns2s+(3.*tauc*ns2s)/(1+omega*omega*MHz2Hz*MHz2Hz*tauc*tauc*ns2s*ns2s))*Kappa;
     191             : 
     192             :   // Optionally add an experimental value (like with RDCs)
     193             :   std::vector<double> exppre;
     194           4 :   exppre.resize( nga.size() );
     195             :   unsigned ntarget=0;
     196          10 :   for(unsigned i=0; i<nga.size(); ++i) {
     197          16 :     if( !parseNumbered( "PREINT", i+1, exppre[i] ) ) {
     198             :       break;
     199             :     }
     200           6 :     ntarget++;
     201             :   }
     202             :   bool addexp=false;
     203           4 :   if(ntarget!=nga.size() && ntarget!=0) {
     204           0 :     error("found wrong number of PREINT values");
     205             :   }
     206           4 :   if(ntarget==nga.size()) {
     207             :     addexp=true;
     208             :   }
     209           4 :   if(getDoScore()&&!addexp) {
     210           0 :     error("with DOSCORE you need to set the PREINT values");
     211             :   }
     212             : 
     213             :   // Create neighbour lists
     214           8 :   nl=Tools::make_unique<NeighborList>(gb_lista,ga_lista,false,true,pbc,getPbc(),comm);
     215             : 
     216             :   // Output details of all contacts
     217             :   unsigned index=0;
     218          16 :   for(unsigned i=0; i<nga.size(); ++i) {
     219          12 :     log.printf("  The %uth PRE is calculated using %u equivalent atoms:\n", i, nga[i]);
     220          12 :     log.printf("    %d", ga_lista[index].serial());
     221          12 :     index++;
     222          16 :     for(unsigned j=1; j<nga[i]; j++) {
     223           4 :       log.printf(" %d", ga_lista[index].serial());
     224           4 :       index++;
     225             :     }
     226          12 :     log.printf("\n");
     227             :   }
     228           4 :   tot_size = index;
     229             : 
     230           4 :   if(pbc) {
     231           4 :     log.printf("  using periodic boundary conditions\n");
     232             :   } else {
     233           0 :     log.printf("  without periodic boundary conditions\n");
     234             :   }
     235             : 
     236           8 :   log << " Bibliography" << plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
     237             : 
     238           4 :   if(!getDoScore()) {
     239           8 :     for(unsigned i=0; i<nga.size(); i++) {
     240             :       std::string num;
     241           6 :       Tools::convert(i,num);
     242           6 :       addComponentWithDerivatives("pre-"+num);
     243          12 :       componentIsNotPeriodic("pre-"+num);
     244             :     }
     245           2 :     if(addexp) {
     246           0 :       for(unsigned i=0; i<nga.size(); i++) {
     247             :         std::string num;
     248           0 :         Tools::convert(i,num);
     249           0 :         addComponent("exp-"+num);
     250           0 :         componentIsNotPeriodic("exp-"+num);
     251           0 :         Value* comp=getPntrToComponent("exp-"+num);
     252           0 :         comp->set(exppre[i]);
     253             :       }
     254             :     }
     255             :   } else {
     256           8 :     for(unsigned i=0; i<nga.size(); i++) {
     257             :       std::string num;
     258           6 :       Tools::convert(i,num);
     259           6 :       addComponent("pre-"+num);
     260          12 :       componentIsNotPeriodic("pre-"+num);
     261             :     }
     262           8 :     for(unsigned i=0; i<nga.size(); i++) {
     263             :       std::string num;
     264           6 :       Tools::convert(i,num);
     265           6 :       addComponent("exp-"+num);
     266           6 :       componentIsNotPeriodic("exp-"+num);
     267           6 :       Value* comp=getPntrToComponent("exp-"+num);
     268           6 :       comp->set(exppre[i]);
     269             :     }
     270             :   }
     271             : 
     272           4 :   requestAtoms(nl->getFullAtomList(), false);
     273           4 :   if(getDoScore()) {
     274           2 :     setParameters(exppre);
     275           2 :     Initialise(nga.size());
     276             :   }
     277           4 :   setDerivatives();
     278           4 :   checkRead();
     279           4 : }
     280             : 
     281         350 : void PRE::calculate() {
     282         350 :   std::vector<Vector> deriv(tot_size, Vector{0,0,0});
     283         350 :   std::vector<double> fact(nga.size(), 0.);
     284             : 
     285             :   // cycle over the number of PRE
     286         350 :   #pragma omp parallel for num_threads(OpenMP::getNumThreads())
     287             :   for(unsigned i=0; i<nga.size(); i++) {
     288             :     Tensor dervir;
     289             :     double pre=0;
     290             :     unsigned index=0;
     291             :     for(unsigned k=0; k<i; k++) {
     292             :       index+=nga[k];
     293             :     }
     294             :     const double c_aver=constant/static_cast<double>(nga[i]);
     295             :     std::string num;
     296             :     Tools::convert(i,num);
     297             :     Value* val=getPntrToComponent("pre-"+num);
     298             :     // cycle over equivalent atoms
     299             :     for(unsigned j=0; j<nga[i]; j++) {
     300             :       // the first atom is always the same (the paramagnetic group)
     301             :       const unsigned i0=nl->getClosePair(index+j).first;
     302             :       const unsigned i1=nl->getClosePair(index+j).second;
     303             : 
     304             :       Vector distance;
     305             :       if(pbc) {
     306             :         distance=pbcDistance(getPosition(i0),getPosition(i1));
     307             :       } else {
     308             :         distance=delta(getPosition(i0),getPosition(i1));
     309             :       }
     310             : 
     311             :       const double r2=distance.modulo2();
     312             :       const double r6=r2*r2*r2;
     313             :       const double r8=r6*r2;
     314             :       const double tmpir6=c_aver/r6;
     315             :       const double tmpir8=-6.*c_aver/r8;
     316             : 
     317             :       pre += tmpir6;
     318             :       deriv[index+j] = -tmpir8*distance;
     319             :       if(!getDoScore()) {
     320             :         dervir   +=  Tensor(distance,deriv[index+j]);
     321             :       }
     322             :     }
     323             :     double tmpratio;
     324             :     if(!doratio) {
     325             :       tmpratio = pre ; //prova a caso per vedere se lui da problemi
     326             :       fact[i] = 1.; //prova a caso per vedere se lui da problemi
     327             :     } else {
     328             :       tmpratio = rtwo[i]*std::exp(-pre*inept) / (rtwo[i]+pre);
     329             :       fact[i] = -tmpratio*(inept+1./(rtwo[i]+pre));
     330             :     }
     331             :     const double ratio = tmpratio;
     332             :     val->set(ratio) ;
     333             :     if(!getDoScore()) {
     334             :       setBoxDerivatives(val, fact[i]*dervir);
     335             :       for(unsigned j=0; j<nga[i]; j++) {
     336             :         const unsigned i0=nl->getClosePair(index+j).first;
     337             :         const unsigned i1=nl->getClosePair(index+j).second;
     338             :         setAtomsDerivatives(val, i0,  fact[i]*deriv[index+j]);
     339             :         setAtomsDerivatives(val, i1, -fact[i]*deriv[index+j]);
     340             :       }
     341             :     } else {
     342             :       setCalcData(i, ratio);
     343             :     }
     344             :   }
     345             : 
     346         350 :   if(getDoScore()) {
     347             :     /* Metainference */
     348         175 :     Tensor dervir;
     349         175 :     double score = getScore();
     350             :     setScore(score);
     351             : 
     352             :     /* calculate final derivatives */
     353         175 :     Value* val=getPntrToComponent("score");
     354         700 :     for(unsigned i=0; i<nga.size(); i++) {
     355             :       unsigned index=0;
     356        1050 :       for(unsigned k=0; k<i; k++) {
     357         525 :         index+=nga[k];
     358             :       }
     359             :       // cycle over equivalent atoms
     360        1225 :       for(unsigned j=0; j<nga[i]; j++) {
     361         700 :         const unsigned i0=nl->getClosePair(index+j).first;
     362         700 :         const unsigned i1=nl->getClosePair(index+j).second;
     363             : 
     364         700 :         Vector distance;
     365         700 :         if(pbc) {
     366         700 :           distance=pbcDistance(getPosition(i0),getPosition(i1));
     367             :         } else {
     368           0 :           distance=delta(getPosition(i0),getPosition(i1));
     369             :         }
     370             : 
     371         700 :         dervir += Tensor(distance,fact[i]*deriv[index+j]*getMetaDer(i));
     372         700 :         setAtomsDerivatives(val, i0,  fact[i]*deriv[index+j]*getMetaDer(i));
     373         700 :         setAtomsDerivatives(val, i1, -fact[i]*deriv[index+j]*getMetaDer(i));
     374             :       }
     375             :     }
     376         175 :     setBoxDerivatives(val, dervir);
     377             :   }
     378         350 : }
     379             : 
     380          20 : void PRE::update() {
     381             :   // write status file
     382          20 :   if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) {
     383           4 :     writeStatus();
     384             :   }
     385          20 : }
     386             : 
     387             : }
     388             : }