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 ANTIBETARMSD
30 : /*
31 : Probe the antiparallel beta sheet content of your protein structure.
32 :
33 : Two protein segments containing three contiguous residues can form an antiparallel beta sheet.
34 : Although if the two segments are part of the same protein chain they must be separated by
35 : a minimum of 2 residues to make room for the turn. This colvar thus generates the set of
36 : all possible six residue sections that could conceivably form an antiparallel 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 antiparallel 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 ideas from the paper cited below. The authors of
44 : this paper use the set of distances from the anti parallel beta sheet configurations to measure
45 : the number of segments that have an configuration that resembles an anti 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 antiparallel 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 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 purely configuration composed of pure anti-parallel beta sheets or the distance between the set of
58 : residues that is closest to an anti-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 antiparallel beta sheet configuration.
65 :
66 : ```plumed
67 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
68 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
69 : ab: ANTIBETARMSD RESIDUES=all STRANDS_CUTOFF=1 R_0=0.1
70 : PRINT ARG=ab 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: ANTIBETARMSD 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 : ab: ANTIBETARMSD RESIDUES=all STRANDS_CUTOFF=1 R_0=0.1 NOPBC
96 : PRINT ARG=ab 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 AntibetaRMSD : public ActionShortcut {
111 : public:
112 : static void registerKeywords( Keywords& keys );
113 : explicit AntibetaRMSD(const ActionOptions&);
114 : };
115 :
116 : PLUMED_REGISTER_ACTION(AntibetaRMSD,"ANTIBETARMSD")
117 :
118 128 : void AntibetaRMSD::registerKeywords( Keywords& keys ) {
119 128 : SecondaryStructureInput::registerKeywords( keys );
120 128 : keys.remove("ATOMS");
121 128 : keys.remove("SEGMENT");
122 128 : keys.remove("STRUCTURE");
123 128 : keys.remove("MASK");
124 128 : keys.remove("ALIGN_STRANDS");
125 256 : keys.setValueDescription("scalar/vector","if LESS_THAN is present the RMSD distance between each residue and the ideal antiparallel beta sheet. If LESS_THAN is not present the number of residue segments where the structure is similar to an anti parallel beta sheet");
126 128 : keys.add("compulsory","STYLE","all","Antiparallel 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 128 : 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 128 : keys.needsAction("DISTANCE");
135 128 : keys.needsAction("CUSTOM");
136 128 : }
137 :
138 24 : AntibetaRMSD::AntibetaRMSD(const ActionOptions&ao):
139 : Action(ao),
140 24 : ActionShortcut(ao) {
141 : // Read in the input and create a string that describes how to compute the less than
142 : std::string ltmap;
143 24 : bool uselessthan=SecondaryStructureInput::readShortcutWords( ltmap, this );
144 : // read in the backbone atoms
145 24 : std::vector<unsigned> chains;
146 24 : std::vector<std::string> all_atoms;
147 48 : SecondaryStructureInput::readBackboneAtoms( this, plumed, "protein", chains, all_atoms );
148 :
149 : bool intra_chain(false), inter_chain(false);
150 : std::string style;
151 48 : parse("STYLE",style);
152 24 : if( Tools::caseInSensStringCompare(style, "all") ) {
153 : intra_chain=true;
154 : inter_chain=true;
155 4 : } else if( Tools::caseInSensStringCompare(style, "inter") ) {
156 : intra_chain=false;
157 : inter_chain=true;
158 1 : } else if( Tools::caseInSensStringCompare(style, "intra") ) {
159 : intra_chain=true;
160 : inter_chain=false;
161 : } else {
162 0 : error( style + " is not a valid directive for the STYLE keyword");
163 : }
164 :
165 24 : double strands_cutoff=-1.0;
166 48 : parse("STRANDS_CUTOFF",strands_cutoff);
167 : std::string scutoff_action;
168 24 : if( strands_cutoff>0 ) {
169 44 : scutoff_action=getShortcutLabel() + "_cut_dists: DISTANCE ";
170 : }
171 :
172 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
173 : std::string seglist;
174 24 : unsigned k=1;
175 24 : if( intra_chain ) {
176 : unsigned nprevious=0;
177 21 : std::vector<unsigned> nlist(30);
178 349 : for(unsigned i=0; i<chains.size(); ++i) {
179 328 : if( chains[i]<40 ) {
180 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");
181 : }
182 : // Loop over all possible triples in each 8 residue segment of protein
183 328 : unsigned nres=chains[i]/5;
184 328 : if( chains[i]%5!=0 ) {
185 0 : error("backbone segment received does not contain a multiple of five residues");
186 : }
187 672 : for(unsigned ires=0; ires<nres-7; ires++) {
188 728 : for(unsigned jres=ires+7; jres<nres; jres++) {
189 6144 : for(unsigned kk=0; kk<15; ++kk) {
190 5760 : nlist[kk]=nprevious + ires*5+kk;
191 5760 : nlist[kk+15]=nprevious + (jres-2)*5+kk;
192 : }
193 : std::string nlstr, num;
194 384 : Tools::convert( nlist[0], nlstr );
195 384 : Tools::convert(k, num);
196 384 : ++k;
197 384 : seglist += " SEGMENT" + num + "=" + nlstr;
198 11520 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
199 11136 : Tools::convert( nlist[kk], nlstr );
200 11136 : seglist += "," + nlstr;
201 : }
202 384 : if( strands_cutoff>0 ) {
203 1536 : scutoff_action += " ATOMS" + num + "=" + all_atoms[nlist[6]] + "," + all_atoms[nlist[21]];
204 : }
205 : }
206 : }
207 328 : nprevious+=chains[i];
208 : }
209 : }
210 24 : if( inter_chain ) {
211 23 : if( chains.size()==1 && !Tools::caseInSensStringCompare(style, "all") ) {
212 0 : error("there is only one chain defined so cannot use inter_chain option");
213 : }
214 23 : std::vector<unsigned> nlist(30);
215 317 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
216 : unsigned iprev=0;
217 2901 : for(unsigned i=0; i<ichain; ++i) {
218 2607 : iprev+=chains[i];
219 : }
220 294 : unsigned inres=chains[ichain]/5;
221 294 : if( chains[ichain]%5!=0 ) {
222 0 : error("backbone segment received does not contain a multiple of five residues");
223 : }
224 294 : if( inres<3 ) {
225 : // A strand contains three residues. Without this the unsigned subtraction in the loop below wraps around
226 1 : continue;
227 : }
228 2049 : for(unsigned ires=0; ires<inres-2; ++ires) {
229 17394 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
230 : unsigned jprev=0;
231 98880 : for(unsigned i=0; i<jchain; ++i) {
232 83242 : jprev+=chains[i];
233 : }
234 15638 : unsigned jnres=chains[jchain]/5;
235 15638 : if( chains[jchain]%5!=0 ) {
236 0 : error("backbone segment received does not contain a multiple of five residues");
237 : }
238 15638 : if( jnres<3 ) {
239 10 : continue;
240 : }
241 109466 : for(unsigned jres=0; jres<jnres-2; ++jres) {
242 1501408 : for(unsigned kk=0; kk<15; ++kk) {
243 1407570 : nlist[kk]=iprev + ires*5+kk;
244 1407570 : nlist[kk+15]=jprev + jres*5+kk;
245 : }
246 : std::string nlstr, num;
247 93838 : Tools::convert( nlist[0], nlstr );
248 93838 : Tools::convert(k, num);
249 93838 : k++;
250 93838 : seglist += " SEGMENT" + num + "=" + nlstr;
251 2815140 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
252 2721302 : Tools::convert( nlist[kk], nlstr );
253 2721302 : seglist += "," + nlstr;
254 : }
255 93838 : if( strands_cutoff>0 ) {
256 375344 : scutoff_action += " ATOMS" + num + "=" + all_atoms[nlist[6]] + "," + all_atoms[nlist[21]];
257 : }
258 : }
259 : }
260 : }
261 : }
262 : }
263 :
264 24 : if( seglist.empty() ) {
265 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");
266 : }
267 :
268 : // Build the reference structure ( in angstroms )
269 24 : std::vector<Vector> reference(30);
270 24 : reference[0]=Vector( 2.263, -3.795, 1.722); // N i
271 24 : reference[1]=Vector( 2.493, -2.426, 2.263); // CA
272 24 : reference[2]=Vector( 3.847, -1.838, 1.761); // CB
273 24 : reference[3]=Vector( 1.301, -1.517, 1.921); // C
274 24 : reference[4]=Vector( 0.852, -1.504, 0.739); // O
275 24 : reference[5]=Vector( 0.818, -0.738, 2.917); // N i+1
276 24 : reference[6]=Vector(-0.299, 0.243, 2.748); // CA
277 24 : reference[7]=Vector(-1.421, -0.076, 3.757); // CB
278 24 : reference[8]=Vector( 0.273, 1.680, 2.854); // C
279 24 : reference[9]=Vector( 0.902, 1.993, 3.888); // O
280 24 : reference[10]=Vector( 0.119, 2.532, 1.813); // N i+2
281 24 : reference[11]=Vector( 0.683, 3.916, 1.680); // CA
282 24 : reference[12]=Vector( 1.580, 3.940, 0.395); // CB
283 24 : reference[13]=Vector(-0.394, 5.011, 1.630); // C
284 24 : reference[14]=Vector(-1.459, 4.814, 0.982); // O
285 24 : reference[15]=Vector(-2.962, 3.559, -1.359); // N j-2
286 24 : reference[16]=Vector(-2.439, 2.526, -2.287); // CA
287 24 : reference[17]=Vector(-1.189, 3.006, -3.087); // CB
288 24 : reference[18]=Vector(-2.081, 1.231, -1.520); // C
289 24 : reference[19]=Vector(-1.524, 1.324, -0.409); // O
290 24 : reference[20]=Vector(-2.326, 0.037, -2.095); // N j-1
291 24 : reference[21]=Vector(-1.858, -1.269, -1.554); // CA
292 24 : reference[22]=Vector(-3.053, -2.199, -1.291); // CB
293 24 : reference[23]=Vector(-0.869, -1.949, -2.512); // C
294 24 : reference[24]=Vector(-1.255, -2.070, -3.710); // O
295 24 : reference[25]=Vector( 0.326, -2.363, -2.072); // N j
296 24 : reference[26]=Vector( 1.405, -2.992, -2.872); // CA
297 24 : reference[27]=Vector( 2.699, -2.129, -2.917); // CB
298 24 : reference[28]=Vector( 1.745, -4.399, -2.330); // C
299 24 : reference[29]=Vector( 1.899, -4.545, -1.102); // O
300 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
301 : std::string ref0, ref1, ref2;
302 24 : Tools::convert( reference[0][0], ref0 );
303 24 : Tools::convert( reference[0][1], ref1 );
304 24 : Tools::convert( reference[0][2], ref2 );
305 72 : std::string structure=" STRUCTURE1=" + ref0 + "," + ref1 + "," + ref2;
306 720 : for(unsigned i=1; i<30; ++i) {
307 2784 : for(unsigned kk=0; kk<3; ++kk) {
308 2088 : Tools::convert( reference[i][kk], ref0 );
309 4176 : structure += "," + ref0;
310 : }
311 : }
312 :
313 24 : std::string nopbcstr="";
314 : bool nopbc;
315 24 : parseFlag("NOPBC",nopbc);
316 24 : if( nopbc ) {
317 : nopbcstr = " NOPBC";
318 : }
319 24 : std::string usegpustr="";
320 : {
321 : bool usegpu;
322 24 : parseFlag("USEGPU",usegpu);
323 24 : if( usegpu ) {
324 : usegpustr = " USEGPU";
325 : }
326 : }
327 : std::string type;
328 24 : parse("TYPE",type);
329 24 : std::string lab = getShortcutLabel() + "_rmsd";
330 24 : if( uselessthan ) {
331 22 : lab = getShortcutLabel();
332 : }
333 24 : std::string atoms="ATOMS=" + all_atoms[0];
334 13385 : for(unsigned i=1; i<all_atoms.size(); ++i) {
335 26722 : atoms += "," + all_atoms[i];
336 : }
337 24 : std::string inputLine = lab+":";
338 24 : if( type=="DRMSD" ) {
339 : inputLine+=" SECONDARY_STRUCTURE_DRMSD BONDLENGTH=0.17";
340 : } else {
341 28 : inputLine+=" SECONDARY_STRUCTURE_RMSD TYPE=" +type;
342 : }
343 24 : inputLine+=" ALIGN_STRANDS " + seglist + structure
344 24 : + " " + atoms + nopbcstr + usegpustr;
345 :
346 24 : if( strands_cutoff>0 ) {
347 : // The distance used for the cutoff has to treat the periodic boundary conditions as the rest of the calculation
348 44 : readInputLine( scutoff_action + nopbcstr );
349 : std::string str_cut;
350 22 : Tools::convert( strands_cutoff, str_cut );
351 44 : readInputLine( getShortcutLabel() + "_cut: CUSTOM ARG=" + getShortcutLabel()
352 22 : + "_cut_dists FUNC=step(" + str_cut + "-x) PERIODIC=NO");
353 22 : inputLine+=" MASK=" + getShortcutLabel() + "_cut";
354 22 : readInputLine(inputLine);
355 22 : if( ltmap.length()>0 ) {
356 : // Create the less than object
357 44 : readInputLine( getShortcutLabel() + "_ltu: LESS_THAN ARG=" + lab
358 66 : + " SWITCH={" + ltmap +"} MASK=" + getShortcutLabel() + "_cut");
359 : // Multiply by the strands cutoff
360 44 : readInputLine( getShortcutLabel()
361 44 : + "_lt: CUSTOM ARG=" + getShortcutLabel() + "_ltu,"
362 66 : + getShortcutLabel() + "_cut"
363 : " FUNC=x*y PERIODIC=NO");
364 : }
365 : } else {
366 2 : readInputLine(inputLine);
367 2 : if( ltmap.length()>0 ) {
368 : // Create the less than object
369 8 : readInputLine( getShortcutLabel() + "_lt: LESS_THAN ARG=" + lab + " SWITCH={" + ltmap +"}");
370 : }
371 : }
372 : // Create the less than object
373 24 : if( ltmap.length()>0 ) {
374 24 : if( uselessthan ) {
375 66 : readInputLine( getShortcutLabel() + "_lessthan: SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
376 : } else {
377 6 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
378 : }
379 : }
380 2920952 : }
381 :
382 : }
383 : }
|
