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