src/kernel/topshake.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_topshake_c = "$Id$";
  30
  31 #include <ctype.h>
  32
  33 #include "sysstuff.h"
  34 #include "physics.h"
  35 #include "macros.h"
  36 #include "readir.h"
  37 #include "typedefs.h"
  38 #include "topshake.h"
  39 #include "toppush.h"
  40 #include "toputil.h"
  41 #include "topdirs.h"
  42 #include "smalloc.h"
  43
  44 static void copy_bond (t_params *pr, int to, int from)
  45 /* copies an entry in a bond list to another position.
  46  * does no allocing or freeing of memory
  47  */
  48 {
  49   memcpy((char*) &(pr->param[to]),(char*) &(pr->param[from]),
  50          (size_t)sizeof(pr->param[from]));
  51 }
  52
  53 static int count_hydrogens (char ***atomname, int nra, atom_id a[])
  54 {
  55   int  i,nh;
  56
  57   if (!atomname)
  58     fatal_error(0,"Cannot call count_hydrogens with no atomname (%s %d)",
  59                 __FILE__,__LINE__);
  60
  61   nh=0;
  62   for (i=0; (i<nra); i++)
  63     if (toupper(**(atomname[a[i]]))=='H')
  64       nh++;
  65
  66   return nh;
  67 }
  68
  69 void make_shake (t_params plist[],t_atoms *atoms,t_atomtype *at,int nshake)
  70 {
  71   char         ***info=atoms->atomname;
  72   t_params     *pr;
  73   t_params     *bonds;
  74   t_param      p,*bond,*ang;
  75   real         b_ij,b_jk;
  76   int          nb,b,i,j,ftype;
  77   bool         bFound;
  78
  79   if (nshake != eshNONE) {
  80     switch (nshake) {
  81     case eshHBONDS:
  82       printf("turning H bonds into constraints...\n");
  83       break;
  84     case eshALLBONDS:
  85       printf("turning all bonds into constraints...\n");
  86       break;
  87     case eshHANGLES:
  88       printf("turning all bonds and H angles into constraints...\n");
  89       break;
  90     case eshALLANGLES:
  91       printf("turning all bonds and angles into constraints...\n");
  92       break;
  93     default:
  94       fatal_error(0,"Invalid option for make_shake (%d)",nshake);
  95     }
  96
  97     /* Add all the angles with hydrogens to the shake list
  98      * and remove them from the bond list
  99      */
 100     for(ftype = 0; (ftype < F_NRE); ftype++) {
 101       if (interaction_function[ftype].flags & IF_BTYPE) {
 102         bonds=&(plist[ftype]);
 103
 104         if ((nshake == eshHANGLES) || (nshake == eshALLANGLES)) {
 105           /* horrible shortcut */
 106           pr = &(plist[F_ANGLES]);
 107           for (i=0; (i < pr->nr); ) {
 108             ang=&(pr->param[i]);
 109 #ifdef DEBUG
 110             printf("Angle: %d-%d-%d\n",ang->AI,ang->AJ,ang->AK);
 111 #endif
 112             if ((nshake == eshALLANGLES) ||
 113                 (count_hydrogens(info,3,ang->a) > 0)) {
 114               /* Can only add hydrogen angle shake, if the two bonds
 115                * are constrained.
 116                * append this angle to the shake list
 117                */
 118               p.AI = ang->AI;
 119               p.AJ = ang->AK;
 120
 121               /* Calculate length of constraint */
 122               bFound=FALSE;
 123               b_ij=b_jk=0.0;
 124               for (j=0; !bFound && (j<bonds->nr); j++) {
 125                 bond=&(bonds->param[j]);
 126                 if (((bond->AI==ang->AI) &&
 127                      (bond->AJ==ang->AJ)) ||
 128                     ((bond->AI==ang->AJ) &&
 129                      (bond->AJ==ang->AI)))
 130                   b_ij=bond->C0;
 131                 if (((bond->AI==ang->AK) &&
 132                      (bond->AJ==ang->AJ)) ||
 133                     ((bond->AI==ang->AJ) &&
 134                      (bond->AJ==ang->AK)))
 135                   b_jk=bond->C0;
 136                 bFound = (b_ij!=0.0) && (b_jk!=0.0);
 137               }
 138               if (!bFound)
 139                 fatal_error(0,"No bond information for bond %s-%s or %s-%s",
 140                             *info[ang->AI],*info[ang->AJ],
 141                             *info[ang->AJ],*info[ang->AK]);
 142               /* apply law of cosines */
 143               p.C0 = sqrt( b_ij*b_ij + b_jk*b_jk -
 144                            2.0*b_ij*b_jk*cos(DEG2RAD*ang->C0) );
 145               p.C1 = p.C0;
 146 #ifdef DEBUG
 147               printf("p: %d, q: %d, dist: %12.5e\n",p.AI,p.AJ,p.C0);
 148 #endif
 149               push_bondnow (&(plist[F_SHAKE]),&p);
 150               /* move the last bond to this position */
 151               copy_bond (pr,i,pr->nr-1);
 152               /* should free memory here!! */
 153               pr->nr--;
 154             }
 155             else
 156               i++;
 157           }
 158         } /* if nshake == ANGLES or ALLANGLES */
 159       } /* if IF_BTYPE */
 160     } /* for ftype */
 161   } /* nshake != NONE */
 162
 163   /* Add all the bonds with hydrogens to the shake list
 164    * and remove them from the bond list
 165    */
 166   if ( (nshake == eshHBONDS)    || (nshake == eshALLBONDS) ||
 167        (nshake == eshALLANGLES) || (nshake == eshHANGLES) ) {
 168     for (ftype=0; (ftype < F_NRE); ftype++) {
 169       if (interaction_function[ftype].flags & IF_BTYPE) {
 170         pr = &(plist[ftype]);
 171         for (i=0; (i < pr->nr); ) {
 172           if ( (nshake != eshHBONDS) ||
 173                (count_hydrogens (info,2,pr->param[i].a) > 0) ) {
 174             /* append this bond to the shake list */
 175             p.AI = pr->param[i].AI;
 176             p.AJ = pr->param[i].AJ;
 177             p.C0 = pr->param[i].C0;
 178             p.C1 = pr->param[i].C2;
 179             push_bondnow (&(plist[F_SHAKE]),&p);
 180
 181             /* move the last bond to this position */
 182             copy_bond (pr,i,pr->nr-1);
 183
 184             /* should free memory here!! */
 185             pr->nr--;
 186           }
 187           else
 188             i++;
 189         }
 190       }
 191     }
 192   }
 193
 194   /* Add all non-connecting shakes to the shake list and throw away
 195      the shakenc list */
 196   for (i=0; i<plist[F_SHAKENC].nr; i++)
 197     push_bondnow(&(plist[F_SHAKE]), &(plist[F_SHAKENC].param[i]));
 198   plist[F_SHAKENC].nr=0;
 199   sfree(plist[F_SHAKENC].param);
 200   plist[F_SHAKENC].param=NULL;
 201 }