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
51 real ellipticity(int nres, t_bb bb[])
57 static const t_ppwstr ppw[] = {
64 { -70.5, -35.8, 0.15 },
71 #define NPPW asize(ppw)
75 real ell, pp2, phi, psi;
78 for (i = 0; (i < nres); i++)
82 for (j = 0; (j < NPPW); j++)
84 pp2 = sqr(phi-ppw[j].phi)+sqr(psi-ppw[j].psi);
96 real ahx_len(int gnx, atom_id index[], rvec x[])
97 /* Assume we have a list of Calpha atoms only! */
101 rvec_sub(x[index[0]], x[index[gnx-1]], dx);
106 real radius(FILE *fp, int nca, atom_id ca_index[], rvec x[])
107 /* Assume we have all the backbone */
113 for (i = 0; (i < nca); i++)
116 dl2 = sqr(x[ai][XX])+sqr(x[ai][YY]);
120 fprintf(fp, " %10g", dl2);
130 return sqrt(dlt/nca);
133 static real rot(rvec x1, rvec x2)
135 real phi1, dphi, cp, sp;
138 phi1 = atan2(x1[YY], x1[XX]);
141 xx = cp*x2[XX]+sp*x2[YY];
142 yy = -sp*x2[XX]+cp*x2[YY];
144 dphi = RAD2DEG*atan2(yy, xx);
149 real twist(int nca, atom_id caindex[], rvec x[])
156 for (i = 1; (i < nca); i++)
160 dphi = rot(x[a0], x[a1]);
172 real ca_phi(int gnx, atom_id index[], rvec x[])
173 /* Assume we have a list of Calpha atoms only! */
176 int i, ai, aj, ak, al, t1, t2, t3;
177 rvec r_ij, r_kj, r_kl, m, n;
186 for (i = 0; (i < gnx-4); i++)
193 dih_angle(x[ai], x[aj], x[ak], x[al], NULL,
194 r_ij, r_kj, r_kl, m, n,
195 &sign, &t1, &t2, &t3);
199 return (phitot/(gnx-4.0));
202 real dip(int nbb, atom_id bbind[], rvec x[], t_atom atom[])
209 for (i = 0; (i < nbb); i++)
213 for (m = 0; (m < DIM); m++)
215 dipje[m] += x[ai][m]*q;
221 real rise(int gnx, atom_id index[], rvec x[])
222 /* Assume we have a list of Calpha atoms only! */
230 for (i = 1; (i < gnx); i++)
238 return (ztot/(gnx-1.0));
241 void av_hblen(FILE *fp3, FILE *fp3a,
242 FILE *fp4, FILE *fp4a,
243 FILE *fp5, FILE *fp5a,
244 real t, int nres, t_bb bb[])
246 int i, n3 = 0, n4 = 0, n5 = 0;
247 real d3 = 0, d4 = 0, d5 = 0;
249 for (i = 0; (i < nres-3); i++)
253 fprintf(fp3a, "%10g", bb[i].d3);
258 fprintf(fp4a, "%10g", bb[i].d4);
264 fprintf(fp5a, "%10g", bb[i].d5);
270 fprintf(fp3, "%10g %10g\n", t, d3/n3);
271 fprintf(fp4, "%10g %10g\n", t, d4/n4);
272 fprintf(fp5, "%10g %10g\n", t, d5/n5);
279 void av_phipsi(FILE *fphi, FILE *fpsi, FILE *fphi2, FILE *fpsi2,
280 real t, int nres, t_bb bb[])
283 real phi = 0, psi = 0;
285 fprintf(fphi2, "%10g", t);
286 fprintf(fpsi2, "%10g", t);
287 for (i = 0; (i < nres); i++)
293 fprintf(fphi2, " %10g", bb[i].phi);
294 fprintf(fpsi2, " %10g", bb[i].psi);
298 fprintf(fphi, "%10g %10g\n", t, (phi/n));
299 fprintf(fpsi, "%10g %10g\n", t, (psi/n));
300 fprintf(fphi2, "\n");
301 fprintf(fpsi2, "\n");
304 static void set_ahcity(int nbb, t_bb bb[])
309 for (n = 0; (n < nbb); n++)
311 pp2 = sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
313 bb[n].bHelix = FALSE;
316 if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
324 t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0,
325 int *nall, atom_id **index,
326 char ***atomname, t_atom atom[],
329 static const char * bb_nm[] = { "N", "H", "CA", "C", "O", "HN" };
330 #define NBB asize(bb_nm)
333 int ai, i, i0, i1, j, k, ri, rnr, gnx, r0, r1;
335 fprintf(stderr, "Please select a group containing the entire backbone\n");
336 rd_index(fn, 1, &gnx, index, &grpname);
338 fprintf(stderr, "Checking group %s\n", grpname);
339 r0 = r1 = atom[(*index)[0]].resind;
340 for (i = 1; (i < gnx); i++)
342 r0 = min(r0, atom[(*index)[i]].resind);
343 r1 = max(r1, atom[(*index)[i]].resind);
346 fprintf(stderr, "There are %d residues\n", rnr);
348 for (i = 0; (i < rnr); i++)
350 bb[i].N = bb[i].H = bb[i].CA = bb[i].C = bb[i].O = -1, bb[i].resno = res0+i;
353 for (i = j = 0; (i < gnx); i++)
356 ri = atom[ai].resind-r0;
357 if (strcmp(*(resinfo[ri].name), "PRO") == 0)
359 if (strcmp(*(atomname[ai]), "CD") == 0)
364 for (k = 0; (k < NBB); k++)
366 if (strcmp(bb_nm[k], *(atomname[ai])) == 0)
379 /* No attempt to address the case where some weird input has both H and HN atoms in the group */
396 for (i0 = 0; (i0 < rnr); i0++)
398 if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
400 (bb[i0].C != -1) && (bb[i0].O != -1))
405 for (i1 = rnr-1; (i1 >= 0); i1--)
407 if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
409 (bb[i1].C != -1) && (bb[i1].O != -1))
423 for (i = i0; (i < i1); i++)
425 bb[i].Cprev = bb[i-1].C;
426 bb[i].Nnext = bb[i+1].N;
428 rnr = max(0, i1-i0+1);
429 fprintf(stderr, "There are %d complete backbone residues (from %d to %d)\n",
430 rnr, bb[i0].resno, bb[i1].resno);
433 gmx_fatal(FARGS, "Zero complete backbone residues were found, cannot proceed");
435 for (i = 0; (i < rnr); i++, i0++)
441 for (i = 0; (i < rnr); i++)
443 ri = atom[bb[i].CA].resind;
444 sprintf(bb[i].label, "%s%d", *(resinfo[ri].name), ri+res0);
448 *nbb = rnr*asize(bb_nm);
453 real pprms(FILE *fp, int nbb, t_bb bb[])
456 real rms, rmst, rms2;
459 for (i = n = 0; (i < nbb); i++)
463 rms = sqrt(bb[i].pprms2);
465 rms2 += bb[i].pprms2;
466 fprintf(fp, "%10g ", rms);
471 rms = sqrt(rms2/n-sqr(rmst/n));
476 void calc_hxprops(int nres, t_bb bb[], rvec x[])
478 int i, ao, an, t1, t2, t3;
479 rvec dx, r_ij, r_kj, r_kl, m, n;
482 for (i = 0; (i < nres); i++)
485 bb[i].d4 = bb[i].d3 = bb[i].d5 = 0;
489 rvec_sub(x[ao], x[an], dx);
495 rvec_sub(x[ao], x[an], dx);
501 rvec_sub(x[ao], x[an], dx);
506 dih_angle(x[bb[i].Cprev], x[bb[i].N], x[bb[i].CA], x[bb[i].C], NULL,
507 r_ij, r_kj, r_kl, m, n,
508 &sign, &t1, &t2, &t3);
510 dih_angle(x[bb[i].N], x[bb[i].CA], x[bb[i].C], x[bb[i].Nnext], NULL,
511 r_ij, r_kj, r_kl, m, n,
512 &sign, &t1, &t2, &t3);
513 bb[i].pprms2 = sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
516 1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
517 4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
518 2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
522 static void check_ahx(int nres, t_bb bb[],
523 int *hstart, int *hend)
525 int h0, h1, h0sav, h1sav;
527 set_ahcity(nres, bb);
528 h0 = h0sav = h1sav = 0;
531 for (; (!bb[h0].bHelix) && (h0 < nres-4); h0++)
535 for (h1 = h0; bb[h1+1].bHelix && (h1 < nres-1); h1++)
541 /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
542 if (h1-h0 > h1sav-h0sav)
555 void do_start_end(int nres, t_bb bb[], int *nbb, atom_id bbindex[],
556 int *nca, atom_id caindex[],
557 gmx_bool bRange, int rStart, int rEnd)
559 int i, j, hstart = 0, hend = 0;
563 for (i = 0; (i < nres); i++)
565 if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
569 if (bb[i].resno == rStart)
573 if (bb[i].resno == rEnd)
581 /* Find start and end of longest helix fragment */
582 check_ahx(nres, bb, &hstart, &hend);
584 fprintf(stderr, "helix from: %d through %d\n",
585 bb[hstart].resno, bb[hend].resno);
587 for (j = 0, i = hstart; (i <= hend); i++)
589 bbindex[j++] = bb[i].N;
590 bbindex[j++] = bb[i].H;
591 bbindex[j++] = bb[i].CA;
592 bbindex[j++] = bb[i].C;
593 bbindex[j++] = bb[i].O;
594 caindex[i-hstart] = bb[i].CA;
597 *nca = (hend-hstart+1);
600 void pr_bb(FILE *fp, int nres, t_bb bb[])
605 fprintf(fp, "%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
606 "AA", "N", "Ca", "C", "O", "Phi", "Psi", "D3", "D4", "D5", "Hx?");
607 for (i = 0; (i < nres); i++)
609 fprintf(fp, "%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
610 bb[i].resno, bb[i].N, bb[i].CA, bb[i].C, bb[i].O,
611 bb[i].phi, bb[i].psi, bb[i].d3, bb[i].d4, bb[i].d5,
612 bb[i].bHelix ? "Yes" : "No");