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