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