3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
43 #include "gmx_fatal.h"
66 void dump_ahx(int nres
,
67 t_bb bb
[],rvec x
[],matrix box
,int teller
)
73 sprintf(buf
,"dump%d.gro",teller
);
75 fprintf(fp
,"Dumping fitted helix frame %d\n",teller
);
76 fprintf(fp
,"%5d\n",nres
*5);
77 for(i
=0; (i
<nres
); i
++) {
78 #define PR(AA) fprintf(fp,"%5d%5s%5s%5d%8.3f%8.3f%8.3f\n",i+1,"GLY",#AA,bb[i].AA,x[bb[i].AA][XX],x[bb[i].AA][YY],x[bb[i].AA][ZZ]); fflush(fp)
87 for(i
=0; (i
<DIM
); i
++)
88 fprintf(fp
,"%10.5f",box
[i
][i
]);
93 void dump_otrj(FILE *otrj
,int natoms
,atom_id all_index
[],rvec x
[],
103 fp
=ffopen("WEDGAMMA10.DAT","w");
108 fprintf(fp
,"%10g\n",fac
);
110 for(i
=0; (i
<natoms
); i
++) {
112 for(m
=0; (m
<DIM
); m
++) {
119 write_gms_ndx(otrj
,natoms
,all_index
,x
,NULL
);
122 int gmx_helix(int argc
,char *argv
[])
124 const char *desc
[] = {
125 "[TT]g_helix[tt] computes all kinds of helix properties. First, the peptide",
126 "is checked to find the longest helical part, as determined by",
127 "hydrogen bonds and [GRK]phi[grk]/[GRK]psi[grk] angles.",
128 "That bit is fitted",
129 "to an ideal helix around the [IT]z[it]-axis and centered around the origin.",
130 "Then the following properties are computed:[PAR]",
131 "[BB]1.[bb] Helix radius (file [TT]radius.xvg[tt]). This is merely the",
132 "RMS deviation in two dimensions for all C[GRK]alpha[grk] atoms.",
133 "it is calced as sqrt((SUM i(x^2(i)+y^2(i)))/N), where N is the number",
134 "of backbone atoms. For an ideal helix the radius is 0.23 nm[BR]",
135 "[BB]2.[bb] Twist (file [TT]twist.xvg[tt]). The average helical angle per",
136 "residue is calculated. For an [GRK]alpha[grk]-helix it is 100 degrees,",
137 "for 3-10 helices it will be smaller, and ",
138 "for 5-helices it will be larger.[BR]",
139 "[BB]3.[bb] Rise per residue (file [TT]rise.xvg[tt]). The helical rise per",
140 "residue is plotted as the difference in [IT]z[it]-coordinate between C[GRK]alpha[grk]",
141 "atoms. For an ideal helix, this is 0.15 nm[BR]",
142 "[BB]4.[bb] Total helix length (file [TT]len-ahx.xvg[tt]). The total length",
144 "helix in nm. This is simply the average rise (see above) times the",
145 "number of helical residues (see below).[BR]",
146 "[BB]5.[bb] Number of helical residues (file [TT]n-ahx.xvg[tt]). The title says",
148 "[BB]6.[bb] Helix dipole, backbone only (file [TT]dip-ahx.xvg[tt]).[BR]",
149 "[BB]7.[bb] RMS deviation from ideal helix, calculated for the C[GRK]alpha[grk]",
150 "atoms only (file [TT]rms-ahx.xvg[tt]).[BR]",
151 "[BB]8.[bb] Average C[GRK]alpha[grk] - C[GRK]alpha[grk] dihedral angle (file [TT]phi-ahx.xvg[tt]).[BR]",
152 "[BB]9.[bb] Average [GRK]phi[grk] and [GRK]psi[grk] angles (file [TT]phipsi.xvg[tt]).[BR]",
153 "[BB]10.[bb] Ellipticity at 222 nm according to Hirst and Brooks.",
156 static const char *ppp
[efhNR
+2] = {
157 NULL
, "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
158 "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5", "CD222", NULL
160 static gmx_bool bCheck
=FALSE
,bFit
=TRUE
,bDBG
=FALSE
,bEV
=FALSE
;
161 static int rStart
=0,rEnd
=0,r0
=1;
163 { "-r0", FALSE
, etINT
, {&r0
},
164 "The first residue number in the sequence" },
165 { "-q", FALSE
, etBOOL
,{&bCheck
},
166 "Check at every step which part of the sequence is helical" },
167 { "-F", FALSE
, etBOOL
,{&bFit
},
168 "Toggle fit to a perfect helix" },
169 { "-db", FALSE
, etBOOL
,{&bDBG
},
170 "Print debug info" },
171 { "-prop", FALSE
, etENUM
, {ppp
},
172 "Select property to weight eigenvectors with. WARNING experimental stuff" },
173 { "-ev", FALSE
, etBOOL
,{&bEV
},
174 "Write a new 'trajectory' file for ED" },
175 { "-ahxstart", FALSE
, etINT
, {&rStart
},
176 "First residue in helix" },
177 { "-ahxend", FALSE
, etINT
, {&rEnd
},
178 "Last residue in helix" }
191 t_xvgrfile xf
[efhNR
] = {
192 { NULL
, NULL
, TRUE
, "radius", "Helix radius", NULL
, "r (nm)" , 0.0 },
193 { NULL
, NULL
, TRUE
, "twist", "Twist per residue", NULL
, "Angle (deg)", 0.0 },
194 { NULL
, NULL
, TRUE
, "rise", "Rise per residue", NULL
, "Rise (nm)", 0.0 },
195 { NULL
, NULL
, FALSE
, "len-ahx", "Length of the Helix", NULL
, "Length (nm)", 0.0 },
196 { NULL
, NULL
, FALSE
, "dip-ahx", "Helix Backbone Dipole", NULL
, "rq (nm e)", 0.0 },
197 { NULL
, NULL
, TRUE
, "rms-ahx", "RMS Deviation from Ideal Helix", NULL
, "RMS (nm)", 0.0 },
198 { NULL
, NULL
, FALSE
, "rmsa-ahx","Average RMSD per Residue", "Residue", "RMS (nm)", 0.0 },
199 { NULL
, NULL
,FALSE
, "cd222", "Ellipticity at 222 nm", NULL
, "nm", 0.0 },
200 { NULL
, NULL
, TRUE
, "pprms", "RMS Distance from \\8a\\4-helix", NULL
, "deg" , 0.0 },
201 { NULL
, NULL
, TRUE
, "caphi", "Average Ca-Ca Dihedral", NULL
, "\\8F\\4(deg)", 0.0 },
202 { NULL
, NULL
, TRUE
, "phi", "Average \\8F\\4 angles", NULL
, "deg" , 0.0 },
203 { NULL
, NULL
, TRUE
, "psi", "Average \\8Y\\4 angles", NULL
, "deg" , 0.0 },
204 { NULL
, NULL
, TRUE
, "hb3", "Average n-n+3 hbond length", NULL
, "nm" , 0.0 },
205 { NULL
, NULL
, TRUE
, "hb4", "Average n-n+4 hbond length", NULL
, "nm" , 0.0 },
206 { NULL
, NULL
, TRUE
, "hb5", "Average n-n+5 hbond length", NULL
, "nm" , 0.0 },
207 { NULL
, NULL
,FALSE
, "JCaHa", "J-Coupling Values", "Residue", "Hz" , 0.0 },
208 { NULL
, NULL
,FALSE
, "helicity","Helicity per Residue", "Residue", "% of time" , 0.0 }
213 char buf
[54],prop
[256];
217 int i
,j
,ai
,m
,nall
,nbb
,nca
,teller
,nSel
=0;
218 atom_id
*bbindex
,*caindex
,*allindex
;
225 gmx_rmpbc_t gpbc
=NULL
;
228 { efTPX
, NULL
, NULL
, ffREAD
},
229 { efNDX
, NULL
, NULL
, ffREAD
},
230 { efTRX
, "-f", NULL
, ffREAD
},
231 { efG87
, "-to", NULL
, ffOPTWR
},
232 { efSTO
, "-cz", "zconf",ffWRITE
},
233 { efSTO
, "-co", "waver",ffWRITE
}
235 #define NFILE asize(fnm)
237 CopyRight(stderr
,argv
[0]);
238 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
239 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
241 bRange
=(opt2parg_bSet("-ahxstart",asize(pa
),pa
) &&
242 opt2parg_bSet("-ahxend",asize(pa
),pa
));
244 top
=read_top(ftp2fn(efTPX
,NFILE
,fnm
),&ePBC
);
246 natoms
=read_first_x(oenv
,&status
,opt2fn("-f",NFILE
,fnm
),&t
,&x
,box
);
248 if (opt2bSet("-to",NFILE
,fnm
)) {
249 otrj
=opt2FILE("-to",NFILE
,fnm
,"w");
251 fprintf(otrj
,"%s Weighted Trajectory: %d atoms, NO box\n",prop
,natoms
);
256 if (natoms
!= top
->atoms
.nr
)
257 gmx_fatal(FARGS
,"Sorry can only run when the number of atoms in the run input file (%d) is equal to the number in the trajectory (%d)",
258 top
->atoms
.nr
,natoms
);
260 bb
=mkbbind(ftp2fn(efNDX
,NFILE
,fnm
),&nres
,&nbb
,r0
,&nall
,&allindex
,
261 top
->atoms
.atomname
,top
->atoms
.atom
,top
->atoms
.resinfo
);
262 snew(bbindex
,natoms
);
265 fprintf(stderr
,"nall=%d\n",nall
);
267 /* Open output files, default x-axis is time */
268 for(i
=0; (i
<efhNR
); i
++) {
269 sprintf(buf
,"%s.xvg",xf
[i
].filenm
);
271 xf
[i
].fp
=xvgropen(buf
,xf
[i
].title
,
272 xf
[i
].xaxis
? xf
[i
].xaxis
: "Time (ps)",
275 sprintf(buf
,"%s.out",xf
[i
].filenm
);
277 xf
[i
].fp2
=ffopen(buf
,"w");
281 /* Read reference frame from tpx file to compute helix length */
282 snew(xref
,top
->atoms
.nr
);
283 read_tpx(ftp2fn(efTPX
,NFILE
,fnm
),
284 NULL
,NULL
,&natoms
,xref
,NULL
,NULL
,NULL
);
285 calc_hxprops(nres
,bb
,xref
,box
);
286 do_start_end(nres
,bb
,xref
,&nbb
,bbindex
,&nca
,caindex
,bRange
,rStart
,rEnd
);
289 fprintf(stderr
,"nca=%d, nbb=%d\n",nca
,nbb
);
290 pr_bb(stdout
,nres
,bb
);
293 gpbc
= gmx_rmpbc_init(&top
->idef
,ePBC
,natoms
,box
);
298 if ((teller
++ % 10) == 0)
299 fprintf(stderr
,"\rt=%.2f",t
);
300 gmx_rmpbc(gpbc
,natoms
,box
,x
);
303 calc_hxprops(nres
,bb
,x
,box
);
305 do_start_end(nres
,bb
,x
,&nbb
,bbindex
,&nca
,caindex
,FALSE
,0,0);
308 rms
=fit_ahx(nres
,bb
,natoms
,nall
,allindex
,x
,nca
,caindex
,box
,bFit
);
311 write_sto_conf(opt2fn("-cz",NFILE
,fnm
),"Helix fitted to Z-Axis",
312 &(top
->atoms
),x
,NULL
,ePBC
,box
);
315 xf
[efhRAD
].val
= radius(xf
[efhRAD
].fp2
,nca
,caindex
,x
);
316 xf
[efhTWIST
].val
= twist(xf
[efhTWIST
].fp2
,nca
,caindex
,x
);
317 xf
[efhRISE
].val
= rise(nca
,caindex
,x
);
318 xf
[efhLEN
].val
= ahx_len(nca
,caindex
,x
,box
);
319 xf
[efhCD222
].val
= ellipticity(nres
,bb
);
320 xf
[efhDIP
].val
= dip(nbb
,bbindex
,x
,top
->atoms
.atom
);
321 xf
[efhRMS
].val
= rms
;
322 xf
[efhCPHI
].val
= ca_phi(nca
,caindex
,x
,box
);
323 xf
[efhPPRMS
].val
= pprms(xf
[efhPPRMS
].fp2
,nres
,bb
);
325 for(j
=0; (j
<=efhCPHI
); j
++)
326 fprintf(xf
[j
].fp
, "%10g %10g\n",t
,xf
[j
].val
);
328 av_phipsi(xf
[efhPHI
].fp
,xf
[efhPSI
].fp
,xf
[efhPHI
].fp2
,xf
[efhPSI
].fp2
,
330 av_hblen(xf
[efhHB3
].fp
,xf
[efhHB3
].fp2
,
331 xf
[efhHB4
].fp
,xf
[efhHB4
].fp2
,
332 xf
[efhHB5
].fp
,xf
[efhHB5
].fp2
,
336 dump_otrj(otrj
,nall
,allindex
,x
,xf
[nSel
].val
,xav
);
338 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
339 fprintf(stderr
,"\n");
341 gmx_rmpbc_done(gpbc
);
348 for(i
=0; (i
<nall
); i
++) {
350 for(m
=0; (m
<DIM
); m
++)
353 write_sto_conf_indexed(opt2fn("-co",NFILE
,fnm
),
354 "Weighted and Averaged conformation",
355 &(top
->atoms
),xav
,NULL
,ePBC
,box
,nall
,allindex
);
358 for(i
=0; (i
<nres
); i
++) {
359 if (bb
[i
].nrms
> 0) {
360 fprintf(xf
[efhRMSA
].fp
,"%10d %10g\n",r0
+i
,bb
[i
].rmsa
/bb
[i
].nrms
);
362 fprintf(xf
[efhAHX
].fp
,"%10d %10g\n",r0
+i
,(bb
[i
].nhx
*100.0)/(real
)teller
);
363 fprintf(xf
[efhJCA
].fp
,"%10d %10g\n",
364 r0
+i
,140.3+(bb
[i
].jcaha
/(double)teller
));
367 for(i
=0; (i
<efhNR
); i
++) {
371 do_view(oenv
,xf
[i
].filenm
,"-nxy");