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