src/gromacs/gmxpreprocess/calch.cpp

   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,2015, 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 <cmath>
  42
  43 #include "gromacs/gmxpreprocess/notset.h"
  44 #include "gromacs/math/units.h"
  45 #include "gromacs/math/utilities.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/utility/fatalerror.h"
  48
  49 #define xAI xa[0]
  50 #define xAJ xa[1]
  51 #define xAK xa[2]
  52 #define xAL xa[3]
  53 #define xH1 xh[0]
  54 #define xH2 xh[1]
  55 #define xH3 xh[2]
  56 #define xH4 xh[3]
  57
  58 /* The source code in this file should be thread-safe.
  59       Please keep it that way. */
  60
  61 static void gen_waterhydrogen(int nh, rvec xa[], rvec xh[], int *l)
  62 {
  63 #define AA 0.081649
  64 #define BB 0.0
  65 #define CC 0.0577350
  66     const  rvec   matrix1[6] = {
  67         { AA,     BB,     CC },
  68         { AA,     BB,     CC },
  69         { AA,     BB,     CC },
  70         { -AA,    BB,     CC },
  71         { -AA,    BB,     CC },
  72         { BB,     AA,    -CC }
  73     };
  74     const  rvec   matrix2[6] = {
  75         { -AA,   BB,   CC },
  76         { BB,    AA,  -CC },
  77         { BB,   -AA,  -CC },
  78         { BB,    AA,  -CC },
  79         { BB,   -AA,  -CC },
  80         { BB,   -AA,  -CC }
  81     };
  82 #undef AA
  83 #undef BB
  84 #undef CC
  85     int        m;
  86
  87     /* This was copied from Gromos */
  88     for (m = 0; (m < DIM); m++)
  89     {
  90         xH1[m] = xAI[m]+matrix1[*l][m];
  91         xH2[m] = xAI[m]+matrix2[*l][m];
  92     }
  93     if (nh > 2)
  94     {
  95         copy_rvec(xAI, xH3);
  96     }
  97     if (nh > 3)
  98     {
  99         copy_rvec(xAI, xH4);
 100     }
 101
 102     *l = (*l+1) % 6;
 103 }
 104
 105 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
 106 {
 107 #define alfaH   (acos(-1/3.0)) /* 109.47 degrees */
 108 #define alfaHpl (2*M_PI/3)     /* 120 degrees */
 109 #define distH   0.1
 110
 111 #define alfaCOM (DEG2RAD*117)
 112 #define alfaCO  (DEG2RAD*121)
 113 #define alfaCOA (DEG2RAD*115)
 114
 115 #define distO   0.123
 116 #define distOA  0.125
 117 #define distOM  0.136
 118
 119     rvec sa, sb, sij;
 120     real s6, rij, ra, rb, xd;
 121     int  d;
 122
 123     s6 = 0.5*std::sqrt(3.e0);
 124
 125     /* common work for constructing one, two or three dihedral hydrogens */
 126     switch (nht)
 127     {
 128         case 2:
 129         case 3:
 130         case 4:
 131         case 8:
 132         case 9:
 133             rij = 0.e0;
 134             for (d = 0; (d < DIM); d++)
 135             {
 136                 xd     = xAJ[d];
 137                 sij[d] = xAI[d]-xd;
 138                 sb[d]  = xd-xAK[d];
 139                 rij   += sqr(sij[d]);
 140             }
 141             rij    = std::sqrt(rij);
 142             sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
 143             sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
 144             sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
 145             ra     = 0.e0;
 146             for (d = 0; (d < DIM); d++)
 147             {
 148                 sij[d] = sij[d]/rij;
 149                 ra    += sqr(sa[d]);
 150             }
 151             ra = std::sqrt(ra);
 152             for (d = 0; (d < DIM); d++)
 153             {
 154                 sa[d] = sa[d]/ra;
 155             }
 156
 157             sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
 158             sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
 159             sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
 160             break;
 161     } /* end switch */
 162
 163     switch (nht)
 164     {
 165         case 1: /* construct one planar hydrogen (peptide,rings) */
 166             rij = 0.e0;
 167             rb  = 0.e0;
 168             for (d = 0; (d < DIM); d++)
 169             {
 170                 sij[d] = xAI[d]-xAJ[d];
 171                 sb[d]  = xAI[d]-xAK[d];
 172                 rij   += sqr(sij[d]);
 173                 rb    += sqr(sb[d]);
 174             }
 175             rij = std::sqrt(rij);
 176             rb  = std::sqrt(rb);
 177             ra  = 0.e0;
 178             for (d = 0; (d < DIM); d++)
 179             {
 180                 sa[d] = sij[d]/rij+sb[d]/rb;
 181                 ra   += sqr(sa[d]);
 182             }
 183             ra = std::sqrt(ra);
 184             for (d = 0; (d < DIM); d++)
 185             {
 186                 xH1[d] = xAI[d]+distH*sa[d]/ra;
 187             }
 188             break;
 189         case 2: /* one single hydrogen, e.g. hydroxyl */
 190             for (d = 0; (d < DIM); d++)
 191             {
 192                 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 193             }
 194             break;
 195         case 3: /* two planar hydrogens, e.g. -NH2 */
 196             for (d = 0; (d < DIM); d++)
 197             {
 198                 xH1[d] = xAI[d]-distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
 199                 xH2[d] = xAI[d]+distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
 200             }
 201             break;
 202         case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
 203             for (d = 0; (d < DIM); d++)
 204             {
 205                 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
 206                 xH2[d] = ( xAI[d]
 207                            - distH*sin(alfaH)*0.5*sb[d]
 208                            + distH*sin(alfaH)*s6*sa[d]
 209                            - distH*cos(alfaH)*sij[d] );
 210                 if (xH3[XX] != NOTSET && xH3[YY] != NOTSET && xH3[ZZ] != NOTSET)
 211                 {
 212                     xH3[d] = ( xAI[d]
 213                                - distH*sin(alfaH)*0.5*sb[d]
 214                                - distH*sin(alfaH)*s6*sa[d]
 215                                - distH*cos(alfaH)*sij[d] );
 216                 }
 217             }
 218             break;
 219         case 5: /* one tetrahedral hydrogen, e.g. C3CH */
 220         {
 221             real center;
 222             rvec dxc;
 223
 224             for (d = 0; (d < DIM); d++)
 225             {
 226                 center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
 227                 dxc[d] = xAI[d]-center;
 228             }
 229             center = norm(dxc);
 230             for (d = 0; (d < DIM); d++)
 231             {
 232                 xH1[d] = xAI[d]+dxc[d]*distH/center;
 233             }
 234             break;
 235         }
 236         case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
 237         {
 238             rvec rBB, rCC1, rCC2, rNN;
 239             real bb, nn;
 240
 241             for (d = 0; (d < DIM); d++)
 242             {
 243                 rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
 244             }
 245             bb = norm(rBB);
 246
 247             rvec_sub(xAI, xAJ, rCC1);
 248             rvec_sub(xAI, xAK, rCC2);
 249             cprod(rCC1, rCC2, rNN);
 250             nn = norm(rNN);
 251
 252             for (d = 0; (d < DIM); d++)
 253             {
 254                 xH1[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
 255                                        sin(alfaH/2.0)*rNN[d]/nn);
 256                 xH2[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
 257                                        sin(alfaH/2.0)*rNN[d]/nn);
 258             }
 259             break;
 260         }
 261         case 7: /* two water hydrogens */
 262             gen_waterhydrogen(2, xa, xh, l);
 263             break;
 264         case 10: /* three water hydrogens */
 265             gen_waterhydrogen(3, xa, xh, l);
 266             break;
 267         case 11: /* four water hydrogens */
 268             gen_waterhydrogen(4, xa, xh, l);
 269             break;
 270         case 8: /* two carboxyl oxygens, -COO- */
 271             for (d = 0; (d < DIM); d++)
 272             {
 273                 xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 274                 xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
 275             }
 276             break;
 277         case 9:          /* carboxyl oxygens and hydrogen, -COOH */
 278         {
 279             rvec xa2[4]; /* i,j,k,l   */
 280
 281             /* first add two oxygens */
 282             for (d = 0; (d < DIM); d++)
 283             {
 284                 xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
 285                 xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
 286             }
 287
 288             /* now use rule 2 to add hydrogen to 2nd oxygen */
 289             copy_rvec(xH2, xa2[0]); /* new i = n' */
 290             copy_rvec(xAI, xa2[1]); /* new j = i  */
 291             copy_rvec(xAJ, xa2[2]); /* new k = j  */
 292             copy_rvec(xAK, xa2[3]); /* new l = k, not used */
 293             calc_h_pos(2, xa2, (xh+2), l);
 294
 295             break;
 296         }
 297         default:
 298             gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);
 299     } /* end switch */
 300 }