Line data    Source code

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