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.
41 #include "gromacs/legacyheaders/macros.h"
43 #include "gromacs/math/utilities.h"
44 #include "gromacs/math/vec.h"
45 #include "gromacs/math/units.h"
47 #include "gromacs/utility/fatalerror.h"
58 /* The source code in this file should be thread-safe.
59 Please keep it that way. */
61 static void gen_waterhydrogen(int nh, rvec xa[], rvec xh[], int *l)
66 const rvec matrix1[6] = {
74 const rvec matrix2[6] = {
88 /* This was copied from Gromos */
89 for (m = 0; (m < DIM); m++)
91 xH1[m] = xAI[m]+matrix1[*l][m];
92 xH2[m] = xAI[m]+matrix2[*l][m];
106 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
108 #define alfaH (acos(-1/3.0)) /* 109.47 degrees */
109 #define alfaHpl (2*M_PI/3) /* 120 degrees */
112 #define alfaCOM (DEG2RAD*117)
113 #define alfaCO (DEG2RAD*121)
114 #define alfaCOA (DEG2RAD*115)
121 real s6, rij, ra, rb, xd;
126 /* common work for constructing one, two or three dihedral hydrogens */
135 for (d = 0; (d < DIM); d++)
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];
147 for (d = 0; (d < DIM); d++)
153 for (d = 0; (d < DIM); d++)
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];
166 case 1: /* construct one planar hydrogen (peptide,rings) */
169 for (d = 0; (d < DIM); d++)
171 sij[d] = xAI[d]-xAJ[d];
172 sb[d] = xAI[d]-xAK[d];
179 for (d = 0; (d < DIM); d++)
181 sa[d] = sij[d]/rij+sb[d]/rb;
185 for (d = 0; (d < DIM); d++)
187 xH1[d] = xAI[d]+distH*sa[d]/ra;
190 case 2: /* one single hydrogen, e.g. hydroxyl */
191 for (d = 0; (d < DIM); d++)
193 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
196 case 3: /* two planar hydrogens, e.g. -NH2 */
197 for (d = 0; (d < DIM); d++)
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];
203 case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
204 for (d = 0; (d < DIM); d++)
206 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[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)
214 - distH*sin(alfaH)*0.5*sb[d]
215 - distH*sin(alfaH)*s6*sa[d]
216 - distH*cos(alfaH)*sij[d] );
220 case 5: /* one tetrahedral hydrogen, e.g. C3CH */
225 for (d = 0; (d < DIM); d++)
227 center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
228 dxc[d] = xAI[d]-center;
231 for (d = 0; (d < DIM); d++)
233 xH1[d] = xAI[d]+dxc[d]*distH/center;
237 case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
239 rvec rBB, rCC1, rCC2, rNN;
242 for (d = 0; (d < DIM); d++)
244 rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
248 rvec_sub(xAI, xAJ, rCC1);
249 rvec_sub(xAI, xAK, rCC2);
250 cprod(rCC1, rCC2, rNN);
253 for (d = 0; (d < DIM); d++)
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);
262 case 7: /* two water hydrogens */
263 gen_waterhydrogen(2, xa, xh, l);
265 case 10: /* three water hydrogens */
266 gen_waterhydrogen(3, xa, xh, l);
268 case 11: /* four water hydrogens */
269 gen_waterhydrogen(4, xa, xh, l);
271 case 8: /* two carboxyl oxygens, -COO- */
272 for (d = 0; (d < DIM); d++)
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];
278 case 9: /* carboxyl oxygens and hydrogen, -COOH */
280 rvec xa2[4]; /* i,j,k,l */
282 /* first add two oxygens */
283 for (d = 0; (d < DIM); d++)
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];
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);
299 gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);