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","scalar","the # PRE");
107 12 : keys.addOutputComponent("exp","PREINT","scalar","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 : }
|
