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 "g_helix computes all kind of helix properties. First, the peptide",
126 "is checked to find the longest helical part. This is determined by",
127 "Hydrogen bonds and Phi/Psi angles.",
128 "That bit is fitted",
129 "to an ideal helix around the Z-axis and centered around the origin.",
130 "Then the following properties are computed:[PAR]",
131 "[BB]1.[bb] Helix radius (file radius.xvg). This is merely the",
132 "RMS deviation in two dimensions for all Calpha 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 twist.xvg). The average helical angle per",
136 "residue is calculated. For alpha helix it is 100 degrees,",
137 "for 3-10 helices it will be smaller,",
138 "for 5-helices it will be larger.[BR]",
139 "[BB]3.[bb] Rise per residue (file rise.xvg). The helical rise per",
140 "residue is plotted as the difference in Z-coordinate between Ca",
141 "atoms. For an ideal helix this is 0.15 nm[BR]",
142 "[BB]4.[bb] Total helix length (file len-ahx.xvg). 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 n-ahx.xvg). The title says",
148 "[BB]6.[bb] Helix Dipole, backbone only (file dip-ahx.xvg).[BR]",
149 "[BB]7.[bb] RMS deviation from ideal helix, calculated for the Calpha",
150 "atoms only (file rms-ahx.xvg).[BR]",
151 "[BB]8.[bb] Average Calpha-Calpha dihedral angle (file phi-ahx.xvg).[BR]",
152 "[BB]9.[bb] Average Phi and Psi angles (file phipsi.xvg).[BR]",
153 "[BB]10.[bb] Ellipticity at 222 nm according to [IT]Hirst and Brooks[it]",
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 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
;
227 { efTPX
, NULL
, NULL
, ffREAD
},
228 { efNDX
, NULL
, NULL
, ffREAD
},
229 { efTRX
, "-f", NULL
, ffREAD
},
230 { efG87
, "-to", NULL
, ffOPTWR
},
231 { efSTO
, "-cz", "zconf",ffWRITE
},
232 { efSTO
, "-co", "waver",ffWRITE
}
234 #define NFILE asize(fnm)
236 CopyRight(stderr
,argv
[0]);
237 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
238 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
240 bRange
=(opt2parg_bSet("-ahxstart",asize(pa
),pa
) &&
241 opt2parg_bSet("-ahxend",asize(pa
),pa
));
243 top
=read_top(ftp2fn(efTPX
,NFILE
,fnm
),&ePBC
);
245 natoms
=read_first_x(oenv
,&status
,opt2fn("-f",NFILE
,fnm
),&t
,&x
,box
);
247 if (opt2bSet("-to",NFILE
,fnm
)) {
248 otrj
=opt2FILE("-to",NFILE
,fnm
,"w");
250 fprintf(otrj
,"%s Weighted Trajectory: %d atoms, NO box\n",prop
,natoms
);
255 if (natoms
!= top
->atoms
.nr
)
256 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)",
257 top
->atoms
.nr
,natoms
);
259 bb
=mkbbind(ftp2fn(efNDX
,NFILE
,fnm
),&nres
,&nbb
,r0
,&nall
,&allindex
,
260 top
->atoms
.atomname
,top
->atoms
.atom
,top
->atoms
.resinfo
);
261 snew(bbindex
,natoms
);
264 fprintf(stderr
,"nall=%d\n",nall
);
266 /* Open output files, default x-axis is time */
267 for(i
=0; (i
<efhNR
); i
++) {
268 sprintf(buf
,"%s.xvg",xf
[i
].filenm
);
270 xf
[i
].fp
=xvgropen(buf
,xf
[i
].title
,
271 xf
[i
].xaxis
? xf
[i
].xaxis
: "Time (ps)",
274 sprintf(buf
,"%s.out",xf
[i
].filenm
);
276 xf
[i
].fp2
=ffopen(buf
,"w");
280 /* Read reference frame from tpx file to compute helix length */
281 snew(xref
,top
->atoms
.nr
);
282 read_tpx(ftp2fn(efTPX
,NFILE
,fnm
),
283 NULL
,NULL
,&natoms
,xref
,NULL
,NULL
,NULL
);
284 calc_hxprops(nres
,bb
,xref
,box
);
285 do_start_end(nres
,bb
,xref
,&nbb
,bbindex
,&nca
,caindex
,bRange
,rStart
,rEnd
);
288 fprintf(stderr
,"nca=%d, nbb=%d\n",nca
,nbb
);
289 pr_bb(stdout
,nres
,bb
);
295 if ((teller
++ % 10) == 0)
296 fprintf(stderr
,"\rt=%.2f",t
);
297 rm_pbc(&(top
->idef
),ePBC
,top
->atoms
.nr
,box
,x
,x
);
299 calc_hxprops(nres
,bb
,x
,box
);
301 do_start_end(nres
,bb
,x
,&nbb
,bbindex
,&nca
,caindex
,FALSE
,0,0);
304 rms
=fit_ahx(nres
,bb
,natoms
,nall
,allindex
,x
,nca
,caindex
,box
,bFit
);
307 write_sto_conf(opt2fn("-cz",NFILE
,fnm
),"Helix fitted to Z-Axis",
308 &(top
->atoms
),x
,NULL
,ePBC
,box
);
311 xf
[efhRAD
].val
= radius(xf
[efhRAD
].fp2
,nca
,caindex
,x
);
312 xf
[efhTWIST
].val
= twist(xf
[efhTWIST
].fp2
,nca
,caindex
,x
);
313 xf
[efhRISE
].val
= rise(nca
,caindex
,x
);
314 xf
[efhLEN
].val
= ahx_len(nca
,caindex
,x
,box
);
315 xf
[efhCD222
].val
= ellipticity(nres
,bb
);
316 xf
[efhDIP
].val
= dip(nbb
,bbindex
,x
,top
->atoms
.atom
);
317 xf
[efhRMS
].val
= rms
;
318 xf
[efhCPHI
].val
= ca_phi(nca
,caindex
,x
,box
);
319 xf
[efhPPRMS
].val
= pprms(xf
[efhPPRMS
].fp2
,nres
,bb
);
321 for(j
=0; (j
<=efhCPHI
); j
++)
322 fprintf(xf
[j
].fp
, "%10g %10g\n",t
,xf
[j
].val
);
324 av_phipsi(xf
[efhPHI
].fp
,xf
[efhPSI
].fp
,xf
[efhPHI
].fp2
,xf
[efhPSI
].fp2
,
326 av_hblen(xf
[efhHB3
].fp
,xf
[efhHB3
].fp2
,
327 xf
[efhHB4
].fp
,xf
[efhHB4
].fp2
,
328 xf
[efhHB5
].fp
,xf
[efhHB5
].fp2
,
332 dump_otrj(otrj
,nall
,allindex
,x
,xf
[nSel
].val
,xav
);
334 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
335 fprintf(stderr
,"\n");
342 for(i
=0; (i
<nall
); i
++) {
344 for(m
=0; (m
<DIM
); m
++)
347 write_sto_conf_indexed(opt2fn("-co",NFILE
,fnm
),
348 "Weighted and Averaged conformation",
349 &(top
->atoms
),xav
,NULL
,ePBC
,box
,nall
,allindex
);
352 for(i
=0; (i
<nres
); i
++) {
353 if (bb
[i
].nrms
> 0) {
354 fprintf(xf
[efhRMSA
].fp
,"%10d %10g\n",r0
+i
,bb
[i
].rmsa
/bb
[i
].nrms
);
356 fprintf(xf
[efhAHX
].fp
,"%10d %10g\n",r0
+i
,(bb
[i
].nhx
*100.0)/(real
)teller
);
357 fprintf(xf
[efhJCA
].fp
,"%10d %10g\n",
358 r0
+i
,140.3+(bb
[i
].jcaha
/(double)teller
));
361 for(i
=0; (i
<efhNR
); i
++) {
365 do_view(oenv
,xf
[i
].filenm
,"-nxy");