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