src/gromacs/gmxpreprocess/calch.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2010,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "calch.h"
  40
  41 #include "gromacs/legacyheaders/macros.h"
  42 #include "gromacs/math/units.h"
  43 #include "gromacs/math/utilities.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/utility/fatalerror.h"
  46
  47 #define xAI xa[0]
  48 #define xAJ xa[1]
  49 #define xAK xa[2]
  50 #define xAL xa[3]
  51 #define xH1 xh[0]
  52 #define xH2 xh[1]
  53 #define xH3 xh[2]
  54 #define xH4 xh[3]
  55
  56 /* The source code in this file should be thread-safe.
  57       Please keep it that way. */
  58
  59 static void gen_waterhydrogen(int nh, rvec xa[], rvec xh[], int *l)
  60 {
  61 #define AA 0.081649
  62 #define BB 0.0
  63 #define CC 0.0577350
  64     const  rvec   matrix1[6] = {
  65         { AA,     BB,     CC },
  66         { AA,     BB,     CC },
  67         { AA,     BB,     CC },
  68         { -AA,    BB,     CC },
  69         { -AA,    BB,     CC },
  70         { BB,     AA,    -CC }
  71     };
  72     const  rvec   matrix2[6] = {
  73         { -AA,   BB,   CC },
  74         { BB,    AA,  -CC },
  75         { BB,   -AA,  -CC },
  76         { BB,    AA,  -CC },
  77         { BB,   -AA,  -CC },
  78         { BB,   -AA,  -CC }
  79     };
  80 #undef AA
  81 #undef BB
  82 #undef CC
  83     int        m;
  84     rvec       kkk;
  85
  86     /* This was copied from Gromos */
  87     for (m = 0; (m < DIM); m++)
  88     {
  89         xH1[m] = xAI[m]+matrix1[*l][m];
  90         xH2[m] = xAI[m]+matrix2[*l][m];
  91     }
  92     if (nh > 2)
  93     {
  94         copy_rvec(xAI, xH3);
  95     }
  96     if (nh > 3)
  97     {
  98         copy_rvec(xAI, xH4);
  99     }
 100
 101     *l = (*l+1) % 6;
 102 }
 103
 104 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
 105 {
 106 #define alfaH   (acos(-1/3.0)) /* 109.47 degrees */
 107 #define alfaHpl (2*M_PI/3)     /* 120 degrees */
 108 #define distH   0.1
 109
 110 #define alfaCOM (DEG2RAD*117)
 111 #define alfaCO  (DEG2RAD*121)
 112 #define alfaCOA (DEG2RAD*115)
 113
 114 #define distO   0.123
 115 #define distOA  0.125
 116 #define distOM  0.136
 117
 118     rvec sa, sb, sij;
 119     real s6, rij, ra, rb, xd;
 120     int  d;
 121
 122     s6 = 0.5*sqrt(3.e0);
 123
 124     /* common work for constructing one, two or three dihedral hydrogens */
 125     switch (nht)
 126     {
 127         case 2:
 128         case 3:
 129         case 4:
 130         case 8:
 131         case 9:
 132             rij = 0.e0;
 133             for (d = 0; (d < DIM); d++)
 134             {
 135                 xd     = xAJ[d];
 136                 sij[d] = xAI[d]-xd;
 137                 sb[d]  = xd-xAK[d];
 138                 rij   += sqr(sij[d]);
 139             }
 140             rij    = sqrt(rij);
 141             sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
 142             sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
 143             sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
 144             ra     = 0.e0;
 145             for (d = 0; (d < DIM); d++)
 146             {
 147                 sij[d] = sij[d]/rij;
 148                 ra    += sqr(sa[d]);
 149             }
 150             ra = sqrt(ra);
 151             for (d = 0; (d < DIM); d++)
 152             {
 153                 sa[d] = sa[d]/ra;
 154             }
 155
 156             sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
 157             sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
 158             sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
 159             break;
 160     } /* end switch */
 161
 162     switch (nht)
 163     {
 164         case 1: /* construct one planar hydrogen (peptide,rings) */
 165             rij = 0.e0;
 166             rb  = 0.e0;
 167             for (d = 0; (d < DIM); d++)
 168             {
 169                 sij[d] = xAI[d]-xAJ[d];
 170                 sb[d]  = xAI[d]-xAK[d];
 171                 rij   += sqr(sij[d]);
 172                 rb    += sqr(sb[d]);
 173             }
 174             rij = sqrt(rij);
 175             rb  = sqrt(rb);
 176             ra  = 0.e0;
 177             for (d = 0; (d < DIM); d++)
 178             {
 179                 sa[d] = sij[d]/rij+sb[d]/rb;
 180                 ra   += sqr(sa[d]);
 181             }
 182             ra = sqrt(ra);
 183             for (d = 0; (d < DIM); d++)
 184             {
 185                 xH1[d] = xAI[d]+distH*sa[d]/ra;
 186             }
 187             break;
 188         case 2: /* one single hydrogen, e.g. hydroxyl */
 189             for (d = 0; (d < DIM); d++)
 190             {
 191                 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 192             }
 193             break;
 194         case 3: /* two planar hydrogens, e.g. -NH2 */
 195             for (d = 0; (d < DIM); d++)
 196             {
 197                 xH1[d] = xAI[d]-distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
 198                 xH2[d] = xAI[d]+distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
 199             }
 200             break;
 201         case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
 202             for (d = 0; (d < DIM); d++)
 203             {
 204                 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 205                 xH2[d] = ( xAI[d]
 206                            - distH*sin(alfaH)*0.5*sb[d]
 207                            + distH*sin(alfaH)*s6*sa[d]
 208                            - distH*cos(alfaH)*sij[d] );
 209                 if (xH3[XX] != NOTSET && xH3[YY] != NOTSET && xH3[ZZ] != NOTSET)
 210                 {
 211                     xH3[d] = ( xAI[d]
 212                                - distH*sin(alfaH)*0.5*sb[d]
 213                                - distH*sin(alfaH)*s6*sa[d]
 214                                - distH*cos(alfaH)*sij[d] );
 215                 }
 216             }
 217             break;
 218         case 5: /* one tetrahedral hydrogen, e.g. C3CH */
 219         {
 220             real center;
 221             rvec dxc;
 222
 223             for (d = 0; (d < DIM); d++)
 224             {
 225                 center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
 226                 dxc[d] = xAI[d]-center;
 227             }
 228             center = norm(dxc);
 229             for (d = 0; (d < DIM); d++)
 230             {
 231                 xH1[d] = xAI[d]+dxc[d]*distH/center;
 232             }
 233             break;
 234         }
 235         case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
 236         {
 237             rvec rBB, rCC1, rCC2, rNN;
 238             real bb, nn;
 239
 240             for (d = 0; (d < DIM); d++)
 241             {
 242                 rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
 243             }
 244             bb = norm(rBB);
 245
 246             rvec_sub(xAI, xAJ, rCC1);
 247             rvec_sub(xAI, xAK, rCC2);
 248             cprod(rCC1, rCC2, rNN);
 249             nn = norm(rNN);
 250
 251             for (d = 0; (d < DIM); d++)
 252             {
 253                 xH1[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
 254                                        sin(alfaH/2.0)*rNN[d]/nn);
 255                 xH2[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
 256                                        sin(alfaH/2.0)*rNN[d]/nn);
 257             }
 258             break;
 259         }
 260         case 7: /* two water hydrogens */
 261             gen_waterhydrogen(2, xa, xh, l);
 262             break;
 263         case 10: /* three water hydrogens */
 264             gen_waterhydrogen(3, xa, xh, l);
 265             break;
 266         case 11: /* four water hydrogens */
 267             gen_waterhydrogen(4, xa, xh, l);
 268             break;
 269         case 8: /* two carboxyl oxygens, -COO- */
 270             for (d = 0; (d < DIM); d++)
 271             {
 272                 xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 273                 xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 274             }
 275             break;
 276         case 9:          /* carboxyl oxygens and hydrogen, -COOH */
 277         {
 278             rvec xa2[4]; /* i,j,k,l   */
 279
 280             /* first add two oxygens */
 281             for (d = 0; (d < DIM); d++)
 282             {
 283                 xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
 284                 xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
 285             }
 286
 287             /* now use rule 2 to add hydrogen to 2nd oxygen */
 288             copy_rvec(xH2, xa2[0]); /* new i = n' */
 289             copy_rvec(xAI, xa2[1]); /* new j = i  */
 290             copy_rvec(xAJ, xa2[2]); /* new k = j  */
 291             copy_rvec(xAK, xa2[3]); /* new l = k, not used */
 292             calc_h_pos(2, xa2, (xh+2), l);
 293
 294             break;
 295         }
 296         default:
 297             gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);
 298     } /* end switch */
 299 }