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