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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_topexcl_c = "$Id$";
  30
  31 #ifdef HAVE_IDENT
  32 #ident "@(#) topexcl.c 1.38 2/2/97"
  33 #endif
  34
  35 #include "assert.h"
  36 #include "sysstuff.h"
  37 #include "smalloc.h"
  38 #include "macros.h"
  39 #include "topexcl.h"
  40 #include "toputil.h"
  41
  42 typedef struct {
  43   int ai,aj;
  44 } sortable;
  45
  46 int bond_sort (const void *a, const void *b)
  47 {
  48   sortable *sa,*sb;
  49
  50   sa = (sortable *) a;
  51   sb = (sortable *) b;
  52
  53   if (sa->ai == sb->ai)
  54     return (sa->aj-sb->aj);
  55   else
  56     return (sa->ai-sb->ai);
  57 }
  58
  59 #ifdef DEBUG
  60 #define prints(str, n, s) __prints(str, n, s)
  61 static void __prints(char *str, int n, sortable *s)
  62 {
  63   int i;
  64
  65   if (debug) {
  66     fprintf(debug,"%s\n",str);
  67     fprintf(debug,"Sortables \n");
  68     for (i=0; (i < n); i++)
  69       fprintf(debug,"%d\t%d\n",s[i].ai,s[i].aj);
  70
  71     fflush(debug);
  72   }
  73 }
  74 #else
  75 #define prints(str, n, s)
  76 #endif
  77
  78 void init_nnb(t_nextnb *nnb, int nr, int nrex)
  79 {
  80   int i;
  81
  82   /* initiate nnb */
  83   nnb->nr   = nr;
  84   nnb->nrex = nrex;
  85
  86   snew(nnb->a,nr);
  87   snew(nnb->nrexcl,nr);
  88   for (i=0; (i < nr); i++) {
  89     snew(nnb->a[i],nrex+1);
  90     snew(nnb->nrexcl[i],nrex+1);
  91   }
  92 }
  93
  94 static void add_nnb (t_nextnb *nnb, int nre, int i, int j)
  95 {
  96   srenew(nnb->a[i][nre],nnb->nrexcl[i][nre]+1);
  97   nnb->a[i][nre][nnb->nrexcl[i][nre]] = j;
  98   nnb->nrexcl[i][nre]++;
  99 }
 100
 101 void done_nnb (t_nextnb *nnb)
 102 {
 103   int i,nre;
 104
 105   for (i=0; (i < nnb->nr); i++) {
 106     for (nre=0; (nre <= nnb->nrex); nre++)
 107       if (nnb->nrexcl[i][nre] > 0)
 108         sfree (nnb->a[i][nre]);
 109     sfree (nnb->nrexcl[i]);
 110     sfree (nnb->a[i]);
 111   }
 112   sfree (nnb->a);
 113   sfree (nnb->nrexcl);
 114   nnb->nr   = 0;
 115   nnb->nrex = 0;
 116 }
 117
 118 #ifdef DEBUG
 119 void __print_nnb(t_nextnb *nnb, char *s)
 120 {
 121   int i,j,k;
 122
 123   if(debug) {
 124     fprintf(debug,"%s\n",s);
 125     fprintf(debug,"nnb->nr: %d\n",nnb->nr);
 126     fprintf(debug,"nnb->nrex: %d\n",nnb->nrex);
 127     for(i=0; (i<nnb->nr); i++)
 128       for(j=0; (j<=nnb->nrex); j++) {
 129         fprintf(debug,"nrexcl[%d][%d]: %d, excl: ",i,j,nnb->nrexcl[i][j]);
 130         for(k=0; (k<nnb->nrexcl[i][j]); k++)
 131           fprintf(debug,"%d, ",nnb->a[i][j][k]);
 132         fprintf(debug,"\n");
 133       }
 134   }
 135 }
 136 #endif
 137
 138 void nnb2excl (t_nextnb *nnb, t_block *excl)
 139 {
 140   int i,j,j_index;
 141   int nre,nrx,nrs,nr_of_sortables;
 142   sortable *s;
 143
 144   srenew(excl->index, nnb->nr+1);
 145   excl->index[0]=0;
 146   for (i=0; (i < nnb->nr); i++) {
 147     /* calculate the total number of exclusions for atom i */
 148     nr_of_sortables=0;
 149     for (nre=0; (nre<=nnb->nrex); nre++)
 150       nr_of_sortables+=nnb->nrexcl[i][nre];
 151
 152     if (debug)
 153       fprintf(debug,"nr_of_sortables: %d\n",nr_of_sortables);
 154     /* make space for sortable array */
 155     snew(s,nr_of_sortables);
 156
 157     /* fill the sortable array and sort it */
 158     nrs=0;
 159     for (nre=0; (nre <= nnb->nrex); nre++)
 160       for (nrx=0; (nrx < nnb->nrexcl[i][nre]); nrx++) {
 161         s[nrs].ai = i;
 162         s[nrs].aj = nnb->a[i][nre][nrx];
 163         nrs++;
 164       }
 165     assert(nrs==nr_of_sortables);
 166     prints("nnb2excl before qsort",nr_of_sortables,s);
 167     qsort ((void *)s,nr_of_sortables,(size_t)sizeof(s[0]),bond_sort);
 168     prints("nnb2excl after qsort",nr_of_sortables,s);
 169
 170     /* remove duplicate entries from the list */
 171     j_index=0;
 172     if (nr_of_sortables > 0) {
 173       for (j=1; (j < nr_of_sortables); j++)
 174         if ((s[j].ai!=s[j-1].ai) || (s[j].aj!=s[j-1].aj))
 175           s[j_index++]=s[j-1];
 176       s[j_index++]=s[j-1];
 177     }
 178     nr_of_sortables=j_index;
 179     prints("after rm-double",j_index,s);
 180
 181     /* make space for arrays */
 182     srenew(excl->a,excl->nra+nr_of_sortables);
 183
 184     /* put the sorted exclusions in the target list */
 185     for (nrs=0; (nrs<nr_of_sortables); nrs++)
 186       excl->a[excl->nra+nrs] = s[nrs].aj;
 187     excl->nra+=nr_of_sortables;
 188     excl->index[i+1]=excl->nra;
 189
 190     /* cleanup temporary space */
 191     sfree (s);
 192   }
 193   if (debug)
 194     print_block(debug,"Exclusions","Atom","Excluded",excl);
 195 }
 196
 197 static void do_gen(int nrbonds,        /* total number of bonds in s    */
 198                    sortable *s,        /* bidirectional list of bonds   */
 199                    t_nextnb *nnb)      /* the tmp storage for excl     */
 200 /* Assume excl is initalised and s[] contains all bonds bidirectional */
 201 {
 202   int i,j,k,n,nb;
 203
 204   /* exclude self */
 205   for(i=0; (i<nnb->nr); i++)
 206     add_nnb(nnb,0,i,i);
 207   print_nnb(nnb,"After exclude self");
 208
 209   /* exclude all the bonded atoms */
 210   if (nnb->nrex > 0)
 211     for (i=0; (i < nrbonds); i++)
 212       add_nnb(nnb,1,s[i].ai,s[i].aj);
 213   print_nnb(nnb,"After exclude bonds");
 214
 215   /* for the nr of exclusions per atom */
 216   for (n=1; (n < nnb->nrex); n++)
 217
 218     /* now for all atoms */
 219     for (i=0; (i < nnb->nr); i++)
 220
 221       /* for all directly bonded atoms of atom i */
 222       for (j=0; (j < nnb->nrexcl[i][1]); j++) {
 223
 224         /* store the 1st neighbour in nb */
 225         nb = nnb->a[i][1][j];
 226
 227         /* store all atoms in nb's n-th list into i's n+1-th list */
 228         for (k=0; (k < nnb->nrexcl[nb][n]); k++)
 229           add_nnb(nnb,n+1,i,nnb->a[nb][n][k]);
 230       }
 231   print_nnb(nnb,"After exclude rest");
 232 }
 233
 234 static void add_b(t_params *bonds, int *nrf, sortable *s)
 235 {
 236   int i;
 237   int ai,aj;
 238
 239   for (i=0; (i < bonds->nr); i++) {
 240     ai = bonds->param[i].AI;
 241     aj = bonds->param[i].AJ;
 242     if ((ai < 0) || (aj < 0))
 243       fatal_error(0,"Impossible atom numbers in bond %d: ai=%d, aj=%d",
 244                   i,ai,aj);
 245     /* Add every bond twice */
 246     s[(*nrf)].ai   = ai;
 247     s[(*nrf)++].aj = aj;
 248     s[(*nrf)].aj   = ai;
 249     s[(*nrf)++].ai = aj;
 250   }
 251 }
 252
 253 void gen_nnb(t_nextnb *nnb,t_params plist[])
 254 {
 255   sortable *s;
 256   int      i,nrbonds,nrf;
 257
 258   nrbonds=0;
 259   for(i=0; (i<F_NRE); i++)
 260     if ( (interaction_function[i].flags & IF_CONNECT) &&
 261          !(interaction_function[i].flags & IF_DUMMY) )
 262       /* we need every bond twice (bidirectional) */
 263       nrbonds += 2*plist[i].nr;
 264
 265   snew(s,nrbonds);
 266
 267   nrf=0;
 268   for(i=0; (i<F_NRE); i++)
 269     if ( (interaction_function[i].flags & IF_CONNECT) &&
 270          !(interaction_function[i].flags & IF_DUMMY) )
 271       add_b(&plist[i],&nrf,s);
 272
 273   /* now sort the bonds */
 274   prints("gen_excl before qsort",nrbonds,s);
 275   qsort((void *) s,nrbonds,(size_t)sizeof(sortable),bond_sort);
 276   prints("gen_excl after qsort",nrbonds,s);
 277
 278   do_gen(nrbonds,s,nnb);
 279   sfree(s);
 280 }
 281
 282 void generate_excl (int nrexcl,int nratoms,t_params plist[],t_block *excl)
 283 {
 284   t_nextnb nnb;
 285
 286   if (nrexcl < 0)
 287     fatal_error(0,"Can't have %d exclusions...",nrexcl);
 288   init_nnb(&nnb,nratoms,nrexcl);
 289   gen_nnb(&nnb,plist);
 290   excl->nr=nratoms;
 291   nnb2excl (&nnb,excl);
 292   done_nnb (&nnb);
 293 }
 294