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 "SecondaryStructureRMSD.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 RMSD 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 following reference \cite pietrucci09jctc. 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 : \f[
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 : \f]
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 \cite pietrucci09jctc. The R_0
54 : parameter must be set by the user - the value used in \cite pietrucci09jctc 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 : Please be aware that for codes like gromacs you must ensure that plumed
64 : reconstructs the chains involved in your CV when you calculate this CV using
65 : anything other than TYPE=DRMSD. For more details as to how to do this see \ref WHOLEMOLECULES.
66 :
67 : \par Examples
68 :
69 : The following input calculates the number of six residue segments of
70 : protein that are in an parallel beta sheet configuration.
71 :
72 : \plumedfile
73 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
74 : MOLINFO STRUCTURE=beta.pdb
75 : pb: PARABETARMSD RESIDUES=all STRANDS_CUTOFF=1
76 : \endplumedfile
77 :
78 : Here the same is done use RMSD instead of DRMSD
79 :
80 : \plumedfile
81 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
82 : MOLINFO STRUCTURE=helix.pdb
83 : WHOLEMOLECULES ENTITY0=1-100
84 : hh: PARABETARMSD RESIDUES=all TYPE=OPTIMAL R_0=0.1 STRANDS_CUTOFF=1
85 : \endplumedfile
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 : class ParabetaRMSD : public ActionShortcut {
91 : public:
92 : static void registerKeywords( Keywords& keys );
93 : explicit ParabetaRMSD(const ActionOptions&);
94 : };
95 :
96 : PLUMED_REGISTER_ACTION(ParabetaRMSD,"PARABETARMSD")
97 :
98 100 : void ParabetaRMSD::registerKeywords( Keywords& keys ) {
99 100 : SecondaryStructureRMSD::registerKeywords( keys );
100 100 : keys.remove("ATOMS");
101 100 : keys.remove("SEGMENT");
102 100 : keys.remove("BONDLENGTH");
103 100 : keys.remove("NO_ACTION_LOG");
104 100 : keys.remove("CUTOFF_ATOMS");
105 100 : keys.remove("STRUCTURE");
106 200 : 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 "
107 : "chain configuration with the appropriate geometry are counted. If STYLE=inter "
108 : "only sheet-like configurations involving two chains are counted, while if STYLE=intra "
109 : "only sheet-like configurations involving a single chain are counted");
110 100 : keys.needsAction("LOWEST");
111 100 : }
112 :
113 18 : ParabetaRMSD::ParabetaRMSD(const ActionOptions&ao):
114 : Action(ao),
115 18 : ActionShortcut(ao) {
116 : // Read in the input and create a string that describes how to compute the less than
117 : std::string ltmap;
118 18 : bool uselessthan=SecondaryStructureRMSD::readShortcutWords( ltmap, this );
119 : // read in the backbone atoms
120 : std::vector<unsigned> chains;
121 : std::string atoms;
122 36 : SecondaryStructureRMSD::readBackboneAtoms( this, plumed, "protein", chains, atoms );
123 :
124 : bool intra_chain(false), inter_chain(false);
125 : std::string seglist;
126 : std::string style;
127 18 : parse("STYLE",style);
128 18 : unsigned jjkk=1;
129 36 : if( Tools::caseInSensStringCompare(style, "all") ) {
130 : intra_chain=true;
131 : inter_chain=true;
132 2 : } else if( Tools::caseInSensStringCompare(style, "inter") ) {
133 : intra_chain=false;
134 : inter_chain=true;
135 0 : } else if( Tools::caseInSensStringCompare(style, "intra") ) {
136 : intra_chain=true;
137 : inter_chain=false;
138 : } else {
139 0 : error( style + " is not a valid directive for the STYLE keyword");
140 : }
141 :
142 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
143 : if( intra_chain ) {
144 : unsigned nprevious=0;
145 17 : std::vector<unsigned> nlist(30);
146 272 : for(unsigned i=0; i<chains.size(); ++i) {
147 255 : if( chains[i]<40 ) {
148 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");
149 : }
150 : // Loop over all possible triples in each 8 residue segment of protein
151 255 : unsigned nres=chains[i]/5;
152 255 : if( chains[i]%5!=0 ) {
153 0 : error("backbone segment received does not contain a multiple of five residues");
154 : }
155 267 : for(unsigned ires=0; ires<nres-8; ires++) {
156 42 : for(unsigned jres=ires+6; jres<nres-2; jres++) {
157 480 : for(unsigned k=0; k<15; ++k) {
158 450 : nlist[k]=nprevious + ires*5+k;
159 450 : nlist[k+15]=nprevious + jres*5+k;
160 : }
161 : std::string nlstr, num;
162 30 : Tools::convert( nlist[0], nlstr );
163 30 : Tools::convert(jjkk, num);
164 30 : jjkk++;
165 60 : seglist += " SEGMENT" + num + "=" + nlstr;
166 900 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
167 870 : Tools::convert( nlist[kk], nlstr );
168 1740 : seglist += "," + nlstr;
169 : }
170 : }
171 : }
172 255 : nprevious+=chains[i];
173 : }
174 : }
175 : // This constructs all conceivable sections of antibeta sheet that form between chains
176 18 : if( inter_chain ) {
177 24 : if( chains.size()==1 && !Tools::caseInSensStringCompare(style, "all") ) {
178 0 : error("there is only one chain defined so cannot use inter_chain option");
179 : }
180 18 : std::vector<unsigned> nlist(30);
181 273 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
182 : unsigned iprev=0;
183 2550 : for(unsigned i=0; i<ichain; ++i) {
184 2295 : iprev+=chains[i];
185 : }
186 255 : unsigned inres=chains[ichain]/5;
187 255 : if( chains[ichain]%5!=0 ) {
188 0 : error("backbone segment received does not contain a multiple of five residues");
189 : }
190 1785 : for(unsigned ires=0; ires<inres-2; ++ires) {
191 15300 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
192 : unsigned jprev=0;
193 87210 : for(unsigned i=0; i<jchain; ++i) {
194 73440 : jprev+=chains[i];
195 : }
196 13770 : unsigned jnres=chains[jchain]/5;
197 13770 : if( chains[jchain]%5!=0 ) {
198 0 : error("backbone segment received does not contain a multiple of five residues");
199 : }
200 96390 : for(unsigned jres=0; jres<jnres-2; ++jres) {
201 1321920 : for(unsigned k=0; k<15; ++k) {
202 1239300 : nlist[k]=iprev + ires*5+k;
203 1239300 : nlist[k+15]=jprev + jres*5+k;
204 : }
205 : std::string nlstr, num;
206 82620 : Tools::convert( nlist[0], nlstr );
207 82620 : Tools::convert(jjkk, num);
208 82620 : jjkk++;
209 165240 : seglist += " SEGMENT" + num + "=" + nlstr;
210 2478600 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
211 2395980 : Tools::convert( nlist[kk], nlstr );
212 4791960 : seglist += "," + nlstr;
213 : }
214 : }
215 : }
216 : }
217 : }
218 : }
219 :
220 : // Build the reference structure ( in angstroms )
221 18 : std::vector<Vector> reference(30);
222 18 : reference[0]=Vector( 1.244, -4.620, -2.127); // N i
223 18 : reference[1]=Vector(-0.016, -4.500, -1.395); // CA
224 18 : reference[2]=Vector( 0.105, -5.089, 0.024); // CB
225 18 : reference[3]=Vector(-0.287, -3.000, -1.301); // C
226 18 : reference[4]=Vector( 0.550, -2.245, -0.822); // O
227 18 : reference[5]=Vector(-1.445, -2.551, -1.779); // N i+1
228 18 : reference[6]=Vector(-1.752, -1.130, -1.677); // CA
229 18 : reference[7]=Vector(-2.113, -0.550, -3.059); // CB
230 18 : reference[8]=Vector(-2.906, -0.961, -0.689); // C
231 18 : reference[9]=Vector(-3.867, -1.738, -0.695); // O
232 18 : reference[10]=Vector(-2.774, 0.034, 0.190); // N i+2
233 18 : reference[11]=Vector(-3.788, 0.331, 1.201); // CA
234 18 : reference[12]=Vector(-3.188, 0.300, 2.624); // CB
235 18 : reference[13]=Vector(-4.294, 1.743, 0.937); // C
236 18 : reference[14]=Vector(-3.503, 2.671, 0.821); // O
237 18 : reference[15]=Vector( 4.746, -2.363, 0.188); // N j
238 18 : reference[16]=Vector( 3.427, -1.839, 0.545); // CA
239 18 : reference[17]=Vector( 3.135, -1.958, 2.074); // CB
240 18 : reference[18]=Vector( 3.346, -0.365, 0.181); // C
241 18 : reference[19]=Vector( 4.237, 0.412, 0.521); // O
242 18 : reference[20]=Vector( 2.261, 0.013, -0.487); // N j+1
243 18 : reference[21]=Vector( 2.024, 1.401, -0.875); // CA
244 18 : reference[22]=Vector( 1.489, 1.514, -2.313); // CB
245 18 : reference[23]=Vector( 0.914, 1.902, 0.044); // C
246 18 : reference[24]=Vector(-0.173, 1.330, 0.052); // O
247 18 : reference[25]=Vector( 1.202, 2.940, 0.828); // N j+2
248 18 : reference[26]=Vector( 0.190, 3.507, 1.718); // CA
249 18 : reference[27]=Vector( 0.772, 3.801, 3.104); // CB
250 18 : reference[28]=Vector(-0.229, 4.791, 1.038); // C
251 18 : reference[29]=Vector( 0.523, 5.771, 0.996); // O
252 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
253 : std::string ref0, ref1, ref2;
254 18 : Tools::convert( reference[0][0], ref0 );
255 18 : Tools::convert( reference[0][1], ref1 );
256 18 : Tools::convert( reference[0][2], ref2 );
257 36 : std::string structure=" STRUCTURE1=" + ref0 + "," + ref1 + "," + ref2;
258 540 : for(unsigned i=1; i<30; ++i) {
259 2088 : for(unsigned k=0; k<3; ++k) {
260 1566 : Tools::convert( reference[i][k], ref0 );
261 3132 : structure += "," + ref0;
262 : }
263 : }
264 :
265 18 : reference[0]=Vector(-1.439, -5.122, -1.144); // N i
266 18 : reference[1]=Vector(-0.816, -3.803, -1.013); // CA
267 18 : reference[2]=Vector( 0.099, -3.509, -2.206); // CB
268 18 : reference[3]=Vector(-1.928, -2.770, -0.952); // C
269 18 : reference[4]=Vector(-2.991, -2.970, -1.551); // O
270 18 : reference[5]=Vector(-1.698, -1.687, -0.215); // N i+1
271 18 : reference[6]=Vector(-2.681, -0.613, -0.143); // CA
272 18 : reference[7]=Vector(-3.323, -0.477, 1.267); // CB
273 18 : reference[8]=Vector(-1.984, 0.681, -0.574); // C
274 18 : reference[9]=Vector(-0.807, 0.921, -0.273); // O
275 18 : reference[10]=Vector(-2.716, 1.492, -1.329); // N i+2
276 18 : reference[11]=Vector(-2.196, 2.731, -1.883); // CA
277 18 : reference[12]=Vector(-2.263, 2.692, -3.418); // CB
278 18 : reference[13]=Vector(-2.989, 3.949, -1.433); // C
279 18 : reference[14]=Vector(-4.214, 3.989, -1.583); // O
280 18 : reference[15]=Vector( 2.464, -4.352, 2.149); // N j
281 18 : reference[16]=Vector( 3.078, -3.170, 1.541); // CA
282 18 : reference[17]=Vector( 3.398, -3.415, 0.060); // CB
283 18 : reference[18]=Vector( 2.080, -2.021, 1.639); // C
284 18 : reference[19]=Vector( 0.938, -2.178, 1.225); // O
285 18 : reference[20]=Vector( 2.525, -0.886, 2.183); // N j+1
286 18 : reference[21]=Vector( 1.692, 0.303, 2.346); // CA
287 18 : reference[22]=Vector( 1.541, 0.665, 3.842); // CB
288 18 : reference[23]=Vector( 2.420, 1.410, 1.608); // C
289 18 : reference[24]=Vector( 3.567, 1.733, 1.937); // O
290 18 : reference[25]=Vector( 1.758, 1.976, 0.600); // N j+2
291 18 : reference[26]=Vector( 2.373, 2.987, -0.238); // CA
292 18 : reference[27]=Vector( 2.367, 2.527, -1.720); // CB
293 18 : reference[28]=Vector( 1.684, 4.331, -0.148); // C
294 18 : reference[29]=Vector( 0.486, 4.430, -0.415); // O
295 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
296 18 : Tools::convert( reference[0][0], ref0 );
297 18 : Tools::convert( reference[0][1], ref1 );
298 18 : Tools::convert( reference[0][2], ref2 );
299 36 : structure +=" STRUCTURE2=" + ref0 + "," + ref1 + "," + ref2;
300 540 : for(unsigned i=1; i<30; ++i) {
301 2088 : for(unsigned k=0; k<3; ++k) {
302 1566 : Tools::convert( reference[i][k], ref0 );
303 3132 : structure += "," + ref0;
304 : }
305 : }
306 :
307 : std::string strands_cutoff;
308 18 : parse("STRANDS_CUTOFF",strands_cutoff);
309 18 : std::string nopbcstr="";
310 : bool nopbc;
311 18 : parseFlag("NOPBC",nopbc);
312 18 : if( nopbc ) {
313 : nopbcstr = " NOPBC";
314 : }
315 18 : if( strands_cutoff.length()>0 ) {
316 32 : strands_cutoff=" CUTOFF_ATOMS=6,21 STRANDS_CUTOFF="+strands_cutoff;
317 : }
318 : std::string type;
319 18 : parse("TYPE",type);
320 18 : std::string lab = getShortcutLabel() + "_low";
321 18 : if( uselessthan ) {
322 17 : lab = getShortcutLabel();
323 : }
324 18 : if( seglist.length()==0 ) {
325 0 : error("no segments to investigate");
326 : }
327 36 : readInputLine( getShortcutLabel() + "_both: SECONDARY_STRUCTURE_RMSD BONDLENGTH=0.17" + seglist + structure + " " + atoms + " TYPE=" + type + strands_cutoff + nopbcstr );
328 18 : if( ltmap.length()>0 ) {
329 : // Create the lowest line
330 36 : readInputLine( lab + ": LOWEST ARG=" + getShortcutLabel() + "_both.struct-1," + getShortcutLabel() + "_both.struct-2" );
331 : // Create the less than object
332 18 : SecondaryStructureRMSD::expandShortcut( uselessthan, getShortcutLabel(), lab, ltmap, this );
333 : }
334 18 : }
335 :
336 : }
337 : }
|
