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 ANTIBETARMSD
30 : /*
31 : Probe the antiparallel beta sheet content of your protein structure.
32 :
33 : Two protein segments containing three contiguous residues can form an antiparallel beta sheet.
34 : Although if the two segments are part of the same protein chain they must be separated by
35 : a minimum of 2 residues to make room for the turn. This colvar thus generates the set of
36 : all possible six residue sections that could conceivably form an antiparallel beta sheet
37 : and calculates the RMSD distance between the configuration in which the residues find themselves
38 : and an idealized antiparallel 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 anti parallel beta sheet configurations to measure
45 : the number of segments that have an configuration that resembles an anti 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 antiparallel 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 purely configuration composed of pure anti-parallel beta sheets or the distance between the set of
58 : residues that is closest to an anti-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 antiparallel beta sheet configuration.
71 :
72 : \plumedfile
73 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
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 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
82 : MOLINFO STRUCTURE=helix.pdb
83 : WHOLEMOLECULES ENTITY0=1-100
84 : hh: ANTIBETARMSD RESIDUES=all TYPE=OPTIMAL R_0=0.1 STRANDS_CUTOFF=1
85 : \endplumedfile
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : class AntibetaRMSD : public ActionShortcut {
90 : public:
91 : static void registerKeywords( Keywords& keys );
92 : explicit AntibetaRMSD(const ActionOptions&);
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(AntibetaRMSD,"ANTIBETARMSD")
96 :
97 101 : void AntibetaRMSD::registerKeywords( Keywords& keys ) {
98 101 : SecondaryStructureRMSD::registerKeywords( keys );
99 101 : keys.remove("ATOMS");
100 101 : keys.remove("SEGMENT");
101 101 : keys.remove("BONDLENGTH");
102 101 : keys.remove("NO_ACTION_LOG");
103 101 : keys.remove("CUTOFF_ATOMS");
104 101 : keys.remove("STRUCTURE");
105 202 : 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 "
106 : "chain configuration with the appropriate geometry are counted. If STYLE=inter "
107 : "only sheet-like configurations involving two chains are counted, while if STYLE=intra "
108 : "only sheet-like configurations involving a single chain are counted");
109 101 : }
110 :
111 17 : AntibetaRMSD::AntibetaRMSD(const ActionOptions&ao):
112 : Action(ao),
113 17 : ActionShortcut(ao) {
114 : // Read in the input and create a string that describes how to compute the less than
115 : std::string ltmap;
116 17 : bool uselessthan=SecondaryStructureRMSD::readShortcutWords( ltmap, this );
117 : // read in the backbone atoms
118 : std::vector<unsigned> chains;
119 : std::string atoms;
120 34 : SecondaryStructureRMSD::readBackboneAtoms( this, plumed, "protein", chains, atoms );
121 :
122 : bool intra_chain(false), inter_chain(false);
123 : std::string style;
124 17 : parse("STYLE",style);
125 34 : if( Tools::caseInSensStringCompare(style, "all") ) {
126 : intra_chain=true;
127 : inter_chain=true;
128 4 : } else if( Tools::caseInSensStringCompare(style, "inter") ) {
129 : intra_chain=false;
130 : inter_chain=true;
131 2 : } else if( Tools::caseInSensStringCompare(style, "intra") ) {
132 : intra_chain=true;
133 : inter_chain=false;
134 : } else {
135 0 : error( style + " is not a valid directive for the STYLE keyword");
136 : }
137 :
138 : // This constructs all conceivable sections of antibeta sheet in the backbone of the chains
139 : std::string seglist;
140 17 : unsigned k=1;
141 17 : if( intra_chain ) {
142 : unsigned nprevious=0;
143 16 : std::vector<unsigned> nlist(30);
144 287 : for(unsigned i=0; i<chains.size(); ++i) {
145 271 : if( chains[i]<40 ) {
146 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");
147 : }
148 : // Loop over all possible triples in each 8 residue segment of protein
149 271 : unsigned nres=chains[i]/5;
150 271 : if( chains[i]%5!=0 ) {
151 0 : error("backbone segment received does not contain a multiple of five residues");
152 : }
153 546 : for(unsigned ires=0; ires<nres-7; ires++) {
154 560 : for(unsigned jres=ires+7; jres<nres; jres++) {
155 4560 : for(unsigned k=0; k<15; ++k) {
156 4275 : nlist[k]=nprevious + ires*5+k;
157 4275 : nlist[k+15]=nprevious + (jres-2)*5+k;
158 : }
159 : std::string nlstr, num;
160 285 : Tools::convert( nlist[0], nlstr );
161 285 : Tools::convert(k, num);
162 285 : k++;
163 570 : seglist += " SEGMENT" + num + "=" + nlstr;
164 8550 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
165 8265 : Tools::convert( nlist[kk], nlstr );
166 16530 : seglist += "," + nlstr;
167 : }
168 : }
169 : }
170 271 : nprevious+=chains[i];
171 : }
172 : }
173 17 : if( inter_chain ) {
174 17 : if( chains.size()==1 && style!="all" ) {
175 0 : error("there is only one chain defined so cannot use inter_chain option");
176 : }
177 16 : std::vector<unsigned> nlist(30);
178 255 : for(unsigned ichain=1; ichain<chains.size(); ++ichain) {
179 : unsigned iprev=0;
180 2382 : for(unsigned i=0; i<ichain; ++i) {
181 2143 : iprev+=chains[i];
182 : }
183 239 : unsigned inres=chains[ichain]/5;
184 239 : if( chains[ichain]%5!=0 ) {
185 0 : error("backbone segment received does not contain a multiple of five residues");
186 : }
187 1668 : for(unsigned ires=0; ires<inres-2; ++ires) {
188 14282 : for(unsigned jchain=0; jchain<ichain; ++jchain) {
189 : unsigned jprev=0;
190 81397 : for(unsigned i=0; i<jchain; ++i) {
191 68544 : jprev+=chains[i];
192 : }
193 12853 : unsigned jnres=chains[jchain]/5;
194 12853 : if( chains[jchain]%5!=0 ) {
195 0 : error("backbone segment received does not contain a multiple of five residues");
196 : }
197 89966 : for(unsigned jres=0; jres<jnres-2; ++jres) {
198 1233808 : for(unsigned k=0; k<15; ++k) {
199 1156695 : nlist[k]=iprev+ ires*5+k;
200 1156695 : nlist[k+15]=jprev+ jres*5+k;
201 : }
202 : std::string nlstr, num;
203 77113 : Tools::convert( nlist[0], nlstr );
204 77113 : Tools::convert(k, num);
205 77113 : k++;
206 154226 : seglist += " SEGMENT" + num + "=" + nlstr;
207 2313390 : for(unsigned kk=1; kk<nlist.size(); ++kk ) {
208 2236277 : Tools::convert( nlist[kk], nlstr );
209 4472554 : seglist += "," + nlstr;
210 : }
211 : }
212 : }
213 : }
214 : }
215 : }
216 :
217 : // Build the reference structure ( in angstroms )
218 17 : std::vector<Vector> reference(30);
219 17 : reference[0]=Vector( 2.263, -3.795, 1.722); // N i
220 17 : reference[1]=Vector( 2.493, -2.426, 2.263); // CA
221 17 : reference[2]=Vector( 3.847, -1.838, 1.761); // CB
222 17 : reference[3]=Vector( 1.301, -1.517, 1.921); // C
223 17 : reference[4]=Vector( 0.852, -1.504, 0.739); // O
224 17 : reference[5]=Vector( 0.818, -0.738, 2.917); // N i+1
225 17 : reference[6]=Vector(-0.299, 0.243, 2.748); // CA
226 17 : reference[7]=Vector(-1.421, -0.076, 3.757); // CB
227 17 : reference[8]=Vector( 0.273, 1.680, 2.854); // C
228 17 : reference[9]=Vector( 0.902, 1.993, 3.888); // O
229 17 : reference[10]=Vector( 0.119, 2.532, 1.813); // N i+2
230 17 : reference[11]=Vector( 0.683, 3.916, 1.680); // CA
231 17 : reference[12]=Vector( 1.580, 3.940, 0.395); // CB
232 17 : reference[13]=Vector(-0.394, 5.011, 1.630); // C
233 17 : reference[14]=Vector(-1.459, 4.814, 0.982); // O
234 17 : reference[15]=Vector(-2.962, 3.559, -1.359); // N j-2
235 17 : reference[16]=Vector(-2.439, 2.526, -2.287); // CA
236 17 : reference[17]=Vector(-1.189, 3.006, -3.087); // CB
237 17 : reference[18]=Vector(-2.081, 1.231, -1.520); // C
238 17 : reference[19]=Vector(-1.524, 1.324, -0.409); // O
239 17 : reference[20]=Vector(-2.326, 0.037, -2.095); // N j-1
240 17 : reference[21]=Vector(-1.858, -1.269, -1.554); // CA
241 17 : reference[22]=Vector(-3.053, -2.199, -1.291); // CB
242 17 : reference[23]=Vector(-0.869, -1.949, -2.512); // C
243 17 : reference[24]=Vector(-1.255, -2.070, -3.710); // O
244 17 : reference[25]=Vector( 0.326, -2.363, -2.072); // N j
245 17 : reference[26]=Vector( 1.405, -2.992, -2.872); // CA
246 17 : reference[27]=Vector( 2.699, -2.129, -2.917); // CB
247 17 : reference[28]=Vector( 1.745, -4.399, -2.330); // C
248 17 : reference[29]=Vector( 1.899, -4.545, -1.102); // O
249 : std::string ref0, ref1, ref2;
250 17 : Tools::convert( reference[0][0], ref0 );
251 17 : Tools::convert( reference[0][1], ref1 );
252 17 : Tools::convert( reference[0][2], ref2 );
253 34 : std::string structure=" STRUCTURE1=" + ref0 + "," + ref1 + "," + ref2;
254 510 : for(unsigned i=1; i<30; ++i) {
255 1972 : for(unsigned k=0; k<3; ++k) {
256 1479 : Tools::convert( reference[i][k], ref0 );
257 2958 : structure += "," + ref0;
258 : }
259 : }
260 :
261 : std::string strands_cutoff;
262 34 : parse("STRANDS_CUTOFF",strands_cutoff);
263 17 : if( strands_cutoff.length()>0 ) {
264 32 : strands_cutoff=" CUTOFF_ATOMS=6,21 STRANDS_CUTOFF="+strands_cutoff;
265 : }
266 : std::string type;
267 17 : parse("TYPE",type);
268 17 : std::string lab = getShortcutLabel() + "_rmsd";
269 17 : if( uselessthan ) {
270 17 : lab = getShortcutLabel();
271 : }
272 17 : std::string nopbcstr="";
273 : bool nopbc;
274 17 : parseFlag("NOPBC",nopbc);
275 17 : if( nopbc ) {
276 : nopbcstr = " NOPBC";
277 : }
278 17 : if( seglist.length()==0 ) {
279 0 : error("no segments to investigate");
280 : }
281 34 : readInputLine( lab + ": SECONDARY_STRUCTURE_RMSD BONDLENGTH=0.17" + seglist + structure + " " + atoms + " TYPE=" + type + strands_cutoff + nopbcstr );
282 : // Create the less than object
283 17 : if( ltmap.length()>0 ) {
284 17 : SecondaryStructureRMSD::expandShortcut( uselessthan, getShortcutLabel(), lab, ltmap, this );
285 : }
286 17 : }
287 :
288 : }
289 : }
|
