2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 #include "gromacs/legacyheaders/macros.h"
41 #include "gromacs/math/utilities.h"
42 #include "gromacs/math/vec.h"
43 #include "gromacs/math/units.h"
45 #include "gromacs/utility/fatalerror.h"
56 /* The source code in this file should be thread-safe.
57 Please keep it that way. */
59 static void gen_waterhydrogen(int nh, rvec xa[], rvec xh[], int *l)
64 const rvec matrix1[6] = {
72 const rvec matrix2[6] = {
86 /* This was copied from Gromos */
87 for (m = 0; (m < DIM); m++)
89 xH1[m] = xAI[m]+matrix1[*l][m];
90 xH2[m] = xAI[m]+matrix2[*l][m];
104 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
106 #define alfaH (acos(-1/3.0)) /* 109.47 degrees */
107 #define alfaHpl (2*M_PI/3) /* 120 degrees */
110 #define alfaCOM (DEG2RAD*117)
111 #define alfaCO (DEG2RAD*121)
112 #define alfaCOA (DEG2RAD*115)
119 real s6, rij, ra, rb, xd;
124 /* common work for constructing one, two or three dihedral hydrogens */
133 for (d = 0; (d < DIM); d++)
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];
145 for (d = 0; (d < DIM); d++)
151 for (d = 0; (d < DIM); d++)
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];
164 case 1: /* construct one planar hydrogen (peptide,rings) */
167 for (d = 0; (d < DIM); d++)
169 sij[d] = xAI[d]-xAJ[d];
170 sb[d] = xAI[d]-xAK[d];
177 for (d = 0; (d < DIM); d++)
179 sa[d] = sij[d]/rij+sb[d]/rb;
183 for (d = 0; (d < DIM); d++)
185 xH1[d] = xAI[d]+distH*sa[d]/ra;
188 case 2: /* one single hydrogen, e.g. hydroxyl */
189 for (d = 0; (d < DIM); d++)
191 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
194 case 3: /* two planar hydrogens, e.g. -NH2 */
195 for (d = 0; (d < DIM); d++)
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];
201 case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
202 for (d = 0; (d < DIM); d++)
204 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[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)
212 - distH*sin(alfaH)*0.5*sb[d]
213 - distH*sin(alfaH)*s6*sa[d]
214 - distH*cos(alfaH)*sij[d] );
218 case 5: /* one tetrahedral hydrogen, e.g. C3CH */
223 for (d = 0; (d < DIM); d++)
225 center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
226 dxc[d] = xAI[d]-center;
229 for (d = 0; (d < DIM); d++)
231 xH1[d] = xAI[d]+dxc[d]*distH/center;
235 case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
237 rvec rBB, rCC1, rCC2, rNN;
240 for (d = 0; (d < DIM); d++)
242 rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
246 rvec_sub(xAI, xAJ, rCC1);
247 rvec_sub(xAI, xAK, rCC2);
248 cprod(rCC1, rCC2, rNN);
251 for (d = 0; (d < DIM); d++)
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);
260 case 7: /* two water hydrogens */
261 gen_waterhydrogen(2, xa, xh, l);
263 case 10: /* three water hydrogens */
264 gen_waterhydrogen(3, xa, xh, l);
266 case 11: /* four water hydrogens */
267 gen_waterhydrogen(4, xa, xh, l);
269 case 8: /* two carboxyl oxygens, -COO- */
270 for (d = 0; (d < DIM); d++)
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];
276 case 9: /* carboxyl oxygens and hydrogen, -COOH */
278 rvec xa2[4]; /* i,j,k,l */
280 /* first add two oxygens */
281 for (d = 0; (d < DIM); d++)
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];
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);
297 gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);