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