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