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 "SecondaryStructureBase.h"
23 : #include "core/ActionShortcut.h"
24 : #include "core/ActionRegister.h"
25 :
26 : namespace PLMD {
27 : namespace secondarystructure {
28 :
29 : //+PLUMEDOC COLVAR PARABETARMSD
30 : /*
31 : Probe the parallel beta sheet content of your protein structure.
32 :
33 : Two protein segments containing three contiguous residues can form a parallel beta sheet.
34 : Although if the two segments are part of the same protein chain they must be separated by
35 : a minimum of 3 residues to make room for the turn. This colvar thus generates the set of
36 : all possible six residue sections that could conceivably form a parallel beta sheet
37 : and calculates the [DRMSD](DRMSD.md) or [RMSD](RMSD.md) distance between the configuration in which the residues find themselves
38 : and an idealized parallel beta sheet structure. These distances can be calculated by either
39 : aligning the instantaneous structure with the reference structure and measuring each
40 : atomic displacement or by calculating differences between the set of inter-atomic
41 : distances in the reference and instantaneous structures.
42 :
43 : This colvar is based on the reference in the bibliography below. The authors of
44 : this paper use the set of distances from the parallel beta sheet configurations to measure
45 : the number of segments whose configuration resembles a parallel beta sheet. This is done by calculating
46 : the following sum of functions of the rmsd distances:
47 :
48 : $$
49 : s = \sum_i \frac{ 1 - \left(\frac{r_i-d_0}{r_0}\right)^n } { 1 - \left(\frac{r_i-d_0}{r_0}\right)^m }
50 : $$
51 :
52 : where the sum runs over all possible segments of parallel beta sheet. By default the
53 : NN, MM and D_0 parameters are set equal to those used in the paper cited below. The R_0
54 : parameter must be set by the user - the value used in the paper cited below was 0.08 nm.
55 :
56 : If you change the function in the above sum you can calculate quantities such as the average
57 : distance from a structure composed of only parallel beta sheets or the distance between the set of
58 : residues that is closest to a parallel beta sheet and the reference configuration. To do these sorts of
59 : calculations you can use the AVERAGE and MIN keywords. In addition you can use the LESS_THAN
60 : keyword if you would like to change the form of the switching function. If you use any of these
61 : options you no longer need to specify NN, R_0, MM and D_0.
62 :
63 : The following input calculates the number of six residue segments of
64 : protein that are in an parallel beta sheet configuration.
65 :
66 : ```plumed
67 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
68 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
69 : pb: PARABETARMSD RESIDUES=all STRANDS_CUTOFF=1 R_0=0.1
70 : PRINT ARG=pb FILE=colvar
71 : ```
72 :
73 : Here the same is done use [RMSD](RMSD.md) instead of [DRMSD](DRMSD.md), which is normally faster.
74 :
75 : ```plumed
76 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
77 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
78 : WHOLEMOLECULES ENTITY0=1-100
79 : hh: PARABETARMSD RESIDUES=all TYPE=OPTIMAL LESS_THAN={RATIONAL R_0=0.1 NN=8 MM=12} STRANDS_CUTOFF=1
80 : PRINT ARG=hh.lessthan FILE=colvar
81 : ```
82 :
83 : __YOUR CALCULATION WILL BE MUCH FASTER IF YOU USE THE `STRANDS_CUTOFF` KEYWORD.__ As you can see from the
84 : expanded version of the inputs above this keyword reduces the computational cost of the calculation by
85 : avoiding calculations of the RMSD values for segments that have the two strands of the beta sheet further apart
86 : than a cutoff.
87 :
88 : ## Periodic boundary conditions
89 :
90 : You can turn off periodic boundary conditions by using the NOPBC flag as shown below:
91 :
92 : ```plumed
93 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
94 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
95 : pb: PARABETARMSD RESIDUES=all STRANDS_CUTOFF=1 R_0=0.1 NOPBC
96 : PRINT ARG=pb FILE=colvar
97 : ```
98 :
99 : If you are using [DRMSD](DRMSD.md) to measure distances and you __don't__ use the NOPBC flag then
100 : all distances in the instantaneous structure are evaluated in a way that takes the periodic boundary conditions
101 : into account. Using the NOPBC flag turns off this treatment.
102 :
103 : If you are using [RMSD](RMSD.md) to measure distances and you __don't__ use the NOPBC flag the instantaneous positions of
104 : each segment for which the RMSD is computed is reconstructed using the procedure outlined in the documentation for [WHOLEMOLECULES](WHOLEMOLECULES.md)
105 : before the RMSD is computed. Using the NOPBC flag turns off this treatment.
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : class ParabetaRMSD : public ActionShortcut {
111 : public:
112 : static void registerKeywords( Keywords& keys );
113 : explicit ParabetaRMSD(const ActionOptions&);
114 : };
115 :
116 : PLUMED_REGISTER_ACTION(ParabetaRMSD,"PARABETARMSD")
117 :
118 129 : void ParabetaRMSD::registerKeywords( Keywords& keys ) {
119 129 : SecondaryStructureInput::registerKeywords( keys );
120 129 : keys.remove("ATOMS");
121 129 : keys.remove("SEGMENT");
122 129 : keys.remove("STRUCTURE");
123 129 : keys.remove("MASK");
124 129 : keys.remove("ALIGN_STRANDS");
125 258 : keys.setValueDescription("scalar/vector","if LESS_THAN is present the RMSD distance between each residue and the ideal parallel beta sheet. If LESS_THAN is not present the number of residue segments where the structure is similar to a parallel beta sheet");
126 129 : keys.add("compulsory","STYLE","all","Parallel beta sheets can either form in a single chain or from a pair of chains. If STYLE=all all "
127 : "chain configuration with the appropriate geometry are counted. If STYLE=inter "
128 : "only sheet-like configurations involving two chains are counted, while if STYLE=intra "
129 : "only sheet-like configurations involving a single chain are counted");
130 129 : keys.add("optional","STRANDS_CUTOFF","If in a segment of protein the two strands are further apart then the calculation "
131 : "of the actual RMSD is skipped as the structure is very far from being beta-sheet like. "
132 : "This keyword speeds up the calculation enormously when you are using the LESS_THAN option. "
133 : "However, if you are using some other option, then this cannot be used");
134 129 : keys.needsAction("LOWEST");
135 129 : keys.needsAction("DISTANCE");
136 129 : keys.needsAction("CUSTOM");
137 129 : }
138 :
139 27 : ParabetaRMSD::ParabetaRMSD(const ActionOptions&ao):
140 : Action(ao),
141 27 : ActionShortcut(ao) {
142 : // Read in the input and create a string that describes how to compute the less than
143 : std::string ltmap;
144 27 : bool uselessthan=SecondaryStructureInput::readShortcutWords( ltmap, this );
145 : // read in the backbone atoms
146 27 : std::vector<unsigned> chains;
147 27 : std::vector<std::string> all_atoms;
148 54 : SecondaryStructureInput::readBackboneAtoms( this, plumed, "protein", chains, all_atoms );
149 :
150 : bool intra_chain(false), inter_chain(false);
151 : std::string style;
152 54 : parse("STYLE",style);
153 27 : if( Tools::caseInSensStringCompare(style, "all") ) {
154 : intra_chain=true;
155 : inter_chain=true;
156 2 : } else if( Tools::caseInSensStringCompare(style, "inter") ) {
157 : intra_chain=false;
158 : inter_chain=true;
159 0 : } else if( Tools::caseInSensStringCompare(style, "intra") ) {
160 : intra_chain=true;
161 : inter_chain=false;
162 : } else {
163 0 : error( style + " is not a valid directive for the STYLE keyword");
164 : }
165 :
166 27 : double strands_cutoff=-1.0;
167 54 : parse("STRANDS_CUTOFF",strands_cutoff);
168 : std::string scutoff_action;
169 27 : if( strands_cutoff>0 ) {
170 46 : scutoff_action=getShortcutLabel() + "_cut_dists: DISTANCE ";
171 : }
172 :
173 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
174 : std::string seglist;
175 27 : unsigned jjkk=1;
176 27 : if( intra_chain ) {
177 : unsigned nprevious=0;
178 25 : std::vector<unsigned> nlist(30);
179 357 : for(unsigned i=0; i<chains.size(); ++i) {
180 332 : if( chains[i]<40 ) {
181 0 : error("segment of backbone is not long enough to form an antiparallel beta hairpin. Each backbone fragment must contain a minimum of 8 residues");
182 : }
183 : // Loop over all possible triples in each 8 residue segment of protein
184 332 : unsigned nres=chains[i]/5;
185 332 : if( chains[i]%5!=0 ) {
186 0 : error("backbone segment received does not contain a multiple of five residues");
187 : }
188 364 : for(unsigned ires=0; ires<nres-8; ires++) {
189 112 : for(unsigned jres=ires+6; jres<nres-2; jres++) {
190 1280 : for(unsigned k=0; k<15; ++k) {
191 1200 : nlist[k]=nprevious + ires*5+k;
192 1200 : nlist[k+15]=nprevious + jres*5+k;
193 : }
194 : std::string nlstr, num;
195 80 : Tools::convert( nlist[0], nlstr );
196 80 : Tools::convert(jjkk, num);
197 80 : jjkk++;
198 80 : seglist += " SEGMENT" + num + "=" + nlstr;
199 2400 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
200 2320 : Tools::convert( nlist[kk], nlstr );
201 2320 : seglist += "," + nlstr;
202 : }
203 80 : if( strands_cutoff>0 ) {
204 160 : scutoff_action += " ATOMS" + num + "=" + all_atoms[nlist[6]] + "," + all_atoms[nlist[21]];
205 : }
206 : }
207 : }
208 332 : nprevious+=chains[i];
209 : }
210 : }
211 : // This constructs all conceivable sections of antibeta sheet that form between chains
212 27 : if( inter_chain ) {
213 27 : if( chains.size()==1 && !Tools::caseInSensStringCompare(style, "all") ) {
214 0 : error("there is only one chain defined so cannot use inter_chain option");
215 : }
216 27 : std::vector<unsigned> nlist(30);
217 353 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
218 : unsigned iprev=0;
219 3237 : for(unsigned i=0; i<ichain; ++i) {
220 2911 : iprev+=chains[i];
221 : }
222 326 : unsigned inres=chains[ichain]/5;
223 326 : if( chains[ichain]%5!=0 ) {
224 0 : error("backbone segment received does not contain a multiple of five residues");
225 : }
226 326 : if( inres<3 ) {
227 : // A strand contains three residues. Without this the unsigned subtraction in the loop below wraps around
228 1 : continue;
229 : }
230 2283 : for(unsigned ires=0; ires<inres-2; ++ires) {
231 19430 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
232 : unsigned jprev=0;
233 110506 : for(unsigned i=0; i<jchain; ++i) {
234 93034 : jprev+=chains[i];
235 : }
236 17472 : unsigned jnres=chains[jchain]/5;
237 17472 : if( chains[jchain]%5!=0 ) {
238 0 : error("backbone segment received does not contain a multiple of five residues");
239 : }
240 17472 : if( jnres<3 ) {
241 10 : continue;
242 : }
243 122314 : for(unsigned jres=0; jres<jnres-2; ++jres) {
244 1677632 : for(unsigned k=0; k<15; ++k) {
245 1572780 : nlist[k]=iprev + ires*5+k;
246 1572780 : nlist[k+15]=jprev + jres*5+k;
247 : }
248 : std::string nlstr, num;
249 104852 : Tools::convert( nlist[0], nlstr );
250 104852 : Tools::convert(jjkk, num);
251 104852 : jjkk++;
252 104852 : seglist += " SEGMENT" + num + "=" + nlstr;
253 3145560 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
254 3040708 : Tools::convert( nlist[kk], nlstr );
255 3040708 : seglist += "," + nlstr;
256 : }
257 104852 : if( strands_cutoff>0 ) {
258 419408 : scutoff_action += " ATOMS" + num + "=" + all_atoms[nlist[6]] + "," + all_atoms[nlist[21]];
259 : }
260 : }
261 : }
262 : }
263 : }
264 : }
265 :
266 27 : if( seglist.empty() ) {
267 0 : error("it is not possible to build any segment of backbone that could form a sheet with the residues and the STYLE that were given");
268 : }
269 :
270 : // Build the reference structure ( in angstroms )
271 27 : std::vector<Vector> reference(30);
272 27 : reference[0]=Vector( 1.244, -4.620, -2.127); // N i
273 27 : reference[1]=Vector(-0.016, -4.500, -1.395); // CA
274 27 : reference[2]=Vector( 0.105, -5.089, 0.024); // CB
275 27 : reference[3]=Vector(-0.287, -3.000, -1.301); // C
276 27 : reference[4]=Vector( 0.550, -2.245, -0.822); // O
277 27 : reference[5]=Vector(-1.445, -2.551, -1.779); // N i+1
278 27 : reference[6]=Vector(-1.752, -1.130, -1.677); // CA
279 27 : reference[7]=Vector(-2.113, -0.550, -3.059); // CB
280 27 : reference[8]=Vector(-2.906, -0.961, -0.689); // C
281 27 : reference[9]=Vector(-3.867, -1.738, -0.695); // O
282 27 : reference[10]=Vector(-2.774, 0.034, 0.190); // N i+2
283 27 : reference[11]=Vector(-3.788, 0.331, 1.201); // CA
284 27 : reference[12]=Vector(-3.188, 0.300, 2.624); // CB
285 27 : reference[13]=Vector(-4.294, 1.743, 0.937); // C
286 27 : reference[14]=Vector(-3.503, 2.671, 0.821); // O
287 27 : reference[15]=Vector( 4.746, -2.363, 0.188); // N j
288 27 : reference[16]=Vector( 3.427, -1.839, 0.545); // CA
289 27 : reference[17]=Vector( 3.135, -1.958, 2.074); // CB
290 27 : reference[18]=Vector( 3.346, -0.365, 0.181); // C
291 27 : reference[19]=Vector( 4.237, 0.412, 0.521); // O
292 27 : reference[20]=Vector( 2.261, 0.013, -0.487); // N j+1
293 27 : reference[21]=Vector( 2.024, 1.401, -0.875); // CA
294 27 : reference[22]=Vector( 1.489, 1.514, -2.313); // CB
295 27 : reference[23]=Vector( 0.914, 1.902, 0.044); // C
296 27 : reference[24]=Vector(-0.173, 1.330, 0.052); // O
297 27 : reference[25]=Vector( 1.202, 2.940, 0.828); // N j+2
298 27 : reference[26]=Vector( 0.190, 3.507, 1.718); // CA
299 27 : reference[27]=Vector( 0.772, 3.801, 3.104); // CB
300 27 : reference[28]=Vector(-0.229, 4.791, 1.038); // C
301 27 : reference[29]=Vector( 0.523, 5.771, 0.996); // O
302 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
303 : std::string ref0, ref1, ref2;
304 27 : Tools::convert( reference[0][0], ref0 );
305 27 : Tools::convert( reference[0][1], ref1 );
306 27 : Tools::convert( reference[0][2], ref2 );
307 81 : std::string structure=" STRUCTURE1=" + ref0 + "," + ref1 + "," + ref2;
308 810 : for(unsigned i=1; i<30; ++i) {
309 3132 : for(unsigned k=0; k<3; ++k) {
310 2349 : Tools::convert( reference[i][k], ref0 );
311 4698 : structure += "," + ref0;
312 : }
313 : }
314 :
315 27 : reference[0]=Vector(-1.439, -5.122, -1.144); // N i
316 27 : reference[1]=Vector(-0.816, -3.803, -1.013); // CA
317 27 : reference[2]=Vector( 0.099, -3.509, -2.206); // CB
318 27 : reference[3]=Vector(-1.928, -2.770, -0.952); // C
319 27 : reference[4]=Vector(-2.991, -2.970, -1.551); // O
320 27 : reference[5]=Vector(-1.698, -1.687, -0.215); // N i+1
321 27 : reference[6]=Vector(-2.681, -0.613, -0.143); // CA
322 27 : reference[7]=Vector(-3.323, -0.477, 1.267); // CB
323 27 : reference[8]=Vector(-1.984, 0.681, -0.574); // C
324 27 : reference[9]=Vector(-0.807, 0.921, -0.273); // O
325 27 : reference[10]=Vector(-2.716, 1.492, -1.329); // N i+2
326 27 : reference[11]=Vector(-2.196, 2.731, -1.883); // CA
327 27 : reference[12]=Vector(-2.263, 2.692, -3.418); // CB
328 27 : reference[13]=Vector(-2.989, 3.949, -1.433); // C
329 27 : reference[14]=Vector(-4.214, 3.989, -1.583); // O
330 27 : reference[15]=Vector( 2.464, -4.352, 2.149); // N j
331 27 : reference[16]=Vector( 3.078, -3.170, 1.541); // CA
332 27 : reference[17]=Vector( 3.398, -3.415, 0.060); // CB
333 27 : reference[18]=Vector( 2.080, -2.021, 1.639); // C
334 27 : reference[19]=Vector( 0.938, -2.178, 1.225); // O
335 27 : reference[20]=Vector( 2.525, -0.886, 2.183); // N j+1
336 27 : reference[21]=Vector( 1.692, 0.303, 2.346); // CA
337 27 : reference[22]=Vector( 1.541, 0.665, 3.842); // CB
338 27 : reference[23]=Vector( 2.420, 1.410, 1.608); // C
339 27 : reference[24]=Vector( 3.567, 1.733, 1.937); // O
340 27 : reference[25]=Vector( 1.758, 1.976, 0.600); // N j+2
341 27 : reference[26]=Vector( 2.373, 2.987, -0.238); // CA
342 27 : reference[27]=Vector( 2.367, 2.527, -1.720); // CB
343 27 : reference[28]=Vector( 1.684, 4.331, -0.148); // C
344 27 : reference[29]=Vector( 0.486, 4.430, -0.415); // O
345 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
346 27 : Tools::convert( reference[0][0], ref0 );
347 27 : Tools::convert( reference[0][1], ref1 );
348 27 : Tools::convert( reference[0][2], ref2 );
349 81 : structure +=" STRUCTURE2=" + ref0 + "," + ref1 + "," + ref2;
350 810 : for(unsigned i=1; i<30; ++i) {
351 3132 : for(unsigned k=0; k<3; ++k) {
352 2349 : Tools::convert( reference[i][k], ref0 );
353 4698 : structure += "," + ref0;
354 : }
355 : }
356 :
357 27 : std::string nopbcstr="";
358 : bool nopbc;
359 27 : parseFlag("NOPBC",nopbc);
360 27 : if( nopbc ) {
361 : nopbcstr = " NOPBC";
362 : }
363 27 : std::string usegpustr="";
364 : {
365 : bool usegpu;
366 27 : parseFlag("USEGPU",usegpu);
367 27 : if( usegpu ) {
368 : usegpustr = " USEGPU";
369 : }
370 : }
371 : std::string type;
372 27 : parse("TYPE",type);
373 27 : std::string lab = getShortcutLabel() + "_low";
374 27 : if( uselessthan ) {
375 22 : lab = getShortcutLabel();
376 : }
377 27 : std::string atoms="ATOMS=" + all_atoms[0];
378 14285 : for(unsigned i=1; i<all_atoms.size(); ++i) {
379 28516 : atoms += "," + all_atoms[i];
380 : }
381 27 : std::string inputLine = getShortcutLabel() + "_both:";
382 27 : if( type=="DRMSD" ) {
383 : inputLine+=" SECONDARY_STRUCTURE_DRMSD BONDLENGTH=0.17";
384 : } else {
385 28 : inputLine+=" SECONDARY_STRUCTURE_RMSD TYPE=" +type;
386 : }
387 27 : inputLine+=" ALIGN_STRANDS " + seglist + structure
388 27 : + " " + atoms + nopbcstr + usegpustr;
389 :
390 27 : if( strands_cutoff>0 ) {
391 : // The distance used for the cutoff has to treat the periodic boundary conditions as the rest of the calculation
392 46 : readInputLine( scutoff_action + nopbcstr );
393 : std::string str_cut;
394 23 : Tools::convert( strands_cutoff, str_cut );
395 46 : readInputLine( getShortcutLabel() + "_cut: CUSTOM ARG=" + getShortcutLabel()
396 23 : + "_cut_dists FUNC=step(" + str_cut + "-x) PERIODIC=NO");
397 23 : inputLine+=" MASK=" + getShortcutLabel() + "_cut";
398 23 : readInputLine(inputLine);
399 23 : if( ltmap.length()>0 ) {
400 : // Create the lowest line
401 46 : readInputLine( lab + ": LOWEST ARG=" + getShortcutLabel() + "_both.struct-1,"
402 69 : + getShortcutLabel() + "_both.struct-2 MASK=" + getShortcutLabel() + "_cut");
403 : // Create the less than object
404 46 : readInputLine( getShortcutLabel() + "_ltu: LESS_THAN ARG=" + lab
405 69 : + " SWITCH={" + ltmap +"} MASK=" + getShortcutLabel() + "_cut");
406 : // Multiply by the strands cutoff
407 46 : readInputLine( getShortcutLabel()
408 46 : + "_lt: CUSTOM ARG=" + getShortcutLabel() + "_ltu,"
409 69 : + getShortcutLabel() + "_cut"
410 : " FUNC=x*y PERIODIC=NO");
411 : }
412 : } else {
413 4 : readInputLine(inputLine);
414 4 : if( ltmap.length()>0 ) {
415 : // Create the lowest line
416 8 : readInputLine( lab + ": LOWEST ARG=" + getShortcutLabel() + "_both.struct-1,"
417 8 : + getShortcutLabel() + "_both.struct-2" );
418 : // Create the less than object
419 16 : readInputLine( getShortcutLabel() + "_lt: LESS_THAN ARG=" + lab + " SWITCH={" + ltmap +"}");
420 : }
421 : }
422 : // Create the less than object
423 27 : if( ltmap.length()>0 ) {
424 27 : if( uselessthan ) {
425 66 : readInputLine( getShortcutLabel() + "_lessthan: SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
426 : } else {
427 15 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
428 : }
429 : }
430 3252933 : }
431 :
432 : }
433 : }
|
