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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Great Red Owns Many ACres of Sand
  31  */
  32 #include "macros.h"
  33 #include "calch.h"
  34 #include "maths.h"
  35 #include "vec.h"
  36 #include "physics.h"
  37
  38 #define xAI xa[0]
  39 #define xAJ xa[1]
  40 #define xAK xa[2]
  41 #define xAL xa[3]
  42 #define xH1 xh[0]
  43 #define xH2 xh[1]
  44 #define xH3 xh[2]
  45
  46 void gen_waterhydrogen(rvec xa[], rvec xh[])
  47 {
  48 #define AA 0.081649
  49 #define BB 0.0
  50 #define CC 0.0577350
  51   const  rvec   matrix1[6] = {
  52     { AA,     BB,     CC },
  53     { AA,     BB,     CC },
  54     { AA,     BB,     CC },
  55     { -AA,    BB,     CC },
  56     { -AA,    BB,     CC },
  57     { BB,     AA,    -CC }
  58   };
  59   const  rvec   matrix2[6] = {
  60     { -AA,   BB,   CC },
  61     { BB,    AA,  -CC },
  62     { BB,   -AA,  -CC },
  63     { BB,    AA,  -CC },
  64     { BB,   -AA,  -CC },
  65     { BB,   -AA,  -CC }
  66   };
  67 #undef AA
  68 #undef BB
  69 #undef CC
  70   static int l=0;
  71   int        m;
  72
  73   /* This was copied from Gromos */
  74   for(m=0; (m<DIM); m++) {
  75     xH1[m]=xAI[m]+matrix1[l][m];
  76     xH2[m]=xAI[m]+matrix2[l][m];
  77   }
  78
  79   l=(l+1) % 6;
  80 }
  81
  82 void calc_h_pos(int nht, rvec xa[], rvec xh[])
  83 {
  84 #define alfaH   (DEG2RAD*109.5)
  85 #define distH   0.1
  86
  87 #define alfaCOM (DEG2RAD*117)
  88 #define alfaCO  (DEG2RAD*121)
  89 #define alfaCOA (DEG2RAD*115)
  90
  91 #define distO   0.123
  92 #define distOA  0.125
  93 #define distOM  0.136
  94
  95   rvec sa,sb,sij;
  96   real s6,rij,ra,rb,xd;
  97   int  d;
  98
  99   s6=0.5*sqrt(3.e0);
 100
 101   /* common work for constructing one, two or three dihedral hydrogens */
 102   switch (nht) {
 103   case 2:
 104   case 3:
 105   case 4:
 106   case 8:
 107   case 9:
 108     rij = 0.e0;
 109     for(d=0; (d<DIM); d++) {
 110       xd     = xAJ[d];
 111       sij[d] = xAI[d]-xd;
 112       sb[d]  = xd-xAK[d];
 113       rij   += sqr(sij[d]);
 114     }
 115     rij = sqrt(rij);
 116     sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
 117     sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
 118     sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
 119     ra = 0.e0;
 120     for(d=0; (d<DIM); d++) {
 121       sij[d] = sij[d]/rij;
 122       ra    += sqr(sa[d]);
 123     }
 124     ra = sqrt(ra);
 125     for(d=0; (d<DIM); d++)
 126       sa[d] = sa[d]/ra;
 127
 128     sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
 129     sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
 130     sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
 131     break;
 132   }/* end switch */
 133
 134   switch (nht) {
 135   case 1: /* construct one planar hydrogen (peptide,rings) */
 136     rij = 0.e0;
 137     rb  = 0.e0;
 138     for(d=0; (d<DIM); d++) {
 139       sij[d] = xAI[d]-xAJ[d];
 140       sb[d]  = xAI[d]-xAK[d];
 141       rij   += sqr(sij[d]);
 142       rb    += sqr(sb[d]);
 143     }
 144     rij = sqrt(rij);
 145     rb  = sqrt(rb);
 146     ra  = 0.e0;
 147     for(d=0; (d<DIM); d++) {
 148       sa[d] = sij[d]/rij+sb[d]/rb;
 149       ra   += sqr(sa[d]);
 150     }
 151     ra = sqrt(ra);
 152     for(d=0; (d<DIM); d++)
 153       xH1[d] = xAI[d]+distH*sa[d]/ra;
 154     break;
 155   case 2: /* one single hydrogen, e.g. hydroxyl */
 156     for(d=0; (d<DIM); d++) {
 157       xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 158     }
 159     break;
 160   case 3: /* two planar hydrogens, e.g. -NH2 */
 161     for(d=0; (d<DIM); d++) {
 162       xH1[d] = xAI[d]-distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 163       xH2[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 164     }
 165     break;
 166   case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
 167     for(d=0; (d<DIM); d++) {
 168       xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 169       xH2[d] = ( xAI[d]
 170                    - distH*sin(alfaH)*0.5*sb[d]
 171                    + distH*sin(alfaH)*s6*sa[d]
 172                    - distH*cos(alfaH)*sij[d] );
 173       if ( xH3[XX]!=NOTSET && xH3[YY]!=NOTSET && xH3[ZZ]!=NOTSET )
 174         xH3[d] = ( xAI[d]
 175                      - distH*sin(alfaH)*0.5*sb[d]
 176                      - distH*sin(alfaH)*s6*sa[d]
 177                      - distH*cos(alfaH)*sij[d] );
 178     }
 179     break;
 180   case 5: { /* one tetrahedral hydrogen, e.g. C3CH */
 181     real center;
 182     rvec dxc;
 183
 184     for(d=0; (d<DIM); d++) {
 185       center=(xAJ[d]+xAK[d]+xAL[d])/3.0;
 186       dxc[d]=xAI[d]-center;
 187     }
 188     center=norm(dxc);
 189     for(d=0; (d<DIM); d++)
 190       xH1[d]=xAI[d]+dxc[d]*distH/center;
 191     break;
 192   }
 193   case 6: { /* two tetrahedral hydrogens, e.g. C-CH2-C */
 194     rvec rBB,rCC1,rCC2,rNN;
 195     real bb,nn;
 196
 197     for(d=0; (d<DIM); d++)
 198       rBB[d]=xAI[d]-0.5*(xAJ[d]+xAK[d]);
 199     bb=norm(rBB);
 200
 201     rvec_sub(xAI,xAJ,rCC1);
 202     rvec_sub(xAI,xAK,rCC2);
 203     oprod(rCC1,rCC2,rNN);
 204     nn=norm(rNN);
 205
 206     for(d=0; (d<DIM); d++) {
 207       xH1[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
 208                                sin(alfaH/2.0)*rNN[d]/nn);
 209       xH2[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
 210                                sin(alfaH/2.0)*rNN[d]/nn);
 211     }
 212     break;
 213   }
 214   case 7: /* two water hydrogens */
 215     gen_waterhydrogen(xa, xh);
 216     break;
 217   case 8: /* two carboxyl oxygens, -COO- */
 218     for(d=0; (d<DIM); d++) {
 219       xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 220       xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 221     }
 222     break;
 223   case 9: { /* carboxyl oxygens and hydrogen, -COOH */
 224     rvec xa2[4]; /* i,j,k,l   */
 225
 226     /* first add two oxygens */
 227     for(d=0; (d<DIM); d++) {
 228       xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
 229       xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
 230     }
 231
 232     /* now use rule 2 to add hydrogen to 2nd oxygen */
 233     copy_rvec(xH2, xa2[0]); /* new i = n' */
 234     copy_rvec(xAI, xa2[1]); /* new j = i  */
 235     copy_rvec(xAJ, xa2[2]); /* new k = j  */
 236     copy_rvec(xAK, xa2[3]); /* new l = k, not used */
 237     calc_h_pos(2, xa2, (xh+2));
 238
 239     break;
 240   }
 241   default:
 242     fatal_error(0,"Invalid argument (%d) for nht in routine genh\n",nht);
 243   } /* end switch */
 244 }