4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
65 void dump_ahx(int nres
,
66 t_bb bb
[],rvec x
[],matrix box
,int teller
)
72 sprintf(buf
,"dump%d.gro",teller
);
74 fprintf(fp
,"Dumping fitted helix frame %d\n",teller
);
75 fprintf(fp
,"%5d\n",nres
*5);
76 for(i
=0; (i
<nres
); i
++) {
77 #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)
86 for(i
=0; (i
<DIM
); i
++)
87 fprintf(fp
,"%10.5f",box
[i
][i
]);
92 int get_prop(char prop
[])
94 static char *ppp
[efhNR
] = {
95 "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
96 "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5", "CD222"
100 fprintf(stderr
,"Select something\n");
101 for(i
=0; (i
<efhNR
); ) {
102 fprintf(stderr
,"%2d: %10s",i
,ppp
[i
]); i
++;
104 fprintf(stderr
," %2d: %10s\n",i
,ppp
[i
]);
108 fprintf(stderr
,"\n");
112 } while ((n
< 0) || (n
>= efhNR
));
120 void dump_otrj(FILE *otrj
,int natoms
,atom_id all_index
[],rvec x
[],
123 static FILE *fp
=NULL
;
124 static real fac0
=1.0;
130 fp
=ffopen("WEDGAMMA10.DAT","w");
135 fprintf(fp
,"%10g\n",fac
);
137 for(i
=0; (i
<natoms
); i
++) {
139 for(m
=0; (m
<DIM
); m
++) {
146 write_gms_ndx(otrj
,natoms
,all_index
,x
,NULL
);
149 int gmx_helix(int argc
,char *argv
[])
151 static char *desc
[] = {
152 "g_helix computes all kind of helix properties. First, the peptide",
153 "is checked to find the longest helical part. This is determined by",
154 "Hydrogen bonds and Phi/Psi angles.",
155 "That bit is fitted",
156 "to an ideal helix around the Z-axis and centered around the origin.",
157 "Then the following properties are computed:[PAR]",
158 "[BB]1.[bb] Helix radius (file radius.xvg). This is merely the",
159 "RMS deviation in two dimensions for all Calpha atoms.",
160 "it is calced as sqrt((SUM i(x^2(i)+y^2(i)))/N), where N is the number",
161 "of backbone atoms. For an ideal helix the radius is 0.23 nm[BR]",
162 "[BB]2.[bb] Twist (file twist.xvg). The average helical angle per",
163 "residue is calculated. For alpha helix it is 100 degrees,",
164 "for 3-10 helices it will be smaller,",
165 "for 5-helices it will be larger.[BR]",
166 "[BB]3.[bb] Rise per residue (file rise.xvg). The helical rise per",
167 "residue is plotted as the difference in Z-coordinate between Ca",
168 "atoms. For an ideal helix this is 0.15 nm[BR]",
169 "[BB]4.[bb] Total helix length (file len-ahx.xvg). The total length",
171 "helix in nm. This is simply the average rise (see above) times the",
172 "number of helical residues (see below).[BR]",
173 "[BB]5.[bb] Number of helical residues (file n-ahx.xvg). The title says",
175 "[BB]6.[bb] Helix Dipole, backbone only (file dip-ahx.xvg).[BR]",
176 "[BB]7.[bb] RMS deviation from ideal helix, calculated for the Calpha",
177 "atoms only (file rms-ahx.xvg).[BR]",
178 "[BB]8.[bb] Average Calpha-Calpha dihedral angle (file phi-ahx.xvg).[BR]",
179 "[BB]9.[bb] Average Phi and Psi angles (file phipsi.xvg).[BR]",
180 "[BB]10.[bb] Ellipticity at 222 nm according to [IT]Hirst and Brooks[it]",
183 static bool bCheck
=FALSE
,bFit
=TRUE
,bDBG
=FALSE
,bEV
=FALSE
;
184 static int rStart
=0,rEnd
=0,r0
=1;
186 { "-r0", FALSE
, etINT
, {&r0
},
187 "The first residue number in the sequence" },
188 { "-q", FALSE
, etBOOL
,{&bCheck
},
189 "Check at every step which part of the sequence is helical" },
190 { "-F", FALSE
, etBOOL
,{&bFit
},
191 "Toggle fit to a perfect helix" },
192 { "-db", FALSE
, etBOOL
,{&bDBG
},
193 "Print debug info" },
194 { "-ev", FALSE
, etBOOL
,{&bEV
},
195 "Write a new 'trajectory' file for ED" },
196 { "-ahxstart", FALSE
, etINT
, {&rStart
},
197 "First residue in helix" },
198 { "-ahxend", FALSE
, etINT
, {&rEnd
},
199 "Last residue in helix" }
212 t_xvgrfile xf
[efhNR
] = {
213 { NULL
, NULL
, TRUE
, "radius", "Helix radius", NULL
, "r (nm)" , 0.0 },
214 { NULL
, NULL
, TRUE
, "twist", "Twist per residue", NULL
, "Angle (deg)", 0.0 },
215 { NULL
, NULL
, TRUE
, "rise", "Rise per residue", NULL
, "Rise (nm)", 0.0 },
216 { NULL
, NULL
, FALSE
, "len-ahx", "Length of the Helix", NULL
, "Length (nm)", 0.0 },
217 { NULL
, NULL
, FALSE
, "dip-ahx", "Helix Backbone Dipole", NULL
, "rq (nm e)", 0.0 },
218 { NULL
, NULL
, TRUE
, "rms-ahx", "RMS Deviation from Ideal Helix", NULL
, "RMS (nm)", 0.0 },
219 { NULL
, NULL
, FALSE
, "rmsa-ahx","Average RMSD per Residue", "Residue", "RMS (nm)", 0.0 },
220 { NULL
, NULL
,FALSE
, "cd222", "Ellipticity at 222 nm", NULL
, "nm", 0.0 },
221 { NULL
, NULL
, TRUE
, "pprms", "RMS Distance from \\8a\\4-helix", NULL
, "deg" , 0.0 },
222 { NULL
, NULL
, TRUE
, "caphi", "Average Ca-Ca Dihedral", NULL
, "\\8F\\4(deg)", 0.0 },
223 { NULL
, NULL
, TRUE
, "phi", "Average \\8F\\4 angles", NULL
, "deg" , 0.0 },
224 { NULL
, NULL
, TRUE
, "psi", "Average \\8Y\\4 angles", NULL
, "deg" , 0.0 },
225 { NULL
, NULL
, TRUE
, "hb3", "Average n-n+3 hbond length", NULL
, "nm" , 0.0 },
226 { NULL
, NULL
, TRUE
, "hb4", "Average n-n+4 hbond length", NULL
, "nm" , 0.0 },
227 { NULL
, NULL
, TRUE
, "hb5", "Average n-n+5 hbond length", NULL
, "nm" , 0.0 },
228 { NULL
, NULL
,FALSE
, "JCaHa", "J-Coupling Values", "Residue", "Hz" , 0.0 },
229 { NULL
, NULL
,FALSE
, "helicity","Helicity per Residue", "Residue", "% of time" , 0.0 }
233 char buf
[54],prop
[256];
235 int natoms
,nre
,nres
,step
;
237 int i
,j
,ai
,m
,nall
,nbb
,nca
,teller
,nSel
=0;
238 atom_id
*bbindex
,*caindex
,*allindex
;
246 { efTPX
, NULL
, NULL
, ffREAD
},
247 { efNDX
, NULL
, NULL
, ffREAD
},
248 { efTRX
, "-f", NULL
, ffREAD
},
249 { efG87
, "-to", NULL
, ffOPTWR
},
250 { efSTO
, "-cz", "zconf",ffWRITE
},
251 { efSTO
, "-co", "waver",ffWRITE
}
253 #define NFILE asize(fnm)
255 CopyRight(stderr
,argv
[0]);
256 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
257 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
);
259 bRange
=(opt2parg_bSet("-ahxstart",asize(pa
),pa
) &&
260 opt2parg_bSet("-ahxend",asize(pa
),pa
));
262 top
=read_top(ftp2fn(efTPX
,NFILE
,fnm
));
264 natoms
=read_first_x(&status
,opt2fn("-f",NFILE
,fnm
),&t
,&x
,box
);
266 if (opt2bSet("-to",NFILE
,fnm
)) {
267 otrj
=opt2FILE("-to",NFILE
,fnm
,"w");
269 fprintf(otrj
,"%s Weighted Trajectory: %d atoms, NO box\n",prop
,natoms
);
274 if (natoms
!= top
->atoms
.nr
)
275 fatal_error(0,"Sorry can only run when the number of atoms in the run input file (%d) is equal to the number in the trajectory (%d)",
276 top
->atoms
.nr
,natoms
);
278 bb
=mkbbind(ftp2fn(efNDX
,NFILE
,fnm
),&nres
,&nbb
,r0
,&nall
,&allindex
,
279 top
->atoms
.atomname
,top
->atoms
.atom
,top
->atoms
.resname
);
280 snew(bbindex
,natoms
);
283 fprintf(stderr
,"nall=%d\n",nall
);
285 /* Open output files, default x-axis is time */
286 for(i
=0; (i
<efhNR
); i
++) {
287 sprintf(buf
,"%s.xvg",xf
[i
].filenm
);
289 xf
[i
].fp
=xvgropen(buf
,xf
[i
].title
,xf
[i
].xaxis
? xf
[i
].xaxis
: "Time (ps)",
292 sprintf(buf
,"%s.out",xf
[i
].filenm
);
294 xf
[i
].fp2
=ffopen(buf
,"w");
298 /* Read reference frame from tpx file to compute helix length */
299 snew(xref
,top
->atoms
.nr
);
300 read_tpx(ftp2fn(efTPX
,NFILE
,fnm
),
301 &step
,&t
,&lambda
,NULL
,NULL
,
302 &natoms
,xref
,NULL
,NULL
,NULL
);
303 calc_hxprops(nres
,bb
,xref
,box
);
304 do_start_end(nres
,bb
,xref
,&nbb
,bbindex
,&nca
,caindex
,bRange
,rStart
,rEnd
);
307 fprintf(stderr
,"nca=%d, nbb=%d\n",nca
,nbb
);
308 pr_bb(stdout
,nres
,bb
);
314 if ((teller
++ % 10) == 0)
315 fprintf(stderr
,"\rt=%.2f",t
);
316 rm_pbc(&(top
->idef
),top
->atoms
.nr
,box
,x
,x
);
318 calc_hxprops(nres
,bb
,x
,box
);
320 do_start_end(nres
,bb
,x
,&nbb
,bbindex
,&nca
,caindex
,FALSE
,0,0);
323 rms
=fit_ahx(nres
,bb
,natoms
,nall
,allindex
,x
,nca
,caindex
,box
,bFit
);
326 write_sto_conf(opt2fn("-cz",NFILE
,fnm
),"Helix fitted to Z-Axis",
327 &(top
->atoms
),x
,NULL
,box
);
330 xf
[efhRAD
].val
= radius(xf
[efhRAD
].fp2
,nca
,caindex
,x
);
331 xf
[efhTWIST
].val
= twist(xf
[efhTWIST
].fp2
,nca
,caindex
,x
);
332 xf
[efhRISE
].val
= rise(nca
,caindex
,x
);
333 xf
[efhLEN
].val
= ahx_len(nca
,caindex
,x
,box
);
334 xf
[efhCD222
].val
= ellipticity(nres
,bb
);
335 xf
[efhDIP
].val
= dip(nbb
,bbindex
,x
,top
->atoms
.atom
);
336 xf
[efhRMS
].val
= rms
;
337 xf
[efhCPHI
].val
= ca_phi(nca
,caindex
,x
,box
);
338 xf
[efhPPRMS
].val
= pprms(xf
[efhPPRMS
].fp2
,nres
,bb
);
340 for(j
=0; (j
<=efhCPHI
); j
++)
341 fprintf(xf
[j
].fp
, "%10g %10g\n",t
,xf
[j
].val
);
343 av_phipsi(xf
[efhPHI
].fp
,xf
[efhPSI
].fp
,xf
[efhPHI
].fp2
,xf
[efhPSI
].fp2
,
345 av_hblen(xf
[efhHB3
].fp
,xf
[efhHB3
].fp2
,
346 xf
[efhHB4
].fp
,xf
[efhHB4
].fp2
,
347 xf
[efhHB5
].fp
,xf
[efhHB5
].fp2
,
351 dump_otrj(otrj
,nall
,allindex
,x
,xf
[nSel
].val
,xav
);
353 } while (read_next_x(status
,&t
,natoms
,x
,box
));
354 fprintf(stderr
,"\n");
361 for(i
=0; (i
<nall
); i
++) {
363 for(m
=0; (m
<DIM
); m
++)
366 write_sto_conf_indexed(opt2fn("-co",NFILE
,fnm
),
367 "Weighted and Averaged conformation",
368 &(top
->atoms
),xav
,NULL
,box
,nall
,allindex
);
371 for(i
=0; (i
<nres
); i
++) {
372 if (bb
[i
].nrms
> 0) {
373 fprintf(xf
[efhRMSA
].fp
,"%10d %10g\n",r0
+i
,bb
[i
].rmsa
/bb
[i
].nrms
);
375 fprintf(xf
[efhAHX
].fp
,"%10d %10g\n",r0
+i
,(bb
[i
].nhx
*100.0)/(real
)teller
);
376 fprintf(xf
[efhJCA
].fp
,"%10d %10g\n",
377 r0
+i
,140.3+(bb
[i
].jcaha
/(double)teller
));
380 for(i
=0; (i
<efhNR
); i
++) {
384 do_view(xf
[i
].filenm
,"-nxy");