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 ALPHARMSD
31 : /*
32 : Probe the alpha helical content of a protein structure.
33 :
34 : Any chain of six contiguous residues in a protein chain can form an alpha helix. This
35 : colvar thus generates the set of all possible six residue sections and calculates
36 : the RMSD distance between the configuration in which the residues find themselves
37 : and an idealized alpha helical structure. These distances can be calculated by either
38 : aligning the instantaneous structure with the reference structure and measuring each
39 : atomic displacement or by calculating differences between the set of interatomic
40 : distances in the reference and instantaneous structures.
41 :
42 : This colvar is based on the following reference \cite pietrucci09jctc. The authors of
43 : this paper use the set of distances from the alpha helix configurations to measure
44 : the number of segments that have an alpha helical configuration. This is done by calculating
45 : the following sum of functions of the rmsd distances:
46 :
47 : \f[
48 : 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 }
49 : \f]
50 :
51 : where the sum runs over all possible segments of alpha helix. By default the
52 : NN, MM and D_0 parameters are set equal to those used in \cite pietrucci09jctc. The R_0
53 : parameter must be set by the user - the value used in \cite pietrucci09jctc was 0.08 nm.
54 :
55 : If you change the function in the above sum you can calculate quantities such as the average
56 : distance from a purely the alpha helical configuration or the distance between the set of
57 : residues that is closest to an alpha helix and the reference configuration. To do these sorts of
58 : calculations you can use the AVERAGE and MIN keywords. In addition you can use the LESS_THAN
59 : keyword if you would like to change the form of the switching function. If you use any of these
60 : options you no longer need to specify NN, R_0, MM and D_0.
61 :
62 : Please be aware that for codes like gromacs you must ensure that plumed
63 : reconstructs the chains involved in your CV when you calculate this CV using
64 : anthing other than TYPE=DRMSD. For more details as to how to do this see \ref WHOLEMOLECULES.
65 :
66 : \par Examples
67 :
68 : The following input calculates the number of six residue segments of
69 : protein that are in an alpha helical configuration.
70 :
71 : \verbatim
72 : MOLINFO STRUCTURE=helix.pdb
73 : ALPHARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
74 : \endverbatim
75 : (see also \ref MOLINFO)
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 6 : class AlphaRMSD : public SecondaryStructureRMSD {
81 : public:
82 : static void registerKeywords( Keywords& keys );
83 : explicit AlphaRMSD(const ActionOptions&);
84 : };
85 :
86 2526 : PLUMED_REGISTER_ACTION(AlphaRMSD,"ALPHARMSD")
87 :
88 4 : void AlphaRMSD::registerKeywords( Keywords& keys ) {
89 4 : SecondaryStructureRMSD::registerKeywords( keys );
90 4 : }
91 :
92 3 : AlphaRMSD::AlphaRMSD(const ActionOptions&ao):
93 : Action(ao),
94 3 : SecondaryStructureRMSD(ao)
95 : {
96 : // read in the backbone atoms
97 3 : std::vector<unsigned> chains; readBackboneAtoms( "protein", chains);
98 :
99 : // This constructs all conceivable sections of alpha helix in the backbone of the chains
100 6 : unsigned nres, nprevious=0; std::vector<unsigned> nlist(30);
101 6 : for(unsigned i=0; i<chains.size(); ++i) {
102 3 : if( chains[i]<30 ) error("segment of backbone defined is not long enough to form an alpha helix. Each backbone fragment must contain a minimum of 6 residues");
103 3 : nres=chains[i]/5;
104 3 : if( chains[i]%5!=0 ) error("backbone segment received does not contain a multiple of five residues");
105 12 : for(unsigned ires=0; ires<nres-5; ires++) {
106 9 : unsigned accum=nprevious + 5*ires;
107 9 : for(unsigned k=0; k<30; ++k) nlist[k] = accum+k;
108 9 : addColvar( nlist );
109 : }
110 3 : nprevious+=chains[i];
111 : }
112 :
113 : // Build the reference structure ( in angstroms )
114 6 : std::vector<Vector> reference(30);
115 3 : reference[0] = Vector( 0.733, 0.519, 5.298 ); // N i
116 3 : reference[1] = Vector( 1.763, 0.810, 4.301 ); // CA
117 3 : reference[2] = Vector( 3.166, 0.543, 4.881 ); // CB
118 3 : reference[3] = Vector( 1.527, -0.045, 3.053 ); // C
119 3 : reference[4] = Vector( 1.646, 0.436, 1.928 ); // O
120 3 : reference[5] = Vector( 1.180, -1.312, 3.254 ); // N i+1
121 3 : reference[6] = Vector( 0.924, -2.203, 2.126 ); // CA
122 3 : reference[7] = Vector( 0.650, -3.626, 2.626 ); // CB
123 3 : reference[8] = Vector(-0.239, -1.711, 1.261 ); // C
124 3 : reference[9] = Vector(-0.190, -1.815, 0.032 ); // O
125 3 : reference[10] = Vector(-1.280, -1.172, 1.891 ); // N i+2
126 3 : reference[11] = Vector(-2.416, -0.661, 1.127 ); // CA
127 3 : reference[12] = Vector(-3.548, -0.217, 2.056 ); // CB
128 3 : reference[13] = Vector(-1.964, 0.529, 0.276 ); // C
129 3 : reference[14] = Vector(-2.364, 0.659, -0.880 ); // O
130 3 : reference[15] = Vector(-1.130, 1.391, 0.856 ); // N i+3
131 3 : reference[16] = Vector(-0.620, 2.565, 0.148 ); // CA
132 3 : reference[17] = Vector( 0.228, 3.439, 1.077 ); // CB
133 3 : reference[18] = Vector( 0.231, 2.129, -1.032 ); // C
134 3 : reference[19] = Vector( 0.179, 2.733, -2.099 ); // O
135 3 : reference[20] = Vector( 1.028, 1.084, -0.833 ); // N i+4
136 3 : reference[21] = Vector( 1.872, 0.593, -1.919 ); // CA
137 3 : reference[22] = Vector( 2.850, -0.462, -1.397 ); // CB
138 3 : reference[23] = Vector( 1.020, 0.020, -3.049 ); // C
139 3 : reference[24] = Vector( 1.317, 0.227, -4.224 ); // O
140 3 : reference[25] = Vector(-0.051, -0.684, -2.696 ); // N i+5
141 3 : reference[26] = Vector(-0.927, -1.261, -3.713 ); // CA
142 3 : reference[27] = Vector(-1.933, -2.219, -3.074 ); // CB
143 3 : reference[28] = Vector(-1.663, -0.171, -4.475 ); // C
144 3 : reference[29] = Vector(-1.916, -0.296, -5.673 ); // O
145 : // Store the secondary structure ( last number makes sure we convert to internal units nm )
146 6 : setSecondaryStructure( reference, 0.17/atoms.getUnits().getLength(), 0.1/atoms.getUnits().getLength() );
147 3 : }
148 :
149 : }
150 2523 : }
|
