src/gromacs/gmxpreprocess/topshake.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <ctype.h>
  41
  42 #include "sysstuff.h"
  43 #include "physics.h"
  44 #include "macros.h"
  45 #include "readir.h"
  46 #include "typedefs.h"
  47 #include "topshake.h"
  48 #include "toppush.h"
  49 #include "toputil.h"
  50 #include "topdirs.h"
  51 #include "smalloc.h"
  52 #include "gmx_fatal.h"
  53
  54 static void copy_bond (t_params *pr, int to, int from)
  55 /* copies an entry in a bond list to another position.
  56  * does no allocing or freeing of memory
  57  */
  58 {
  59   /*memcpy((char*) &(pr->param[to]),(char*) &(pr->param[from]),
  60     (size_t)sizeof(pr->param[from]));*/
  61   int i;
  62
  63   if (to != from) {
  64     range_check(to,0,pr->nr);
  65     range_check(from,0,pr->nr);
  66     for(i=0; (i<MAXATOMLIST); i++)
  67       pr->param[to].a[i] = pr->param[from].a[i];
  68     for(i=0; (i<MAXFORCEPARAM); i++)
  69       pr->param[to].c[i] = pr->param[from].c[i];
  70     for(i=0; (i<MAXSLEN); i++)
  71       pr->param[to].s[i] = pr->param[from].s[i];
  72   }
  73 }
  74
  75 static int count_hydrogens (char ***atomname, int nra, atom_id a[])
  76 {
  77   int  i,nh;
  78
  79   if (!atomname)
  80     gmx_fatal(FARGS,"Cannot call count_hydrogens with no atomname (%s %d)",
  81                 __FILE__,__LINE__);
  82
  83   nh=0;
  84   for (i=0; (i<nra); i++)
  85     if (toupper(**(atomname[a[i]]))=='H')
  86       nh++;
  87
  88   return nh;
  89 }
  90
  91 void make_shake (t_params plist[],t_atoms *atoms,gpp_atomtype_t at,int nshake)
  92 {
  93   char         ***info=atoms->atomname;
  94   t_params     *pr;
  95   t_params     *bonds;
  96   t_param      p,*bond,*ang;
  97   real         b_ij,b_jk;
  98   int          nb,b,i,j,ftype,ftype_a;
  99   gmx_bool         bFound;
 100
 101   if (nshake != eshNONE) {
 102     switch (nshake) {
 103     case eshHBONDS:
 104       printf("turning H bonds into constraints...\n");
 105       break;
 106     case eshALLBONDS:
 107       printf("turning all bonds into constraints...\n");
 108       break;
 109     case eshHANGLES:
 110       printf("turning all bonds and H angles into constraints...\n");
 111       break;
 112     case eshALLANGLES:
 113       printf("turning all bonds and angles into constraints...\n");
 114       break;
 115     default:
 116       gmx_fatal(FARGS,"Invalid option for make_shake (%d)",nshake);
 117     }
 118
 119     if ((nshake == eshHANGLES) || (nshake == eshALLANGLES)) {
 120       /* Add all the angles with hydrogens to the shake list
 121        * and remove them from the bond list
 122        */
 123       for(ftype = 0; (ftype < F_NRE); ftype++) {
 124         if (interaction_function[ftype].flags & IF_BTYPE) {
 125           bonds=&(plist[ftype]);
 126
 127           for(ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) {
 128             if (interaction_function[ftype_a].flags & IF_ATYPE) {
 129               pr = &(plist[ftype_a]);
 130
 131               for (i=0; (i < pr->nr); ) {
 132                 int numhydrogens;
 133
 134                 ang=&(pr->param[i]);
 135 #ifdef DEBUG
 136                 printf("Angle: %d-%d-%d\n",ang->AI,ang->AJ,ang->AK);
 137 #endif
 138                 numhydrogens = count_hydrogens(info,3,ang->a);
 139                 if ((nshake == eshALLANGLES) ||
 140                     (numhydrogens > 1) ||
 141                     (numhydrogens == 1 && toupper(**(info[ang->a[1]]))=='O')) {
 142                   /* Can only add hydrogen angle shake, if the two bonds
 143                    * are constrained.
 144                    * append this angle to the shake list
 145                    */
 146                   p.AI = ang->AI;
 147                   p.AJ = ang->AK;
 148
 149                   /* Calculate length of constraint */
 150                   bFound=FALSE;
 151                   b_ij=b_jk=0.0;
 152                   for (j=0; !bFound && (j<bonds->nr); j++) {
 153                     bond=&(bonds->param[j]);
 154                     if (((bond->AI==ang->AI) &&
 155                          (bond->AJ==ang->AJ)) ||
 156                         ((bond->AI==ang->AJ) &&
 157                          (bond->AJ==ang->AI)))
 158                       b_ij=bond->C0;
 159                     if (((bond->AI==ang->AK) &&
 160                          (bond->AJ==ang->AJ)) ||
 161                         ((bond->AI==ang->AJ) &&
 162                          (bond->AJ==ang->AK)))
 163                       b_jk=bond->C0;
 164                     bFound = (b_ij!=0.0) && (b_jk!=0.0);
 165                   }
 166                   if (bFound) {
 167                     /* apply law of cosines */
 168                     p.C0 = sqrt( b_ij*b_ij + b_jk*b_jk -
 169                                  2.0*b_ij*b_jk*cos(DEG2RAD*ang->C0) );
 170                     p.C1 = p.C0;
 171 #ifdef DEBUG
 172                     printf("p: %d, q: %d, dist: %12.5e\n",p.AI,p.AJ,p.C0);
 173 #endif
 174                     add_param_to_list (&(plist[F_CONSTR]),&p);
 175                     /* move the last bond to this position */
 176                     copy_bond (pr,i,pr->nr-1);
 177                     /* should free memory here!! */
 178                     pr->nr--;
 179                   }
 180                 }
 181                 else
 182                   i++;
 183               }
 184             } /* if IF_ATYPE */
 185           } /* for ftype_A */
 186         } /* if IF_BTYPE */
 187       } /* for ftype */
 188     } /* if shake angles */
 189
 190     /* Add all the bonds with hydrogens to the shake list
 191      * and remove them from the bond list
 192      */
 193     for (ftype=0; (ftype < F_NRE); ftype++) {
 194       if (interaction_function[ftype].flags & IF_BTYPE) {
 195         pr = &(plist[ftype]);
 196         j = 0;
 197         for(i=0; i<pr->nr; i++) {
 198           if ( (nshake != eshHBONDS) ||
 199                (count_hydrogens (info,2,pr->param[i].a) > 0) ) {
 200             /* append this bond to the shake list */
 201             p.AI = pr->param[i].AI;
 202             p.AJ = pr->param[i].AJ;
 203             p.C0 = pr->param[i].C0;
 204             p.C1 = pr->param[i].C2;
 205             add_param_to_list (&(plist[F_CONSTR]),&p);
 206           } else {
 207               copy_bond(pr,j++,i);
 208           }
 209         }
 210         pr->nr = j;
 211       }
 212     }
 213   }
 214 }