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 #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)
76 else if ((ac=(p1->AI-p2->AI))!=0)
79 return (p1->AK-p2->AK);
82 static int pcomp(const void *a1, const void *a2)
89 if ((pc=(p1->AI-p2->AI))!=0)
92 return (p1->AJ-p2->AJ);
95 static int dcomp(const void *d1, const void *d2)
102 /* First sort by J & K (the two central) atoms */
103 if ((dc=(p1->AJ-p2->AJ))!=0)
107 else if ((dc=(p1->AK-p2->AK))!=0)
111 /* Then make sure to put rtp dihedrals before generated ones */
112 else if (was_dihedral_set_in_rtp(p1) &&
113 !was_dihedral_set_in_rtp(p2))
117 else if (!was_dihedral_set_in_rtp(p1) &&
118 was_dihedral_set_in_rtp(p2))
122 /* Finally, sort by I and J (two outer) atoms */
123 else if ((dc=(p1->AI-p2->AI))!=0)
129 return (p1->AL-p2->AL);
134 static gmx_bool is_dihedral_on_same_bond(t_param *p1, t_param *p2)
136 if (((p1->AJ==p2->AJ) && (p1->AK==p2->AK)) ||
137 ((p1->AJ==p2->AK) && (p1->AK==p2->AJ)))
144 static gmx_bool preq(t_param *p1, t_param *p2)
146 if ((p1->AI==p2->AI) && (p1->AJ==p2->AJ))
152 static void rm2par(t_param p[], int *np, peq eq)
163 for(i=1; (i<(*np)); i++)
164 if (!eq(&p[i],&p[i-1]))
166 /* Index now holds pointers to all the non-equal params,
167 * this only works when p is sorted of course
169 for(i=0; (i<nind); i++) {
170 for(j=0; (j<MAXATOMLIST); j++)
171 p[i].a[j]=p[index[i]].a[j];
172 for(j=0; (j<MAXFORCEPARAM); j++)
173 p[i].c[j]=p[index[i]].c[j];
174 if (p[index[i]].a[0] == p[index[i]].a[1]) {
177 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
178 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
179 p[i].a[0],p[i].a[1],p[i].a[2],p[i].a[3]);
181 } else if (index[i] > i) {
182 /* Copy the string only if it comes from somewhere else
183 * otherwise we will end up copying a random (newly freed) pointer.
184 * Since the index is sorted we only have to test for index[i] > i.
186 strcpy(p[i].s,p[index[i]].s);
194 static void cppar(t_param p[], int np, t_params plist[], int ftype)
205 for(i=0; (i<np); i++) {
206 for(j=0; (j<nral); j++)
207 ps->param[ps->nr].a[j] = p[i].a[j];
208 for(j=0; (j<nrfp); j++)
209 ps->param[ps->nr].c[j] = p[i].c[j];
210 for(j=0; (j<MAXSLEN); j++)
211 ps->param[ps->nr].s[j] = p[i].s[j];
216 static void cpparam(t_param *dest, t_param *src)
220 for(j=0; (j<MAXATOMLIST); j++)
221 dest->a[j] = src->a[j];
222 for(j=0; (j<MAXFORCEPARAM); j++)
223 dest->c[j] = src->c[j];
224 for(j=0; (j<MAXSLEN); j++)
225 dest->s[j] = src->s[j];
228 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
234 for(j=0; (j<MAXFORCEPARAM); j++)
243 static int int_comp(const void *a,const void *b)
245 return (*(int *)a) - (*(int *)b);
248 static int atom_id_comp(const void *a,const void *b)
250 return (*(atom_id *)a) - (*(atom_id *)b);
253 static int eq_imp(atom_id a1[],atom_id a2[])
258 for(j=0; (j<4); j++) {
262 qsort(b1,4,(size_t)sizeof(b1[0]),int_comp);
263 qsort(b2,4,(size_t)sizeof(b2[0]),int_comp);
272 static int idcomp(const void *a,const void *b)
279 if ((d=(pa->a[0]-pb->a[0])) != 0)
281 else if ((d=(pa->a[3]-pb->a[3])) != 0)
283 else if ((d=(pa->a[1]-pb->a[1])) != 0)
286 return (int) (pa->a[2]-pb->a[2]);
289 static void sort_id(int nr,t_param ps[])
293 /* First swap order of atoms around if necessary */
294 for(i=0; (i<nr); i++) {
295 if (ps[i].a[3] < ps[i].a[0]) {
296 tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
297 tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
302 qsort(ps,nr,(size_t)sizeof(ps[0]),idcomp);
305 static int n_hydro(atom_id a[],char ***atomname)
310 for(i=0; (i<4); i+=3) {
311 aname=*atomname[a[i]];
312 c0=toupper(aname[0]);
315 else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9')) {
316 c1=toupper(aname[1]);
324 /* Clean up the dihedrals (both generated and read from the .rtp
326 static void clean_dih(t_param *dih, int *ndih,t_param improper[],int nimproper,
327 t_atoms *atoms,gmx_bool bKeepAllGeneratedDihedrals,
328 gmx_bool bRemoveDihedralIfWithImproper)
333 /* Construct the list of the indices of the dihedrals
334 * (i.e. generated or read) that might be kept. */
335 snew(index, *ndih+1);
336 if (bKeepAllGeneratedDihedrals)
338 fprintf(stderr,"Keeping all generated dihedrals\n");
340 for(i = 0; i < nind; i++)
349 /* Check if generated dihedral i should be removed. The
350 * dihedrals have been sorted by dcomp() above, so all those
351 * on the same two central atoms are together, with those from
352 * the .rtp file preceding those that were automatically
353 * generated. We remove the latter if the former exist. */
354 for(i = 0; i < *ndih; i++)
356 /* Keep the dihedrals that were defined in the .rtp file,
357 * and the dihedrals that were generated and different
358 * from the last one (whether it was generated or not). */
359 if (was_dihedral_set_in_rtp(&dih[i]) ||
361 !is_dihedral_on_same_bond(&dih[i],&dih[i-1]))
370 for(i=0; i<nind; i++)
372 gmx_bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
373 gmx_bool bKeep = TRUE;
374 if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
376 /* Remove the dihedral if there is an improper on the same
378 for(j = 0; j < nimproper && bKeep; j++)
380 bKeep = !is_dihedral_on_same_bond(&dih[index[i]],&improper[j]);
386 /* If we don't want all dihedrals, we want to select the
387 * ones with the fewest hydrogens. Note that any generated
388 * dihedrals on the same bond as an .rtp dihedral may have
389 * been already pruned above in the construction of
390 * index[]. However, their parameters are still present,
391 * and l is looping over this dihedral and all of its
392 * pruned siblings. */
393 int bestl = index[i];
394 if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
396 /* Minimum number of hydrogens for i and l atoms */
400 is_dihedral_on_same_bond(&dih[index[i]],&dih[l]));
403 int nh = n_hydro(dih[l].a,atoms->atomname);
417 cpparam(&dih[k],&dih[bestl]);
423 for (i = k; i < *ndih; i++)
432 static int get_impropers(t_atoms *atoms,t_hackblock hb[],t_param **improper,
433 gmx_bool bAllowMissing)
436 t_rbondeds *impropers;
437 t_rbonded *hbimproper;
438 int nimproper,i,j,k,r,start,ninc,nalloc;
439 atom_id ai[MAXATOMLIST];
444 snew(*improper,nalloc);
446 /* Add all the impropers from the residue database to the list. */
450 for(i=0; (i<atoms->nres); i++) {
451 impropers=&hb[i].rb[ebtsIDIHS];
452 for(j=0; (j<impropers->nb); j++) {
454 for(k=0; (k<4) && !bStop; k++) {
455 ai[k] = search_atom(impropers->b[j].a[k],start,
457 "improper",bAllowMissing);
458 if (ai[k] == NO_ATID)
462 if (nimproper == nalloc) {
464 srenew(*improper,nalloc);
467 set_p(&((*improper)[nimproper]),ai,NULL,impropers->b[j].s);
471 while ((start<atoms->nr) && (atoms->atom[start].resind == i))
479 static int nb_dist(t_nextnb *nnb,int ai,int aj)
489 nrexcl=nnb->nrexcl[ai];
490 for(nre=1; (nre < nnb->nrex); nre++) {
492 for(nrx=0; (nrx < nrexcl[nre]); nrx++) {
493 if ((aj == a[nrx]) && (NRE == -1))
500 gmx_bool is_hydro(t_atoms *atoms,int ai)
502 return ((*(atoms->atomname[ai]))[0] == 'H');
505 static void get_atomnames_min(int n,char **anm,
506 int resind,t_atoms *atoms,atom_id *a)
510 /* Assume ascending residue numbering */
512 if (atoms->atom[a[m]].resind < resind)
514 else if (atoms->atom[a[m]].resind > resind)
518 strcat(anm[m],*(atoms->atomname[a[m]]));
522 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
523 gmx_bool bAllowMissing)
526 atom_id a,astart,i1,i2,itmp;
532 for(a=0; a<atoms->nr; a++) {
533 r = atoms->atom[a].resind;
534 if (a==atoms->nr-1 || atoms->atom[a+1].resind != r) {
535 hbexcl = &hb[r].rb[ebtsEXCLS];
537 for(e=0; e<hbexcl->nb; e++) {
538 anm = hbexcl->b[e].a[0];
539 i1 = search_atom(anm,astart,atoms,
540 "exclusion",bAllowMissing);
541 anm = hbexcl->b[e].a[1];
542 i2 = search_atom(anm,astart,atoms,
543 "exclusion",bAllowMissing);
544 if (i1!=NO_ATID && i2!=NO_ATID) {
550 srenew(excls[i1].e,excls[i1].nr+1);
551 excls[i1].e[excls[i1].nr] = i2;
560 for(a=0; a<atoms->nr; a++)
562 qsort(excls[a].e,excls[a].nr,(size_t)sizeof(atom_id),atom_id_comp);
565 static void remove_excl(t_excls *excls, int remove)
569 for(i=remove+1; i<excls->nr; i++)
570 excls->e[i-1] = excls->e[i];
575 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
577 int i,j,j1,k,k1,l,l1,m,n,e;
581 /* extract all i-j-k-l neighbours from nnb struct */
582 for(i=0; (i<nnb->nr); i++) {
583 /* For all particles */
586 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
587 /* For all first neighbours */
590 for(e=0; e<excl->nr; e++)
591 if (excl->e[e] == j1)
595 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
596 /* For all first neighbours of j1 */
599 for(e=0; e<excl->nr; e++)
600 if (excl->e[e] == k1)
604 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
605 /* For all first neighbours of k1 */
608 for(e=0; e<excl->nr; e++)
609 if (excl->e[e] == l1)
617 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
619 int i,j,j1,k,k1,l,l1,m,n,e,N;
622 for(N=1; (N<min(nrexcl,nnb->nrex)); N++) {
623 /* extract all i-j-k-l neighbours from nnb struct */
624 for(i=0; (i<nnb->nr); i++) {
625 /* For all particles */
628 excl->nr += nnb->nrexcl[i][N];
629 srenew(excl->e,excl->nr);
630 for(j=0; (j<nnb->nrexcl[i][N]); j++)
631 /* For all first neighbours */
632 if (nnb->a[i][N][j] != i)
633 excl->e[n++] = nnb->a[i][N][j];
638 /* Generate pairs, angles and dihedrals from .rtp settings */
639 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
640 t_params plist[], t_excls excls[], t_hackblock hb[],
641 gmx_bool bAllowMissing)
643 t_param *ang,*dih,*pai,*improper;
644 t_rbondeds *hbang, *hbdih;
646 int res,minres,maxres;
647 int i,j,j1,k,k1,l,l1,m,n,i1,i2;
648 int ninc,maxang,maxdih,maxpai;
649 int nang,ndih,npai,nimproper,nbd;
651 gmx_bool bFound,bExcl;
654 /* These are the angles, dihedrals and pairs that we generate
655 * from the bonds. The ones that are already there from the rtp file
662 maxang = maxdih = maxpai = ninc;
672 gen_excls(atoms,excls,hb,bAllowMissing);
674 /* Extract all i-j-k-l neighbours from nnb struct to generate all
675 * angles and dihedrals. */
676 for(i=0; (i<nnb->nr); i++)
677 /* For all particles */
678 for(j=0; (j<nnb->nrexcl[i][1]); j++) {
679 /* For all first neighbours */
681 for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
682 /* For all first neighbours of j1 */
685 /* Generate every angle only once */
687 if (nang == maxang) {
696 set_p_string(&(ang[nang]),"");
698 minres = atoms->atom[ang[nang].a[0]].resind;
701 minres = min(minres,atoms->atom[ang[nang].a[m]].resind);
702 maxres = max(maxres,atoms->atom[ang[nang].a[m]].resind);
704 res = 2*minres-maxres;
706 res += maxres-minres;
707 get_atomnames_min(3,anm,res,atoms,ang[nang].a);
708 hbang=&hb[res].rb[ebtsANGLES];
709 for(l=0; (l<hbang->nb); l++) {
710 if (strcmp(anm[1],hbang->b[l].AJ)==0) {
714 ((strcmp(anm[m],hbang->b[l].AI)==0) &&
715 (strcmp(anm[2-m],hbang->b[l].AK)==0)));
717 set_p_string(&(ang[nang]),hbang->b[l].s);
721 } while (res < maxres);
725 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
728 for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
729 /* For all first neighbours of k1 */
731 if ((l1 != i) && (l1 != j1)) {
732 if (ndih == maxdih) {
740 for (m=0; m<MAXFORCEPARAM; m++)
741 dih[ndih].c[m]=NOTSET;
742 set_p_string(&(dih[ndih]),"");
745 minres = atoms->atom[dih[ndih].a[0]].resind;
748 minres = min(minres,atoms->atom[dih[ndih].a[m]].resind);
749 maxres = max(maxres,atoms->atom[dih[ndih].a[m]].resind);
751 res = 2*minres-maxres;
753 res += maxres-minres;
754 get_atomnames_min(4,anm,res,atoms,dih[ndih].a);
755 hbdih=&hb[res].rb[ebtsPDIHS];
756 for(n=0; (n<hbdih->nb); n++) {
760 ((strcmp(anm[3*m], hbdih->b[n].AI)==0) &&
761 (strcmp(anm[1+m], hbdih->b[n].AJ)==0) &&
762 (strcmp(anm[2-m], hbdih->b[n].AK)==0) &&
763 (strcmp(anm[3-3*m],hbdih->b[n].AL)==0)));
765 set_p_string(&dih[ndih],hbdih->b[n].s);
767 /* Set the last parameter to be able to see
768 if the dihedral was in the rtp list.
770 dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
773 /* Set the next direct in case the rtp contains
774 multiple entries for this dihedral.
776 if (ndih == maxdih) {
784 for (m=0; m<MAXFORCEPARAM; m++)
785 dih[ndih].c[m]=NOTSET;
788 } while (res < maxres);
791 if (ndih == maxdih) {
799 for (m=0; m<MAXFORCEPARAM; m++)
800 dih[ndih].c[m]=NOTSET;
801 set_p_string(&(dih[ndih]),"");
805 nbd=nb_dist(nnb,i,l1);
807 fprintf(debug,"Distance (%d-%d) = %d\n",i+1,l1+1,nbd);
812 for(m=0; m<excls[i1].nr; m++)
813 bExcl = bExcl || excls[i1].e[m]==i2;
815 if (rtp[0].bGenerateHH14Interactions ||
816 !(is_hydro(atoms,i1) && is_hydro(atoms,i2))) {
817 if (npai == maxpai) {
825 set_p_string(&(pai[npai]),"");
837 /* Sort angles with respect to j-i-k (middle atom first) */
839 qsort(ang,nang,(size_t)sizeof(ang[0]),acomp);
841 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
843 qsort(dih,ndih,(size_t)sizeof(dih[0]),dcomp);
847 qsort(pai,npai,(size_t)sizeof(pai[0]),pcomp);
849 /* Remove doubles, could occur in 6-rings, such as phenyls,
850 maybe one does not want this when fudgeQQ < 1.
852 fprintf(stderr,"Before cleaning: %d pairs\n",npai);
853 rm2par(pai,&npai,preq);
856 /* Get the impropers from the database */
857 nimproper = get_impropers(atoms,hb,&improper,bAllowMissing);
859 /* Sort the impropers */
860 sort_id(nimproper,improper);
863 fprintf(stderr,"Before cleaning: %d dihedrals\n",ndih);
864 clean_dih(dih,&ndih,improper,nimproper,atoms,
865 rtp[0].bKeepAllGeneratedDihedrals,
866 rtp[0].bRemoveDihedralIfWithImproper);
869 /* Now we have unique lists of angles and dihedrals
870 * Copy them into the destination struct
872 cppar(ang, nang, plist,F_ANGLES);
873 cppar(dih, ndih, plist,F_PDIHS);
874 cppar(improper,nimproper,plist,F_IDIHS);
875 cppar(pai, npai, plist,F_LJ14);
877 /* Remove all exclusions which are within nrexcl */
878 clean_excls(nnb,rtp[0].nrexcl,excls);