src/gmxlib/calch.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  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_calch_c = "$Id$";
  30
  31 #include "macros.h"
  32 #include "calch.h"
  33 #include "maths.h"
  34 #include "vec.h"
  35 #include "physics.h"
  36
  37 #define xAI xa[0]
  38 #define xAJ xa[1]
  39 #define xAK xa[2]
  40 #define xAL xa[3]
  41 #define xH1 xh[0]
  42 #define xH2 xh[1]
  43 #define xH3 xh[2]
  44
  45 void gen_waterhydrogen(rvec xa[], rvec xh[])
  46 {
  47 #define AA 0.081649
  48 #define BB 0.0
  49 #define CC 0.0577350
  50   const  rvec   matrix1[6] = {
  51     { AA,     BB,     CC },
  52     { AA,     BB,     CC },
  53     { AA,     BB,     CC },
  54     { -AA,    BB,     CC },
  55     { -AA,    BB,     CC },
  56     { BB,     AA,    -CC }
  57   };
  58   const  rvec   matrix2[6] = {
  59     { -AA,   BB,   CC },
  60     { BB,    AA,  -CC },
  61     { BB,   -AA,  -CC },
  62     { BB,    AA,  -CC },
  63     { BB,   -AA,  -CC },
  64     { BB,   -AA,  -CC }
  65   };
  66 #undef AA
  67 #undef BB
  68 #undef CC
  69   static int l=0;
  70   int        m;
  71
  72   /* This was copied from Gromos */
  73   for(m=0; (m<DIM); m++) {
  74     xH1[m]=xAI[m]+matrix1[l][m];
  75     xH2[m]=xAI[m]+matrix2[l][m];
  76   }
  77
  78   l=(l+1) % 6;
  79 }
  80
  81 void calc_h_pos(int nht, rvec xa[], rvec xh[])
  82 {
  83 #define alfaH   (DEG2RAD*109.5)
  84 #define distH   0.1
  85
  86 #define alfaCOM (DEG2RAD*117)
  87 #define alfaCO  (DEG2RAD*121)
  88 #define alfaCOA (DEG2RAD*115)
  89
  90 #define distO   0.123
  91 #define distOA  0.125
  92 #define distOM  0.136
  93
  94   rvec sa,sb,sij;
  95   real s6,rij,ra,rb,xd;
  96   int  d;
  97
  98   s6=0.5*sqrt(3.e0);
  99
 100   /* common work for constructing one, two or three dihedral hydrogens */
 101   switch (nht) {
 102   case 2:
 103   case 3:
 104   case 4:
 105   case 8:
 106   case 9:
 107     rij = 0.e0;
 108     for(d=0; (d<DIM); d++) {
 109       xd     = xAJ[d];
 110       sij[d] = xAI[d]-xd;
 111       sb[d]  = xd-xAK[d];
 112       rij   += sqr(sij[d]);
 113     }
 114     rij = sqrt(rij);
 115     sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
 116     sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
 117     sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
 118     ra = 0.e0;
 119     for(d=0; (d<DIM); d++) {
 120       sij[d] = sij[d]/rij;
 121       ra    += sqr(sa[d]);
 122     }
 123     ra = sqrt(ra);
 124     for(d=0; (d<DIM); d++)
 125       sa[d] = sa[d]/ra;
 126
 127     sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
 128     sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
 129     sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
 130     break;
 131   }/* end switch */
 132
 133   switch (nht) {
 134   case 1: /* construct one planar hydrogen (peptide,rings) */
 135     rij = 0.e0;
 136     rb  = 0.e0;
 137     for(d=0; (d<DIM); d++) {
 138       sij[d] = xAI[d]-xAJ[d];
 139       sb[d]  = xAI[d]-xAK[d];
 140       rij   += sqr(sij[d]);
 141       rb    += sqr(sb[d]);
 142     }
 143     rij = sqrt(rij);
 144     rb  = sqrt(rb);
 145     ra  = 0.e0;
 146     for(d=0; (d<DIM); d++) {
 147       sa[d] = sij[d]/rij+sb[d]/rb;
 148       ra   += sqr(sa[d]);
 149     }
 150     ra = sqrt(ra);
 151     for(d=0; (d<DIM); d++)
 152       xH1[d] = xAI[d]+distH*sa[d]/ra;
 153     break;
 154   case 2: /* one single hydrogen, e.g. hydroxyl */
 155     for(d=0; (d<DIM); d++) {
 156       xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 157     }
 158     break;
 159   case 3: /* two planar hydrogens, e.g. -NH2 */
 160     for(d=0; (d<DIM); d++) {
 161       xH1[d] = xAI[d]-distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 162       xH2[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 163     }
 164     break;
 165   case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
 166     for(d=0; (d<DIM); d++) {
 167       xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 168       xH2[d] = ( xAI[d]
 169                    - distH*sin(alfaH)*0.5*sb[d]
 170                    + distH*sin(alfaH)*s6*sa[d]
 171                    - distH*cos(alfaH)*sij[d] );
 172       if ( xH3[XX]!=NOTSET && xH3[YY]!=NOTSET && xH3[ZZ]!=NOTSET )
 173         xH3[d] = ( xAI[d]
 174                      - distH*sin(alfaH)*0.5*sb[d]
 175                      - distH*sin(alfaH)*s6*sa[d]
 176                      - distH*cos(alfaH)*sij[d] );
 177     }
 178     break;
 179   case 5: { /* one tetrahedral hydrogen, e.g. C3CH */
 180     real center;
 181     rvec dxc;
 182
 183     for(d=0; (d<DIM); d++) {
 184       center=(xAJ[d]+xAK[d]+xAL[d])/3.0;
 185       dxc[d]=xAI[d]-center;
 186     }
 187     center=norm(dxc);
 188     for(d=0; (d<DIM); d++)
 189       xH1[d]=xAI[d]+dxc[d]*distH/center;
 190     break;
 191   }
 192   case 6: { /* two tetrahedral hydrogens, e.g. C-CH2-C */
 193     rvec rBB,rCC1,rCC2,rNN;
 194     real bb,nn;
 195
 196     for(d=0; (d<DIM); d++)
 197       rBB[d]=xAI[d]-0.5*(xAJ[d]+xAK[d]);
 198     bb=norm(rBB);
 199
 200     rvec_sub(xAI,xAJ,rCC1);
 201     rvec_sub(xAI,xAK,rCC2);
 202     oprod(rCC1,rCC2,rNN);
 203     nn=norm(rNN);
 204
 205     for(d=0; (d<DIM); d++) {
 206       xH1[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
 207                                sin(alfaH/2.0)*rNN[d]/nn);
 208       xH2[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
 209                                sin(alfaH/2.0)*rNN[d]/nn);
 210     }
 211     break;
 212   }
 213   case 7: /* two water hydrogens */
 214     gen_waterhydrogen(xa, xh);
 215     break;
 216   case 8: /* two carboxyl oxygens, -COO- */
 217     for(d=0; (d<DIM); d++) {
 218       xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 219       xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 220     }
 221     break;
 222   case 9: { /* carboxyl oxygens and hydrogen, -COOH */
 223     rvec xa2[4]; /* i,j,k,l   */
 224
 225     /* first add two oxygens */
 226     for(d=0; (d<DIM); d++) {
 227       xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
 228       xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
 229     }
 230
 231     /* now use rule 2 to add hydrogen to 2nd oxygen */
 232     copy_rvec(xH2, xa2[0]); /* new i = n' */
 233     copy_rvec(xAI, xa2[1]); /* new j = i  */
 234     copy_rvec(xAJ, xa2[2]); /* new k = j  */
 235     copy_rvec(xAK, xa2[3]); /* new l = k, not used */
 236     calc_h_pos(2, xa2, &xH3);
 237
 238     break;
 239   }
 240   default:
 241     fatal_error(0,"Invalid argument (%d) for nht in routine genh\n",nht);
 242   } /* end switch */
 243 }