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,2015, 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! */
50 #include "gromacs/fileio/confio.h"
51 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
52 #include "gromacs/gmxpreprocess/notset.h"
53 #include "gromacs/gmxpreprocess/pgutil.h"
54 #include "gromacs/gmxpreprocess/resall.h"
55 #include "gromacs/gmxpreprocess/topio.h"
56 #include "gromacs/gmxpreprocess/toputil.h"
57 #include "gromacs/legacyheaders/types/ifunc.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/utility/cstringutil.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/smalloc.h"
63 #define DIHEDRAL_WAS_SET_IN_RTP 0
64 static gmx_bool
was_dihedral_set_in_rtp(t_param
*dih
)
66 return dih
->c
[MAXFORCEPARAM
-1] == DIHEDRAL_WAS_SET_IN_RTP
;
69 typedef gmx_bool (*peq
)(t_param
*p1
, t_param
*p2
);
71 static int acomp(const void *a1
, const void *a2
)
78 if ((ac
= (p1
->aj()-p2
->aj())) != 0)
82 else if ((ac
= (p1
->ai()-p2
->ai())) != 0)
88 return (p1
->ak()-p2
->ak());
92 static int pcomp(const void *a1
, const void *a2
)
99 if ((pc
= (p1
->ai()-p2
->ai())) != 0)
105 return (p1
->aj()-p2
->aj());
109 static int dcomp(const void *d1
, const void *d2
)
116 /* First sort by J & K (the two central) atoms */
117 if ((dc
= (p1
->aj()-p2
->aj())) != 0)
121 else if ((dc
= (p1
->ak()-p2
->ak())) != 0)
125 /* Then make sure to put rtp dihedrals before generated ones */
126 else if (was_dihedral_set_in_rtp(p1
) &&
127 !was_dihedral_set_in_rtp(p2
))
131 else if (!was_dihedral_set_in_rtp(p1
) &&
132 was_dihedral_set_in_rtp(p2
))
136 /* Finally, sort by I and J (two outer) atoms */
137 else if ((dc
= (p1
->ai()-p2
->ai())) != 0)
143 return (p1
->al()-p2
->al());
148 static gmx_bool
is_dihedral_on_same_bond(t_param
*p1
, t_param
*p2
)
150 if (((p1
->aj() == p2
->aj()) && (p1
->ak() == p2
->ak())) ||
151 ((p1
->aj() == p2
->ak()) && (p1
->ak() == p2
->aj())))
162 static gmx_bool
preq(t_param
*p1
, t_param
*p2
)
164 if ((p1
->ai() == p2
->ai()) && (p1
->aj() == p2
->aj()))
174 static void rm2par(t_param p
[], int *np
, peq eq
)
187 for (i
= 1; (i
< (*np
)); i
++)
189 if (!eq(&p
[i
], &p
[i
-1]))
194 /* Index now holds pointers to all the non-equal params,
195 * this only works when p is sorted of course
197 for (i
= 0; (i
< nind
); i
++)
199 for (j
= 0; (j
< MAXATOMLIST
); j
++)
201 p
[i
].a
[j
] = p
[index
[i
]].a
[j
];
203 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
205 p
[i
].c
[j
] = p
[index
[i
]].c
[j
];
207 if (p
[index
[i
]].a
[0] == p
[index
[i
]].a
[1])
212 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
213 "a[0] %d a[1] %d a[2] %d a[3] %d\n",
214 p
[i
].a
[0], p
[i
].a
[1], p
[i
].a
[2], p
[i
].a
[3]);
218 else if (index
[i
] > i
)
220 /* Copy the string only if it comes from somewhere else
221 * otherwise we will end up copying a random (newly freed) pointer.
222 * Since the index is sorted we only have to test for index[i] > i.
224 strcpy(p
[i
].s
, p
[index
[i
]].s
);
232 static void cppar(t_param p
[], int np
, t_params plist
[], int ftype
)
234 int i
, j
, nral
, nrfp
;
243 for (i
= 0; (i
< np
); i
++)
245 for (j
= 0; (j
< nral
); j
++)
247 ps
->param
[ps
->nr
].a
[j
] = p
[i
].a
[j
];
249 for (j
= 0; (j
< nrfp
); j
++)
251 ps
->param
[ps
->nr
].c
[j
] = p
[i
].c
[j
];
253 for (j
= 0; (j
< MAXSLEN
); j
++)
255 ps
->param
[ps
->nr
].s
[j
] = p
[i
].s
[j
];
261 static void cpparam(t_param
*dest
, t_param
*src
)
265 for (j
= 0; (j
< MAXATOMLIST
); j
++)
267 dest
->a
[j
] = src
->a
[j
];
269 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
271 dest
->c
[j
] = src
->c
[j
];
273 for (j
= 0; (j
< MAXSLEN
); j
++)
275 dest
->s
[j
] = src
->s
[j
];
279 static void set_p(t_param
*p
, atom_id ai
[4], real
*c
, char *s
)
283 for (j
= 0; (j
< 4); j
++)
287 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
302 static int int_comp(const void *a
, const void *b
)
304 return (*(int *)a
) - (*(int *)b
);
307 static int atom_id_comp(const void *a
, const void *b
)
309 return (*(atom_id
*)a
) - (*(atom_id
*)b
);
312 static int eq_imp(atom_id a1
[], atom_id a2
[])
317 for (j
= 0; (j
< 4); j
++)
322 qsort(b1
, 4, (size_t)sizeof(b1
[0]), int_comp
);
323 qsort(b2
, 4, (size_t)sizeof(b2
[0]), int_comp
);
325 for (j
= 0; (j
< 4); j
++)
336 static int idcomp(const void *a
, const void *b
)
343 if ((d
= (pa
->a
[0]-pb
->a
[0])) != 0)
347 else if ((d
= (pa
->a
[3]-pb
->a
[3])) != 0)
351 else if ((d
= (pa
->a
[1]-pb
->a
[1])) != 0)
357 return (int) (pa
->a
[2]-pb
->a
[2]);
361 static void sort_id(int nr
, t_param ps
[])
365 /* First swap order of atoms around if necessary */
366 for (i
= 0; (i
< nr
); i
++)
368 if (ps
[i
].a
[3] < ps
[i
].a
[0])
370 tmp
= ps
[i
].a
[3]; ps
[i
].a
[3] = ps
[i
].a
[0]; ps
[i
].a
[0] = tmp
;
371 tmp
= ps
[i
].a
[2]; ps
[i
].a
[2] = ps
[i
].a
[1]; ps
[i
].a
[1] = tmp
;
377 qsort(ps
, nr
, (size_t)sizeof(ps
[0]), idcomp
);
381 static int n_hydro(atom_id a
[], char ***atomname
)
386 for (i
= 0; (i
< 4); i
+= 3)
388 aname
= *atomname
[a
[i
]];
389 c0
= toupper(aname
[0]);
394 else if (((int)strlen(aname
) > 1) && (c0
>= '0') && (c0
<= '9'))
396 c1
= toupper(aname
[1]);
406 /* Clean up the dihedrals (both generated and read from the .rtp
408 static void clean_dih(t_param
*dih
, int *ndih
, t_param improper
[], int nimproper
,
409 t_atoms
*atoms
, gmx_bool bKeepAllGeneratedDihedrals
,
410 gmx_bool bRemoveDihedralIfWithImproper
)
415 /* Construct the list of the indices of the dihedrals
416 * (i.e. generated or read) that might be kept. */
417 snew(index
, *ndih
+1);
418 if (bKeepAllGeneratedDihedrals
)
420 fprintf(stderr
, "Keeping all generated dihedrals\n");
422 for (i
= 0; i
< nind
; i
++)
431 /* Check if generated dihedral i should be removed. The
432 * dihedrals have been sorted by dcomp() above, so all those
433 * on the same two central atoms are together, with those from
434 * the .rtp file preceding those that were automatically
435 * generated. We remove the latter if the former exist. */
436 for (i
= 0; i
< *ndih
; i
++)
438 /* Keep the dihedrals that were defined in the .rtp file,
439 * and the dihedrals that were generated and different
440 * from the last one (whether it was generated or not). */
441 if (was_dihedral_set_in_rtp(&dih
[i
]) ||
443 !is_dihedral_on_same_bond(&dih
[i
], &dih
[i
-1]))
452 for (i
= 0; i
< nind
; i
++)
454 gmx_bool bWasSetInRTP
= was_dihedral_set_in_rtp(&dih
[index
[i
]]);
455 gmx_bool bKeep
= TRUE
;
456 if (!bWasSetInRTP
&& bRemoveDihedralIfWithImproper
)
458 /* Remove the dihedral if there is an improper on the same
460 for (j
= 0; j
< nimproper
&& bKeep
; j
++)
462 bKeep
= !is_dihedral_on_same_bond(&dih
[index
[i
]], &improper
[j
]);
468 /* If we don't want all dihedrals, we want to select the
469 * ones with the fewest hydrogens. Note that any generated
470 * dihedrals on the same bond as an .rtp dihedral may have
471 * been already pruned above in the construction of
472 * index[]. However, their parameters are still present,
473 * and l is looping over this dihedral and all of its
474 * pruned siblings. */
475 int bestl
= index
[i
];
476 if (!bKeepAllGeneratedDihedrals
&& !bWasSetInRTP
)
478 /* Minimum number of hydrogens for i and l atoms */
482 is_dihedral_on_same_bond(&dih
[index
[i
]], &dih
[l
]));
485 int nh
= n_hydro(dih
[l
].a
, atoms
->atomname
);
499 cpparam(&dih
[k
], &dih
[bestl
]);
505 for (i
= k
; i
< *ndih
; i
++)
507 strcpy(dih
[i
].s
, "");
514 static int get_impropers(t_atoms
*atoms
, t_hackblock hb
[], t_param
**improper
,
515 gmx_bool bAllowMissing
)
517 t_rbondeds
*impropers
;
518 int nimproper
, i
, j
, k
, start
, ninc
, nalloc
;
519 atom_id ai
[MAXATOMLIST
];
524 snew(*improper
, nalloc
);
526 /* Add all the impropers from the residue database to the list. */
531 for (i
= 0; (i
< atoms
->nres
); i
++)
533 impropers
= &hb
[i
].rb
[ebtsIDIHS
];
534 for (j
= 0; (j
< impropers
->nb
); j
++)
537 for (k
= 0; (k
< 4) && !bStop
; k
++)
539 ai
[k
] = search_atom(impropers
->b
[j
].a
[k
], start
,
541 "improper", bAllowMissing
);
542 if (ai
[k
] == NO_ATID
)
549 if (nimproper
== nalloc
)
552 srenew(*improper
, nalloc
);
555 set_p(&((*improper
)[nimproper
]), ai
, NULL
, impropers
->b
[j
].s
);
559 while ((start
< atoms
->nr
) && (atoms
->atom
[start
].resind
== i
))
569 static int nb_dist(t_nextnb
*nnb
, int ai
, int aj
)
581 nrexcl
= nnb
->nrexcl
[ai
];
582 for (nre
= 1; (nre
< nnb
->nrex
); nre
++)
585 for (nrx
= 0; (nrx
< nrexcl
[nre
]); nrx
++)
587 if ((aj
== a
[nrx
]) && (NRE
== -1))
596 gmx_bool
is_hydro(t_atoms
*atoms
, int ai
)
598 return ((*(atoms
->atomname
[ai
]))[0] == 'H');
601 static void get_atomnames_min(int n
, char **anm
,
602 int resind
, t_atoms
*atoms
, atom_id
*a
)
606 /* Assume ascending residue numbering */
607 for (m
= 0; m
< n
; m
++)
609 if (atoms
->atom
[a
[m
]].resind
< resind
)
613 else if (atoms
->atom
[a
[m
]].resind
> resind
)
621 strcat(anm
[m
], *(atoms
->atomname
[a
[m
]]));
625 static void gen_excls(t_atoms
*atoms
, t_excls
*excls
, t_hackblock hb
[],
626 gmx_bool bAllowMissing
)
629 atom_id a
, astart
, i1
, i2
, itmp
;
635 for (a
= 0; a
< atoms
->nr
; a
++)
637 r
= atoms
->atom
[a
].resind
;
638 if (a
== atoms
->nr
-1 || atoms
->atom
[a
+1].resind
!= r
)
640 hbexcl
= &hb
[r
].rb
[ebtsEXCLS
];
642 for (e
= 0; e
< hbexcl
->nb
; e
++)
644 anm
= hbexcl
->b
[e
].a
[0];
645 i1
= search_atom(anm
, astart
, atoms
,
646 "exclusion", bAllowMissing
);
647 anm
= hbexcl
->b
[e
].a
[1];
648 i2
= search_atom(anm
, astart
, atoms
,
649 "exclusion", bAllowMissing
);
650 if (i1
!= NO_ATID
&& i2
!= NO_ATID
)
658 srenew(excls
[i1
].e
, excls
[i1
].nr
+1);
659 excls
[i1
].e
[excls
[i1
].nr
] = i2
;
668 for (a
= 0; a
< atoms
->nr
; a
++)
672 qsort(excls
[a
].e
, excls
[a
].nr
, (size_t)sizeof(atom_id
), atom_id_comp
);
677 static void remove_excl(t_excls
*excls
, int remove
)
681 for (i
= remove
+1; i
< excls
->nr
; i
++)
683 excls
->e
[i
-1] = excls
->e
[i
];
689 void clean_excls(t_nextnb
*nnb
, int nrexcl
, t_excls excls
[])
691 int i
, j
, j1
, k
, k1
, l
, l1
, e
;
696 /* extract all i-j-k-l neighbours from nnb struct */
697 for (i
= 0; (i
< nnb
->nr
); i
++)
699 /* For all particles */
702 for (j
= 0; (j
< nnb
->nrexcl
[i
][1]); j
++)
704 /* For all first neighbours */
705 j1
= nnb
->a
[i
][1][j
];
707 for (e
= 0; e
< excl
->nr
; e
++)
709 if (excl
->e
[e
] == j1
)
711 remove_excl(excl
, e
);
717 for (k
= 0; (k
< nnb
->nrexcl
[j1
][1]); k
++)
719 /* For all first neighbours of j1 */
720 k1
= nnb
->a
[j1
][1][k
];
722 for (e
= 0; e
< excl
->nr
; e
++)
724 if (excl
->e
[e
] == k1
)
726 remove_excl(excl
, e
);
732 for (l
= 0; (l
< nnb
->nrexcl
[k1
][1]); l
++)
734 /* For all first neighbours of k1 */
735 l1
= nnb
->a
[k1
][1][l
];
737 for (e
= 0; e
< excl
->nr
; e
++)
739 if (excl
->e
[e
] == l1
)
741 remove_excl(excl
, e
);
753 void generate_excls(t_nextnb
*nnb
, int nrexcl
, t_excls excls
[])
758 for (N
= 1; (N
< std::min(nrexcl
, nnb
->nrex
)); N
++)
760 /* extract all i-j-k-l neighbours from nnb struct */
761 for (i
= 0; (i
< nnb
->nr
); i
++)
763 /* For all particles */
766 excl
->nr
+= nnb
->nrexcl
[i
][N
];
767 srenew(excl
->e
, excl
->nr
);
768 for (j
= 0; (j
< nnb
->nrexcl
[i
][N
]); j
++)
770 /* For all first neighbours */
771 if (nnb
->a
[i
][N
][j
] != i
)
773 excl
->e
[n
++] = nnb
->a
[i
][N
][j
];
780 /* Generate pairs, angles and dihedrals from .rtp settings */
781 void gen_pad(t_nextnb
*nnb
, t_atoms
*atoms
, t_restp rtp
[],
782 t_params plist
[], t_excls excls
[], t_hackblock hb
[],
783 gmx_bool bAllowMissing
)
785 t_param
*ang
, *dih
, *pai
, *improper
;
786 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
;
796 /* These are the angles, dihedrals and pairs that we generate
797 * from the bonds. The ones that are already there from the rtp file
804 maxang
= maxdih
= maxpai
= ninc
;
810 for (i
= 0; i
< 4; i
++)
817 gen_excls(atoms
, excls
, hb
, bAllowMissing
);
818 /* mark all entries as not matched yet */
819 for (i
= 0; i
< atoms
->nres
; i
++)
821 for (j
= 0; j
< ebtsNR
; j
++)
823 for (k
= 0; k
< hb
[i
].rb
[j
].nb
; k
++)
825 hb
[i
].rb
[j
].b
[k
].match
= FALSE
;
831 /* Extract all i-j-k-l neighbours from nnb struct to generate all
832 * angles and dihedrals. */
833 for (i
= 0; (i
< nnb
->nr
); i
++)
835 /* For all particles */
836 for (j
= 0; (j
< nnb
->nrexcl
[i
][1]); j
++)
838 /* For all first neighbours */
839 j1
= nnb
->a
[i
][1][j
];
840 for (k
= 0; (k
< nnb
->nrexcl
[j1
][1]); k
++)
842 /* For all first neighbours of j1 */
843 k1
= nnb
->a
[j1
][1][k
];
846 /* Generate every angle only once */
857 ang
[nang
].c0() = NOTSET
;
858 ang
[nang
].c1() = NOTSET
;
859 set_p_string(&(ang
[nang
]), "");
862 minres
= atoms
->atom
[ang
[nang
].a
[0]].resind
;
864 for (m
= 1; m
< 3; m
++)
866 minres
= std::min(minres
, atoms
->atom
[ang
[nang
].a
[m
]].resind
);
867 maxres
= std::max(maxres
, atoms
->atom
[ang
[nang
].a
[m
]].resind
);
869 res
= 2*minres
-maxres
;
872 res
+= maxres
-minres
;
873 get_atomnames_min(3, anm
, res
, atoms
, ang
[nang
].a
);
874 hbang
= &hb
[res
].rb
[ebtsANGLES
];
875 for (l
= 0; (l
< hbang
->nb
); l
++)
877 if (strcmp(anm
[1], hbang
->b
[l
].aj()) == 0)
880 for (m
= 0; m
< 3; m
+= 2)
883 ((strcmp(anm
[m
], hbang
->b
[l
].ai()) == 0) &&
884 (strcmp(anm
[2-m
], hbang
->b
[l
].ak()) == 0)));
888 set_p_string(&(ang
[nang
]), hbang
->b
[l
].s
);
889 /* Mark that we found a match for this entry */
890 hbang
->b
[l
].match
= TRUE
;
895 while (res
< maxres
);
899 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
903 for (l
= 0; (l
< nnb
->nrexcl
[k1
][1]); l
++)
905 /* For all first neighbours of k1 */
906 l1
= nnb
->a
[k1
][1][l
];
907 if ((l1
!= i
) && (l1
!= j1
))
918 for (m
= 0; m
< MAXFORCEPARAM
; m
++)
920 dih
[ndih
].c
[m
] = NOTSET
;
922 set_p_string(&(dih
[ndih
]), "");
926 minres
= atoms
->atom
[dih
[ndih
].a
[0]].resind
;
928 for (m
= 1; m
< 4; m
++)
930 minres
= std::min(minres
, atoms
->atom
[dih
[ndih
].a
[m
]].resind
);
931 maxres
= std::max(maxres
, atoms
->atom
[dih
[ndih
].a
[m
]].resind
);
933 res
= 2*minres
-maxres
;
936 res
+= maxres
-minres
;
937 get_atomnames_min(4, anm
, res
, atoms
, dih
[ndih
].a
);
938 hbdih
= &hb
[res
].rb
[ebtsPDIHS
];
939 for (n
= 0; (n
< hbdih
->nb
); n
++)
942 for (m
= 0; m
< 2; m
++)
945 ((strcmp(anm
[3*m
], hbdih
->b
[n
].ai()) == 0) &&
946 (strcmp(anm
[1+m
], hbdih
->b
[n
].aj()) == 0) &&
947 (strcmp(anm
[2-m
], hbdih
->b
[n
].ak()) == 0) &&
948 (strcmp(anm
[3-3*m
], hbdih
->b
[n
].al()) == 0)));
952 set_p_string(&dih
[ndih
], hbdih
->b
[n
].s
);
953 /* Mark that we found a match for this entry */
954 hbdih
->b
[n
].match
= TRUE
;
956 /* Set the last parameter to be able to see
957 if the dihedral was in the rtp list.
959 dih
[ndih
].c
[MAXFORCEPARAM
-1] = DIHEDRAL_WAS_SET_IN_RTP
;
962 /* Set the next direct in case the rtp contains
963 multiple entries for this dihedral.
974 for (m
= 0; m
< MAXFORCEPARAM
; m
++)
976 dih
[ndih
].c
[m
] = NOTSET
;
981 while (res
< maxres
);
994 for (m
= 0; m
< MAXFORCEPARAM
; m
++)
996 dih
[ndih
].c
[m
] = NOTSET
;
998 set_p_string(&(dih
[ndih
]), "");
1002 nbd
= nb_dist(nnb
, i
, l1
);
1005 fprintf(debug
, "Distance (%d-%d) = %d\n", i
+1, l1
+1, nbd
);
1009 i1
= std::min(i
, l1
);
1010 i2
= std::max(i
, l1
);
1012 for (m
= 0; m
< excls
[i1
].nr
; m
++)
1014 bExcl
= bExcl
|| excls
[i1
].e
[m
] == i2
;
1018 if (rtp
[0].bGenerateHH14Interactions
||
1019 !(is_hydro(atoms
, i1
) && is_hydro(atoms
, i2
)))
1024 srenew(pai
, maxpai
);
1026 pai
[npai
].ai() = i1
;
1027 pai
[npai
].aj() = i2
;
1028 pai
[npai
].c0() = NOTSET
;
1029 pai
[npai
].c1() = NOTSET
;
1030 set_p_string(&(pai
[npai
]), "");
1045 /* The above approach is great in that we double-check that e.g. an angle
1046 * really corresponds to three atoms connected by bonds, but this is not
1047 * generally true. Go through the angle and dihedral hackblocks to add
1048 * entries that we have not yet marked as matched when going through bonds.
1050 for (i
= 0; i
< atoms
->nres
; i
++)
1052 /* Add remaining angles from hackblock */
1053 hbang
= &hb
[i
].rb
[ebtsANGLES
];
1054 for (j
= 0; j
< hbang
->nb
; j
++)
1056 if (hbang
->b
[j
].match
== TRUE
)
1058 /* We already used this entry, continue to the next */
1061 /* Hm - entry not used, let's see if we can find all atoms */
1065 srenew(ang
, maxang
);
1068 for (k
= 0; k
< 3 && bFound
; k
++)
1070 p
= hbang
->b
[j
].a
[k
];
1077 else if (p
[0] == '+')
1082 ang
[nang
].a
[k
] = search_res_atom(p
, res
, atoms
, "angle", TRUE
);
1083 bFound
= (ang
[nang
].a
[k
] != NO_ATID
);
1085 ang
[nang
].c0() = NOTSET
;
1086 ang
[nang
].c1() = NOTSET
;
1090 set_p_string(&(ang
[nang
]), hbang
->b
[j
].s
);
1091 hbang
->b
[j
].match
= TRUE
;
1092 /* Incrementing nang means we save this angle */
1097 /* Add remaining dihedrals from hackblock */
1098 hbdih
= &hb
[i
].rb
[ebtsPDIHS
];
1099 for (j
= 0; j
< hbdih
->nb
; j
++)
1101 if (hbdih
->b
[j
].match
== TRUE
)
1103 /* We already used this entry, continue to the next */
1106 /* Hm - entry not used, let's see if we can find all atoms */
1110 srenew(dih
, maxdih
);
1113 for (k
= 0; k
< 4 && bFound
; k
++)
1115 p
= hbdih
->b
[j
].a
[k
];
1122 else if (p
[0] == '+')
1127 dih
[ndih
].a
[k
] = search_res_atom(p
, res
, atoms
, "dihedral", TRUE
);
1128 bFound
= (dih
[ndih
].a
[k
] != NO_ATID
);
1130 for (m
= 0; m
< MAXFORCEPARAM
; m
++)
1132 dih
[ndih
].c
[m
] = NOTSET
;
1137 set_p_string(&(dih
[ndih
]), hbdih
->b
[j
].s
);
1138 hbdih
->b
[j
].match
= TRUE
;
1139 /* Incrementing ndih means we save this dihedral */
1146 /* Sort angles with respect to j-i-k (middle atom first) */
1149 qsort(ang
, nang
, (size_t)sizeof(ang
[0]), acomp
);
1152 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
1155 qsort(dih
, ndih
, (size_t)sizeof(dih
[0]), dcomp
);
1158 /* Sort the pairs */
1161 qsort(pai
, npai
, (size_t)sizeof(pai
[0]), pcomp
);
1165 /* Remove doubles, could occur in 6-rings, such as phenyls,
1166 maybe one does not want this when fudgeQQ < 1.
1168 fprintf(stderr
, "Before cleaning: %d pairs\n", npai
);
1169 rm2par(pai
, &npai
, preq
);
1172 /* Get the impropers from the database */
1173 nimproper
= get_impropers(atoms
, hb
, &improper
, bAllowMissing
);
1175 /* Sort the impropers */
1176 sort_id(nimproper
, improper
);
1180 fprintf(stderr
, "Before cleaning: %d dihedrals\n", ndih
);
1181 clean_dih(dih
, &ndih
, improper
, nimproper
, atoms
,
1182 rtp
[0].bKeepAllGeneratedDihedrals
,
1183 rtp
[0].bRemoveDihedralIfWithImproper
);
1186 /* Now we have unique lists of angles and dihedrals
1187 * Copy them into the destination struct
1189 cppar(ang
, nang
, plist
, F_ANGLES
);
1190 cppar(dih
, ndih
, plist
, F_PDIHS
);
1191 cppar(improper
, nimproper
, plist
, F_IDIHS
);
1192 cppar(pai
, npai
, plist
, F_LJ14
);
1194 /* Remove all exclusions which are within nrexcl */
1195 clean_excls(nnb
, rtp
[0].nrexcl
, excls
);