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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35 /* This file is completely threadsafe - keep it that way! */
46 #include "gromacs/fileio/confio.h"
51 #include "gpp_nextnb.h"
54 #include "gmx_fatal.h"
59 #define DIHEDRAL_WAS_SET_IN_RTP 0
60 static gmx_bool was_dihedral_set_in_rtp(t_param *dih)
62 return dih->c[MAXFORCEPARAM-1] == DIHEDRAL_WAS_SET_IN_RTP;
65 typedef gmx_bool (*peq)(t_param *p1, t_param *p2);
67 static int acomp(const void *a1, const void *a2)
74 if ((ac = (p1->AJ-p2->AJ)) != 0)
78 else if ((ac = (p1->AI-p2->AI)) != 0)
84 return (p1->AK-p2->AK);
88 static int pcomp(const void *a1, const void *a2)
95 if ((pc = (p1->AI-p2->AI)) != 0)
101 return (p1->AJ-p2->AJ);
105 static int dcomp(const void *d1, const void *d2)
112 /* First sort by J & K (the two central) atoms */
113 if ((dc = (p1->AJ-p2->AJ)) != 0)
117 else if ((dc = (p1->AK-p2->AK)) != 0)
121 /* Then make sure to put rtp dihedrals before generated ones */
122 else if (was_dihedral_set_in_rtp(p1) &&
123 !was_dihedral_set_in_rtp(p2))
127 else if (!was_dihedral_set_in_rtp(p1) &&
128 was_dihedral_set_in_rtp(p2))
132 /* Finally, sort by I and J (two outer) atoms */
133 else if ((dc = (p1->AI-p2->AI)) != 0)
139 return (p1->AL-p2->AL);
144 static gmx_bool is_dihedral_on_same_bond(t_param *p1, t_param *p2)
146 if (((p1->AJ == p2->AJ) && (p1->AK == p2->AK)) ||
147 ((p1->AJ == p2->AK) && (p1->AK == p2->AJ)))
158 static gmx_bool preq(t_param *p1, t_param *p2)
160 if ((p1->AI == p2->AI) && (p1->AJ == p2->AJ))
170 static void rm2par(t_param p[], int *np, peq eq)
183 for (i = 1; (i < (*np)); i++)
185 if (!eq(&p[i], &p[i-1]))
190 /* Index now holds pointers to all the non-equal params,
191 * this only works when p is sorted of course
193 for (i = 0; (i < nind); i++)
195 for (j = 0; (j < MAXATOMLIST); j++)
197 p[i].a[j] = p[index[i]].a[j];
199 for (j = 0; (j < MAXFORCEPARAM); j++)
201 p[i].c[j] = p[index[i]].c[j];
203 if (p[index[i]].a[0] == p[index[i]].a[1])
208 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
209 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
210 p[i].a[0], p[i].a[1], p[i].a[2], p[i].a[3]);
214 else if (index[i] > i)
216 /* Copy the string only if it comes from somewhere else
217 * otherwise we will end up copying a random (newly freed) pointer.
218 * Since the index is sorted we only have to test for index[i] > i.
220 strcpy(p[i].s, p[index[i]].s);
228 static void cppar(t_param p[], int np, t_params plist[], int ftype)
230 int i, j, nral, nrfp;
239 for (i = 0; (i < np); i++)
241 for (j = 0; (j < nral); j++)
243 ps->param[ps->nr].a[j] = p[i].a[j];
245 for (j = 0; (j < nrfp); j++)
247 ps->param[ps->nr].c[j] = p[i].c[j];
249 for (j = 0; (j < MAXSLEN); j++)
251 ps->param[ps->nr].s[j] = p[i].s[j];
257 static void cpparam(t_param *dest, t_param *src)
261 for (j = 0; (j < MAXATOMLIST); j++)
263 dest->a[j] = src->a[j];
265 for (j = 0; (j < MAXFORCEPARAM); j++)
267 dest->c[j] = src->c[j];
269 for (j = 0; (j < MAXSLEN); j++)
271 dest->s[j] = src->s[j];
275 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
279 for (j = 0; (j < 4); j++)
283 for (j = 0; (j < MAXFORCEPARAM); j++)
298 static int int_comp(const void *a, const void *b)
300 return (*(int *)a) - (*(int *)b);
303 static int atom_id_comp(const void *a, const void *b)
305 return (*(atom_id *)a) - (*(atom_id *)b);
308 static int eq_imp(atom_id a1[], atom_id a2[])
313 for (j = 0; (j < 4); j++)
318 qsort(b1, 4, (size_t)sizeof(b1[0]), int_comp);
319 qsort(b2, 4, (size_t)sizeof(b2[0]), int_comp);
321 for (j = 0; (j < 4); j++)
332 static int idcomp(const void *a, const void *b)
339 if ((d = (pa->a[0]-pb->a[0])) != 0)
343 else if ((d = (pa->a[3]-pb->a[3])) != 0)
347 else if ((d = (pa->a[1]-pb->a[1])) != 0)
353 return (int) (pa->a[2]-pb->a[2]);
357 static void sort_id(int nr, t_param ps[])
361 /* First swap order of atoms around if necessary */
362 for (i = 0; (i < nr); i++)
364 if (ps[i].a[3] < ps[i].a[0])
366 tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
367 tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
373 qsort(ps, nr, (size_t)sizeof(ps[0]), idcomp);
377 static int n_hydro(atom_id a[], char ***atomname)
382 for (i = 0; (i < 4); i += 3)
384 aname = *atomname[a[i]];
385 c0 = toupper(aname[0]);
390 else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9'))
392 c1 = toupper(aname[1]);
402 /* Clean up the dihedrals (both generated and read from the .rtp
404 static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper,
405 t_atoms *atoms, gmx_bool bKeepAllGeneratedDihedrals,
406 gmx_bool bRemoveDihedralIfWithImproper)
411 /* Construct the list of the indices of the dihedrals
412 * (i.e. generated or read) that might be kept. */
413 snew(index, *ndih+1);
414 if (bKeepAllGeneratedDihedrals)
416 fprintf(stderr, "Keeping all generated dihedrals\n");
418 for (i = 0; i < nind; i++)
427 /* Check if generated dihedral i should be removed. The
428 * dihedrals have been sorted by dcomp() above, so all those
429 * on the same two central atoms are together, with those from
430 * the .rtp file preceding those that were automatically
431 * generated. We remove the latter if the former exist. */
432 for (i = 0; i < *ndih; i++)
434 /* Keep the dihedrals that were defined in the .rtp file,
435 * and the dihedrals that were generated and different
436 * from the last one (whether it was generated or not). */
437 if (was_dihedral_set_in_rtp(&dih[i]) ||
439 !is_dihedral_on_same_bond(&dih[i], &dih[i-1]))
448 for (i = 0; i < nind; i++)
450 gmx_bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
451 gmx_bool bKeep = TRUE;
452 if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
454 /* Remove the dihedral if there is an improper on the same
456 for (j = 0; j < nimproper && bKeep; j++)
458 bKeep = !is_dihedral_on_same_bond(&dih[index[i]], &improper[j]);
464 /* If we don't want all dihedrals, we want to select the
465 * ones with the fewest hydrogens. Note that any generated
466 * dihedrals on the same bond as an .rtp dihedral may have
467 * been already pruned above in the construction of
468 * index[]. However, their parameters are still present,
469 * and l is looping over this dihedral and all of its
470 * pruned siblings. */
471 int bestl = index[i];
472 if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
474 /* Minimum number of hydrogens for i and l atoms */
478 is_dihedral_on_same_bond(&dih[index[i]], &dih[l]));
481 int nh = n_hydro(dih[l].a, atoms->atomname);
495 cpparam(&dih[k], &dih[bestl]);
501 for (i = k; i < *ndih; i++)
503 strcpy(dih[i].s, "");
510 static int get_impropers(t_atoms *atoms, t_hackblock hb[], t_param **improper,
511 gmx_bool bAllowMissing)
514 t_rbondeds *impropers;
515 t_rbonded *hbimproper;
516 int nimproper, i, j, k, r, start, ninc, nalloc;
517 atom_id ai[MAXATOMLIST];
522 snew(*improper, nalloc);
524 /* Add all the impropers from the residue database to the list. */
529 for (i = 0; (i < atoms->nres); i++)
531 impropers = &hb[i].rb[ebtsIDIHS];
532 for (j = 0; (j < impropers->nb); j++)
535 for (k = 0; (k < 4) && !bStop; k++)
537 ai[k] = search_atom(impropers->b[j].a[k], start,
539 "improper", bAllowMissing);
540 if (ai[k] == NO_ATID)
547 if (nimproper == nalloc)
550 srenew(*improper, nalloc);
553 set_p(&((*improper)[nimproper]), ai, NULL, impropers->b[j].s);
557 while ((start < atoms->nr) && (atoms->atom[start].resind == i))
567 static int nb_dist(t_nextnb *nnb, int ai, int aj)
579 nrexcl = nnb->nrexcl[ai];
580 for (nre = 1; (nre < nnb->nrex); nre++)
583 for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
585 if ((aj == a[nrx]) && (NRE == -1))
594 gmx_bool is_hydro(t_atoms *atoms, int ai)
596 return ((*(atoms->atomname[ai]))[0] == 'H');
599 static void get_atomnames_min(int n, char **anm,
600 int resind, t_atoms *atoms, atom_id *a)
604 /* Assume ascending residue numbering */
605 for (m = 0; m < n; m++)
607 if (atoms->atom[a[m]].resind < resind)
611 else if (atoms->atom[a[m]].resind > resind)
619 strcat(anm[m], *(atoms->atomname[a[m]]));
623 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
624 gmx_bool bAllowMissing)
627 atom_id a, astart, i1, i2, itmp;
633 for (a = 0; a < atoms->nr; a++)
635 r = atoms->atom[a].resind;
636 if (a == atoms->nr-1 || atoms->atom[a+1].resind != r)
638 hbexcl = &hb[r].rb[ebtsEXCLS];
640 for (e = 0; e < hbexcl->nb; e++)
642 anm = hbexcl->b[e].a[0];
643 i1 = search_atom(anm, astart, atoms,
644 "exclusion", bAllowMissing);
645 anm = hbexcl->b[e].a[1];
646 i2 = search_atom(anm, astart, atoms,
647 "exclusion", bAllowMissing);
648 if (i1 != NO_ATID && i2 != NO_ATID)
656 srenew(excls[i1].e, excls[i1].nr+1);
657 excls[i1].e[excls[i1].nr] = i2;
666 for (a = 0; a < atoms->nr; a++)
670 qsort(excls[a].e, excls[a].nr, (size_t)sizeof(atom_id), atom_id_comp);
675 static void remove_excl(t_excls *excls, int remove)
679 for (i = remove+1; i < excls->nr; i++)
681 excls->e[i-1] = excls->e[i];
687 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
689 int i, j, j1, k, k1, l, l1, m, n, e;
694 /* extract all i-j-k-l neighbours from nnb struct */
695 for (i = 0; (i < nnb->nr); i++)
697 /* For all particles */
700 for (j = 0; (j < nnb->nrexcl[i][1]); j++)
702 /* For all first neighbours */
703 j1 = nnb->a[i][1][j];
705 for (e = 0; e < excl->nr; e++)
707 if (excl->e[e] == j1)
709 remove_excl(excl, e);
715 for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
717 /* For all first neighbours of j1 */
718 k1 = nnb->a[j1][1][k];
720 for (e = 0; e < excl->nr; e++)
722 if (excl->e[e] == k1)
724 remove_excl(excl, e);
730 for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
732 /* For all first neighbours of k1 */
733 l1 = nnb->a[k1][1][l];
735 for (e = 0; e < excl->nr; e++)
737 if (excl->e[e] == l1)
739 remove_excl(excl, e);
751 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
753 int i, j, j1, k, k1, l, l1, m, n, e, N;
756 for (N = 1; (N < min(nrexcl, nnb->nrex)); N++)
758 /* extract all i-j-k-l neighbours from nnb struct */
759 for (i = 0; (i < nnb->nr); i++)
761 /* For all particles */
764 excl->nr += nnb->nrexcl[i][N];
765 srenew(excl->e, excl->nr);
766 for (j = 0; (j < nnb->nrexcl[i][N]); j++)
768 /* For all first neighbours */
769 if (nnb->a[i][N][j] != i)
771 excl->e[n++] = nnb->a[i][N][j];
778 /* Generate pairs, angles and dihedrals from .rtp settings */
779 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
780 t_params plist[], t_excls excls[], t_hackblock hb[],
781 gmx_bool bAllowMissing)
783 t_param *ang, *dih, *pai, *improper;
784 t_rbondeds *hbang, *hbdih;
786 int res, minres, maxres;
787 int i, j, j1, k, k1, l, l1, m, n, i1, i2;
788 int ninc, maxang, maxdih, maxpai;
789 int nang, ndih, npai, nimproper, nbd;
791 gmx_bool bFound, bExcl;
794 /* These are the angles, dihedrals and pairs that we generate
795 * from the bonds. The ones that are already there from the rtp file
802 maxang = maxdih = maxpai = ninc;
808 for (i = 0; i < 4; i++)
815 gen_excls(atoms, excls, hb, bAllowMissing);
818 /* Extract all i-j-k-l neighbours from nnb struct to generate all
819 * angles and dihedrals. */
820 for (i = 0; (i < nnb->nr); i++)
822 /* For all particles */
823 for (j = 0; (j < nnb->nrexcl[i][1]); j++)
825 /* For all first neighbours */
826 j1 = nnb->a[i][1][j];
827 for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
829 /* For all first neighbours of j1 */
830 k1 = nnb->a[j1][1][k];
833 /* Generate every angle only once */
844 ang[nang].C0 = NOTSET;
845 ang[nang].C1 = NOTSET;
846 set_p_string(&(ang[nang]), "");
849 minres = atoms->atom[ang[nang].a[0]].resind;
851 for (m = 1; m < 3; m++)
853 minres = min(minres, atoms->atom[ang[nang].a[m]].resind);
854 maxres = max(maxres, atoms->atom[ang[nang].a[m]].resind);
856 res = 2*minres-maxres;
859 res += maxres-minres;
860 get_atomnames_min(3, anm, res, atoms, ang[nang].a);
861 hbang = &hb[res].rb[ebtsANGLES];
862 for (l = 0; (l < hbang->nb); l++)
864 if (strcmp(anm[1], hbang->b[l].AJ) == 0)
867 for (m = 0; m < 3; m += 2)
870 ((strcmp(anm[m], hbang->b[l].AI) == 0) &&
871 (strcmp(anm[2-m], hbang->b[l].AK) == 0)));
875 set_p_string(&(ang[nang]), hbang->b[l].s);
880 while (res < maxres);
884 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
888 for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
890 /* For all first neighbours of k1 */
891 l1 = nnb->a[k1][1][l];
892 if ((l1 != i) && (l1 != j1))
903 for (m = 0; m < MAXFORCEPARAM; m++)
905 dih[ndih].c[m] = NOTSET;
907 set_p_string(&(dih[ndih]), "");
911 minres = atoms->atom[dih[ndih].a[0]].resind;
913 for (m = 1; m < 4; m++)
915 minres = min(minres, atoms->atom[dih[ndih].a[m]].resind);
916 maxres = max(maxres, atoms->atom[dih[ndih].a[m]].resind);
918 res = 2*minres-maxres;
921 res += maxres-minres;
922 get_atomnames_min(4, anm, res, atoms, dih[ndih].a);
923 hbdih = &hb[res].rb[ebtsPDIHS];
924 for (n = 0; (n < hbdih->nb); n++)
927 for (m = 0; m < 2; m++)
930 ((strcmp(anm[3*m], hbdih->b[n].AI) == 0) &&
931 (strcmp(anm[1+m], hbdih->b[n].AJ) == 0) &&
932 (strcmp(anm[2-m], hbdih->b[n].AK) == 0) &&
933 (strcmp(anm[3-3*m], hbdih->b[n].AL) == 0)));
937 set_p_string(&dih[ndih], hbdih->b[n].s);
939 /* Set the last parameter to be able to see
940 if the dihedral was in the rtp list.
942 dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
945 /* Set the next direct in case the rtp contains
946 multiple entries for this dihedral.
957 for (m = 0; m < MAXFORCEPARAM; m++)
959 dih[ndih].c[m] = NOTSET;
964 while (res < maxres);
977 for (m = 0; m < MAXFORCEPARAM; m++)
979 dih[ndih].c[m] = NOTSET;
981 set_p_string(&(dih[ndih]), "");
985 nbd = nb_dist(nnb, i, l1);
988 fprintf(debug, "Distance (%d-%d) = %d\n", i+1, l1+1, nbd);
995 for (m = 0; m < excls[i1].nr; m++)
997 bExcl = bExcl || excls[i1].e[m] == i2;
1001 if (rtp[0].bGenerateHH14Interactions ||
1002 !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
1007 srenew(pai, maxpai);
1011 pai[npai].C0 = NOTSET;
1012 pai[npai].C1 = NOTSET;
1013 set_p_string(&(pai[npai]), "");
1026 /* Sort angles with respect to j-i-k (middle atom first) */
1029 qsort(ang, nang, (size_t)sizeof(ang[0]), acomp);
1032 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
1035 qsort(dih, ndih, (size_t)sizeof(dih[0]), dcomp);
1038 /* Sort the pairs */
1041 qsort(pai, npai, (size_t)sizeof(pai[0]), pcomp);
1045 /* Remove doubles, could occur in 6-rings, such as phenyls,
1046 maybe one does not want this when fudgeQQ < 1.
1048 fprintf(stderr, "Before cleaning: %d pairs\n", npai);
1049 rm2par(pai, &npai, preq);
1052 /* Get the impropers from the database */
1053 nimproper = get_impropers(atoms, hb, &improper, bAllowMissing);
1055 /* Sort the impropers */
1056 sort_id(nimproper, improper);
1060 fprintf(stderr, "Before cleaning: %d dihedrals\n", ndih);
1061 clean_dih(dih, &ndih, improper, nimproper, atoms,
1062 rtp[0].bKeepAllGeneratedDihedrals,
1063 rtp[0].bRemoveDihedralIfWithImproper);
1066 /* Now we have unique lists of angles and dihedrals
1067 * Copy them into the destination struct
1069 cppar(ang, nang, plist, F_ANGLES);
1070 cppar(dih, ndih, plist, F_PDIHS);
1071 cppar(improper, nimproper, plist, F_IDIHS);
1072 cppar(pai, npai, plist, F_LJ14);
1074 /* Remove all exclusions which are within nrexcl */
1075 clean_excls(nnb, rtp[0].nrexcl, excls);