src/local/init_sh.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_init_sh_c = "$Id$";
  30
  31 #include "init_sh.h"
  32 #include "smalloc.h"
  33 #include "assert.h"
  34 #include "names.h"
  35
  36 static void pr_shell(FILE *log,int ns,t_shell s[])
  37 {
  38   int i;
  39
  40   fprintf(log,"SHELL DATA\n");
  41   fprintf(log,"%5s  %8s  %5s  %5s  %5s\n",
  42           "Shell","Force k","Nucl1","Nucl2","Nucl3");
  43   for(i=0; (i<ns); i++) {
  44     fprintf(log,"%5d  %8.3f  %5d",s[i].shell,1.0/s[i].k_1,s[i].nucl1);
  45     if (s[i].nnucl == 2)
  46       fprintf(log,"  %5d\n",s[i].nucl2);
  47     else if (s[i].nnucl == 3)
  48       fprintf(log,"  %5d  %5d\n",s[i].nucl2,s[i].nucl3);
  49     else
  50       fprintf(log,"\n");
  51   }
  52 }
  53
  54 t_shell *init_shells(FILE *log,int start,int homenr,
  55                      t_idef *idef,t_mdatoms *md,int *nshell)
  56 {
  57   t_shell     *shell=NULL;
  58   int         *shell_index;
  59   int         n[eptNR],ns,nsi;
  60   int         i,j,type,ftype,nra;
  61   int         pt1,pt2,a1,a2;
  62   bool        bS1,bS2;
  63   t_iatom     *ia;
  64
  65   for(i=0; (i<eptNR); i++)
  66     n[i]=0;
  67   snew(shell_index,homenr);
  68   nsi = 0;
  69   for(i=start; (i<start+homenr); i++) {
  70     n[md->ptype[i]]++;
  71     if (md->ptype[i] == eptShell)
  72       shell_index[i-start] = nsi++;
  73   }
  74   assert(nsi == n[eptShell]);
  75   for(i=0; (i<eptNR); i++)
  76     if (n[i]!=0)
  77       fprintf(log,"There are: %d %s\n",n[i],ptype_str[i]);
  78
  79   ns      = n[eptShell];
  80   *nshell = ns;
  81   if (ns > 0) {
  82     snew(shell,ns);
  83
  84     /* Initiate the shell structures */
  85     for(i=0; (i<ns); i++) {
  86       shell[i].shell=NO_ATID;
  87       shell[i].nucl1=NO_ATID;
  88       shell[i].nucl2=NO_ATID;
  89       shell[i].nucl3=NO_ATID;
  90       shell[i].nnucl=0;
  91       shell[i].k_1=0;
  92       shell[i].k=0;
  93     }
  94
  95     /* Now fill the structures */
  96     ns=0;
  97     ia=idef->il[F_BONDS].iatoms;
  98     for(i=0; (i<idef->il[F_BONDS].nr); ) {
  99       type  = ia[0];
 100       ftype = idef->functype[type];
 101       nra   = interaction_function[ftype].nratoms;
 102
 103       /* Check whether we have a bond */
 104
 105       if (md->ptype[ia[1]] == eptShell) {
 106         a1 = ia[1];
 107         a2 = ia[2];
 108       }
 109       else if (md->ptype[ia[2]] == eptShell) {
 110         a1 = ia[2];
 111         a2 = ia[1];
 112       }
 113       else {
 114         i  += nra+1;
 115         ia += nra+1;
 116         continue;
 117       }
 118       /* Check whether one of the particles is a shell... */
 119       nsi = shell_index[a1-start];
 120       if ((nsi < 0) || (nsi >= *nshell))
 121         fatal_error(0,"nsi is %d should be within 0 - %d. a1 = %d",
 122                     nsi,*nshell,a1);
 123       if (shell[nsi].shell == NO_ATID) {
 124         shell[nsi].shell = a1;
 125         ns ++;
 126       }
 127       else if (shell[nsi].shell != a1)
 128         fatal_error(0,"What is this?");
 129
 130       if      (shell[nsi].nucl1 == NO_ATID)
 131         shell[nsi].nucl1 = a2;
 132       else if (shell[nsi].nucl2 == NO_ATID)
 133         shell[nsi].nucl2 = a2;
 134       else if (shell[nsi].nucl3 == NO_ATID)
 135         shell[nsi].nucl3 = a2;
 136       else {
 137         pr_shell(log,ns,shell);
 138         fatal_error(0,"Can not handle more than three bonds per shell\n");
 139       }
 140       shell[nsi].k    += idef->iparams[type].harmonic.krA;
 141       shell[nsi].nnucl++;
 142
 143       ia += nra+1;
 144       i  += nra+1;
 145     }
 146     ia = idef->il[F_WPOL].iatoms;
 147     for(i=0; (i<idef->il[F_WPOL].nr); ) {
 148       type  = ia[0];
 149       ftype = idef->functype[type];
 150       nra   = interaction_function[ftype].nratoms;
 151
 152       a1    = ia[1]+4;  /* Shell */
 153       a2    = ia[1]+3;  /* Dummy */
 154
 155       /* Check whether one of the particles is a shell... */
 156       nsi = shell_index[a1-start];
 157       if ((nsi < 0) || (nsi >= *nshell))
 158         fatal_error(0,"nsi is %d should be within 0 - %d. a1 = %d",
 159                     nsi,*nshell,a1);
 160       if (shell[nsi].shell == NO_ATID) {
 161         shell[nsi].shell = a1;
 162         ns ++;
 163       }
 164       else if (shell[nsi].shell != a1)
 165         fatal_error(0,"What is this? shell=%d, a1=%d",shell[nsi].shell,a1);
 166
 167       shell[nsi].nucl1 = a2;
 168       shell[nsi].k     = (idef->iparams[type].wpol.kx+
 169                           idef->iparams[type].wpol.ky+
 170                           idef->iparams[type].wpol.kz)/3.0;
 171       shell[nsi].nnucl++;
 172
 173       ia += nra+1;
 174       i  += nra+1;
 175     }
 176     /* Verify whether it's all correct */
 177     assert(ns == *nshell);
 178
 179     for(i=0; (i<ns); i++)
 180       shell[i].k_1 = 1.0/shell[i].k;
 181
 182     if (debug)
 183       pr_shell(debug,ns,shell);
 184   }
 185   return shell;
 186 }
 187