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