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             : {
      73             : private:
      74             :   bool             pbc;
      75             :   bool             doratio;
      76             :   double           constant;
      77             :   double           inept;
      78             :   std::vector<double>   rtwo;
      79             :   std::vector<unsigned> nga;
      80             :   std::unique_ptr<NeighborList> nl;
      81             :   unsigned         tot_size;
      82             : public:
      83             :   static void registerKeywords( Keywords& keys );
      84             :   explicit PRE(const ActionOptions&);
      85             :   void calculate() override;
      86             :   void update() override;
      87             : };
      88             : 
      89             : PLUMED_REGISTER_ACTION(PRE,"PRE")
      90             : 
      91           6 : void PRE::registerKeywords( Keywords& keys ) {
      92           6 :   MetainferenceBase::registerKeywords(keys);
      93          12 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
      94          12 :   keys.addFlag("NORATIO",false,"Set to TRUE if you want to compute PRE without Intensity Ratio");
      95          12 :   keys.add("compulsory","INEPT","is the INEPT time (in ms).");
      96          12 :   keys.add("compulsory","TAUC","is the correlation time (in ns) for this electron-nuclear interaction.");
      97          12 :   keys.add("compulsory","OMEGA","is the Larmor frequency of the nuclear spin (in MHz).");
      98          12 :   keys.add("atoms","SPINLABEL","The atom to be used as the paramagnetic center.");
      99          12 :   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          12 :   keys.reset_style("GROUPA","atoms");
     103          12 :   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          12 :   keys.add("numbered","PREINT","Add an experimental value for each PRE.");
     106          12 :   keys.addOutputComponent("pre","default","the # PRE");
     107          12 :   keys.addOutputComponent("exp","PREINT","the # PRE experimental intensity");
     108           6 : }
     109             : 
     110           4 : PRE::PRE(const ActionOptions&ao):
     111             :   PLUMED_METAINF_INIT(ao),
     112           4 :   pbc(true),
     113           4 :   doratio(true)
     114             : {
     115           4 :   bool nopbc=!pbc;
     116           4 :   parseFlag("NOPBC",nopbc);
     117           4 :   pbc=!nopbc;
     118             : 
     119           4 :   bool noratio=!doratio;
     120           4 :   parseFlag("NORATIO",noratio);
     121           4 :   doratio=!noratio;
     122             : 
     123             :   std::vector<AtomNumber> atom;
     124           8 :   parseAtomList("SPINLABEL",atom);
     125           4 :   if(atom.size()!=1) error("Number of specified atom should be 1");
     126             : 
     127             :   // Read in the atoms
     128             :   std::vector<AtomNumber> t, ga_lista, gb_lista;
     129          12 :   for(int i=1;; ++i ) {
     130          32 :     parseAtomList("GROUPA", i, t );
     131          16 :     if( t.empty() ) break;
     132          28 :     for(unsigned j=0; j<t.size(); j++) {ga_lista.push_back(t[j]); gb_lista.push_back(atom[0]);}
     133          12 :     nga.push_back(t.size());
     134          12 :     t.resize(0);
     135          12 :   }
     136             : 
     137             :   // Read in reference values
     138           4 :   rtwo.resize( nga.size() );
     139           4 :   if(doratio) {
     140             :     unsigned ntarget=0;
     141           4 :     for(unsigned i=0; i<nga.size(); ++i) {
     142           8 :       if( !parseNumbered( "RTWO", i+1, rtwo[i] ) ) break;
     143           0 :       ntarget++;
     144             :     }
     145           4 :     if( ntarget==0 ) {
     146           4 :       parse("RTWO",rtwo[0]);
     147          12 :       for(unsigned i=1; i<nga.size(); ++i) rtwo[i]=rtwo[0];
     148           0 :     } else if( ntarget!=nga.size() ) error("found wrong number of RTWO values");
     149             :   }
     150             : 
     151           4 :   double tauc=0.;
     152           4 :   parse("TAUC",tauc);
     153           4 :   if(tauc==0.) error("TAUC must be set");
     154             : 
     155           4 :   double omega=0.;
     156           4 :   parse("OMEGA",omega);
     157           4 :   if(omega==0.) error("OMEGA must be set");
     158             : 
     159           4 :   inept=0.;
     160           4 :   if(doratio) {
     161           4 :     parse("INEPT",inept);
     162           4 :     if(inept==0.) error("INEPT must be set");
     163           4 :     inept *= 0.001; // ms2s
     164             :   }
     165             : 
     166             :   const double ns2s   = 0.000000001;
     167             :   const double MHz2Hz = 1000000.;
     168             :   const double Kappa  = 12300000000.00; // this is 1/15*S*(S+1)*gamma^2*g^2*beta^2
     169             :   // where gamma is the nuclear gyromagnetic ratio,
     170             :   // g is the electronic g factor, and beta is the Bohr magneton
     171             :   // in nm^6/s^2
     172           4 :   constant = (4.*tauc*ns2s+(3.*tauc*ns2s)/(1+omega*omega*MHz2Hz*MHz2Hz*tauc*tauc*ns2s*ns2s))*Kappa;
     173             : 
     174             :   // Optionally add an experimental value (like with RDCs)
     175             :   std::vector<double> exppre;
     176           4 :   exppre.resize( nga.size() );
     177             :   unsigned ntarget=0;
     178          10 :   for(unsigned i=0; i<nga.size(); ++i) {
     179          16 :     if( !parseNumbered( "PREINT", i+1, exppre[i] ) ) break;
     180           6 :     ntarget++;
     181             :   }
     182             :   bool addexp=false;
     183           4 :   if(ntarget!=nga.size() && ntarget!=0) error("found wrong number of PREINT values");
     184           4 :   if(ntarget==nga.size()) addexp=true;
     185           4 :   if(getDoScore()&&!addexp) error("with DOSCORE you need to set the PREINT values");
     186             : 
     187             :   // Create neighbour lists
     188           8 :   nl=Tools::make_unique<NeighborList>(gb_lista,ga_lista,false,true,pbc,getPbc(),comm);
     189             : 
     190             :   // Output details of all contacts
     191             :   unsigned index=0;
     192          16 :   for(unsigned i=0; i<nga.size(); ++i) {
     193          12 :     log.printf("  The %uth PRE is calculated using %u equivalent atoms:\n", i, nga[i]);
     194          12 :     log.printf("    %d", ga_lista[index].serial());
     195          12 :     index++;
     196          16 :     for(unsigned j=1; j<nga[i]; j++) {
     197           4 :       log.printf(" %d", ga_lista[index].serial());
     198           4 :       index++;
     199             :     }
     200          12 :     log.printf("\n");
     201             :   }
     202           4 :   tot_size = index;
     203             : 
     204           4 :   if(pbc)      log.printf("  using periodic boundary conditions\n");
     205           0 :   else         log.printf("  without periodic boundary conditions\n");
     206             : 
     207           8 :   log << " Bibliography" << plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
     208             : 
     209           4 :   if(!getDoScore()) {
     210           8 :     for(unsigned i=0; i<nga.size(); i++) {
     211           6 :       std::string num; Tools::convert(i,num);
     212          12 :       addComponentWithDerivatives("pre-"+num);
     213          12 :       componentIsNotPeriodic("pre-"+num);
     214             :     }
     215           2 :     if(addexp) {
     216           0 :       for(unsigned i=0; i<nga.size(); i++) {
     217           0 :         std::string num; Tools::convert(i,num);
     218           0 :         addComponent("exp-"+num);
     219           0 :         componentIsNotPeriodic("exp-"+num);
     220           0 :         Value* comp=getPntrToComponent("exp-"+num);
     221           0 :         comp->set(exppre[i]);
     222             :       }
     223             :     }
     224             :   } else {
     225           8 :     for(unsigned i=0; i<nga.size(); i++) {
     226           6 :       std::string num; Tools::convert(i,num);
     227          12 :       addComponent("pre-"+num);
     228          12 :       componentIsNotPeriodic("pre-"+num);
     229             :     }
     230           8 :     for(unsigned i=0; i<nga.size(); i++) {
     231           6 :       std::string num; Tools::convert(i,num);
     232          12 :       addComponent("exp-"+num);
     233           6 :       componentIsNotPeriodic("exp-"+num);
     234           6 :       Value* comp=getPntrToComponent("exp-"+num);
     235           6 :       comp->set(exppre[i]);
     236             :     }
     237             :   }
     238             : 
     239           4 :   requestAtoms(nl->getFullAtomList(), false);
     240           4 :   if(getDoScore()) {
     241           2 :     setParameters(exppre);
     242           2 :     Initialise(nga.size());
     243             :   }
     244           4 :   setDerivatives();
     245           4 :   checkRead();
     246           4 : }
     247             : 
     248         350 : void PRE::calculate()
     249             : {
     250         350 :   std::vector<Vector> deriv(tot_size, Vector{0,0,0});
     251         350 :   std::vector<double> fact(nga.size(), 0.);
     252             : 
     253             :   // cycle over the number of PRE
     254         350 :   #pragma omp parallel for num_threads(OpenMP::getNumThreads())
     255             :   for(unsigned i=0; i<nga.size(); i++) {
     256             :     Tensor dervir;
     257             :     double pre=0;
     258             :     unsigned index=0;
     259             :     for(unsigned k=0; k<i; k++) index+=nga[k];
     260             :     const double c_aver=constant/static_cast<double>(nga[i]);
     261             :     std::string num; Tools::convert(i,num);
     262             :     Value* val=getPntrToComponent("pre-"+num);
     263             :     // cycle over equivalent atoms
     264             :     for(unsigned j=0; j<nga[i]; j++) {
     265             :       // the first atom is always the same (the paramagnetic group)
     266             :       const unsigned i0=nl->getClosePair(index+j).first;
     267             :       const unsigned i1=nl->getClosePair(index+j).second;
     268             : 
     269             :       Vector distance;
     270             :       if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
     271             :       else    distance=delta(getPosition(i0),getPosition(i1));
     272             : 
     273             :       const double r2=distance.modulo2();
     274             :       const double r6=r2*r2*r2;
     275             :       const double r8=r6*r2;
     276             :       const double tmpir6=c_aver/r6;
     277             :       const double tmpir8=-6.*c_aver/r8;
     278             : 
     279             :       pre += tmpir6;
     280             :       deriv[index+j] = -tmpir8*distance;
     281             :       if(!getDoScore()) dervir   +=  Tensor(distance,deriv[index+j]);
     282             :     }
     283             :     double tmpratio;
     284             :     if(!doratio) {
     285             :       tmpratio = pre ; //prova a caso per vedere se lui da problemi
     286             :       fact[i] = 1.; //prova a caso per vedere se lui da problemi
     287             :     } else {
     288             :       tmpratio = rtwo[i]*std::exp(-pre*inept) / (rtwo[i]+pre);
     289             :       fact[i] = -tmpratio*(inept+1./(rtwo[i]+pre));
     290             :     }
     291             :     const double ratio = tmpratio;
     292             :     val->set(ratio) ;
     293             :     if(!getDoScore()) {
     294             :       setBoxDerivatives(val, fact[i]*dervir);
     295             :       for(unsigned j=0; j<nga[i]; j++) {
     296             :         const unsigned i0=nl->getClosePair(index+j).first;
     297             :         const unsigned i1=nl->getClosePair(index+j).second;
     298             :         setAtomsDerivatives(val, i0,  fact[i]*deriv[index+j]);
     299             :         setAtomsDerivatives(val, i1, -fact[i]*deriv[index+j]);
     300             :       }
     301             :     } else setCalcData(i, ratio);
     302             :   }
     303             : 
     304         350 :   if(getDoScore()) {
     305             :     /* Metainference */
     306         175 :     Tensor dervir;
     307         175 :     double score = getScore();
     308         175 :     setScore(score);
     309             : 
     310             :     /* calculate final derivatives */
     311         175 :     Value* val=getPntrToComponent("score");
     312         700 :     for(unsigned i=0; i<nga.size(); i++) {
     313             :       unsigned index=0;
     314        1050 :       for(unsigned k=0; k<i; k++) index+=nga[k];
     315             :       // cycle over equivalent atoms
     316        1225 :       for(unsigned j=0; j<nga[i]; j++) {
     317         700 :         const unsigned i0=nl->getClosePair(index+j).first;
     318         700 :         const unsigned i1=nl->getClosePair(index+j).second;
     319             : 
     320         700 :         Vector distance;
     321         700 :         if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
     322           0 :         else    distance=delta(getPosition(i0),getPosition(i1));
     323             : 
     324         700 :         dervir += Tensor(distance,fact[i]*deriv[index+j]*getMetaDer(i));
     325         700 :         setAtomsDerivatives(val, i0,  fact[i]*deriv[index+j]*getMetaDer(i));
     326         700 :         setAtomsDerivatives(val, i1, -fact[i]*deriv[index+j]*getMetaDer(i));
     327             :       }
     328             :     }
     329         175 :     setBoxDerivatives(val, dervir);
     330             :   }
     331         350 : }
     332             : 
     333          20 : void PRE::update() {
     334             :   // write status file
     335          20 :   if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
     336          20 : }
     337             : 
     338             : }
     339             : }