Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 "Colvar.h"
23 : #include "ActionRegister.h"
24 : #include "tools/NeighborList.h"
25 : #include "tools/OpenMP.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace colvar {
34 :
35 : //+PLUMEDOC COLVAR PRE
36 : /*
37 : Calculates the Paramegnetic Resonance Enhancement intensity ratio between two atoms.
38 : The reference atom for the spin label is added with SPINLABEL, the affected atom(s)
39 : are give as numbered GROUPA1, GROUPA2, ...
40 : The additional parameters needed for the calculation are given as INEPT, the inept
41 : time, TAUC the correlation time, OMEGA, the larmor frequency and RTWO for the relaxation
42 : time.
43 :
44 : \par Examples
45 :
46 : In the following example five PRE intensities are calculated using the distance between the
47 : oxigen of the spin label and the backbone hydrogens. Omega is the NMR frequency, RTWO the
48 : R2 for the hydrogens, INEPT of 8 ms for the experiment and a TAUC of 1.21 ns
49 :
50 : \verbatim
51 : PRE ...
52 : LABEL=HN_pre
53 : INEPT=8
54 : TAUC=1.21
55 : OMEGA=900
56 : SPINLABEL=1818
57 : GROUPA1=86 RTWO1=0.0120272827
58 : GROUPA2=177 RTWO2=0.0263953158
59 : GROUPA3=285 RTWO3=0.0058899829
60 : GROUPA4=335 RTWO4=0.0102072646
61 : GROUPA5=451 RTWO5=0.0086341843
62 : ... PRE
63 :
64 : PRINT ARG=HN_pre.* FILE=PRE.dat STRIDE=1
65 :
66 : \endverbatim
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : class PRE : public Colvar {
72 : private:
73 : bool pbc;
74 : double constant, inept;
75 : vector<double> rtwo;
76 : vector<unsigned> nga;
77 : NeighborList *nl;
78 : public:
79 : static void registerKeywords( Keywords& keys );
80 : explicit PRE(const ActionOptions&);
81 : ~PRE();
82 : virtual void calculate();
83 : };
84 :
85 2525 : PLUMED_REGISTER_ACTION(PRE,"PRE")
86 :
87 3 : void PRE::registerKeywords( Keywords& keys ) {
88 3 : Colvar::registerKeywords( keys );
89 3 : componentsAreNotOptional(keys);
90 3 : useCustomisableComponents(keys);
91 3 : keys.add("compulsory","INEPT","is the INEPT time (in ms).");
92 3 : keys.add("compulsory","TAUC","is the correlation time (in ns) for this electron-nuclear interaction.");
93 3 : keys.add("compulsory","OMEGA","is the Larmor frequency of the nuclear spin (in MHz).");
94 3 : keys.add("atoms","SPINLABEL","The atom to be used as the paramagnetic center.");
95 : keys.add("numbered","GROUPA","the atoms involved in each of the contacts you wish to calculate. "
96 : "Keywords like GROUPA1, GROUPA2, GROUPA3,... should be listed and one contact will be "
97 3 : "calculated for each ATOM keyword you specify.");
98 3 : keys.reset_style("GROUPA","atoms");
99 : keys.add("numbered","RTWO","The relaxation of the atom/atoms in the corresponding GROUPA of atoms. "
100 3 : "Keywords like RTWO1, RTWO2, RTWO3,... should be listed.");
101 3 : keys.addFlag("ADDEXP",false,"Set to TRUE if you want to have fixed components with the experimetnal values.");
102 3 : keys.add("numbered","PREINT","Add an experimental value for each PRE.");
103 3 : keys.addOutputComponent("pre","default","the # PRE");
104 3 : keys.addOutputComponent("exp","ADDEXP","the # PRE experimental intensity");
105 3 : }
106 :
107 2 : PRE::PRE(const ActionOptions&ao):
108 : PLUMED_COLVAR_INIT(ao),
109 2 : pbc(true)
110 : {
111 2 : bool nopbc=!pbc;
112 2 : parseFlag("NOPBC",nopbc);
113 2 : pbc=!nopbc;
114 :
115 2 : vector<AtomNumber> atom;
116 2 : parseAtomList("SPINLABEL",atom);
117 2 : if(atom.size()!=1) error("Number of specified atom should be 1");
118 :
119 : // Read in the atoms
120 4 : vector<AtomNumber> t, ga_lista, gb_lista;
121 8 : for(int i=1;; ++i ) {
122 8 : parseAtomList("GROUPA", i, t );
123 8 : if( t.empty() ) break;
124 6 : for(unsigned j=0; j<t.size(); j++) {ga_lista.push_back(t[j]); gb_lista.push_back(atom[0]);}
125 6 : nga.push_back(t.size());
126 6 : t.resize(0);
127 6 : }
128 :
129 : // Read in reference values
130 2 : rtwo.resize( nga.size() );
131 2 : unsigned ntarget=0;
132 2 : for(unsigned i=0; i<nga.size(); ++i) {
133 2 : if( !parseNumbered( "RTWO", i+1, rtwo[i] ) ) break;
134 0 : ntarget++;
135 : }
136 2 : if( ntarget==0 ) {
137 2 : parse("RTWO",rtwo[0]);
138 2 : for(unsigned i=1; i<nga.size(); ++i) rtwo[i]=rtwo[0];
139 0 : } else if( ntarget!=nga.size() ) error("found wrong number of RTWO values");
140 :
141 2 : double tauc=0.;
142 2 : parse("TAUC",tauc);
143 2 : if(tauc==0.) error("TAUC must be set");
144 :
145 2 : double omega=0.;
146 2 : parse("OMEGA",omega);
147 2 : if(omega==0.) error("OMEGA must be set");
148 :
149 2 : inept=0.;
150 2 : parse("INEPT",inept);
151 2 : if(inept==0.) error("INEPT must be set");
152 2 : inept *= 0.001; // ms2s
153 :
154 2 : const double ns2s = 0.000000001;
155 2 : const double MHz2Hz = 1000000.;
156 2 : const double Kappa = 12300000000.00; // this is 1/15*S*(S+1)*gamma^2*g^2*beta^2
157 : // where gamma is the nuclear gyromagnetic ratio,
158 : // g is the electronic g factor, and beta is the Bohr magneton
159 : // in nm^6/s^2
160 2 : constant = (4.*tauc*ns2s+(3.*tauc*ns2s)/(1+omega*omega*MHz2Hz*MHz2Hz*tauc*tauc*ns2s*ns2s))*Kappa;
161 :
162 2 : bool addexp=false;
163 2 : parseFlag("ADDEXP",addexp);
164 :
165 4 : vector<double> exppre;
166 2 : if(addexp) {
167 0 : exppre.resize( nga.size() );
168 0 : unsigned ntarget=0;
169 :
170 0 : for(unsigned i=0; i<nga.size(); ++i) {
171 0 : if( !parseNumbered( "PREINT", i+1, exppre[i] ) ) break;
172 0 : ntarget++;
173 : }
174 0 : if( ntarget!=nga.size() ) error("found wrong number of PREINT values");
175 : }
176 :
177 : // Create neighbour lists
178 2 : nl= new NeighborList(gb_lista,ga_lista,true,pbc,getPbc());
179 :
180 : // Ouput details of all contacts
181 2 : unsigned index=0;
182 8 : for(unsigned i=0; i<nga.size(); ++i) {
183 6 : log.printf(" The %uth PRE is calculated using %u equivalent atoms:\n", i, nga[i]);
184 6 : log.printf(" %d", ga_lista[index].serial());
185 6 : index++;
186 8 : for(unsigned j=1; j<nga[i]; j++) {
187 2 : log.printf(" %d", ga_lista[index].serial());
188 2 : index++;
189 : }
190 6 : log.printf("\n");
191 : }
192 :
193 2 : if(pbc) log.printf(" using periodic boundary conditions\n");
194 0 : else log.printf(" without periodic boundary conditions\n");
195 :
196 8 : for(unsigned i=0; i<nga.size(); i++) {
197 6 : string num; Tools::convert(i,num);
198 6 : addComponentWithDerivatives("pre_"+num);
199 6 : componentIsNotPeriodic("pre_"+num);
200 6 : }
201 :
202 2 : if(addexp) {
203 0 : for(unsigned i=0; i<nga.size(); i++) {
204 0 : string num; Tools::convert(i,num);
205 0 : addComponent("exp_"+num);
206 0 : componentIsNotPeriodic("exp_"+num);
207 0 : Value* comp=getPntrToComponent("exp_"+num);
208 0 : comp->set(exppre[i]);
209 0 : }
210 : }
211 :
212 2 : requestAtoms(nl->getFullAtomList());
213 4 : checkRead();
214 2 : }
215 :
216 8 : PRE::~PRE() {
217 2 : delete nl;
218 6 : }
219 :
220 175 : void PRE::calculate()
221 : {
222 : // cycle over the number of PRE
223 523 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
224 348 : for(unsigned i=0; i<nga.size(); i++) {
225 525 : vector<Vector> deriv;
226 524 : Tensor dervir;
227 525 : double pre=0;
228 525 : unsigned index=0;
229 1050 : for(unsigned k=0; k<i; k++) index+=nga[k];
230 : // cycle over equivalent atoms
231 1225 : for(unsigned j=0; j<nga[i]; j++) {
232 700 : const double c_aver=constant/((double)nga[i]);
233 : // the first atom is always the same (the paramagnetic group)
234 700 : const unsigned i0=nl->getClosePair(index+j).first;
235 700 : const unsigned i1=nl->getClosePair(index+j).second;
236 :
237 700 : Vector distance;
238 1396 : if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
239 0 : else distance=delta(getPosition(i0),getPosition(i1));
240 :
241 700 : const double r2=distance.modulo2();
242 698 : const double r6=r2*r2*r2;
243 698 : const double r8=r6*r2;
244 698 : const double tmpir6=c_aver/r6;
245 698 : const double tmpir8=-6.*c_aver/r8;
246 :
247 698 : pre += tmpir6;
248 :
249 698 : Vector tmpv = -tmpir8*distance;
250 699 : deriv.push_back(tmpv);
251 700 : dervir += Tensor(distance,tmpv);
252 : }
253 525 : const double ratio = rtwo[i]*exp(-pre*inept) / (rtwo[i]+pre);
254 525 : const double fact = -ratio*(inept+1./(rtwo[i]+pre));
255 :
256 525 : Value* val=getPntrToComponent(i);
257 524 : val->set(ratio);
258 524 : setBoxDerivatives(val, fact*dervir);
259 :
260 1225 : for(unsigned j=0; j<nga[i]; j++) {
261 700 : const unsigned i0=nl->getClosePair(index+j).first;
262 700 : const unsigned i1=nl->getClosePair(index+j).second;
263 700 : setAtomsDerivatives(val, i0, fact*deriv[j]);
264 700 : setAtomsDerivatives(val, i1, -fact*deriv[j]);
265 525 : }
266 : }
267 175 : }
268 :
269 : }
270 2523 : }
|
