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 ALPHARMSD
30 : /*
31 : Probe the alpha helical content of a protein structure.
32 :
33 : Any chain of six contiguous residues in a protein chain can form an alpha helix. This
34 : colvar thus generates the set of all possible six residue sections and calculates
35 : the [DRMSD](DRMSD.md) or [RMSD](RMSD.md) distance between the configuration in which the residues find themselves
36 : and an idealized alpha helical structure. These distances can be calculated by either
37 : aligning the instantaneous structure with the reference structure and measuring each
38 : atomic displacement or by calculating differences between the set of inter-atomic
39 : distances in the reference and instantaneous structures.
40 :
41 : In the original paper on this method (see below) the authors used the set of distances from the alpha helix configurations to measure
42 : the number of segments that have an alpha helical configuration. This is done by calculating
43 : the following sum of functions of the rmsd distances:
44 :
45 : $$
46 : 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 }
47 : $$
48 :
49 : where the sum runs over all possible segments of alpha helix. By default the
50 : NN, MM and D_0 parameters are set equal to those used in the paper cited below. The R_0
51 : parameter must be set by the user - the value used in the paper cited below was 0.08 nm.
52 :
53 : If you change the function in the above sum you can calculate quantities such as the average
54 : distance from a purely the alpha helical configuration or the distance between the set of
55 : residues that is closest to an alpha helix and the reference configuration. To do these sorts of
56 : calculations you can use the AVERAGE and MIN keywords. In addition you can use the LESS_THAN
57 : keyword if you would like to change the form of the switching function. If you use any of these
58 : options you no longer need to specify NN, R_0, MM and D_0.
59 :
60 : The following input calculates the number of six residue segments of
61 : protein that are in an alpha helical configuration.
62 :
63 : ```plumed
64 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
65 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
66 : alpha: ALPHARMSD RESIDUES=all R_0=0.1
67 : PRINT ARG=alpha FILE=colvar
68 : ```
69 :
70 : Here the same is done use [RMSD](RMSD.md) instead of [DRMSD](DRMSD.md), which is normally faster.
71 :
72 : ```plumed
73 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
74 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
75 : WHOLEMOLECULES ENTITY0=1-100
76 : alpha: ALPHARMSD RESIDUES=all TYPE=OPTIMAL LESS_THAN={RATIONAL R_0=0.1 NN=8 MM=12}
77 : PRINT ARG=alpha.lessthan FILE=colvar
78 : ```
79 :
80 : ## Periodic boundary conditions
81 :
82 : You can turn off periodic boundary conditions by using the NOPBC flag as shown below:
83 :
84 : ```plumed
85 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
86 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
87 : alpha: ALPHARMSD RESIDUES=all R_0=0.1 NOPBC
88 : PRINT ARG=alpha FILE=colvar
89 : ```
90 :
91 : If you are using [DRMSD](DRMSD.md) to measure distances and you __don't__ use the NOPBC flag then
92 : all distances in the instaneous structure are evaluated in a way that takes the periodic boundary conditions
93 : into account. Using the NOPBC flag turns off this treatment.
94 :
95 : If you are using [RMSD](RMSD.md) to measure distances and you __don't__ use the NOPBC flag the the instaneous positions of
96 : each segment for which the RMSD is computed is reconstructed using the procedure outlined in the documentation for [WHOLEMOLECULES](WHOLEMOLECULES.md)
97 : before the RMSD is computed. Using the NOPBC flag turns off this treatment.
98 :
99 : */
100 : //+ENDPLUMEDOC
101 :
102 : class AlphaRMSD : public ActionShortcut {
103 : public:
104 : static void registerKeywords( Keywords& keys );
105 : explicit AlphaRMSD(const ActionOptions&);
106 : };
107 :
108 : PLUMED_REGISTER_ACTION(AlphaRMSD,"ALPHARMSD")
109 :
110 8 : void AlphaRMSD::registerKeywords( Keywords& keys ) {
111 8 : SecondaryStructureBase<Vector>::registerKeywords( keys );
112 8 : keys.remove("ATOMS");
113 8 : keys.remove("SEGMENT");
114 8 : keys.remove("STRUCTURE");
115 8 : keys.remove("MASK");
116 8 : keys.remove("ALIGN_STRANDS");
117 16 : keys.setValueDescription("scalar/vector","if LESS_THAN is present the RMSD distance between each residue and the ideal alpha helix. If LESS_THAN is not present the number of residue segments where the structure is similar to an alpha helix");
118 8 : }
119 :
120 3 : AlphaRMSD::AlphaRMSD(const ActionOptions&ao):
121 : Action(ao),
122 3 : ActionShortcut(ao) {
123 : // Read in the input and create a string that describes how to compute the less than
124 : std::string ltmap;
125 3 : bool uselessthan=SecondaryStructureBase<Vector>::readShortcutWords( ltmap, this );
126 : // read in the backbone atoms
127 : std::vector<unsigned> chains;
128 : std::vector<std::string> all_atoms;
129 3 : SecondaryStructureBase<Vector>::readBackboneAtoms( this, plumed, "protein", chains, all_atoms );
130 :
131 : // This constructs all conceivable sections of alpha helix in the backbone of the chains
132 3 : unsigned nprevious=0, segno=1;
133 : std::string seglist;
134 6 : for(unsigned i=0; i<chains.size(); ++i) {
135 3 : if( chains[i]<30 ) {
136 0 : error("segment of backbone defined is not long enough to form an alpha helix. Each backbone fragment must contain a minimum of 6 residues");
137 : }
138 3 : unsigned nres=chains[i]/5;
139 3 : if( chains[i]%5!=0 ) {
140 0 : error("backbone segment received does not contain a multiple of five residues");
141 : }
142 18 : for(unsigned ires=0; ires<nres-5; ires++) {
143 15 : unsigned accum=nprevious + 5*ires;
144 : std::string strval, num;
145 15 : Tools::convert( segno, num );
146 15 : Tools::convert( accum, strval );
147 30 : seglist += " SEGMENT" + num + "=" + strval;
148 450 : for(unsigned k=1; k<30; ++k) {
149 435 : Tools::convert( accum+k, strval );
150 870 : seglist += "," + strval;
151 : }
152 15 : segno++;
153 : }
154 3 : nprevious+=chains[i];
155 : }
156 :
157 : // Build the reference structure ( in angstroms )
158 3 : std::vector<Vector> reference(30);
159 3 : reference[0] = Vector( 0.733, 0.519, 5.298 ); // N i
160 3 : reference[1] = Vector( 1.763, 0.810, 4.301 ); // CA
161 3 : reference[2] = Vector( 3.166, 0.543, 4.881 ); // CB
162 3 : reference[3] = Vector( 1.527, -0.045, 3.053 ); // C
163 3 : reference[4] = Vector( 1.646, 0.436, 1.928 ); // O
164 3 : reference[5] = Vector( 1.180, -1.312, 3.254 ); // N i+1
165 3 : reference[6] = Vector( 0.924, -2.203, 2.126 ); // CA
166 3 : reference[7] = Vector( 0.650, -3.626, 2.626 ); // CB
167 3 : reference[8] = Vector(-0.239, -1.711, 1.261 ); // C
168 3 : reference[9] = Vector(-0.190, -1.815, 0.032 ); // O
169 3 : reference[10] = Vector(-1.280, -1.172, 1.891 ); // N i+2
170 3 : reference[11] = Vector(-2.416, -0.661, 1.127 ); // CA
171 3 : reference[12] = Vector(-3.548, -0.217, 2.056 ); // CB
172 3 : reference[13] = Vector(-1.964, 0.529, 0.276 ); // C
173 3 : reference[14] = Vector(-2.364, 0.659, -0.880 ); // O
174 3 : reference[15] = Vector(-1.130, 1.391, 0.856 ); // N i+3
175 3 : reference[16] = Vector(-0.620, 2.565, 0.148 ); // CA
176 3 : reference[17] = Vector( 0.228, 3.439, 1.077 ); // CB
177 3 : reference[18] = Vector( 0.231, 2.129, -1.032 ); // C
178 3 : reference[19] = Vector( 0.179, 2.733, -2.099 ); // O
179 3 : reference[20] = Vector( 1.028, 1.084, -0.833 ); // N i+4
180 3 : reference[21] = Vector( 1.872, 0.593, -1.919 ); // CA
181 3 : reference[22] = Vector( 2.850, -0.462, -1.397 ); // CB
182 3 : reference[23] = Vector( 1.020, 0.020, -3.049 ); // C
183 3 : reference[24] = Vector( 1.317, 0.227, -4.224 ); // O
184 3 : reference[25] = Vector(-0.051, -0.684, -2.696 ); // N i+5
185 3 : reference[26] = Vector(-0.927, -1.261, -3.713 ); // CA
186 3 : reference[27] = Vector(-1.933, -2.219, -3.074 ); // CB
187 3 : reference[28] = Vector(-1.663, -0.171, -4.475 ); // C
188 3 : reference[29] = Vector(-1.916, -0.296, -5.673 ); // O
189 : std::string ref0, ref1, ref2;
190 3 : Tools::convert( reference[0][0], ref0 );
191 3 : Tools::convert( reference[0][1], ref1 );
192 3 : Tools::convert( reference[0][2], ref2 );
193 6 : std::string structure=" STRUCTURE1=" + ref0 + "," + ref1 + "," + ref2;
194 90 : for(unsigned i=1; i<30; ++i) {
195 348 : for(unsigned k=0; k<3; ++k) {
196 261 : Tools::convert( reference[i][k], ref0 );
197 522 : structure += "," + ref0;
198 : }
199 : }
200 :
201 : std::string type;
202 3 : parse("TYPE",type);
203 3 : std::string lab = getShortcutLabel() + "_rmsd";
204 3 : if( uselessthan ) {
205 1 : lab = getShortcutLabel();
206 : }
207 3 : std::string nopbcstr="";
208 : bool nopbc;
209 3 : parseFlag("NOPBC",nopbc);
210 3 : if( nopbc ) {
211 : nopbcstr = " NOPBC";
212 : }
213 3 : std::string usegpustr="";
214 : {
215 : bool usegpu;
216 3 : parseFlag("USEGPU",usegpu);
217 3 : if( usegpu ) {
218 : usegpustr = " USEGPU";
219 : }
220 : }
221 :
222 3 : std::string atoms="ATOMS=" + all_atoms[0];
223 150 : for(unsigned i=1; i<all_atoms.size(); ++i) {
224 294 : atoms += "," + all_atoms[i];
225 : }
226 3 : if( type=="DRMSD" ) {
227 6 : readInputLine( lab + ": SECONDARY_STRUCTURE_DRMSD BONDLENGTH=0.17" + seglist + structure + " " + atoms + nopbcstr + usegpustr);
228 : } else {
229 0 : readInputLine( lab + ": SECONDARY_STRUCTURE_RMSD " + seglist + structure + " " + atoms + " TYPE=" + type + nopbcstr + usegpustr);
230 : }
231 : // Create the less than object
232 3 : if( ltmap.length()>0 ) {
233 6 : readInputLine( getShortcutLabel() + "_lt: LESS_THAN ARG=" + lab + " SWITCH={" + ltmap +"}");
234 3 : if( uselessthan ) {
235 2 : readInputLine( getShortcutLabel() + "_lessthan: SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
236 : } else {
237 4 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_lt PERIODIC=NO");
238 : }
239 : }
240 6 : }
241 :
242 : }
243 : }
|
