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! */
51 #include "gpp_nextnb.h"
54 #include "gmx_fatal.h"
59 typedef gmx_bool (*peq)(t_param *p1, t_param *p2);
61 static int acomp(const void *a1, const void *a2)
68 if ((ac=(p1->AJ-p2->AJ))!=0)
70 else if ((ac=(p1->AI-p2->AI))!=0)
73 return (p1->AK-p2->AK);
76 static int pcomp(const void *a1, const void *a2)
83 if ((pc=(p1->AI-p2->AI))!=0)
86 return (p1->AJ-p2->AJ);
89 static int dcomp(const void *d1, const void *d2)
96 /* First sort by J & K (the two central) atoms */
97 if ((dc=(p1->AJ-p2->AJ))!=0)
99 else if ((dc=(p1->AK-p2->AK))!=0)
101 /* Then make sure to put rtp dihedrals before generated ones */
102 else if (p1->c[MAXFORCEPARAM-1]==0 && p2->c[MAXFORCEPARAM-1]==NOTSET)
104 else if (p1->c[MAXFORCEPARAM-1]==NOTSET && p2->c[MAXFORCEPARAM-1]==0)
106 /* Finally, sort by I and J (two outer) atoms */
107 else if ((dc=(p1->AI-p2->AI))!=0)
110 return (p1->AL-p2->AL);
114 static gmx_bool aeq(t_param *p1, t_param *p2)
118 else if (((p1->AI==p2->AI) && (p1->AK==p2->AK)) ||
119 ((p1->AI==p2->AK) && (p1->AK==p2->AI)))
127 static gmx_bool deq(t_param *p1, t_param *p2)
129 if (((p1->AJ==p2->AJ) && (p1->AK==p2->AK)) ||
130 ((p1->AJ==p2->AK) && (p1->AK==p2->AJ)))
137 static gmx_bool remove_dih(t_param *p, int i, int np)
138 /* check if dihedral p[i] should be removed */
143 if (p[i].c[MAXFORCEPARAM-1]==NOTSET) {
145 bRem = deq(&p[i],&p[i-1]);
148 /* also remove p[i] if there is a dihedral on the same bond
149 which has parameters set */
151 while (!bRem && (j<np) && deq(&p[i],&p[j])) {
152 bRem = (p[j].c[MAXFORCEPARAM-1] != NOTSET);
161 static gmx_bool preq(t_param *p1, t_param *p2)
163 if ((p1->AI==p2->AI) && (p1->AJ==p2->AJ))
169 static void rm2par(t_param p[], int *np, peq eq)
180 for(i=1; (i<(*np)); i++)
181 if (!eq(&p[i],&p[i-1]))
183 /* Index now holds pointers to all the non-equal params,
184 * this only works when p is sorted of course
186 for(i=0; (i<nind); i++) {
187 for(j=0; (j<MAXATOMLIST); j++)
188 p[i].a[j]=p[index[i]].a[j];
189 for(j=0; (j<MAXFORCEPARAM); j++)
190 p[i].c[j]=p[index[i]].c[j];
191 if (p[index[i]].a[0] == p[index[i]].a[1]) {
194 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
195 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
196 p[i].a[0],p[i].a[1],p[i].a[2],p[i].a[3]);
198 } else if (index[i] > i) {
199 /* Copy the string only if it comes from somewhere else
200 * otherwise we will end up copying a random (newly freed) pointer.
201 * Since the index is sorted we only have to test for index[i] > i.
203 strcpy(p[i].s,p[index[i]].s);
211 static void cppar(t_param p[], int np, t_params plist[], int ftype)
222 for(i=0; (i<np); i++) {
223 for(j=0; (j<nral); j++)
224 ps->param[ps->nr].a[j] = p[i].a[j];
225 for(j=0; (j<nrfp); j++)
226 ps->param[ps->nr].c[j] = p[i].c[j];
227 for(j=0; (j<MAXSLEN); j++)
228 ps->param[ps->nr].s[j] = p[i].s[j];
233 static void cpparam(t_param *dest, t_param *src)
237 for(j=0; (j<MAXATOMLIST); j++)
238 dest->a[j] = src->a[j];
239 for(j=0; (j<MAXFORCEPARAM); j++)
240 dest->c[j] = src->c[j];
241 for(j=0; (j<MAXSLEN); j++)
242 dest->s[j] = src->s[j];
245 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
251 for(j=0; (j<MAXFORCEPARAM); j++)
260 static int int_comp(const void *a,const void *b)
262 return (*(int *)a) - (*(int *)b);
265 static int atom_id_comp(const void *a,const void *b)
267 return (*(atom_id *)a) - (*(atom_id *)b);
270 static int eq_imp(atom_id a1[],atom_id a2[])
275 for(j=0; (j<4); j++) {
279 qsort(b1,4,(size_t)sizeof(b1[0]),int_comp);
280 qsort(b2,4,(size_t)sizeof(b2[0]),int_comp);
289 static gmx_bool ideq(t_param *p1, t_param *p2)
291 return eq_imp(p1->a,p2->a);
294 static int idcomp(const void *a,const void *b)
301 if ((d=(pa->a[0]-pb->a[0])) != 0)
303 else if ((d=(pa->a[3]-pb->a[3])) != 0)
305 else if ((d=(pa->a[1]-pb->a[1])) != 0)
308 return (int) (pa->a[2]-pb->a[2]);
311 static void sort_id(int nr,t_param ps[])
315 /* First swap order of atoms around if necessary */
316 for(i=0; (i<nr); i++) {
317 if (ps[i].a[3] < ps[i].a[0]) {
318 tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
319 tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
324 qsort(ps,nr,(size_t)sizeof(ps[0]),idcomp);
327 static void dump_param(FILE *fp,char *title,int n,t_param ps[])
331 fprintf(fp,"%s: %d entries\n",title,n);
332 for(i=0; (i<n); i++) {
333 fprintf(fp,"%3d: A=[ ",i);
334 for(j=0; (j<MAXATOMLIST); j++)
335 fprintf(fp," %5d",ps[i].a[j]);
337 for(j=0; (j<MAXFORCEPARAM); j++)
338 fprintf(fp," %10.5e",ps[i].c[j]);
343 static int n_hydro(atom_id a[],char ***atomname)
348 for(i=0; (i<4); i+=3) {
349 aname=*atomname[a[i]];
350 c0=toupper(aname[0]);
353 else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9')) {
354 c1=toupper(aname[1]);
362 static void clean_dih(t_param *dih, int *ndih,t_param idih[],int nidih,
363 t_atoms *atoms,gmx_bool bAlldih, gmx_bool bRemoveDih)
367 gmx_bool bIsSet,bKeep;
372 fprintf(stderr,"Keeping all generated dihedrals\n");
374 for(i=0; i<nind; i++)
378 /* Make an index of all dihedrals over each bond */
380 for(i=0; i<*ndih; i++)
381 if (!remove_dih(dih,i,*ndih))
386 /* if we don't want all dihedrals, we need to select the ones with the
391 for(i=0; i<nind; i++) {
392 bIsSet = (dih[index[i]].c[MAXFORCEPARAM-1] != NOTSET);
394 if (!bIsSet && bRemoveDih)
395 /* remove the dihedral if there is an improper on the same bond */
396 for(j=0; (j<nidih) && bKeep; j++)
397 bKeep = !deq(&dih[index[i]],&idih[j]);
400 /* Now select the "fittest" dihedral:
401 * the one with as few hydrogens as possible
404 /* Best choice to get dihedral from */
406 if (!bAlldih && !bIsSet) {
407 /* Minimum number of hydrogens for i and l atoms */
409 for(l=index[i]; (l<index[i+1]) && deq(&dih[index[i]],&dih[l]); l++) {
410 if ((nh=n_hydro(dih[l].a,atoms->atomname)) < minh) {
419 cpparam(&(dih[k]),&dih[bestl]);
424 for (i=k; i<*ndih; i++)
431 static int get_impropers(t_atoms *atoms,t_hackblock hb[],t_param **idih,
432 gmx_bool bAllowMissing)
437 int nidih,i,j,k,r,start,ninc,nalloc;
438 atom_id ai[MAXATOMLIST];
445 /* Add all the impropers from the residue database to the list. */
449 for(i=0; (i<atoms->nres); i++) {
450 idihs=&hb[i].rb[ebtsIDIHS];
451 for(j=0; (j<idihs->nb); j++) {
453 for(k=0; (k<4) && !bStop; k++) {
454 ai[k] = search_atom(idihs->b[j].a[k],start,
455 atoms->nr,atoms->atom,atoms->atomname,
456 "improper",bAllowMissing);
457 if (ai[k] == NO_ATID)
461 if (nidih == nalloc) {
463 srenew(*idih,nalloc);
466 set_p(&((*idih)[nidih]),ai,NULL,idihs->b[j].s);
470 while ((start<atoms->nr) && (atoms->atom[start].resind == i))
478 static int nb_dist(t_nextnb *nnb,int ai,int aj)
488 nrexcl=nnb->nrexcl[ai];
489 for(nre=1; (nre < nnb->nrex); nre++) {
491 for(nrx=0; (nrx < nrexcl[nre]); nrx++) {
492 if ((aj == a[nrx]) && (NRE == -1))
499 gmx_bool is_hydro(t_atoms *atoms,int ai)
501 return ((*(atoms->atomname[ai]))[0] == 'H');
504 static void get_atomnames_min(int n,char **anm,
505 int resind,t_atoms *atoms,atom_id *a)
509 /* Assume ascending residue numbering */
511 if (atoms->atom[a[m]].resind < resind)
513 else if (atoms->atom[a[m]].resind > resind)
517 strcat(anm[m],*(atoms->atomname[a[m]]));
521 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
522 gmx_bool bAllowMissing)
525 atom_id a,astart,i1,i2,itmp;
531 for(a=0; a<atoms->nr; a++) {
532 r = atoms->atom[a].resind;
533 if (a==atoms->nr-1 || atoms->atom[a+1].resind != r) {
534 hbexcl = &hb[r].rb[ebtsEXCLS];
536 for(e=0; e<hbexcl->nb; e++) {
537 anm = hbexcl->b[e].a[0];
538 i1 = search_atom(anm,astart,atoms->nr,atoms->atom,atoms->atomname,
539 "exclusion",bAllowMissing);
540 anm = hbexcl->b[e].a[1];
541 i2 = search_atom(anm,astart,atoms->nr,atoms->atom,atoms->atomname,
542 "exclusion",bAllowMissing);
543 if (i1!=NO_ATID && i2!=NO_ATID) {
549 srenew(excls[i1].e,excls[i1].nr+1);
550 excls[i1].e[excls[i1].nr] = i2;
559 for(a=0; a<atoms->nr; a++)
561 qsort(excls[a].e,excls[a].nr,(size_t)sizeof(atom_id),atom_id_comp);
564 static void remove_excl(t_excls *excls, int remove)
568 for(i=remove+1; i<excls->nr; i++)
569 excls->e[i-1] = excls->e[i];
574 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
576 int i,j,j1,k,k1,l,l1,m,n,e;
580 /* extract all i-j-k-l neighbours from nnb struct */
581 for(i=0; (i<nnb->nr); i++) {
582 /* For all particles */
585 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
586 /* For all first neighbours */
589 for(e=0; e<excl->nr; e++)
590 if (excl->e[e] == j1)
594 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
595 /* For all first neighbours of j1 */
598 for(e=0; e<excl->nr; e++)
599 if (excl->e[e] == k1)
603 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
604 /* For all first neighbours of k1 */
607 for(e=0; e<excl->nr; e++)
608 if (excl->e[e] == l1)
616 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
618 int i,j,j1,k,k1,l,l1,m,n,e,N;
621 for(N=1; (N<min(nrexcl,nnb->nrex)); N++) {
622 /* extract all i-j-k-l neighbours from nnb struct */
623 for(i=0; (i<nnb->nr); i++) {
624 /* For all particles */
627 excl->nr += nnb->nrexcl[i][N];
628 srenew(excl->e,excl->nr);
629 for(j=0; (j<nnb->nrexcl[i][N]); j++)
630 /* For all first neighbours */
631 if (nnb->a[i][N][j] != i)
632 excl->e[n++] = nnb->a[i][N][j];
637 void gen_pad(t_nextnb *nnb, t_atoms *atoms, int nrexcl, gmx_bool bH14,
638 t_params plist[], t_excls excls[], t_hackblock hb[],
639 gmx_bool bAlldih, gmx_bool bRemoveDih, gmx_bool bAllowMissing)
641 t_param *ang,*dih,*pai,*idih;
642 t_rbondeds *hbang, *hbdih;
644 int res,minres,maxres;
645 int i,j,j1,k,k1,l,l1,m,n,i1,i2;
646 int ninc,maxang,maxdih,maxpai;
647 int nang,ndih,npai,nidih,nbd;
649 gmx_bool bFound,bExcl;
652 /* These are the angles, dihedrals and pairs that we generate
653 * from the bonds. The ones that are already there from the rtp file
660 maxang = maxdih = maxpai = ninc;
670 gen_excls(atoms,excls,hb,bAllowMissing);
672 /* extract all i-j-k-l neighbours from nnb struct */
673 for(i=0; (i<nnb->nr); i++)
674 /* For all particles */
675 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
676 /* For all first neighbours */
678 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
679 /* For all first neighbours of j1 */
682 /* Generate every angle only once */
684 if (nang == maxang) {
693 set_p_string(&(ang[nang]),"");
695 minres = atoms->atom[ang[nang].a[0]].resind;
698 minres = min(minres,atoms->atom[ang[nang].a[m]].resind);
699 maxres = max(maxres,atoms->atom[ang[nang].a[m]].resind);
701 res = 2*minres-maxres;
703 res += maxres-minres;
704 get_atomnames_min(3,anm,res,atoms,ang[nang].a);
705 hbang=&hb[res].rb[ebtsANGLES];
706 for(l=0; (l<hbang->nb); l++) {
707 if (strcmp(anm[1],hbang->b[l].AJ)==0) {
711 ((strcmp(anm[m],hbang->b[l].AI)==0) &&
712 (strcmp(anm[2-m],hbang->b[l].AK)==0)));
714 set_p_string(&(ang[nang]),hbang->b[l].s);
718 } while (res < maxres);
722 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
725 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
726 /* For all first neighbours of k1 */
728 if ((l1 != i) && (l1 != j1)) {
729 if (ndih == maxdih) {
737 for (m=0; m<MAXFORCEPARAM; m++)
738 dih[ndih].c[m]=NOTSET;
739 set_p_string(&(dih[ndih]),"");
742 minres = atoms->atom[dih[ndih].a[0]].resind;
745 minres = min(minres,atoms->atom[dih[ndih].a[m]].resind);
746 maxres = max(maxres,atoms->atom[dih[ndih].a[m]].resind);
748 res = 2*minres-maxres;
750 res += maxres-minres;
751 get_atomnames_min(4,anm,res,atoms,dih[ndih].a);
752 hbdih=&hb[res].rb[ebtsPDIHS];
753 for(n=0; (n<hbdih->nb); n++) {
757 ((strcmp(anm[3*m], hbdih->b[n].AI)==0) &&
758 (strcmp(anm[1+m], hbdih->b[n].AJ)==0) &&
759 (strcmp(anm[2-m], hbdih->b[n].AK)==0) &&
760 (strcmp(anm[3-3*m],hbdih->b[n].AL)==0)));
762 set_p_string(&dih[ndih],hbdih->b[n].s);
764 /* Set the last parameter to be able to see
765 if the dihedral was in the rtp list.
767 dih[ndih].c[MAXFORCEPARAM-1] = 0;
770 /* Set the next direct in case the rtp contains
771 multiple entries for this dihedral.
773 if (ndih == maxdih) {
781 for (m=0; m<MAXFORCEPARAM; m++)
782 dih[ndih].c[m]=NOTSET;
785 } while (res < maxres);
788 if (ndih == maxdih) {
796 for (m=0; m<MAXFORCEPARAM; m++)
797 dih[ndih].c[m]=NOTSET;
798 set_p_string(&(dih[ndih]),"");
802 nbd=nb_dist(nnb,i,l1);
804 fprintf(debug,"Distance (%d-%d) = %d\n",i+1,l1+1,nbd);
809 for(m=0; m<excls[i1].nr; m++)
810 bExcl = bExcl || excls[i1].e[m]==i2;
812 if (bH14 || !(is_hydro(atoms,i1) && is_hydro(atoms,i2))) {
813 if (npai == maxpai) {
821 set_p_string(&(pai[npai]),"");
833 /* Sort angles with respect to j-i-k (middle atom first) */
835 qsort(ang,nang,(size_t)sizeof(ang[0]),acomp);
837 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
839 qsort(dih,ndih,(size_t)sizeof(dih[0]),dcomp);
843 qsort(pai,npai,(size_t)sizeof(pai[0]),pcomp);
845 /* Remove doubles, could occur in 6-rings, such as phenyls,
846 maybe one does not want this when fudgeQQ < 1.
848 fprintf(stderr,"Before cleaning: %d pairs\n",npai);
849 rm2par(pai,&npai,preq);
852 /* Get the impropers from the database */
853 nidih = get_impropers(atoms,hb,&idih,bAllowMissing);
855 /* Sort the impropers */
859 /* Remove dihedrals which are impropers
860 and when bAlldih is not set remove multiple dihedrals over one bond.
862 fprintf(stderr,"Before cleaning: %d dihedrals\n",ndih);
863 clean_dih(dih,&ndih,idih,nidih,atoms,bAlldih,bRemoveDih);
866 /* Now we have unique lists of angles and dihedrals
867 * Copy them into the destination struct
869 cppar(ang, nang, plist,F_ANGLES);
870 cppar(dih, ndih, plist,F_PDIHS);
871 cppar(idih,nidih,plist,F_IDIHS);
872 cppar(pai, npai, plist,F_LJ14);
874 /* Remove all exclusions which are within nrexcl */
875 clean_excls(nnb,nrexcl,excls);