src/kernel/topexcl.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37 #include "assert.h"
  38 #include "sysstuff.h"
  39 #include "smalloc.h"
  40 #include "macros.h"
  41 #include "topexcl.h"
  42 #include "toputil.h"
  43
  44 typedef struct {
  45   int ai,aj;
  46 } sortable;
  47
  48 int bond_sort (const void *a, const void *b)
  49 {
  50   sortable *sa,*sb;
  51
  52   sa = (sortable *) a;
  53   sb = (sortable *) b;
  54
  55   if (sa->ai == sb->ai)
  56     return (sa->aj-sb->aj);
  57   else
  58     return (sa->ai-sb->ai);
  59 }
  60
  61 #ifdef DEBUG
  62 #define prints(str, n, s) __prints(str, n, s)
  63 static void __prints(char *str, int n, sortable *s)
  64 {
  65   int i;
  66
  67   if (debug) {
  68     fprintf(debug,"%s\n",str);
  69     fprintf(debug,"Sortables \n");
  70     for (i=0; (i < n); i++)
  71       fprintf(debug,"%d\t%d\n",s[i].ai,s[i].aj);
  72
  73     fflush(debug);
  74   }
  75 }
  76 #else
  77 #define prints(str, n, s)
  78 #endif
  79
  80 void init_nnb(t_nextnb *nnb, int nr, int nrex)
  81 {
  82   int i;
  83
  84   /* initiate nnb */
  85   nnb->nr   = nr;
  86   nnb->nrex = nrex;
  87
  88   snew(nnb->a,nr);
  89   snew(nnb->nrexcl,nr);
  90   for (i=0; (i < nr); i++) {
  91     snew(nnb->a[i],nrex+1);
  92     snew(nnb->nrexcl[i],nrex+1);
  93   }
  94 }
  95
  96 static void add_nnb (t_nextnb *nnb, int nre, int i, int j)
  97 {
  98   srenew(nnb->a[i][nre],nnb->nrexcl[i][nre]+1);
  99   nnb->a[i][nre][nnb->nrexcl[i][nre]] = j;
 100   nnb->nrexcl[i][nre]++;
 101 }
 102
 103 void done_nnb (t_nextnb *nnb)
 104 {
 105   int i,nre;
 106
 107   for (i=0; (i < nnb->nr); i++) {
 108     for (nre=0; (nre <= nnb->nrex); nre++)
 109       if (nnb->nrexcl[i][nre] > 0)
 110         sfree (nnb->a[i][nre]);
 111     sfree (nnb->nrexcl[i]);
 112     sfree (nnb->a[i]);
 113   }
 114   sfree (nnb->a);
 115   sfree (nnb->nrexcl);
 116   nnb->nr   = 0;
 117   nnb->nrex = 0;
 118 }
 119
 120 #ifdef DEBUG
 121 void __print_nnb(t_nextnb *nnb, char *s)
 122 {
 123   int i,j,k;
 124
 125   if(debug) {
 126     fprintf(debug,"%s\n",s);
 127     fprintf(debug,"nnb->nr: %d\n",nnb->nr);
 128     fprintf(debug,"nnb->nrex: %d\n",nnb->nrex);
 129     for(i=0; (i<nnb->nr); i++)
 130       for(j=0; (j<=nnb->nrex); j++) {
 131         fprintf(debug,"nrexcl[%d][%d]: %d, excl: ",i,j,nnb->nrexcl[i][j]);
 132         for(k=0; (k<nnb->nrexcl[i][j]); k++)
 133           fprintf(debug,"%d, ",nnb->a[i][j][k]);
 134         fprintf(debug,"\n");
 135       }
 136   }
 137 }
 138 #endif
 139
 140 void nnb2excl (t_nextnb *nnb, t_block *excl)
 141 {
 142   int i,j,j_index;
 143   int nre,nrx,nrs,nr_of_sortables;
 144   sortable *s;
 145
 146   srenew(excl->index, nnb->nr+1);
 147   excl->index[0]=0;
 148   for (i=0; (i < nnb->nr); i++) {
 149     /* calculate the total number of exclusions for atom i */
 150     nr_of_sortables=0;
 151     for (nre=0; (nre<=nnb->nrex); nre++)
 152       nr_of_sortables+=nnb->nrexcl[i][nre];
 153
 154     if (debug)
 155       fprintf(debug,"nr_of_sortables: %d\n",nr_of_sortables);
 156     /* make space for sortable array */
 157     snew(s,nr_of_sortables);
 158
 159     /* fill the sortable array and sort it */
 160     nrs=0;
 161     for (nre=0; (nre <= nnb->nrex); nre++)
 162       for (nrx=0; (nrx < nnb->nrexcl[i][nre]); nrx++) {
 163         s[nrs].ai = i;
 164         s[nrs].aj = nnb->a[i][nre][nrx];
 165         nrs++;
 166       }
 167     assert(nrs==nr_of_sortables);
 168     prints("nnb2excl before qsort",nr_of_sortables,s);
 169     if (nr_of_sortables > 1) {
 170       qsort ((void *)s,nr_of_sortables,(size_t)sizeof(s[0]),bond_sort);
 171       prints("nnb2excl after qsort",nr_of_sortables,s);
 172     }
 173
 174     /* remove duplicate entries from the list */
 175     j_index=0;
 176     if (nr_of_sortables > 0) {
 177       for (j=1; (j < nr_of_sortables); j++)
 178         if ((s[j].ai!=s[j-1].ai) || (s[j].aj!=s[j-1].aj))
 179           s[j_index++]=s[j-1];
 180       s[j_index++]=s[j-1];
 181     }
 182     nr_of_sortables=j_index;
 183     prints("after rm-double",j_index,s);
 184
 185     /* make space for arrays */
 186     srenew(excl->a,excl->nra+nr_of_sortables);
 187
 188     /* put the sorted exclusions in the target list */
 189     for (nrs=0; (nrs<nr_of_sortables); nrs++)
 190       excl->a[excl->nra+nrs] = s[nrs].aj;
 191     excl->nra+=nr_of_sortables;
 192     excl->index[i+1]=excl->nra;
 193
 194     /* cleanup temporary space */
 195     sfree (s);
 196   }
 197   if (debug)
 198     print_block(debug,"Exclusions","Atom","Excluded",excl);
 199 }
 200
 201 static void do_gen(int nrbonds,        /* total number of bonds in s    */
 202                    sortable *s,        /* bidirectional list of bonds   */
 203                    t_nextnb *nnb)      /* the tmp storage for excl     */
 204 /* Assume excl is initalised and s[] contains all bonds bidirectional */
 205 {
 206   int i,j,k,n,nb;
 207
 208   /* exclude self */
 209   for(i=0; (i<nnb->nr); i++)
 210     add_nnb(nnb,0,i,i);
 211   print_nnb(nnb,"After exclude self");
 212
 213   /* exclude all the bonded atoms */
 214   if (nnb->nrex > 0)
 215     for (i=0; (i < nrbonds); i++)
 216       add_nnb(nnb,1,s[i].ai,s[i].aj);
 217   print_nnb(nnb,"After exclude bonds");
 218
 219   /* for the nr of exclusions per atom */
 220   for (n=1; (n < nnb->nrex); n++)
 221
 222     /* now for all atoms */
 223     for (i=0; (i < nnb->nr); i++)
 224
 225       /* for all directly bonded atoms of atom i */
 226       for (j=0; (j < nnb->nrexcl[i][1]); j++) {
 227
 228         /* store the 1st neighbour in nb */
 229         nb = nnb->a[i][1][j];
 230
 231         /* store all atoms in nb's n-th list into i's n+1-th list */
 232         for (k=0; (k < nnb->nrexcl[nb][n]); k++)
 233           add_nnb(nnb,n+1,i,nnb->a[nb][n][k]);
 234       }
 235   print_nnb(nnb,"After exclude rest");
 236 }
 237
 238 static void add_b(t_params *bonds, int *nrf, sortable *s)
 239 {
 240   int i;
 241   int ai,aj;
 242
 243   for (i=0; (i < bonds->nr); i++) {
 244     ai = bonds->param[i].AI;
 245     aj = bonds->param[i].AJ;
 246     if ((ai < 0) || (aj < 0))
 247       fatal_error(0,"Impossible atom numbers in bond %d: ai=%d, aj=%d",
 248                   i,ai,aj);
 249     /* Add every bond twice */
 250     s[(*nrf)].ai   = ai;
 251     s[(*nrf)++].aj = aj;
 252     s[(*nrf)].aj   = ai;
 253     s[(*nrf)++].ai = aj;
 254   }
 255 }
 256
 257 void gen_nnb(t_nextnb *nnb,t_params plist[])
 258 {
 259   sortable *s;
 260   int      i,nrbonds,nrf;
 261
 262   nrbonds=0;
 263   for(i=0; (i<F_NRE); i++)
 264     if (IS_CHEMBOND(i))
 265       /* we need every bond twice (bidirectional) */
 266       nrbonds += 2*plist[i].nr;
 267
 268   snew(s,nrbonds);
 269
 270   nrf=0;
 271   for(i=0; (i<F_NRE); i++)
 272     if (IS_CHEMBOND(i))
 273       add_b(&plist[i],&nrf,s);
 274
 275   /* now sort the bonds */
 276   prints("gen_excl before qsort",nrbonds,s);
 277   if (nrbonds > 1) {
 278     qsort((void *) s,nrbonds,(size_t)sizeof(sortable),bond_sort);
 279     prints("gen_excl after qsort",nrbonds,s);
 280   }
 281
 282   do_gen(nrbonds,s,nnb);
 283   sfree(s);
 284 }
 285
 286 void generate_excl (int nrexcl,int nratoms,t_params plist[],t_block *excl)
 287 {
 288   t_nextnb nnb;
 289
 290   if (nrexcl < 0)
 291     fatal_error(0,"Can't have %d exclusions...",nrexcl);
 292   init_nnb(&nnb,nratoms,nrexcl);
 293   gen_nnb(&nnb,plist);
 294   excl->nr=nratoms;
 295   nnb2excl (&nnb,excl);
 296   done_nnb (&nnb);
 297 }
 298