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,2015, 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.
43 #include "gromacs/gmxpreprocess/notset.h"
44 #include "gromacs/math/units.h"
45 #include "gromacs/math/utilities.h"
46 #include "gromacs/math/vec.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] = {
87 /* This was copied from Gromos */
88 for (m
= 0; (m
< DIM
); m
++)
90 xH1
[m
] = xAI
[m
]+matrix1
[*l
][m
];
91 xH2
[m
] = xAI
[m
]+matrix2
[*l
][m
];
105 void calc_h_pos(int nht
, rvec xa
[], rvec xh
[], int *l
)
107 #define alfaH (acos(-1/3.0)) /* 109.47 degrees */
108 #define alfaHpl (2*M_PI/3) /* 120 degrees */
111 #define alfaCOM (DEG2RAD*117)
112 #define alfaCO (DEG2RAD*121)
113 #define alfaCOA (DEG2RAD*115)
120 real s6
, rij
, ra
, rb
, xd
;
123 s6
= 0.5*std::sqrt(3.e0
);
125 /* common work for constructing one, two or three dihedral hydrogens */
134 for (d
= 0; (d
< DIM
); d
++)
141 rij
= std::sqrt(rij
);
142 sa
[XX
] = sij
[YY
]*sb
[ZZ
]-sij
[ZZ
]*sb
[YY
];
143 sa
[YY
] = sij
[ZZ
]*sb
[XX
]-sij
[XX
]*sb
[ZZ
];
144 sa
[ZZ
] = sij
[XX
]*sb
[YY
]-sij
[YY
]*sb
[XX
];
146 for (d
= 0; (d
< DIM
); d
++)
152 for (d
= 0; (d
< DIM
); d
++)
157 sb
[XX
] = sa
[YY
]*sij
[ZZ
]-sa
[ZZ
]*sij
[YY
];
158 sb
[YY
] = sa
[ZZ
]*sij
[XX
]-sa
[XX
]*sij
[ZZ
];
159 sb
[ZZ
] = sa
[XX
]*sij
[YY
]-sa
[YY
]*sij
[XX
];
165 case 1: /* construct one planar hydrogen (peptide,rings) */
168 for (d
= 0; (d
< DIM
); d
++)
170 sij
[d
] = xAI
[d
]-xAJ
[d
];
171 sb
[d
] = xAI
[d
]-xAK
[d
];
175 rij
= std::sqrt(rij
);
178 for (d
= 0; (d
< DIM
); d
++)
180 sa
[d
] = sij
[d
]/rij
+sb
[d
]/rb
;
184 for (d
= 0; (d
< DIM
); d
++)
186 xH1
[d
] = xAI
[d
]+distH
*sa
[d
]/ra
;
189 case 2: /* one single hydrogen, e.g. hydroxyl */
190 for (d
= 0; (d
< DIM
); d
++)
192 xH1
[d
] = xAI
[d
]+distH
*sin(alfaH
)*sb
[d
]-distH
*cos(alfaH
)*sij
[d
];
195 case 3: /* two planar hydrogens, e.g. -NH2 */
196 for (d
= 0; (d
< DIM
); d
++)
198 xH1
[d
] = xAI
[d
]-distH
*sin(alfaHpl
)*sb
[d
]-distH
*cos(alfaHpl
)*sij
[d
];
199 xH2
[d
] = xAI
[d
]+distH
*sin(alfaHpl
)*sb
[d
]-distH
*cos(alfaHpl
)*sij
[d
];
202 case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
203 for (d
= 0; (d
< DIM
); d
++)
205 xH1
[d
] = xAI
[d
]+distH
*sin(alfaH
)*sb
[d
]-distH
*cos(alfaH
)*sij
[d
];
207 - distH
*sin(alfaH
)*0.5*sb
[d
]
208 + distH
*sin(alfaH
)*s6
*sa
[d
]
209 - distH
*cos(alfaH
)*sij
[d
] );
210 if (xH3
[XX
] != NOTSET
&& xH3
[YY
] != NOTSET
&& xH3
[ZZ
] != NOTSET
)
213 - distH
*sin(alfaH
)*0.5*sb
[d
]
214 - distH
*sin(alfaH
)*s6
*sa
[d
]
215 - distH
*cos(alfaH
)*sij
[d
] );
219 case 5: /* one tetrahedral hydrogen, e.g. C3CH */
224 for (d
= 0; (d
< DIM
); d
++)
226 center
= (xAJ
[d
]+xAK
[d
]+xAL
[d
])/3.0;
227 dxc
[d
] = xAI
[d
]-center
;
230 for (d
= 0; (d
< DIM
); d
++)
232 xH1
[d
] = xAI
[d
]+dxc
[d
]*distH
/center
;
236 case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
238 rvec rBB
, rCC1
, rCC2
, rNN
;
241 for (d
= 0; (d
< DIM
); d
++)
243 rBB
[d
] = xAI
[d
]-0.5*(xAJ
[d
]+xAK
[d
]);
247 rvec_sub(xAI
, xAJ
, rCC1
);
248 rvec_sub(xAI
, xAK
, rCC2
);
249 cprod(rCC1
, rCC2
, rNN
);
252 for (d
= 0; (d
< DIM
); d
++)
254 xH1
[d
] = xAI
[d
]+distH
*(cos(alfaH
/2.0)*rBB
[d
]/bb
+
255 sin(alfaH
/2.0)*rNN
[d
]/nn
);
256 xH2
[d
] = xAI
[d
]+distH
*(cos(alfaH
/2.0)*rBB
[d
]/bb
-
257 sin(alfaH
/2.0)*rNN
[d
]/nn
);
261 case 7: /* two water hydrogens */
262 gen_waterhydrogen(2, xa
, xh
, l
);
264 case 10: /* three water hydrogens */
265 gen_waterhydrogen(3, xa
, xh
, l
);
267 case 11: /* four water hydrogens */
268 gen_waterhydrogen(4, xa
, xh
, l
);
270 case 8: /* two carboxyl oxygens, -COO- */
271 for (d
= 0; (d
< DIM
); d
++)
273 xH1
[d
] = xAI
[d
]-distOM
*sin(alfaCOM
)*sb
[d
]-distOM
*cos(alfaCOM
)*sij
[d
];
274 xH2
[d
] = xAI
[d
]+distOM
*sin(alfaCOM
)*sb
[d
]-distOM
*cos(alfaCOM
)*sij
[d
];
277 case 9: /* carboxyl oxygens and hydrogen, -COOH */
279 rvec xa2
[4]; /* i,j,k,l */
281 /* first add two oxygens */
282 for (d
= 0; (d
< DIM
); d
++)
284 xH1
[d
] = xAI
[d
]-distO
*sin(alfaCO
)*sb
[d
]-distO
*cos(alfaCO
)*sij
[d
];
285 xH2
[d
] = xAI
[d
]+distOA
*sin(alfaCOA
)*sb
[d
]-distOA
*cos(alfaCOA
)*sij
[d
];
288 /* now use rule 2 to add hydrogen to 2nd oxygen */
289 copy_rvec(xH2
, xa2
[0]); /* new i = n' */
290 copy_rvec(xAI
, xa2
[1]); /* new j = i */
291 copy_rvec(xAJ
, xa2
[2]); /* new k = j */
292 copy_rvec(xAK
, xa2
[3]); /* new l = k, not used */
293 calc_h_pos(2, xa2
, (xh
+2), l
);
298 gmx_fatal(FARGS
, "Invalid argument (%d) for nht in routine genh\n", nht
);