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