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