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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
57 /* The source code in this file should be thread-safe.
58 Please keep it that way. */
60 static void gen_waterhydrogen(int nh,rvec xa[], rvec xh[],int *l)
65 const rvec matrix1[6] = {
73 const rvec matrix2[6] = {
87 /* This was copied from Gromos */
88 for(m=0; (m<DIM); m++) {
89 xH1[m]=xAI[m]+matrix1[*l][m];
90 xH2[m]=xAI[m]+matrix2[*l][m];
100 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
102 #define alfaH (acos(-1/3.0)) /* 109.47 degrees */
103 #define alfaHpl (2*M_PI/3) /* 120 degrees */
106 #define alfaCOM (DEG2RAD*117)
107 #define alfaCO (DEG2RAD*121)
108 #define alfaCOA (DEG2RAD*115)
115 real s6,rij,ra,rb,xd;
120 /* common work for constructing one, two or three dihedral hydrogens */
128 for(d=0; (d<DIM); d++) {
135 sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
136 sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
137 sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
139 for(d=0; (d<DIM); d++) {
144 for(d=0; (d<DIM); d++)
147 sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
148 sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
149 sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
154 case 1: /* construct one planar hydrogen (peptide,rings) */
157 for(d=0; (d<DIM); d++) {
158 sij[d] = xAI[d]-xAJ[d];
159 sb[d] = xAI[d]-xAK[d];
166 for(d=0; (d<DIM); d++) {
167 sa[d] = sij[d]/rij+sb[d]/rb;
171 for(d=0; (d<DIM); d++)
172 xH1[d] = xAI[d]+distH*sa[d]/ra;
174 case 2: /* one single hydrogen, e.g. hydroxyl */
175 for(d=0; (d<DIM); d++) {
176 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
179 case 3: /* two planar hydrogens, e.g. -NH2 */
180 for(d=0; (d<DIM); d++) {
181 xH1[d] = xAI[d]-distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
182 xH2[d] = xAI[d]+distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
185 case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
186 for(d=0; (d<DIM); d++) {
187 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
189 - distH*sin(alfaH)*0.5*sb[d]
190 + distH*sin(alfaH)*s6*sa[d]
191 - distH*cos(alfaH)*sij[d] );
192 if ( xH3[XX]!=NOTSET && xH3[YY]!=NOTSET && xH3[ZZ]!=NOTSET )
194 - distH*sin(alfaH)*0.5*sb[d]
195 - distH*sin(alfaH)*s6*sa[d]
196 - distH*cos(alfaH)*sij[d] );
199 case 5: { /* one tetrahedral hydrogen, e.g. C3CH */
203 for(d=0; (d<DIM); d++) {
204 center=(xAJ[d]+xAK[d]+xAL[d])/3.0;
205 dxc[d]=xAI[d]-center;
208 for(d=0; (d<DIM); d++)
209 xH1[d]=xAI[d]+dxc[d]*distH/center;
212 case 6: { /* two tetrahedral hydrogens, e.g. C-CH2-C */
213 rvec rBB,rCC1,rCC2,rNN;
216 for(d=0; (d<DIM); d++)
217 rBB[d]=xAI[d]-0.5*(xAJ[d]+xAK[d]);
220 rvec_sub(xAI,xAJ,rCC1);
221 rvec_sub(xAI,xAK,rCC2);
222 cprod(rCC1,rCC2,rNN);
225 for(d=0; (d<DIM); d++) {
226 xH1[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
227 sin(alfaH/2.0)*rNN[d]/nn);
228 xH2[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
229 sin(alfaH/2.0)*rNN[d]/nn);
233 case 7: /* two water hydrogens */
234 gen_waterhydrogen(2, xa, xh, l);
236 case 10: /* three water hydrogens */
237 gen_waterhydrogen(3, xa, xh, l);
239 case 11: /* four water hydrogens */
240 gen_waterhydrogen(4, xa, xh, l);
242 case 8: /* two carboxyl oxygens, -COO- */
243 for(d=0; (d<DIM); d++) {
244 xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
245 xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
248 case 9: { /* carboxyl oxygens and hydrogen, -COOH */
249 rvec xa2[4]; /* i,j,k,l */
251 /* first add two oxygens */
252 for(d=0; (d<DIM); d++) {
253 xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
254 xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
257 /* now use rule 2 to add hydrogen to 2nd oxygen */
258 copy_rvec(xH2, xa2[0]); /* new i = n' */
259 copy_rvec(xAI, xa2[1]); /* new j = i */
260 copy_rvec(xAJ, xa2[2]); /* new k = j */
261 copy_rvec(xAK, xa2[3]); /* new l = k, not used */
262 calc_h_pos(2, xa2, (xh+2), l);
267 gmx_fatal(FARGS,"Invalid argument (%d) for nht in routine genh\n",nht);