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