Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "Colvar.h"
23 : #include "tools/Communicator.h"
24 : #include "tools/NeighborList.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/SwitchingFunction.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR CONTACTMAP
32 : /*
33 : Calculate the distances between a number of pairs of atoms and transform each distance by a switching function.
34 :
35 : This CV calculates a series of distances between pairs of atoms and transforms them by a switching function.
36 : The transformed distance can be compared with a reference value in order to calculate the squared distance
37 : between two contact maps. Each distance can also be weighted for a given value. CONTACTMAP can be used together
38 : with [FUNCPATHMSD](FUNCPATHMSD.md) to define a path in contactmap space as was done in the 2008 paper by Bonomi that
39 : is cited below.
40 :
41 : The individual contact map distances related to each contact can be accessed as components
42 : named `cm.contact-1`, `cm.contact-2`, etc, assuming that the label of the CONTACTMAP is `cm`.
43 :
44 : ## Examples
45 :
46 : The following example calculates switching functions based on the distances between atoms
47 : 1 and 2, 3 and 4 and 4 and 5. The values of these three switching functions are then output
48 : to a file named colvar.
49 :
50 : ```plumed
51 : f1: CONTACTMAP ATOMS1=1,2 ATOMS2=3,4 ATOMS3=4,5 ATOMS4=5,6 SWITCH={RATIONAL R_0=1.5}
52 : PRINT ARG=f1.* FILE=colvar
53 : ```
54 :
55 : The following example calculates the difference of the current contact map with respect
56 : to a reference provided. In this case REFERENCE is the fraction of contact that is formed
57 : (i.e. the distance between two atoms transformed with the SWITCH), while R_0 is the contact
58 : distance. WEIGHT gives the relative weight of each contact to the final distance measure.
59 :
60 : ```plumed
61 : cmap: CONTACTMAP ...
62 : ATOMS1=1,2 REFERENCE1=0.1 WEIGHT1=0.5
63 : ATOMS2=3,4 REFERENCE2=0.5 WEIGHT2=1.0
64 : ATOMS3=4,5 REFERENCE3=0.25 WEIGHT3=1.0
65 : ATOMS4=5,6 REFERENCE4=0.0 WEIGHT4=0.5
66 : SWITCH={RATIONAL R_0=1.5}
67 : CMDIST
68 : ...
69 :
70 : PRINT ARG=cmap FILE=colvar
71 : ```
72 :
73 : The next example calculates calculates fraction of native contacts (Q)
74 : for Trp-cage mini-protein. R_0 is the distance at which the switch function is guaranteed to
75 : be 1.0 – it doesn't really matter for Q and should be something very small, like 1 A.
76 : REF is the reference distance for the contact, e.g. the distance from a crystal structure.
77 : LAMBDA is the tolerance for the distance – if set to 1.0, the contact would have to have exactly
78 : the reference value to be formed; instead for lambda values of 1.5–1.8 are usually used to allow some slack.
79 : BETA is the softness of the switch function, default is 50nm.
80 : WEIGHT is the 1/(number of contacts) giving equal weight to each contact.
81 :
82 : When using native contact Q switch function, please cite the 2013 paper from Best from the bibliography below
83 :
84 : ```plumed
85 : # The full (much-longer) example available in regtest/basic/rt72/
86 :
87 : cmap: CONTACTMAP ...
88 : ATOMS1=1,67 SWITCH1={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.4059} WEIGHT1=0.003597
89 : ATOMS2=1,68 SWITCH2={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.4039} WEIGHT2=0.003597
90 : ATOMS3=1,69 SWITCH3={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.3215} WEIGHT3=0.003597
91 : ATOMS4=5,61 SWITCH4={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.4277} WEIGHT4=0.003597
92 : ATOMS5=5,67 SWITCH5={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.3851} WEIGHT5=0.003597
93 : ATOMS6=5,68 SWITCH6={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.3811} WEIGHT6=0.003597
94 : ATOMS7=5,69 SWITCH7={Q R_0=0.01 BETA=50.0 LAMBDA=1.5 REF=0.3133} WEIGHT7=0.003597
95 : SUM
96 : ...
97 :
98 : PRINT ARG=cmap FILE=colvar
99 : ```
100 :
101 : */
102 : //+ENDPLUMEDOC
103 :
104 : class ContactMap : public Colvar {
105 : private:
106 : bool pbc, serial, docomp, dosum, docmdist;
107 : std::unique_ptr<NeighborList> nl;
108 : std::vector<SwitchingFunction> sfs;
109 : std::vector<double> reference, weight;
110 : public:
111 : static void registerKeywords( Keywords& keys );
112 : explicit ContactMap(const ActionOptions&);
113 : // active methods:
114 : void calculate() override;
115 0 : void checkFieldsAllowed() override {}
116 : };
117 :
118 : PLUMED_REGISTER_ACTION(ContactMap,"CONTACTMAP")
119 :
120 21 : void ContactMap::registerKeywords( Keywords& keys ) {
121 21 : Colvar::registerKeywords( keys );
122 21 : keys.add("numbered","ATOMS","the atoms involved in each of the contacts you wish to calculate. "
123 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one contact will be "
124 : "calculated for each ATOM keyword you specify.");
125 42 : keys.reset_style("ATOMS","atoms");
126 21 : keys.add("numbered","SWITCH","The switching functions to use for each of the contacts in your map. "
127 : "You can either specify a global switching function using SWITCH or one "
128 : "switching function for each contact. Details of the various switching "
129 : "functions you can use are provided in the documentation for LESS_THAN.");
130 21 : keys.add("numbered","REFERENCE","A reference value for a given contact, by default is 0.0 "
131 : "You can either specify a global reference value using REFERENCE or one "
132 : "reference value for each contact.");
133 21 : keys.add("numbered","WEIGHT","A weight value for a given contact, by default is 1.0 "
134 : "You can either specify a global weight value using WEIGHT or one "
135 : "weight value for each contact.");
136 42 : keys.reset_style("SWITCH","compulsory");
137 42 : keys.linkActionInDocs("SWITCH","LESS_THAN");
138 21 : keys.addFlag("SUM",false,"calculate the sum of all the contacts in the input");
139 21 : keys.addFlag("CMDIST",false,"calculate the distance with respect to the provided reference contact map");
140 21 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
141 42 : keys.addOutputComponent("contact","default","scalar","By not using SUM or CMDIST each contact will be stored in a component");
142 42 : keys.setValueDescription("scalar","the sum of all the switching function on all the distances");
143 21 : keys.addDOI("10.1021/ja803652f");
144 21 : keys.addDOI("10.1073/pnas.1311599110");
145 21 : }
146 :
147 17 : ContactMap::ContactMap(const ActionOptions&ao):
148 : PLUMED_COLVAR_INIT(ao),
149 17 : pbc(true),
150 17 : serial(false),
151 17 : docomp(true),
152 17 : dosum(false),
153 17 : docmdist(false) {
154 17 : parseFlag("SERIAL",serial);
155 17 : parseFlag("SUM",dosum);
156 17 : parseFlag("CMDIST",docmdist);
157 17 : if(docmdist==true&&dosum==true) {
158 0 : error("You cannot use SUM and CMDIST together");
159 : }
160 17 : bool nopbc=!pbc;
161 17 : parseFlag("NOPBC",nopbc);
162 17 : pbc=!nopbc;
163 :
164 : // Read in the atoms
165 17 : std::vector<AtomNumber> t, ga_lista, gb_lista;
166 17 : for(int i=1;; ++i ) {
167 2364 : parseAtomList("ATOMS", i, t );
168 1182 : if( t.empty() ) {
169 : break;
170 : }
171 :
172 1165 : if( t.size()!=2 ) {
173 : std::string ss;
174 0 : Tools::convert(i,ss);
175 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
176 : }
177 1165 : ga_lista.push_back(t[0]);
178 1165 : gb_lista.push_back(t[1]);
179 1165 : t.resize(0);
180 :
181 : // Add a value for this contact
182 : std::string num;
183 1165 : Tools::convert(i,num);
184 1165 : if(!dosum&&!docmdist) {
185 10 : addComponentWithDerivatives("contact-"+num);
186 10 : componentIsNotPeriodic("contact-"+num);
187 : }
188 1165 : }
189 17 : if( ga_lista.empty() ) {
190 1 : error("no contacts have been specified, use the ATOMS1, ATOMS2, ... keywords");
191 : }
192 : // Create neighbour lists
193 32 : nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,true,pbc,getPbc(),comm);
194 :
195 : // Read in switching functions
196 : std::string errors;
197 16 : sfs.resize( ga_lista.size() );
198 : unsigned nswitch=0;
199 1180 : for(unsigned i=0; i<ga_lista.size(); ++i) {
200 : std::string num, sw1;
201 1165 : Tools::convert(i+1, num);
202 2330 : if( !parseNumbered( "SWITCH", i+1, sw1 ) ) {
203 : break;
204 : }
205 1164 : nswitch++;
206 1164 : sfs[i].set(sw1,errors);
207 1164 : if( errors.length()!=0 ) {
208 0 : error("problem reading SWITCH" + num + " keyword : " + errors );
209 : }
210 : }
211 16 : if( nswitch==0 ) {
212 : std::string sw;
213 2 : parse("SWITCH",sw);
214 1 : if(sw.length()==0) {
215 0 : error("no switching function specified use SWITCH keyword");
216 : }
217 2 : for(unsigned i=0; i<ga_lista.size(); ++i) {
218 1 : sfs[i].set(sw,errors);
219 1 : if( errors.length()!=0 ) {
220 0 : error("problem reading SWITCH keyword : " + errors );
221 : }
222 : }
223 15 : } else if( nswitch!=sfs.size() ) {
224 : std::string num;
225 0 : Tools::convert(nswitch+1, num);
226 0 : error("missing SWITCH" + num + " keyword");
227 : }
228 :
229 : // Read in reference values
230 : nswitch=0;
231 16 : reference.resize(ga_lista.size(), 0.);
232 21 : for(unsigned i=0; i<ga_lista.size(); ++i) {
233 40 : if( !parseNumbered( "REFERENCE", i+1, reference[i] ) ) {
234 : break;
235 : }
236 5 : nswitch++;
237 : }
238 16 : if( nswitch==0 ) {
239 : // a single value (REFERENCE) for all the contacts
240 : std::string global_ref;
241 30 : parse("REFERENCE",global_ref);
242 15 : if( global_ref.length()>0 ) {
243 1 : double ref0=0.0;
244 1 : if( !Tools::convertNoexcept(global_ref,ref0) ) {
245 0 : error("cannot understand the value " + global_ref + " given to the REFERENCE keyword");
246 : }
247 2 : for(unsigned i=0; i<ga_lista.size(); ++i) {
248 1 : reference[i]=ref0;
249 : }
250 1 : nswitch = ga_lista.size();
251 : }
252 : }
253 15 : if(nswitch == 0 && docmdist) {
254 0 : error("with CMDIST one must use REFERENCE to setup the reference contact map");
255 : }
256 :
257 : // Read in weight values
258 : nswitch=0;
259 16 : weight.resize(ga_lista.size(), 1.0);
260 1175 : for(unsigned i=0; i<ga_lista.size(); ++i) {
261 2322 : if( !parseNumbered( "WEIGHT", i+1, weight[i] ) ) {
262 : break;
263 : }
264 1159 : nswitch++;
265 : }
266 16 : if( nswitch==0 ) {
267 2 : parse("WEIGHT",weight[0]);
268 6 : for(unsigned i=1; i<ga_lista.size(); ++i) {
269 4 : weight[i]=weight[0];
270 : }
271 : nswitch = ga_lista.size();
272 : }
273 :
274 : // Output details of all contacts
275 1181 : for(unsigned i=0; i<sfs.size(); ++i) {
276 1165 : log.printf(" The %uth contact is calculated from atoms : %d %d. Inflection point of switching function is at %s. Reference contact value is %f\n",
277 2330 : i+1, ga_lista[i].serial(), gb_lista[i].serial(), ( sfs[i].description() ).c_str(), reference[i] );
278 : }
279 :
280 16 : if(dosum) {
281 13 : addValueWithDerivatives();
282 13 : setNotPeriodic();
283 13 : log.printf(" colvar is sum of all contacts in contact map\n");
284 : }
285 16 : if(docmdist) {
286 2 : addValueWithDerivatives();
287 2 : setNotPeriodic();
288 2 : log.printf(" colvar is distance between the contact map matrix and the provided reference matrix\n");
289 : }
290 :
291 16 : if(dosum || docmdist) {
292 15 : docomp=false;
293 : } else {
294 1 : serial=true;
295 1 : docomp=true;
296 : }
297 :
298 : // Set up if it is just a list of contacts
299 16 : requestAtoms(nl->getFullAtomList());
300 16 : checkRead();
301 25 : }
302 :
303 6027 : void ContactMap::calculate() {
304 :
305 6027 : double ncoord=0.;
306 6027 : Tensor virial;
307 6027 : std::vector<Vector> deriv(getNumberOfAtoms());
308 :
309 : unsigned stride;
310 : unsigned rank;
311 6027 : if(serial) {
312 : // when using components the parallelisation do not work
313 : stride=1;
314 : rank=0;
315 : } else {
316 6022 : stride=comm.Get_size();
317 6022 : rank=comm.Get_rank();
318 : }
319 :
320 : // sum over close pairs
321 445235 : for(unsigned i=rank; i<nl->size(); i+=stride) {
322 439208 : Vector distance;
323 439208 : unsigned i0=nl->getClosePair(i).first;
324 439208 : unsigned i1=nl->getClosePair(i).second;
325 439208 : if(pbc) {
326 439208 : distance=pbcDistance(getPosition(i0),getPosition(i1));
327 : } else {
328 0 : distance=delta(getPosition(i0),getPosition(i1));
329 : }
330 :
331 439208 : double dfunc=0.;
332 439208 : double coord = weight[i]*(sfs[i].calculate(distance.modulo(), dfunc) - reference[i]);
333 439208 : Vector tmpder = weight[i]*dfunc*distance;
334 439208 : Tensor tmpvir = weight[i]*dfunc*Tensor(distance,distance);
335 439208 : if(!docmdist) {
336 439179 : deriv[i0] -= tmpder;
337 439179 : deriv[i1] += tmpder;
338 : virial -= tmpvir;
339 439179 : ncoord += coord;
340 : } else {
341 29 : tmpder *= 2.*coord;
342 : tmpvir *= 2.*coord;
343 29 : deriv[i0] -= tmpder;
344 29 : deriv[i1] += tmpder;
345 : virial -= tmpvir;
346 29 : ncoord += coord*coord;
347 : }
348 :
349 439208 : if(docomp) {
350 25 : Value* val=getPntrToComponent( i );
351 25 : setAtomsDerivatives( val, i0, deriv[i0] );
352 25 : setAtomsDerivatives( val, i1, deriv[i1] );
353 25 : setBoxDerivatives( val, -tmpvir );
354 : val->set(coord);
355 : }
356 : }
357 :
358 6027 : if(!serial) {
359 6022 : comm.Sum(&ncoord,1);
360 6022 : if(!deriv.empty()) {
361 6022 : comm.Sum(&deriv[0][0],3*deriv.size());
362 : }
363 6022 : comm.Sum(&virial[0][0],9);
364 : }
365 :
366 6027 : if( !docomp ) {
367 884388 : for(unsigned i=0; i<deriv.size(); ++i) {
368 878366 : setAtomsDerivatives(i,deriv[i]);
369 : }
370 6022 : setValue (ncoord);
371 6022 : setBoxDerivatives (virial);
372 : }
373 6027 : }
374 :
375 : }
376 : }
|
