Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 : #include <memory>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace isdb {
35 :
36 : //+PLUMEDOC ISDB_COLVAR NOE
37 : /*
38 : Calculates NOE intensities as sums of 1/r^6, also averaging over multiple equivalent atoms
39 : or ambiguous NOE.
40 :
41 : Each NOE is defined by two groups containing the same number of atoms, distances are
42 : calculated in pairs, transformed in 1/r^6, summed and saved as components.
43 :
44 : \f[
45 : NOE() = (\frac{1}{N_{eq}}\sum_j^{N_{eq}} (\frac{1}{r_j^6}))
46 : \f]
47 :
48 : NOE can be used to calculate a Metainference score over one or more replicas using the intrinsic implementation
49 : of \ref METAINFERENCE that is activated by DOSCORE.
50 :
51 : \par Examples
52 : In the following examples three noes are defined, the first is calculated based on the distances
53 : of atom 1-2 and 3-2; the second is defined by the distance 5-7 and the third by the distances
54 : 4-15,4-16,8-15,8-16. \ref METAINFERENCE is activated using DOSCORE.
55 :
56 : \plumedfile
57 : NOE ...
58 : GROUPA1=1,3 GROUPB1=2,2
59 : GROUPA2=5 GROUPB2=7
60 : GROUPA3=4,4,8,8 GROUPB3=15,16,15,16
61 : DOSCORE
62 : LABEL=noes
63 : ... NOE
64 :
65 : PRINT ARG=noes.* FILE=colvar
66 : \endplumedfile
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 33 : class NOE :
72 : public MetainferenceBase
73 : {
74 : private:
75 : bool pbc;
76 : vector<unsigned> nga;
77 : std::unique_ptr<NeighborList> nl;
78 : unsigned tot_size;
79 : public:
80 : static void registerKeywords( Keywords& keys );
81 : explicit NOE(const ActionOptions&);
82 : void calculate() override;
83 : void update() override;
84 : };
85 :
86 7854 : PLUMED_REGISTER_ACTION(NOE,"NOE")
87 :
88 12 : void NOE::registerKeywords( Keywords& keys ) {
89 12 : componentsAreNotOptional(keys);
90 12 : MetainferenceBase::registerKeywords(keys);
91 36 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
92 : keys.add("numbered","GROUPA","the atoms involved in each of the contacts you wish to calculate. "
93 : "Keywords like GROUPA1, GROUPA2, GROUPA3,... should be listed and one contact will be "
94 48 : "calculated for each ATOM keyword you specify.");
95 : keys.add("numbered","GROUPB","the atoms involved in each of the contacts you wish to calculate. "
96 : "Keywords like GROUPB1, GROUPB2, GROUPB3,... should be listed and one contact will be "
97 48 : "calculated for each ATOM keyword you specify.");
98 36 : keys.reset_style("GROUPA","atoms");
99 36 : keys.reset_style("GROUPB","atoms");
100 48 : keys.add("numbered","NOEDIST","Add an experimental value for each NOE.");
101 48 : keys.addOutputComponent("noe","default","the # NOE");
102 48 : keys.addOutputComponent("exp","NOEDIST","the # NOE experimental distance");
103 12 : }
104 :
105 11 : NOE::NOE(const ActionOptions&ao):
106 : PLUMED_METAINF_INIT(ao),
107 22 : pbc(true)
108 : {
109 11 : bool nopbc=!pbc;
110 22 : parseFlag("NOPBC",nopbc);
111 11 : pbc=!nopbc;
112 :
113 : // Read in the atoms
114 : vector<AtomNumber> t, ga_lista, gb_lista;
115 22 : for(int i=1;; ++i ) {
116 66 : parseAtomList("GROUPA", i, t );
117 33 : if( t.empty() ) break;
118 88 : for(unsigned j=0; j<t.size(); j++) ga_lista.push_back(t[j]);
119 44 : nga.push_back(t.size());
120 22 : t.resize(0);
121 22 : }
122 : vector<unsigned> ngb;
123 22 : for(int i=1;; ++i ) {
124 66 : parseAtomList("GROUPB", i, t );
125 33 : if( t.empty() ) break;
126 88 : for(unsigned j=0; j<t.size(); j++) gb_lista.push_back(t[j]);
127 44 : ngb.push_back(t.size());
128 66 : if(ngb[i-1]!=nga[i-1]) error("The same number of atoms is expected for the same GROUPA-GROUPB couple");
129 22 : t.resize(0);
130 22 : }
131 11 : if(nga.size()!=ngb.size()) error("There should be the same number of GROUPA and GROUPB keywords");
132 : // Create neighbour lists
133 22 : nl.reset( new NeighborList(ga_lista,gb_lista,true,pbc,getPbc()) );
134 :
135 : // Optionally add an experimental value (like with RDCs)
136 : vector<double> noedist;
137 11 : noedist.resize( nga.size() );
138 : unsigned ntarget=0;
139 58 : for(unsigned i=0; i<nga.size(); ++i) {
140 40 : if( !parseNumbered( "NOEDIST", i+1, noedist[i] ) ) break;
141 18 : ntarget++;
142 : }
143 : bool addexp=false;
144 22 : if(ntarget!=nga.size() && ntarget!=0) error("found wrong number of NOEDIST values");
145 11 : if(ntarget==nga.size()) addexp=true;
146 11 : if(getDoScore()&&!addexp) error("with DOSCORE you need to set the NOEDIST values");
147 :
148 : // Ouput details of all contacts
149 : unsigned index=0;
150 55 : for(unsigned i=0; i<nga.size(); ++i) {
151 22 : log.printf(" The %uth NOE is calculated using %u equivalent couples of atoms\n", i, nga[i]);
152 88 : for(unsigned j=0; j<nga[i]; j++) {
153 66 : log.printf(" couple %u is %d %d.\n", j, ga_lista[index].serial(), gb_lista[index].serial() );
154 33 : index++;
155 : }
156 : }
157 11 : tot_size = index;
158 :
159 11 : if(pbc) log.printf(" using periodic boundary conditions\n");
160 0 : else log.printf(" without periodic boundary conditions\n");
161 :
162 33 : log << " Bibliography" << plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
163 :
164 11 : if(!getDoScore()) {
165 35 : for(unsigned i=0; i<nga.size(); i++) {
166 14 : string num; Tools::convert(i,num);
167 28 : addComponentWithDerivatives("noe-"+num);
168 28 : componentIsNotPeriodic("noe-"+num);
169 : }
170 7 : if(addexp) {
171 25 : for(unsigned i=0; i<nga.size(); i++) {
172 10 : string num; Tools::convert(i,num);
173 20 : addComponent("exp-"+num);
174 20 : componentIsNotPeriodic("exp-"+num);
175 20 : Value* comp=getPntrToComponent("exp-"+num);
176 10 : comp->set(noedist[i]);
177 : }
178 : }
179 : } else {
180 20 : for(unsigned i=0; i<nga.size(); i++) {
181 8 : string num; Tools::convert(i,num);
182 16 : addComponent("noe-"+num);
183 16 : componentIsNotPeriodic("noe-"+num);
184 : }
185 20 : for(unsigned i=0; i<nga.size(); i++) {
186 8 : string num; Tools::convert(i,num);
187 16 : addComponent("exp-"+num);
188 16 : componentIsNotPeriodic("exp-"+num);
189 16 : Value* comp=getPntrToComponent("exp-"+num);
190 8 : comp->set(noedist[i]);
191 : }
192 : }
193 :
194 11 : requestAtoms(nl->getFullAtomList(), false);
195 11 : if(getDoScore()) {
196 4 : setParameters(noedist);
197 4 : Initialise(nga.size());
198 : }
199 11 : setDerivatives();
200 11 : checkRead();
201 11 : }
202 :
203 456 : void NOE::calculate()
204 : {
205 456 : const unsigned ngasz=nga.size();
206 456 : vector<Vector> deriv(tot_size, Vector{0,0,0});
207 :
208 1366 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
209 : for(unsigned i=0; i<ngasz; i++) {
210 906 : Tensor dervir;
211 : double noe=0;
212 : unsigned index=0;
213 1821 : for(unsigned k=0; k<i; k++) index+=nga[k];
214 909 : string num; Tools::convert(i,num);
215 1813 : Value* val=getPntrToComponent("noe-"+num);
216 : // cycle over equivalent atoms
217 4552 : for(unsigned j=0; j<nga[i]; j++) {
218 2728 : const unsigned i0=nl->getClosePair(index+j).first;
219 1366 : const unsigned i1=nl->getClosePair(index+j).second;
220 :
221 1368 : Vector distance;
222 5468 : if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
223 0 : else distance=delta(getPosition(i0),getPosition(i1));
224 :
225 1366 : const double ir2=1./distance.modulo2();
226 1368 : const double ir6=ir2*ir2*ir2;
227 1368 : const double ir8=6*ir6*ir2;
228 :
229 1368 : noe += ir6;
230 2736 : deriv[index+j] = ir8*distance;
231 1367 : if(!getDoScore()) {
232 2446 : dervir += Tensor(distance, deriv[index+j]);
233 2448 : setAtomsDerivatives(val, i0, deriv[index+j]);
234 2446 : setAtomsDerivatives(val, i1, -deriv[index+j]);
235 : }
236 : }
237 : val->set(noe);
238 911 : if(!getDoScore()) {
239 816 : setBoxDerivatives(val, dervir);
240 : } else setCalcData(i, noe);
241 : }
242 :
243 456 : if(getDoScore()) {
244 : /* Metainference */
245 48 : Tensor dervir;
246 48 : double score = getScore();
247 : setScore(score);
248 :
249 : /* calculate final derivatives */
250 96 : Value* val=getPntrToComponent("score");
251 144 : for(unsigned i=0; i<ngasz; i++) {
252 : unsigned index=0;
253 96 : for(unsigned k=0; k<i; k++) index+=nga[k];
254 : // cycle over equivalent atoms
255 384 : for(unsigned j=0; j<nga[i]; j++) {
256 288 : const unsigned i0=nl->getClosePair(index+j).first;
257 144 : const unsigned i1=nl->getClosePair(index+j).second;
258 :
259 144 : Vector distance;
260 432 : if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
261 0 : else distance=delta(getPosition(i0),getPosition(i1));
262 :
263 288 : dervir += Tensor(distance,deriv[index+j]*getMetaDer(i));
264 144 : setAtomsDerivatives(val, i0, deriv[index+j]*getMetaDer(i));
265 144 : setAtomsDerivatives(val, i1, -deriv[index+j]*getMetaDer(i));
266 : }
267 : }
268 48 : setBoxDerivatives(val, dervir);
269 : }
270 456 : }
271 :
272 132 : void NOE::update() {
273 : // write status file
274 264 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
275 132 : }
276 :
277 : }
278 5874 : }
|
