Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 ANTIBETARMSD
31 : /*
32 : Probe the antiparallel beta sheet content of your protein structure.
33 :
34 : Two protein segments containing three continguous residues can form an antiparallel beta sheet.
35 : Although if the two segments are part of the same protein chain they must be separated by
36 : a minimum of 2 residues to make room for the turn. This colvar thus generates the set of
37 : all possible six residue sections that could conceivably form an antiparallel beta sheet
38 : and calculates the RMSD distance between the configuration in which the residues find themselves
39 : and an idealized antiparallel 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 anti parallel beta sheet configurations to measure
46 : the number of segments that have an configuration that resemebles an anti paralel 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 antiparallel 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 purely configuration composed of pure anti-parallel beta sheets or the distance between the set of
59 : residues that is closest to an anti-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 antiparallel beta sheet configuration.
72 :
73 : \plumedfile
74 : MOLINFO STRUCTURE=beta.pdb
75 : ab: ANTIBETARMSD RESIDUES=all STRANDS_CUTOFF=1
76 : \endplumedfile
77 :
78 : Here the same is done use RMSD instead of DRMSD
79 :
80 : \plumedfile
81 : MOLINFO STRUCTURE=helix.pdb
82 : WHOLEMOLECULES ENTITY0=1-100
83 : hh: ANTIBETARMSD RESIDUES=all TYPE=OPTIMAL R_0=0.1 STRANDS_CUTOFF=1
84 : \endplumedfile
85 : */
86 : //+ENDPLUMEDOC
87 :
88 24 : class AntibetaRMSD : public SecondaryStructureRMSD {
89 : public:
90 : static void registerKeywords( Keywords& keys );
91 : explicit AntibetaRMSD(const ActionOptions&);
92 : };
93 :
94 6464 : PLUMED_REGISTER_ACTION(AntibetaRMSD,"ANTIBETARMSD")
95 :
96 13 : void AntibetaRMSD::registerKeywords( Keywords& keys ) {
97 13 : SecondaryStructureRMSD::registerKeywords( keys );
98 65 : 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 "
99 : "chain configuration with the appropriate geometry are counted. If STYLE=inter "
100 : "only sheet-like configurations involving two chains are counted, while if STYLE=intra "
101 : "only sheet-like configurations involving a single chain are counted");
102 26 : keys.use("STRANDS_CUTOFF");
103 13 : }
104 :
105 12 : AntibetaRMSD::AntibetaRMSD(const ActionOptions&ao):
106 : Action(ao),
107 12 : SecondaryStructureRMSD(ao)
108 : {
109 : // read in the backbone atoms
110 24 : std::vector<unsigned> chains; readBackboneAtoms( "protein", chains );
111 :
112 : bool intra_chain(false), inter_chain(false);
113 24 : std::string style; parse("STYLE",style);
114 12 : if( style=="all" ) {
115 : intra_chain=true; inter_chain=true;
116 2 : } else if( style=="inter") {
117 : intra_chain=false; inter_chain=true;
118 0 : } else if( style=="intra") {
119 : intra_chain=true; inter_chain=false;
120 : } else {
121 0 : error( style + " is not a valid directive for the STYLE keyword");
122 : }
123 :
124 : // Align the atoms based on the positions of these two atoms
125 12 : setAtomsFromStrands( 6, 21 );
126 :
127 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
128 12 : if( intra_chain ) {
129 10 : unsigned nprevious=0; std::vector<unsigned> nlist(30);
130 509 : for(unsigned i=0; i<chains.size(); ++i) {
131 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");
132 : // Loop over all possible triples in each 8 residue segment of protein
133 163 : unsigned nres=chains[i]/5;
134 163 : if( chains[i]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
135 497 : for(unsigned ires=0; ires<nres-7; ires++) {
136 344 : for(unsigned jres=ires+7; jres<nres; jres++) {
137 5487 : for(unsigned k=0; k<15; ++k) {
138 5310 : nlist[k]=nprevious + ires*5+k;
139 5310 : nlist[k+15]=nprevious + (jres-2)*5+k;
140 : }
141 177 : addColvar( nlist );
142 : }
143 : }
144 163 : nprevious+=chains[i];
145 : }
146 : }
147 12 : if( inter_chain ) {
148 13 : if( chains.size()==1 && style!="all" ) error("there is only one chain defined so cannot use inter_chain option");
149 12 : std::vector<unsigned> nlist(30);
150 489 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
151 2913 : unsigned iprev=0; for(unsigned i=0; i<ichain; ++i) iprev+=chains[i];
152 155 : unsigned inres=chains[ichain]/5;
153 155 : if( chains[ichain]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
154 1995 : for(unsigned ires=0; ires<inres-2; ++ires) {
155 17448 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
156 96392 : unsigned jprev=0; for(unsigned i=0; i<jchain; ++i) jprev+=chains[i];
157 16528 : unsigned jnres=chains[jchain]/5;
158 8264 : if( chains[jchain]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
159 107412 : for(unsigned jres=0; jres<jnres-2; ++jres) {
160 1536794 : for(unsigned k=0; k<15; ++k) {
161 1487220 : nlist[k]=iprev+ ires*5+k;
162 1487220 : nlist[k+15]=jprev+ jres*5+k;
163 : }
164 49574 : addColvar( nlist );
165 : }
166 : }
167 : }
168 : }
169 : }
170 :
171 : // Build the reference structure ( in angstroms )
172 12 : std::vector<Vector> reference(30);
173 12 : reference[0]=Vector( 2.263, -3.795, 1.722); // N i
174 12 : reference[1]=Vector( 2.493, -2.426, 2.263); // CA
175 12 : reference[2]=Vector( 3.847, -1.838, 1.761); // CB
176 12 : reference[3]=Vector( 1.301, -1.517, 1.921); // C
177 12 : reference[4]=Vector( 0.852, -1.504, 0.739); // O
178 12 : reference[5]=Vector( 0.818, -0.738, 2.917); // N i+1
179 12 : reference[6]=Vector(-0.299, 0.243, 2.748); // CA
180 12 : reference[7]=Vector(-1.421, -0.076, 3.757); // CB
181 12 : reference[8]=Vector( 0.273, 1.680, 2.854); // C
182 12 : reference[9]=Vector( 0.902, 1.993, 3.888); // O
183 12 : reference[10]=Vector( 0.119, 2.532, 1.813); // N i+2
184 12 : reference[11]=Vector( 0.683, 3.916, 1.680); // CA
185 12 : reference[12]=Vector( 1.580, 3.940, 0.395); // CB
186 12 : reference[13]=Vector(-0.394, 5.011, 1.630); // C
187 12 : reference[14]=Vector(-1.459, 4.814, 0.982); // O
188 12 : reference[15]=Vector(-2.962, 3.559, -1.359); // N j-2
189 12 : reference[16]=Vector(-2.439, 2.526, -2.287); // CA
190 12 : reference[17]=Vector(-1.189, 3.006, -3.087); // CB
191 12 : reference[18]=Vector(-2.081, 1.231, -1.520); // C
192 12 : reference[19]=Vector(-1.524, 1.324, -0.409); // O
193 12 : reference[20]=Vector(-2.326, 0.037, -2.095); // N j-1
194 12 : reference[21]=Vector(-1.858, -1.269, -1.554); // CA
195 12 : reference[22]=Vector(-3.053, -2.199, -1.291); // CB
196 12 : reference[23]=Vector(-0.869, -1.949, -2.512); // C
197 12 : reference[24]=Vector(-1.255, -2.070, -3.710); // O
198 12 : reference[25]=Vector( 0.326, -2.363, -2.072); // N j
199 12 : reference[26]=Vector( 1.405, -2.992, -2.872); // CA
200 12 : reference[27]=Vector( 2.699, -2.129, -2.917); // CB
201 12 : reference[28]=Vector( 1.745, -4.399, -2.330); // C
202 12 : reference[29]=Vector( 1.899, -4.545, -1.102); // O
203 :
204 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
205 12 : setSecondaryStructure( reference, 0.17/atoms.getUnits().getLength(), 0.1/atoms.getUnits().getLength() );
206 12 : }
207 :
208 : }
209 4839 : }
|
