2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 /* This file is completely threadsafe - keep it that way! */
48 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/utility/cstringutil.h"
50 #include "gromacs/fileio/confio.h"
51 #include "gromacs/math/vec.h"
54 #include "gpp_nextnb.h"
57 #include "gromacs/utility/fatalerror.h"
62 #define DIHEDRAL_WAS_SET_IN_RTP 0
63 static gmx_bool was_dihedral_set_in_rtp(t_param *dih)
65 return dih->c[MAXFORCEPARAM-1] == DIHEDRAL_WAS_SET_IN_RTP;
68 typedef gmx_bool (*peq)(t_param *p1, t_param *p2);
70 static int acomp(const void *a1, const void *a2)
77 if ((ac = (p1->AJ-p2->AJ)) != 0)
81 else if ((ac = (p1->AI-p2->AI)) != 0)
87 return (p1->AK-p2->AK);
91 static int pcomp(const void *a1, const void *a2)
98 if ((pc = (p1->AI-p2->AI)) != 0)
104 return (p1->AJ-p2->AJ);
108 static int dcomp(const void *d1, const void *d2)
115 /* First sort by J & K (the two central) atoms */
116 if ((dc = (p1->AJ-p2->AJ)) != 0)
120 else if ((dc = (p1->AK-p2->AK)) != 0)
124 /* Then make sure to put rtp dihedrals before generated ones */
125 else if (was_dihedral_set_in_rtp(p1) &&
126 !was_dihedral_set_in_rtp(p2))
130 else if (!was_dihedral_set_in_rtp(p1) &&
131 was_dihedral_set_in_rtp(p2))
135 /* Finally, sort by I and J (two outer) atoms */
136 else if ((dc = (p1->AI-p2->AI)) != 0)
142 return (p1->AL-p2->AL);
147 static gmx_bool is_dihedral_on_same_bond(t_param *p1, t_param *p2)
149 if (((p1->AJ == p2->AJ) && (p1->AK == p2->AK)) ||
150 ((p1->AJ == p2->AK) && (p1->AK == p2->AJ)))
161 static gmx_bool preq(t_param *p1, t_param *p2)
163 if ((p1->AI == p2->AI) && (p1->AJ == p2->AJ))
173 static void rm2par(t_param p[], int *np, peq eq)
186 for (i = 1; (i < (*np)); i++)
188 if (!eq(&p[i], &p[i-1]))
193 /* Index now holds pointers to all the non-equal params,
194 * this only works when p is sorted of course
196 for (i = 0; (i < nind); i++)
198 for (j = 0; (j < MAXATOMLIST); j++)
200 p[i].a[j] = p[index[i]].a[j];
202 for (j = 0; (j < MAXFORCEPARAM); j++)
204 p[i].c[j] = p[index[i]].c[j];
206 if (p[index[i]].a[0] == p[index[i]].a[1])
211 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
212 "a[0] %d a[1] %d a[2] %d a[3] %d\n",
213 p[i].a[0], p[i].a[1], p[i].a[2], p[i].a[3]);
217 else if (index[i] > i)
219 /* Copy the string only if it comes from somewhere else
220 * otherwise we will end up copying a random (newly freed) pointer.
221 * Since the index is sorted we only have to test for index[i] > i.
223 strcpy(p[i].s, p[index[i]].s);
231 static void cppar(t_param p[], int np, t_params plist[], int ftype)
233 int i, j, nral, nrfp;
242 for (i = 0; (i < np); i++)
244 for (j = 0; (j < nral); j++)
246 ps->param[ps->nr].a[j] = p[i].a[j];
248 for (j = 0; (j < nrfp); j++)
250 ps->param[ps->nr].c[j] = p[i].c[j];
252 for (j = 0; (j < MAXSLEN); j++)
254 ps->param[ps->nr].s[j] = p[i].s[j];
260 static void cpparam(t_param *dest, t_param *src)
264 for (j = 0; (j < MAXATOMLIST); j++)
266 dest->a[j] = src->a[j];
268 for (j = 0; (j < MAXFORCEPARAM); j++)
270 dest->c[j] = src->c[j];
272 for (j = 0; (j < MAXSLEN); j++)
274 dest->s[j] = src->s[j];
278 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
282 for (j = 0; (j < 4); j++)
286 for (j = 0; (j < MAXFORCEPARAM); j++)
301 static int int_comp(const void *a, const void *b)
303 return (*(int *)a) - (*(int *)b);
306 static int atom_id_comp(const void *a, const void *b)
308 return (*(atom_id *)a) - (*(atom_id *)b);
311 static int eq_imp(atom_id a1[], atom_id a2[])
316 for (j = 0; (j < 4); j++)
321 qsort(b1, 4, (size_t)sizeof(b1[0]), int_comp);
322 qsort(b2, 4, (size_t)sizeof(b2[0]), int_comp);
324 for (j = 0; (j < 4); j++)
335 static int idcomp(const void *a, const void *b)
342 if ((d = (pa->a[0]-pb->a[0])) != 0)
346 else if ((d = (pa->a[3]-pb->a[3])) != 0)
350 else if ((d = (pa->a[1]-pb->a[1])) != 0)
356 return (int) (pa->a[2]-pb->a[2]);
360 static void sort_id(int nr, t_param ps[])
364 /* First swap order of atoms around if necessary */
365 for (i = 0; (i < nr); i++)
367 if (ps[i].a[3] < ps[i].a[0])
369 tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
370 tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
376 qsort(ps, nr, (size_t)sizeof(ps[0]), idcomp);
380 static int n_hydro(atom_id a[], char ***atomname)
385 for (i = 0; (i < 4); i += 3)
387 aname = *atomname[a[i]];
388 c0 = toupper(aname[0]);
393 else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9'))
395 c1 = toupper(aname[1]);
405 /* Clean up the dihedrals (both generated and read from the .rtp
407 static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper,
408 t_atoms *atoms, gmx_bool bKeepAllGeneratedDihedrals,
409 gmx_bool bRemoveDihedralIfWithImproper)
414 /* Construct the list of the indices of the dihedrals
415 * (i.e. generated or read) that might be kept. */
416 snew(index, *ndih+1);
417 if (bKeepAllGeneratedDihedrals)
419 fprintf(stderr, "Keeping all generated dihedrals\n");
421 for (i = 0; i < nind; i++)
430 /* Check if generated dihedral i should be removed. The
431 * dihedrals have been sorted by dcomp() above, so all those
432 * on the same two central atoms are together, with those from
433 * the .rtp file preceding those that were automatically
434 * generated. We remove the latter if the former exist. */
435 for (i = 0; i < *ndih; i++)
437 /* Keep the dihedrals that were defined in the .rtp file,
438 * and the dihedrals that were generated and different
439 * from the last one (whether it was generated or not). */
440 if (was_dihedral_set_in_rtp(&dih[i]) ||
442 !is_dihedral_on_same_bond(&dih[i], &dih[i-1]))
451 for (i = 0; i < nind; i++)
453 gmx_bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
454 gmx_bool bKeep = TRUE;
455 if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
457 /* Remove the dihedral if there is an improper on the same
459 for (j = 0; j < nimproper && bKeep; j++)
461 bKeep = !is_dihedral_on_same_bond(&dih[index[i]], &improper[j]);
467 /* If we don't want all dihedrals, we want to select the
468 * ones with the fewest hydrogens. Note that any generated
469 * dihedrals on the same bond as an .rtp dihedral may have
470 * been already pruned above in the construction of
471 * index[]. However, their parameters are still present,
472 * and l is looping over this dihedral and all of its
473 * pruned siblings. */
474 int bestl = index[i];
475 if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
477 /* Minimum number of hydrogens for i and l atoms */
481 is_dihedral_on_same_bond(&dih[index[i]], &dih[l]));
484 int nh = n_hydro(dih[l].a, atoms->atomname);
498 cpparam(&dih[k], &dih[bestl]);
504 for (i = k; i < *ndih; i++)
506 strcpy(dih[i].s, "");
513 static int get_impropers(t_atoms *atoms, t_hackblock hb[], t_param **improper,
514 gmx_bool bAllowMissing)
517 t_rbondeds *impropers;
518 t_rbonded *hbimproper;
519 int nimproper, i, j, k, r, start, ninc, nalloc;
520 atom_id ai[MAXATOMLIST];
525 snew(*improper, nalloc);
527 /* Add all the impropers from the residue database to the list. */
532 for (i = 0; (i < atoms->nres); i++)
534 impropers = &hb[i].rb[ebtsIDIHS];
535 for (j = 0; (j < impropers->nb); j++)
538 for (k = 0; (k < 4) && !bStop; k++)
540 ai[k] = search_atom(impropers->b[j].a[k], start,
542 "improper", bAllowMissing);
543 if (ai[k] == NO_ATID)
550 if (nimproper == nalloc)
553 srenew(*improper, nalloc);
556 set_p(&((*improper)[nimproper]), ai, NULL, impropers->b[j].s);
560 while ((start < atoms->nr) && (atoms->atom[start].resind == i))
570 static int nb_dist(t_nextnb *nnb, int ai, int aj)
582 nrexcl = nnb->nrexcl[ai];
583 for (nre = 1; (nre < nnb->nrex); nre++)
586 for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
588 if ((aj == a[nrx]) && (NRE == -1))
597 gmx_bool is_hydro(t_atoms *atoms, int ai)
599 return ((*(atoms->atomname[ai]))[0] == 'H');
602 static void get_atomnames_min(int n, char **anm,
603 int resind, t_atoms *atoms, atom_id *a)
607 /* Assume ascending residue numbering */
608 for (m = 0; m < n; m++)
610 if (atoms->atom[a[m]].resind < resind)
614 else if (atoms->atom[a[m]].resind > resind)
622 strcat(anm[m], *(atoms->atomname[a[m]]));
626 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
627 gmx_bool bAllowMissing)
630 atom_id a, astart, i1, i2, itmp;
636 for (a = 0; a < atoms->nr; a++)
638 r = atoms->atom[a].resind;
639 if (a == atoms->nr-1 || atoms->atom[a+1].resind != r)
641 hbexcl = &hb[r].rb[ebtsEXCLS];
643 for (e = 0; e < hbexcl->nb; e++)
645 anm = hbexcl->b[e].a[0];
646 i1 = search_atom(anm, astart, atoms,
647 "exclusion", bAllowMissing);
648 anm = hbexcl->b[e].a[1];
649 i2 = search_atom(anm, astart, atoms,
650 "exclusion", bAllowMissing);
651 if (i1 != NO_ATID && i2 != NO_ATID)
659 srenew(excls[i1].e, excls[i1].nr+1);
660 excls[i1].e[excls[i1].nr] = i2;
669 for (a = 0; a < atoms->nr; a++)
673 qsort(excls[a].e, excls[a].nr, (size_t)sizeof(atom_id), atom_id_comp);
678 static void remove_excl(t_excls *excls, int remove)
682 for (i = remove+1; i < excls->nr; i++)
684 excls->e[i-1] = excls->e[i];
690 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
692 int i, j, j1, k, k1, l, l1, m, n, e;
697 /* extract all i-j-k-l neighbours from nnb struct */
698 for (i = 0; (i < nnb->nr); i++)
700 /* For all particles */
703 for (j = 0; (j < nnb->nrexcl[i][1]); j++)
705 /* For all first neighbours */
706 j1 = nnb->a[i][1][j];
708 for (e = 0; e < excl->nr; e++)
710 if (excl->e[e] == j1)
712 remove_excl(excl, e);
718 for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
720 /* For all first neighbours of j1 */
721 k1 = nnb->a[j1][1][k];
723 for (e = 0; e < excl->nr; e++)
725 if (excl->e[e] == k1)
727 remove_excl(excl, e);
733 for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
735 /* For all first neighbours of k1 */
736 l1 = nnb->a[k1][1][l];
738 for (e = 0; e < excl->nr; e++)
740 if (excl->e[e] == l1)
742 remove_excl(excl, e);
754 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
756 int i, j, j1, k, k1, l, l1, m, n, e, N;
759 for (N = 1; (N < min(nrexcl, nnb->nrex)); N++)
761 /* extract all i-j-k-l neighbours from nnb struct */
762 for (i = 0; (i < nnb->nr); i++)
764 /* For all particles */
767 excl->nr += nnb->nrexcl[i][N];
768 srenew(excl->e, excl->nr);
769 for (j = 0; (j < nnb->nrexcl[i][N]); j++)
771 /* For all first neighbours */
772 if (nnb->a[i][N][j] != i)
774 excl->e[n++] = nnb->a[i][N][j];
781 /* Generate pairs, angles and dihedrals from .rtp settings */
782 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
783 t_params plist[], t_excls excls[], t_hackblock hb[],
784 gmx_bool bAllowMissing)
786 t_param *ang, *dih, *pai, *improper;
787 t_rbondeds *hbang, *hbdih;
789 int res, minres, maxres;
790 int i, j, j1, k, k1, l, l1, m, n, i1, i2;
791 int ninc, maxang, maxdih, maxpai;
792 int nang, ndih, npai, nimproper, nbd;
794 gmx_bool bFound, bExcl;
797 /* These are the angles, dihedrals and pairs that we generate
798 * from the bonds. The ones that are already there from the rtp file
805 maxang = maxdih = maxpai = ninc;
811 for (i = 0; i < 4; i++)
818 gen_excls(atoms, excls, hb, bAllowMissing);
821 /* Extract all i-j-k-l neighbours from nnb struct to generate all
822 * angles and dihedrals. */
823 for (i = 0; (i < nnb->nr); i++)
825 /* For all particles */
826 for (j = 0; (j < nnb->nrexcl[i][1]); j++)
828 /* For all first neighbours */
829 j1 = nnb->a[i][1][j];
830 for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
832 /* For all first neighbours of j1 */
833 k1 = nnb->a[j1][1][k];
836 /* Generate every angle only once */
847 ang[nang].C0 = NOTSET;
848 ang[nang].C1 = NOTSET;
849 set_p_string(&(ang[nang]), "");
852 minres = atoms->atom[ang[nang].a[0]].resind;
854 for (m = 1; m < 3; m++)
856 minres = min(minres, atoms->atom[ang[nang].a[m]].resind);
857 maxres = max(maxres, atoms->atom[ang[nang].a[m]].resind);
859 res = 2*minres-maxres;
862 res += maxres-minres;
863 get_atomnames_min(3, anm, res, atoms, ang[nang].a);
864 hbang = &hb[res].rb[ebtsANGLES];
865 for (l = 0; (l < hbang->nb); l++)
867 if (strcmp(anm[1], hbang->b[l].AJ) == 0)
870 for (m = 0; m < 3; m += 2)
873 ((strcmp(anm[m], hbang->b[l].AI) == 0) &&
874 (strcmp(anm[2-m], hbang->b[l].AK) == 0)));
878 set_p_string(&(ang[nang]), hbang->b[l].s);
883 while (res < maxres);
887 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
891 for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
893 /* For all first neighbours of k1 */
894 l1 = nnb->a[k1][1][l];
895 if ((l1 != i) && (l1 != j1))
906 for (m = 0; m < MAXFORCEPARAM; m++)
908 dih[ndih].c[m] = NOTSET;
910 set_p_string(&(dih[ndih]), "");
914 minres = atoms->atom[dih[ndih].a[0]].resind;
916 for (m = 1; m < 4; m++)
918 minres = min(minres, atoms->atom[dih[ndih].a[m]].resind);
919 maxres = max(maxres, atoms->atom[dih[ndih].a[m]].resind);
921 res = 2*minres-maxres;
924 res += maxres-minres;
925 get_atomnames_min(4, anm, res, atoms, dih[ndih].a);
926 hbdih = &hb[res].rb[ebtsPDIHS];
927 for (n = 0; (n < hbdih->nb); n++)
930 for (m = 0; m < 2; m++)
933 ((strcmp(anm[3*m], hbdih->b[n].AI) == 0) &&
934 (strcmp(anm[1+m], hbdih->b[n].AJ) == 0) &&
935 (strcmp(anm[2-m], hbdih->b[n].AK) == 0) &&
936 (strcmp(anm[3-3*m], hbdih->b[n].AL) == 0)));
940 set_p_string(&dih[ndih], hbdih->b[n].s);
942 /* Set the last parameter to be able to see
943 if the dihedral was in the rtp list.
945 dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
948 /* Set the next direct in case the rtp contains
949 multiple entries for this dihedral.
960 for (m = 0; m < MAXFORCEPARAM; m++)
962 dih[ndih].c[m] = NOTSET;
967 while (res < maxres);
980 for (m = 0; m < MAXFORCEPARAM; m++)
982 dih[ndih].c[m] = NOTSET;
984 set_p_string(&(dih[ndih]), "");
988 nbd = nb_dist(nnb, i, l1);
991 fprintf(debug, "Distance (%d-%d) = %d\n", i+1, l1+1, nbd);
998 for (m = 0; m < excls[i1].nr; m++)
1000 bExcl = bExcl || excls[i1].e[m] == i2;
1004 if (rtp[0].bGenerateHH14Interactions ||
1005 !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
1010 srenew(pai, maxpai);
1014 pai[npai].C0 = NOTSET;
1015 pai[npai].C1 = NOTSET;
1016 set_p_string(&(pai[npai]), "");
1029 /* Sort angles with respect to j-i-k (middle atom first) */
1032 qsort(ang, nang, (size_t)sizeof(ang[0]), acomp);
1035 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
1038 qsort(dih, ndih, (size_t)sizeof(dih[0]), dcomp);
1041 /* Sort the pairs */
1044 qsort(pai, npai, (size_t)sizeof(pai[0]), pcomp);
1048 /* Remove doubles, could occur in 6-rings, such as phenyls,
1049 maybe one does not want this when fudgeQQ < 1.
1051 fprintf(stderr, "Before cleaning: %d pairs\n", npai);
1052 rm2par(pai, &npai, preq);
1055 /* Get the impropers from the database */
1056 nimproper = get_impropers(atoms, hb, &improper, bAllowMissing);
1058 /* Sort the impropers */
1059 sort_id(nimproper, improper);
1063 fprintf(stderr, "Before cleaning: %d dihedrals\n", ndih);
1064 clean_dih(dih, &ndih, improper, nimproper, atoms,
1065 rtp[0].bKeepAllGeneratedDihedrals,
1066 rtp[0].bRemoveDihedralIfWithImproper);
1069 /* Now we have unique lists of angles and dihedrals
1070 * Copy them into the destination struct
1072 cppar(ang, nang, plist, F_ANGLES);
1073 cppar(dih, ndih, plist, F_PDIHS);
1074 cppar(improper, nimproper, plist, F_IDIHS);
1075 cppar(pai, npai, plist, F_LJ14);
1077 /* Remove all exclusions which are within nrexcl */
1078 clean_excls(nnb, rtp[0].nrexcl, excls);