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