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 "SecondaryStructureRMSD.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/Atoms.h"
26 :
27 : namespace PLMD {
28 : namespace secondarystructure {
29 :
30 : //+PLUMEDOC COLVAR PARABETARMSD
31 : /*
32 : Probe the parallel beta sheet content of your protein structure.
33 :
34 : Two protein segments containing three continguous residues can form a parallel beta sheet.
35 : Although if the two segments are part of the same protein chain they must be separated by
36 : a minimum of 3 residues to make room for the turn. This colvar thus generates the set of
37 : all possible six residue sections that could conceivably form a parallel beta sheet
38 : and calculates the RMSD distance between the configuration in which the residues find themselves
39 : and an idealized parallel beta sheet structure. These distances can be calculated by either
40 : aligning the instantaneous structure with the reference structure and measuring each
41 : atomic displacement or by calculating differences between the set of interatomic
42 : distances in the reference and instantaneous structures.
43 :
44 : This colvar is based on the following reference \cite pietrucci09jctc. The authors of
45 : this paper use the set of distances from the parallel beta sheet configurations to measure
46 : the number of segments whose configuration resembles a parallel beta sheet. This is done by calculating
47 : the following sum of functions of the rmsd distances:
48 :
49 : \f[
50 : 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 }
51 : \f]
52 :
53 : where the sum runs over all possible segments of parallel beta sheet. By default the
54 : NN, MM and D_0 parameters are set equal to those used in \cite pietrucci09jctc. The R_0
55 : parameter must be set by the user - the value used in \cite pietrucci09jctc was 0.08 nm.
56 :
57 : If you change the function in the above sum you can calculate quantities such as the average
58 : distance from a structure composed of only parallel beta sheets or the distance between the set of
59 : residues that is closest to a parallel beta sheet and the reference configuration. To do these sorts of
60 : calculations you can use the AVERAGE and MIN keywords. In addition you can use the LESS_THAN
61 : keyword if you would like to change the form of the switching function. If you use any of these
62 : options you no longer need to specify NN, R_0, MM and D_0.
63 :
64 : Please be aware that for codes like gromacs you must ensure that plumed
65 : reconstructs the chains involved in your CV when you calculate this CV using
66 : anthing other than TYPE=DRMSD. For more details as to how to do this see \ref WHOLEMOLECULES.
67 :
68 : \par Examples
69 :
70 : The following input calculates the number of six residue segments of
71 : protein that are in an parallel beta sheet configuration.
72 :
73 : \verbatim
74 : MOLINFO STRUCTURE=helix.pdb
75 : PARABETARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
76 : \endverbatim
77 : (see also \ref MOLINFO)
78 :
79 : */
80 : //+ENDPLUMEDOC
81 :
82 20 : class ParabetaRMSD : public SecondaryStructureRMSD {
83 : public:
84 : static void registerKeywords( Keywords& keys );
85 : explicit ParabetaRMSD(const ActionOptions&);
86 : };
87 :
88 2533 : PLUMED_REGISTER_ACTION(ParabetaRMSD,"PARABETARMSD")
89 :
90 11 : void ParabetaRMSD::registerKeywords( Keywords& keys ) {
91 11 : SecondaryStructureRMSD::registerKeywords( keys );
92 : 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 "
93 : "chain configuration with the appropriate geometry are counted. If STYLE=inter "
94 : "only sheet-like configurations involving two chains are counted, while if STYLE=intra "
95 11 : "only sheet-like configurations involving a single chain are counted");
96 11 : keys.use("STRANDS_CUTOFF");
97 11 : }
98 :
99 10 : ParabetaRMSD::ParabetaRMSD(const ActionOptions&ao):
100 : Action(ao),
101 10 : SecondaryStructureRMSD(ao)
102 : {
103 : // read in the backbone atoms
104 10 : std::vector<unsigned> chains; readBackboneAtoms( "protein", chains );
105 :
106 10 : bool intra_chain(false), inter_chain(false);
107 20 : std::string style; parse("STYLE",style);
108 10 : if( style=="all" ) {
109 10 : intra_chain=true; inter_chain=true;
110 0 : } else if( style=="inter") {
111 0 : intra_chain=false; inter_chain=true;
112 0 : } else if( style=="intra") {
113 0 : intra_chain=true; inter_chain=false;
114 : } else {
115 0 : error( style + " is not a valid directive for the STYLE keyword");
116 : }
117 :
118 : // Align the atoms based on the positions of these two atoms
119 10 : setAtomsFromStrands( 6, 21 );
120 :
121 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
122 10 : if( intra_chain ) {
123 10 : unsigned nres, nprevious=0; std::vector<unsigned> nlist(30);
124 173 : for(unsigned i=0; i<chains.size(); ++i) {
125 163 : if( chains[i]<40 ) error("segment of backbone is not long enough to form an antiparallel beta hairpin. Each backbone fragment must contain a minimum of 8 residues");
126 : // Loop over all possible triples in each 8 residue segment of protein
127 163 : nres=chains[i]/5;
128 163 : if( chains[i]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
129 167 : for(unsigned ires=0; ires<nres-8; ires++) {
130 14 : for(unsigned jres=ires+6; jres<nres-2; jres++) {
131 160 : for(unsigned k=0; k<15; ++k) {
132 150 : nlist[k]=nprevious + ires*5+k;
133 150 : nlist[k+15]=nprevious + jres*5+k;
134 : }
135 10 : addColvar( nlist );
136 : }
137 : }
138 163 : nprevious+=chains[i];
139 10 : }
140 : }
141 : // This constructs all conceivable sections of antibeta sheet that form between chains
142 10 : if( inter_chain ) {
143 10 : if( chains.size()==1 && style!="all" ) error("there is only one chain defined so cannot use inter_chain option");
144 10 : unsigned iprev,jprev,inres,jnres; std::vector<unsigned> nlist(30);
145 163 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
146 153 : iprev=0; for(unsigned i=0; i<ichain; ++i) iprev+=chains[i];
147 153 : inres=chains[ichain]/5;
148 153 : if( chains[ichain]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
149 1071 : for(unsigned ires=0; ires<inres-2; ++ires) {
150 9180 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
151 8262 : jprev=0; for(unsigned i=0; i<jchain; ++i) jprev+=chains[i];
152 8262 : jnres=chains[jchain]/5;
153 8262 : if( chains[jchain]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
154 57834 : for(unsigned jres=0; jres<jnres-2; ++jres) {
155 793152 : for(unsigned k=0; k<15; ++k) {
156 743580 : nlist[k]=iprev + ires*5+k;
157 743580 : nlist[k+15]=jprev + jres*5+k;
158 : }
159 49572 : addColvar( nlist );
160 : }
161 : }
162 : }
163 10 : }
164 : }
165 :
166 : // Build the reference structure ( in angstroms )
167 20 : std::vector<Vector> reference(30);
168 10 : reference[0]=Vector( 1.244, -4.620, -2.127); // N i
169 10 : reference[1]=Vector(-0.016, -4.500, -1.395); // CA
170 10 : reference[2]=Vector( 0.105, -5.089, 0.024); // CB
171 10 : reference[3]=Vector(-0.287, -3.000, -1.301); // C
172 10 : reference[4]=Vector( 0.550, -2.245, -0.822); // O
173 10 : reference[5]=Vector(-1.445, -2.551, -1.779); // N i+1
174 10 : reference[6]=Vector(-1.752, -1.130, -1.677); // CA
175 10 : reference[7]=Vector(-2.113, -0.550, -3.059); // CB
176 10 : reference[8]=Vector(-2.906, -0.961, -0.689); // C
177 10 : reference[9]=Vector(-3.867, -1.738, -0.695); // O
178 10 : reference[10]=Vector(-2.774, 0.034, 0.190); // N i+2
179 10 : reference[11]=Vector(-3.788, 0.331, 1.201); // CA
180 10 : reference[12]=Vector(-3.188, 0.300, 2.624); // CB
181 10 : reference[13]=Vector(-4.294, 1.743, 0.937); // C
182 10 : reference[14]=Vector(-3.503, 2.671, 0.821); // O
183 10 : reference[15]=Vector( 4.746, -2.363, 0.188); // N j
184 10 : reference[16]=Vector( 3.427, -1.839, 0.545); // CA
185 10 : reference[17]=Vector( 3.135, -1.958, 2.074); // CB
186 10 : reference[18]=Vector( 3.346, -0.365, 0.181); // C
187 10 : reference[19]=Vector( 4.237, 0.412, 0.521); // O
188 10 : reference[20]=Vector( 2.261, 0.013, -0.487); // N j+1
189 10 : reference[21]=Vector( 2.024, 1.401, -0.875); // CA
190 10 : reference[22]=Vector( 1.489, 1.514, -2.313); // CB
191 10 : reference[23]=Vector( 0.914, 1.902, 0.044); // C
192 10 : reference[24]=Vector(-0.173, 1.330, 0.052); // O
193 10 : reference[25]=Vector( 1.202, 2.940, 0.828); // N j+2
194 10 : reference[26]=Vector( 0.190, 3.507, 1.718); // CA
195 10 : reference[27]=Vector( 0.772, 3.801, 3.104); // CB
196 10 : reference[28]=Vector(-0.229, 4.791, 1.038); // C
197 10 : reference[29]=Vector( 0.523, 5.771, 0.996); // O
198 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
199 10 : setSecondaryStructure( reference, 0.17/atoms.getUnits().getLength(), 0.1/atoms.getUnits().getLength() );
200 :
201 10 : reference[0]=Vector(-1.439, -5.122, -1.144); // N i
202 10 : reference[1]=Vector(-0.816, -3.803, -1.013); // CA
203 10 : reference[2]=Vector( 0.099, -3.509, -2.206); // CB
204 10 : reference[3]=Vector(-1.928, -2.770, -0.952); // C
205 10 : reference[4]=Vector(-2.991, -2.970, -1.551); // O
206 10 : reference[5]=Vector(-1.698, -1.687, -0.215); // N i+1
207 10 : reference[6]=Vector(-2.681, -0.613, -0.143); // CA
208 10 : reference[7]=Vector(-3.323, -0.477, 1.267); // CB
209 10 : reference[8]=Vector(-1.984, 0.681, -0.574); // C
210 10 : reference[9]=Vector(-0.807, 0.921, -0.273); // O
211 10 : reference[10]=Vector(-2.716, 1.492, -1.329); // N i+2
212 10 : reference[11]=Vector(-2.196, 2.731, -1.883); // CA
213 10 : reference[12]=Vector(-2.263, 2.692, -3.418); // CB
214 10 : reference[13]=Vector(-2.989, 3.949, -1.433); // C
215 10 : reference[14]=Vector(-4.214, 3.989, -1.583); // O
216 10 : reference[15]=Vector( 2.464, -4.352, 2.149); // N j
217 10 : reference[16]=Vector( 3.078, -3.170, 1.541); // CA
218 10 : reference[17]=Vector( 3.398, -3.415, 0.060); // CB
219 10 : reference[18]=Vector( 2.080, -2.021, 1.639); // C
220 10 : reference[19]=Vector( 0.938, -2.178, 1.225); // O
221 10 : reference[20]=Vector( 2.525, -0.886, 2.183); // N j+1
222 10 : reference[21]=Vector( 1.692, 0.303, 2.346); // CA
223 10 : reference[22]=Vector( 1.541, 0.665, 3.842); // CB
224 10 : reference[23]=Vector( 2.420, 1.410, 1.608); // C
225 10 : reference[24]=Vector( 3.567, 1.733, 1.937); // O
226 10 : reference[25]=Vector( 1.758, 1.976, 0.600); // N j+2
227 10 : reference[26]=Vector( 2.373, 2.987, -0.238); // CA
228 10 : reference[27]=Vector( 2.367, 2.527, -1.720); // CB
229 10 : reference[28]=Vector( 1.684, 4.331, -0.148); // C
230 10 : reference[29]=Vector( 0.486, 4.430, -0.415); // O
231 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
232 20 : setSecondaryStructure( reference, 0.17/atoms.getUnits().getLength(), 0.1/atoms.getUnits().getLength() );
233 10 : }
234 :
235 : }
236 2523 : }
|
