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