3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
54 /* The source code in this file should be thread-safe.
55 Please keep it that way. */
57 static void gen_waterhydrogen(int nh,rvec xa[], rvec xh[],int *l)
62 const rvec matrix1[6] = {
70 const rvec matrix2[6] = {
84 /* This was copied from Gromos */
85 for(m=0; (m<DIM); m++) {
86 xH1[m]=xAI[m]+matrix1[*l][m];
87 xH2[m]=xAI[m]+matrix2[*l][m];
97 void calc_h_pos(int nht, rvec xa[], rvec xh[], int *l)
99 #define alfaH (acos(-1/3.0)) /* 109.47 degrees */
100 #define alfaHpl (2*M_PI/3) /* 120 degrees */
103 #define alfaCOM (DEG2RAD*117)
104 #define alfaCO (DEG2RAD*121)
105 #define alfaCOA (DEG2RAD*115)
112 real s6,rij,ra,rb,xd;
117 /* common work for constructing one, two or three dihedral hydrogens */
125 for(d=0; (d<DIM); d++) {
132 sa[XX] = sij[YY]*sb[ZZ]-sij[ZZ]*sb[YY];
133 sa[YY] = sij[ZZ]*sb[XX]-sij[XX]*sb[ZZ];
134 sa[ZZ] = sij[XX]*sb[YY]-sij[YY]*sb[XX];
136 for(d=0; (d<DIM); d++) {
141 for(d=0; (d<DIM); d++)
144 sb[XX] = sa[YY]*sij[ZZ]-sa[ZZ]*sij[YY];
145 sb[YY] = sa[ZZ]*sij[XX]-sa[XX]*sij[ZZ];
146 sb[ZZ] = sa[XX]*sij[YY]-sa[YY]*sij[XX];
151 case 1: /* construct one planar hydrogen (peptide,rings) */
154 for(d=0; (d<DIM); d++) {
155 sij[d] = xAI[d]-xAJ[d];
156 sb[d] = xAI[d]-xAK[d];
163 for(d=0; (d<DIM); d++) {
164 sa[d] = sij[d]/rij+sb[d]/rb;
168 for(d=0; (d<DIM); d++)
169 xH1[d] = xAI[d]+distH*sa[d]/ra;
171 case 2: /* one single hydrogen, e.g. hydroxyl */
172 for(d=0; (d<DIM); d++) {
173 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
176 case 3: /* two planar hydrogens, e.g. -NH2 */
177 for(d=0; (d<DIM); d++) {
178 xH1[d] = xAI[d]-distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
179 xH2[d] = xAI[d]+distH*sin(alfaHpl)*sb[d]-distH*cos(alfaHpl)*sij[d];
182 case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
183 for(d=0; (d<DIM); d++) {
184 xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
186 - distH*sin(alfaH)*0.5*sb[d]
187 + distH*sin(alfaH)*s6*sa[d]
188 - distH*cos(alfaH)*sij[d] );
189 if ( xH3[XX]!=NOTSET && xH3[YY]!=NOTSET && xH3[ZZ]!=NOTSET )
191 - distH*sin(alfaH)*0.5*sb[d]
192 - distH*sin(alfaH)*s6*sa[d]
193 - distH*cos(alfaH)*sij[d] );
196 case 5: { /* one tetrahedral hydrogen, e.g. C3CH */
200 for(d=0; (d<DIM); d++) {
201 center=(xAJ[d]+xAK[d]+xAL[d])/3.0;
202 dxc[d]=xAI[d]-center;
205 for(d=0; (d<DIM); d++)
206 xH1[d]=xAI[d]+dxc[d]*distH/center;
209 case 6: { /* two tetrahedral hydrogens, e.g. C-CH2-C */
210 rvec rBB,rCC1,rCC2,rNN;
213 for(d=0; (d<DIM); d++)
214 rBB[d]=xAI[d]-0.5*(xAJ[d]+xAK[d]);
217 rvec_sub(xAI,xAJ,rCC1);
218 rvec_sub(xAI,xAK,rCC2);
219 cprod(rCC1,rCC2,rNN);
222 for(d=0; (d<DIM); d++) {
223 xH1[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
224 sin(alfaH/2.0)*rNN[d]/nn);
225 xH2[d]=xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
226 sin(alfaH/2.0)*rNN[d]/nn);
230 case 7: /* two water hydrogens */
231 gen_waterhydrogen(2, xa, xh, l);
233 case 10: /* three water hydrogens */
234 gen_waterhydrogen(3, xa, xh, l);
236 case 11: /* four water hydrogens */
237 gen_waterhydrogen(4, xa, xh, l);
239 case 8: /* two carboxyl oxygens, -COO- */
240 for(d=0; (d<DIM); d++) {
241 xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
242 xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
245 case 9: { /* carboxyl oxygens and hydrogen, -COOH */
246 rvec xa2[4]; /* i,j,k,l */
248 /* first add two oxygens */
249 for(d=0; (d<DIM); d++) {
250 xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
251 xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
254 /* now use rule 2 to add hydrogen to 2nd oxygen */
255 copy_rvec(xH2, xa2[0]); /* new i = n' */
256 copy_rvec(xAI, xa2[1]); /* new j = i */
257 copy_rvec(xAJ, xa2[2]); /* new k = j */
258 copy_rvec(xAK, xa2[3]); /* new l = k, not used */
259 calc_h_pos(2, xa2, (xh+2), l);
264 gmx_fatal(FARGS,"Invalid argument (%d) for nht in routine genh\n",nht);