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);
113 static gmx_bool deq(t_param *p1, t_param *p2)
115 if (((p1->AJ==p2->AJ) && (p1->AK==p2->AK)) ||
116 ((p1->AJ==p2->AK) && (p1->AK==p2->AJ)))
123 static gmx_bool remove_dih(t_param *p, int i, int np)
124 /* check if dihedral p[i] should be removed */
129 if (p[i].c[MAXFORCEPARAM-1]==NOTSET) {
131 bRem = deq(&p[i],&p[i-1]);
134 /* also remove p[i] if there is a dihedral on the same bond
135 which has parameters set */
137 while (!bRem && (j<np) && deq(&p[i],&p[j])) {
138 bRem = (p[j].c[MAXFORCEPARAM-1] != NOTSET);
147 static gmx_bool preq(t_param *p1, t_param *p2)
149 if ((p1->AI==p2->AI) && (p1->AJ==p2->AJ))
155 static void rm2par(t_param p[], int *np, peq eq)
166 for(i=1; (i<(*np)); i++)
167 if (!eq(&p[i],&p[i-1]))
169 /* Index now holds pointers to all the non-equal params,
170 * this only works when p is sorted of course
172 for(i=0; (i<nind); i++) {
173 for(j=0; (j<MAXATOMLIST); j++)
174 p[i].a[j]=p[index[i]].a[j];
175 for(j=0; (j<MAXFORCEPARAM); j++)
176 p[i].c[j]=p[index[i]].c[j];
177 if (p[index[i]].a[0] == p[index[i]].a[1]) {
180 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
181 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
182 p[i].a[0],p[i].a[1],p[i].a[2],p[i].a[3]);
184 } else if (index[i] > i) {
185 /* Copy the string only if it comes from somewhere else
186 * otherwise we will end up copying a random (newly freed) pointer.
187 * Since the index is sorted we only have to test for index[i] > i.
189 strcpy(p[i].s,p[index[i]].s);
197 static void cppar(t_param p[], int np, t_params plist[], int ftype)
208 for(i=0; (i<np); i++) {
209 for(j=0; (j<nral); j++)
210 ps->param[ps->nr].a[j] = p[i].a[j];
211 for(j=0; (j<nrfp); j++)
212 ps->param[ps->nr].c[j] = p[i].c[j];
213 for(j=0; (j<MAXSLEN); j++)
214 ps->param[ps->nr].s[j] = p[i].s[j];
219 static void cpparam(t_param *dest, t_param *src)
223 for(j=0; (j<MAXATOMLIST); j++)
224 dest->a[j] = src->a[j];
225 for(j=0; (j<MAXFORCEPARAM); j++)
226 dest->c[j] = src->c[j];
227 for(j=0; (j<MAXSLEN); j++)
228 dest->s[j] = src->s[j];
231 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
237 for(j=0; (j<MAXFORCEPARAM); j++)
246 static int int_comp(const void *a,const void *b)
248 return (*(int *)a) - (*(int *)b);
251 static int atom_id_comp(const void *a,const void *b)
253 return (*(atom_id *)a) - (*(atom_id *)b);
256 static int eq_imp(atom_id a1[],atom_id a2[])
261 for(j=0; (j<4); j++) {
265 qsort(b1,4,(size_t)sizeof(b1[0]),int_comp);
266 qsort(b2,4,(size_t)sizeof(b2[0]),int_comp);
275 static int idcomp(const void *a,const void *b)
282 if ((d=(pa->a[0]-pb->a[0])) != 0)
284 else if ((d=(pa->a[3]-pb->a[3])) != 0)
286 else if ((d=(pa->a[1]-pb->a[1])) != 0)
289 return (int) (pa->a[2]-pb->a[2]);
292 static void sort_id(int nr,t_param ps[])
296 /* First swap order of atoms around if necessary */
297 for(i=0; (i<nr); i++) {
298 if (ps[i].a[3] < ps[i].a[0]) {
299 tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
300 tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
305 qsort(ps,nr,(size_t)sizeof(ps[0]),idcomp);
308 static int n_hydro(atom_id a[],char ***atomname)
313 for(i=0; (i<4); i+=3) {
314 aname=*atomname[a[i]];
315 c0=toupper(aname[0]);
318 else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9')) {
319 c1=toupper(aname[1]);
327 static void clean_dih(t_param *dih, int *ndih,t_param idih[],int nidih,
328 t_atoms *atoms,gmx_bool bAlldih, gmx_bool bRemoveDih)
332 gmx_bool bIsSet,bKeep;
337 fprintf(stderr,"Keeping all generated dihedrals\n");
339 for(i=0; i<nind; i++)
343 /* Make an index of all dihedrals over each bond */
345 for(i=0; i<*ndih; i++)
346 if (!remove_dih(dih,i,*ndih))
351 /* if we don't want all dihedrals, we need to select the ones with the
356 for(i=0; i<nind; i++) {
357 bIsSet = (dih[index[i]].c[MAXFORCEPARAM-1] != NOTSET);
359 if (!bIsSet && bRemoveDih)
360 /* remove the dihedral if there is an improper on the same bond */
361 for(j=0; (j<nidih) && bKeep; j++)
362 bKeep = !deq(&dih[index[i]],&idih[j]);
365 /* Now select the "fittest" dihedral:
366 * the one with as few hydrogens as possible
369 /* Best choice to get dihedral from */
371 if (!bAlldih && !bIsSet) {
372 /* Minimum number of hydrogens for i and l atoms */
374 for(l=index[i]; (l<index[i+1]) && deq(&dih[index[i]],&dih[l]); l++) {
375 if ((nh=n_hydro(dih[l].a,atoms->atomname)) < minh) {
384 cpparam(&(dih[k]),&dih[bestl]);
389 for (i=k; i<*ndih; i++)
396 static int get_impropers(t_atoms *atoms,t_hackblock hb[],t_param **idih,
397 gmx_bool bAllowMissing)
402 int nidih,i,j,k,r,start,ninc,nalloc;
403 atom_id ai[MAXATOMLIST];
410 /* Add all the impropers from the residue database to the list. */
414 for(i=0; (i<atoms->nres); i++) {
415 idihs=&hb[i].rb[ebtsIDIHS];
416 for(j=0; (j<idihs->nb); j++) {
418 for(k=0; (k<4) && !bStop; k++) {
419 ai[k] = search_atom(idihs->b[j].a[k],start,
421 "improper",bAllowMissing);
422 if (ai[k] == NO_ATID)
426 if (nidih == nalloc) {
428 srenew(*idih,nalloc);
431 set_p(&((*idih)[nidih]),ai,NULL,idihs->b[j].s);
435 while ((start<atoms->nr) && (atoms->atom[start].resind == i))
443 static int nb_dist(t_nextnb *nnb,int ai,int aj)
453 nrexcl=nnb->nrexcl[ai];
454 for(nre=1; (nre < nnb->nrex); nre++) {
456 for(nrx=0; (nrx < nrexcl[nre]); nrx++) {
457 if ((aj == a[nrx]) && (NRE == -1))
464 gmx_bool is_hydro(t_atoms *atoms,int ai)
466 return ((*(atoms->atomname[ai]))[0] == 'H');
469 static void get_atomnames_min(int n,char **anm,
470 int resind,t_atoms *atoms,atom_id *a)
474 /* Assume ascending residue numbering */
476 if (atoms->atom[a[m]].resind < resind)
478 else if (atoms->atom[a[m]].resind > resind)
482 strcat(anm[m],*(atoms->atomname[a[m]]));
486 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
487 gmx_bool bAllowMissing)
490 atom_id a,astart,i1,i2,itmp;
496 for(a=0; a<atoms->nr; a++) {
497 r = atoms->atom[a].resind;
498 if (a==atoms->nr-1 || atoms->atom[a+1].resind != r) {
499 hbexcl = &hb[r].rb[ebtsEXCLS];
501 for(e=0; e<hbexcl->nb; e++) {
502 anm = hbexcl->b[e].a[0];
503 i1 = search_atom(anm,astart,atoms,
504 "exclusion",bAllowMissing);
505 anm = hbexcl->b[e].a[1];
506 i2 = search_atom(anm,astart,atoms,
507 "exclusion",bAllowMissing);
508 if (i1!=NO_ATID && i2!=NO_ATID) {
514 srenew(excls[i1].e,excls[i1].nr+1);
515 excls[i1].e[excls[i1].nr] = i2;
524 for(a=0; a<atoms->nr; a++)
526 qsort(excls[a].e,excls[a].nr,(size_t)sizeof(atom_id),atom_id_comp);
529 static void remove_excl(t_excls *excls, int remove)
533 for(i=remove+1; i<excls->nr; i++)
534 excls->e[i-1] = excls->e[i];
539 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
541 int i,j,j1,k,k1,l,l1,m,n,e;
545 /* extract all i-j-k-l neighbours from nnb struct */
546 for(i=0; (i<nnb->nr); i++) {
547 /* For all particles */
550 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
551 /* For all first neighbours */
554 for(e=0; e<excl->nr; e++)
555 if (excl->e[e] == j1)
559 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
560 /* For all first neighbours of j1 */
563 for(e=0; e<excl->nr; e++)
564 if (excl->e[e] == k1)
568 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
569 /* For all first neighbours of k1 */
572 for(e=0; e<excl->nr; e++)
573 if (excl->e[e] == l1)
581 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
583 int i,j,j1,k,k1,l,l1,m,n,e,N;
586 for(N=1; (N<min(nrexcl,nnb->nrex)); N++) {
587 /* extract all i-j-k-l neighbours from nnb struct */
588 for(i=0; (i<nnb->nr); i++) {
589 /* For all particles */
592 excl->nr += nnb->nrexcl[i][N];
593 srenew(excl->e,excl->nr);
594 for(j=0; (j<nnb->nrexcl[i][N]); j++)
595 /* For all first neighbours */
596 if (nnb->a[i][N][j] != i)
597 excl->e[n++] = nnb->a[i][N][j];
602 void gen_pad(t_nextnb *nnb, t_atoms *atoms, int nrexcl, gmx_bool bH14,
603 t_params plist[], t_excls excls[], t_hackblock hb[],
604 gmx_bool bAlldih, gmx_bool bRemoveDih, gmx_bool bAllowMissing)
606 t_param *ang,*dih,*pai,*idih;
607 t_rbondeds *hbang, *hbdih;
609 int res,minres,maxres;
610 int i,j,j1,k,k1,l,l1,m,n,i1,i2;
611 int ninc,maxang,maxdih,maxpai;
612 int nang,ndih,npai,nidih,nbd;
614 gmx_bool bFound,bExcl;
617 /* These are the angles, dihedrals and pairs that we generate
618 * from the bonds. The ones that are already there from the rtp file
625 maxang = maxdih = maxpai = ninc;
635 gen_excls(atoms,excls,hb,bAllowMissing);
637 /* extract all i-j-k-l neighbours from nnb struct */
638 for(i=0; (i<nnb->nr); i++)
639 /* For all particles */
640 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
641 /* For all first neighbours */
643 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
644 /* For all first neighbours of j1 */
647 /* Generate every angle only once */
649 if (nang == maxang) {
658 set_p_string(&(ang[nang]),"");
660 minres = atoms->atom[ang[nang].a[0]].resind;
663 minres = min(minres,atoms->atom[ang[nang].a[m]].resind);
664 maxres = max(maxres,atoms->atom[ang[nang].a[m]].resind);
666 res = 2*minres-maxres;
668 res += maxres-minres;
669 get_atomnames_min(3,anm,res,atoms,ang[nang].a);
670 hbang=&hb[res].rb[ebtsANGLES];
671 for(l=0; (l<hbang->nb); l++) {
672 if (strcmp(anm[1],hbang->b[l].AJ)==0) {
676 ((strcmp(anm[m],hbang->b[l].AI)==0) &&
677 (strcmp(anm[2-m],hbang->b[l].AK)==0)));
679 set_p_string(&(ang[nang]),hbang->b[l].s);
683 } while (res < maxres);
687 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
690 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
691 /* For all first neighbours of k1 */
693 if ((l1 != i) && (l1 != j1)) {
694 if (ndih == maxdih) {
702 for (m=0; m<MAXFORCEPARAM; m++)
703 dih[ndih].c[m]=NOTSET;
704 set_p_string(&(dih[ndih]),"");
707 minres = atoms->atom[dih[ndih].a[0]].resind;
710 minres = min(minres,atoms->atom[dih[ndih].a[m]].resind);
711 maxres = max(maxres,atoms->atom[dih[ndih].a[m]].resind);
713 res = 2*minres-maxres;
715 res += maxres-minres;
716 get_atomnames_min(4,anm,res,atoms,dih[ndih].a);
717 hbdih=&hb[res].rb[ebtsPDIHS];
718 for(n=0; (n<hbdih->nb); n++) {
722 ((strcmp(anm[3*m], hbdih->b[n].AI)==0) &&
723 (strcmp(anm[1+m], hbdih->b[n].AJ)==0) &&
724 (strcmp(anm[2-m], hbdih->b[n].AK)==0) &&
725 (strcmp(anm[3-3*m],hbdih->b[n].AL)==0)));
727 set_p_string(&dih[ndih],hbdih->b[n].s);
729 /* Set the last parameter to be able to see
730 if the dihedral was in the rtp list.
732 dih[ndih].c[MAXFORCEPARAM-1] = 0;
735 /* Set the next direct in case the rtp contains
736 multiple entries for this dihedral.
738 if (ndih == maxdih) {
746 for (m=0; m<MAXFORCEPARAM; m++)
747 dih[ndih].c[m]=NOTSET;
750 } while (res < maxres);
753 if (ndih == maxdih) {
761 for (m=0; m<MAXFORCEPARAM; m++)
762 dih[ndih].c[m]=NOTSET;
763 set_p_string(&(dih[ndih]),"");
767 nbd=nb_dist(nnb,i,l1);
769 fprintf(debug,"Distance (%d-%d) = %d\n",i+1,l1+1,nbd);
774 for(m=0; m<excls[i1].nr; m++)
775 bExcl = bExcl || excls[i1].e[m]==i2;
777 if (bH14 || !(is_hydro(atoms,i1) && is_hydro(atoms,i2))) {
778 if (npai == maxpai) {
786 set_p_string(&(pai[npai]),"");
798 /* Sort angles with respect to j-i-k (middle atom first) */
800 qsort(ang,nang,(size_t)sizeof(ang[0]),acomp);
802 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
804 qsort(dih,ndih,(size_t)sizeof(dih[0]),dcomp);
808 qsort(pai,npai,(size_t)sizeof(pai[0]),pcomp);
810 /* Remove doubles, could occur in 6-rings, such as phenyls,
811 maybe one does not want this when fudgeQQ < 1.
813 fprintf(stderr,"Before cleaning: %d pairs\n",npai);
814 rm2par(pai,&npai,preq);
817 /* Get the impropers from the database */
818 nidih = get_impropers(atoms,hb,&idih,bAllowMissing);
820 /* Sort the impropers */
824 /* Remove dihedrals which are impropers
825 and when bAlldih is not set remove multiple dihedrals over one bond.
827 fprintf(stderr,"Before cleaning: %d dihedrals\n",ndih);
828 clean_dih(dih,&ndih,idih,nidih,atoms,bAlldih,bRemoveDih);
831 /* Now we have unique lists of angles and dihedrals
832 * Copy them into the destination struct
834 cppar(ang, nang, plist,F_ANGLES);
835 cppar(dih, ndih, plist,F_PDIHS);
836 cppar(idih,nidih,plist,F_IDIHS);
837 cppar(pai, npai, plist,F_LJ14);
839 /* Remove all exclusions which are within nrexcl */
840 clean_excls(nnb,nrexcl,excls);