Bug Summary

File:gromacs/gmxpreprocess/gen_vsite.c
Location:line 740, column 5
Description:Value stored to 'xCE1' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <stdio.h>
43#include <stdlib.h>
44#include <string.h>
45
46#include "gen_vsite.h"
47#include "resall.h"
48#include "add_par.h"
49#include "gromacs/math/vec.h"
50#include "toputil.h"
51#include "physics.h"
52#include "index.h"
53#include "names.h"
54#include "gromacs/utility/futil.h"
55#include "gpp_atomtype.h"
56#include "fflibutil.h"
57#include "macros.h"
58
59#include "gromacs/utility/cstringutil.h"
60#include "gromacs/utility/fatalerror.h"
61#include "gromacs/utility/smalloc.h"
62
63#define MAXNAME32 32
64#define OPENDIR'[' '[' /* starting sign for directive */
65#define CLOSEDIR']' ']' /* ending sign for directive */
66
67typedef struct {
68 char atomtype[MAXNAME32]; /* Type for the XH3/XH2 atom */
69 gmx_bool isplanar; /* If true, the atomtype above and the three connected
70 * ones are in a planar geometry. The two next entries
71 * are undefined in that case
72 */
73 int nhydrogens; /* number of connected hydrogens */
74 char nextheavytype[MAXNAME32]; /* Type for the heavy atom bonded to XH2/XH3 */
75 char dummymass[MAXNAME32]; /* The type of MNH* or MCH3* dummy mass to use */
76} t_vsiteconf;
77
78
79/* Structure to represent average bond and angles values in vsite aromatic
80 * residues. Note that these are NOT necessarily the bonds and angles from the
81 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
82 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
83 */
84typedef struct {
85 char resname[MAXNAME32];
86 int nbonds;
87 int nangles;
88 struct vsitetop_bond {
89 char atom1[MAXNAME32];
90 char atom2[MAXNAME32];
91 float value;
92 } *bond; /* list of bonds */
93 struct vsitetop_angle {
94 char atom1[MAXNAME32];
95 char atom2[MAXNAME32];
96 char atom3[MAXNAME32];
97 float value;
98 } *angle; /* list of angles */
99} t_vsitetop;
100
101
102enum {
103 DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
104 DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
105};
106
107typedef char t_dirname[STRLEN4096];
108
109static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
110 "CH3",
111 "NH3",
112 "NH2",
113 "PHE",
114 "TYR",
115 "TRP",
116 "HISA",
117 "HISB",
118 "HISH"
119};
120
121static int ddb_name2dir(char *name)
122{
123 /* Translate a directive name to the number of the directive.
124 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
125 */
126
127 int i, index;
128
129 index = -1;
130
131 for (i = 0; i < DDB_DIR_NR && index < 0; i++)
132 {
133 if (!gmx_strcasecmp(name, ddb_dirnames[i]))
134 {
135 index = i;
136 }
137 }
138
139 return index;
140}
141
142
143static void read_vsite_database(const char *ddbname,
144 t_vsiteconf **pvsiteconflist, int *nvsiteconf,
145 t_vsitetop **pvsitetoplist, int *nvsitetop)
146{
147 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
148 * and aromatic vsite parameters by reading them from a ff???.vsd file.
149 *
150 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
151 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
152 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
153 * the type of the next heavy atom it is bonded to, and the third field the type
154 * of dummy mass that will be used for this group.
155 *
156 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
157 * case the second field should just be the word planar.
158 */
159
160 FILE *ddb;
161 char dirstr[STRLEN4096];
162 char pline[STRLEN4096];
163 int i, j, n, k, nvsite, ntop, curdir, prevdir;
164 t_vsiteconf *vsiteconflist;
165 t_vsitetop *vsitetoplist;
166 char *ch;
167 char s1[MAXNAME32], s2[MAXNAME32], s3[MAXNAME32], s4[MAXNAME32];
168
169 ddb = libopen(ddbname);
170
171 nvsite = *nvsiteconf;
172 vsiteconflist = *pvsiteconflist;
173 ntop = *nvsitetop;
174 vsitetoplist = *pvsitetoplist;
175
176 curdir = -1;
177
178 snew(vsiteconflist, 1)(vsiteconflist) = save_calloc("vsiteconflist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 178, (1), sizeof(*(vsiteconflist)));
179 snew(vsitetoplist, 1)(vsitetoplist) = save_calloc("vsitetoplist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 179, (1), sizeof(*(vsitetoplist)));
180
181 while (fgets2(pline, STRLEN4096-2, ddb) != NULL((void*)0))
182 {
183 strip_comment(pline);
184 trim(pline);
185 if (strlen(pline) > 0)
186 {
187 if (pline[0] == OPENDIR'[')
188 {
189 strncpy(dirstr, pline+1, STRLEN-2)__builtin_strncpy (dirstr, pline+1, 4096 -2);
190 if ((ch = strchr (dirstr, CLOSEDIR)(__extension__ (__builtin_constant_p (']') && !__builtin_constant_p
(dirstr) && (']') == '\0' ? (char *) __rawmemchr (dirstr
, ']') : __builtin_strchr (dirstr, ']')))) != NULL((void*)0))
191 {
192 (*ch) = 0;
193 }
194 trim (dirstr);
195
196 if (!gmx_strcasecmp(dirstr, "HID") ||
197 !gmx_strcasecmp(dirstr, "HISD"))
198 {
199 sprintf(dirstr, "HISA");
200 }
201 else if (!gmx_strcasecmp(dirstr, "HIE") ||
202 !gmx_strcasecmp(dirstr, "HISE"))
203 {
204 sprintf(dirstr, "HISB");
205 }
206 else if (!gmx_strcasecmp(dirstr, "HIP"))
207 {
208 sprintf(dirstr, "HISH");
209 }
210
211 curdir = ddb_name2dir(dirstr);
212 if (curdir < 0)
213 {
214 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 214, "Invalid directive %s in vsite database %s",
215 dirstr, ddbname);
216 }
217 }
218 else
219 {
220 switch (curdir)
221 {
222 case -1:
223 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 223, "First entry in vsite database must be a directive.\n");
224 break;
225 case DDB_CH3:
226 case DDB_NH3:
227 case DDB_NH2:
228 n = sscanf(pline, "%s%s%s", s1, s2, s3);
229 if (n < 3 && !gmx_strcasecmp(s2, "planar"))
230 {
231 srenew(vsiteconflist, nvsite+1)(vsiteconflist) = save_realloc("vsiteconflist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 231, (vsiteconflist), (nvsite+1), sizeof(*(vsiteconflist)));
232 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1)__builtin_strncpy (vsiteconflist[nvsite].atomtype, s1, 32 -1);
233 vsiteconflist[nvsite].isplanar = TRUE1;
234 vsiteconflist[nvsite].nextheavytype[0] = 0;
235 vsiteconflist[nvsite].dummymass[0] = 0;
236 vsiteconflist[nvsite].nhydrogens = 2;
237 nvsite++;
238 }
239 else if (n == 3)
240 {
241 srenew(vsiteconflist, (nvsite+1))(vsiteconflist) = save_realloc("vsiteconflist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 241, (vsiteconflist), ((nvsite+1)), sizeof(*(vsiteconflist)
));
242 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1)__builtin_strncpy (vsiteconflist[nvsite].atomtype, s1, 32 -1);
243 vsiteconflist[nvsite].isplanar = FALSE0;
244 strncpy(vsiteconflist[nvsite].nextheavytype, s2, MAXNAME-1)__builtin_strncpy (vsiteconflist[nvsite].nextheavytype, s2, 32
-1);
245 strncpy(vsiteconflist[nvsite].dummymass, s3, MAXNAME-1)__builtin_strncpy (vsiteconflist[nvsite].dummymass, s3, 32 -1
);
246 if (curdir == DDB_NH2)
247 {
248 vsiteconflist[nvsite].nhydrogens = 2;
249 }
250 else
251 {
252 vsiteconflist[nvsite].nhydrogens = 3;
253 }
254 nvsite++;
255 }
256 else
257 {
258 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 258, "Not enough directives in vsite database line: %s\n", pline);
259 }
260 break;
261 case DDB_PHE:
262 case DDB_TYR:
263 case DDB_TRP:
264 case DDB_HISA:
265 case DDB_HISB:
266 case DDB_HISH:
267 i = 0;
268 while ((i < ntop) && gmx_strcasecmp(dirstr, vsitetoplist[i].resname))
269 {
270 i++;
271 }
272 /* Allocate a new topology entry if this is a new residue */
273 if (i == ntop)
274 {
275 srenew(vsitetoplist, ntop+1)(vsitetoplist) = save_realloc("vsitetoplist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 275, (vsitetoplist), (ntop+1), sizeof(*(vsitetoplist)));
276 ntop++; /* i still points to current vsite topology entry */
277 strncpy(vsitetoplist[i].resname, dirstr, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].resname, dirstr, 32 -1);
278 vsitetoplist[i].nbonds = vsitetoplist[i].nangles = 0;
279 snew(vsitetoplist[i].bond, 1)(vsitetoplist[i].bond) = save_calloc("vsitetoplist[i].bond", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 279, (1), sizeof(*(vsitetoplist[i].bond)));
280 snew(vsitetoplist[i].angle, 1)(vsitetoplist[i].angle) = save_calloc("vsitetoplist[i].angle"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 280, (1), sizeof(*(vsitetoplist[i].angle)));
281 }
282 n = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
283 if (n == 3)
284 {
285 /* bond */
286 k = vsitetoplist[i].nbonds++;
287 srenew(vsitetoplist[i].bond, k+1)(vsitetoplist[i].bond) = save_realloc("vsitetoplist[i].bond",
"/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 287, (vsitetoplist[i].bond), (k+1), sizeof(*(vsitetoplist[i
].bond)));
288 strncpy(vsitetoplist[i].bond[k].atom1, s1, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].bond[k].atom1, s1, 32 -1);
289 strncpy(vsitetoplist[i].bond[k].atom2, s2, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].bond[k].atom2, s2, 32 -1);
290 vsitetoplist[i].bond[k].value = strtod(s3, NULL((void*)0));
291 }
292 else if (n == 4)
293 {
294 /* angle */
295 k = vsitetoplist[i].nangles++;
296 srenew(vsitetoplist[i].angle, k+1)(vsitetoplist[i].angle) = save_realloc("vsitetoplist[i].angle"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 296, (vsitetoplist[i].angle), (k+1), sizeof(*(vsitetoplist[
i].angle)));
297 strncpy(vsitetoplist[i].angle[k].atom1, s1, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].angle[k].atom1, s1, 32 -1);
298 strncpy(vsitetoplist[i].angle[k].atom2, s2, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].angle[k].atom2, s2, 32 -1);
299 strncpy(vsitetoplist[i].angle[k].atom3, s3, MAXNAME-1)__builtin_strncpy (vsitetoplist[i].angle[k].atom3, s3, 32 -1);
300 vsitetoplist[i].angle[k].value = strtod(s4, NULL((void*)0));
301 }
302 else
303 {
304 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 304, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
305 }
306 break;
307 default:
308 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 308, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
309 }
310 }
311 }
312 }
313
314 *pvsiteconflist = vsiteconflist;
315 *pvsitetoplist = vsitetoplist;
316 *nvsiteconf = nvsite;
317 *nvsitetop = ntop;
318
319 gmx_ffclose(ddb);
320}
321
322static int nitrogen_is_planar(t_vsiteconf vsiteconflist[], int nvsiteconf, char atomtype[])
323{
324 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
325 * and -1 if not found.
326 */
327 int i, res;
328 gmx_bool found = FALSE0;
329 for (i = 0; i < nvsiteconf && !found; i++)
330 {
331 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atomtype) && (vsiteconflist[i].nhydrogens == 2));
332 }
333 if (found)
334 {
335 res = (vsiteconflist[i-1].isplanar == TRUE1);
336 }
337 else
338 {
339 res = -1;
340 }
341
342 return res;
343}
344
345static char *get_dummymass_name(t_vsiteconf vsiteconflist[], int nvsiteconf, char atom[], char nextheavy[])
346{
347 /* Return the dummy mass name if found, or NULL if not set in ddb database */
348 int i;
349 gmx_bool found = FALSE0;
350 for (i = 0; i < nvsiteconf && !found; i++)
351 {
352 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atom) &&
353 !gmx_strcasecmp(vsiteconflist[i].nextheavytype, nextheavy));
354 }
355 if (found)
356 {
357 return vsiteconflist[i-1].dummymass;
358 }
359 else
360 {
361 return NULL((void*)0);
362 }
363}
364
365
366
367static real get_ddb_bond(t_vsitetop *vsitetop, int nvsitetop,
368 const char res[],
369 const char atom1[], const char atom2[])
370{
371 int i, j;
372
373 i = 0;
374 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
375 {
376 i++;
377 }
378 if (i == nvsitetop)
379 {
380 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 380, "No vsite information for residue %s found in vsite database.\n", res);
381 }
382 j = 0;
383 while (j < vsitetop[i].nbonds &&
384 ( strcmp(atom1, vsitetop[i].bond[j].atom1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom1) && __builtin_constant_p (vsitetop[i].bond[j]
.atom1) && (__s1_len = strlen (atom1), __s2_len = strlen
(vsitetop[i].bond[j].atom1), (!((size_t)(const void *)((atom1
) + 1) - (size_t)(const void *)(atom1) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].bond[j]
.atom1) + 1) - (size_t)(const void *)(vsitetop[i].bond[j].atom1
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom1, vsitetop
[i].bond[j].atom1) : (__builtin_constant_p (atom1) &&
((size_t)(const void *)((atom1) + 1) - (size_t)(const void *
)(atom1) == 1) && (__s1_len = strlen (atom1), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].bond[j].atom1) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom1) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom1) == 1) ? __builtin_strcmp
(atom1, vsitetop[i].bond[j].atom1) : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(vsitetop[i].bond[j].atom1); int __result = (((const unsigned
char *) (const char *) (atom1))[0] - __s2[0]); if (__s1_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (atom1))[1] - __s2[1]); if (__s1_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (atom1))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (atom1))[3] - __s2[3]); } } __result; }))) :
(__builtin_constant_p (vsitetop[i].bond[j].atom1) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom1) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom1) == 1) &&
(__s2_len = strlen (vsitetop[i].bond[j].atom1), __s2_len <
4) ? (__builtin_constant_p (atom1) && ((size_t)(const
void *)((atom1) + 1) - (size_t)(const void *)(atom1) == 1) ?
__builtin_strcmp (atom1, vsitetop[i].bond[j].atom1) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (atom1); int __result = (((const unsigned char *) (const
char *) (vsitetop[i].bond[j].atom1))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (vsitetop[i].bond[j].atom1))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (vsitetop[i].bond
[j].atom1))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (vsitetop
[i].bond[j].atom1))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom1, vsitetop[i].bond[j].atom1)))); }) || strcmp(atom2, vsitetop[i].bond[j].atom2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom2) && __builtin_constant_p (vsitetop[i].bond[j]
.atom2) && (__s1_len = strlen (atom2), __s2_len = strlen
(vsitetop[i].bond[j].atom2), (!((size_t)(const void *)((atom2
) + 1) - (size_t)(const void *)(atom2) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].bond[j]
.atom2) + 1) - (size_t)(const void *)(vsitetop[i].bond[j].atom2
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom2, vsitetop
[i].bond[j].atom2) : (__builtin_constant_p (atom2) &&
((size_t)(const void *)((atom2) + 1) - (size_t)(const void *
)(atom2) == 1) && (__s1_len = strlen (atom2), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].bond[j].atom2) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom2) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom2) == 1) ? __builtin_strcmp
(atom2, vsitetop[i].bond[j].atom2) : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(vsitetop[i].bond[j].atom2); int __result = (((const unsigned
char *) (const char *) (atom2))[0] - __s2[0]); if (__s1_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (atom2))[1] - __s2[1]); if (__s1_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (atom2))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (atom2))[3] - __s2[3]); } } __result; }))) :
(__builtin_constant_p (vsitetop[i].bond[j].atom2) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom2) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom2) == 1) &&
(__s2_len = strlen (vsitetop[i].bond[j].atom2), __s2_len <
4) ? (__builtin_constant_p (atom2) && ((size_t)(const
void *)((atom2) + 1) - (size_t)(const void *)(atom2) == 1) ?
__builtin_strcmp (atom2, vsitetop[i].bond[j].atom2) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (atom2); int __result = (((const unsigned char *) (const
char *) (vsitetop[i].bond[j].atom2))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (vsitetop[i].bond[j].atom2))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (vsitetop[i].bond
[j].atom2))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (vsitetop
[i].bond[j].atom2))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom2, vsitetop[i].bond[j].atom2)))); })) &&
385 ( strcmp(atom2, vsitetop[i].bond[j].atom1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom2) && __builtin_constant_p (vsitetop[i].bond[j]
.atom1) && (__s1_len = strlen (atom2), __s2_len = strlen
(vsitetop[i].bond[j].atom1), (!((size_t)(const void *)((atom2
) + 1) - (size_t)(const void *)(atom2) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].bond[j]
.atom1) + 1) - (size_t)(const void *)(vsitetop[i].bond[j].atom1
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom2, vsitetop
[i].bond[j].atom1) : (__builtin_constant_p (atom2) &&
((size_t)(const void *)((atom2) + 1) - (size_t)(const void *
)(atom2) == 1) && (__s1_len = strlen (atom2), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].bond[j].atom1) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom1) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom1) == 1) ? __builtin_strcmp
(atom2, vsitetop[i].bond[j].atom1) : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(vsitetop[i].bond[j].atom1); int __result = (((const unsigned
char *) (const char *) (atom2))[0] - __s2[0]); if (__s1_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (atom2))[1] - __s2[1]); if (__s1_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (atom2))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (atom2))[3] - __s2[3]); } } __result; }))) :
(__builtin_constant_p (vsitetop[i].bond[j].atom1) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom1) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom1) == 1) &&
(__s2_len = strlen (vsitetop[i].bond[j].atom1), __s2_len <
4) ? (__builtin_constant_p (atom2) && ((size_t)(const
void *)((atom2) + 1) - (size_t)(const void *)(atom2) == 1) ?
__builtin_strcmp (atom2, vsitetop[i].bond[j].atom1) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (atom2); int __result = (((const unsigned char *) (const
char *) (vsitetop[i].bond[j].atom1))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (vsitetop[i].bond[j].atom1))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (vsitetop[i].bond
[j].atom1))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (vsitetop
[i].bond[j].atom1))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom2, vsitetop[i].bond[j].atom1)))); }) || strcmp(atom1, vsitetop[i].bond[j].atom2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom1) && __builtin_constant_p (vsitetop[i].bond[j]
.atom2) && (__s1_len = strlen (atom1), __s2_len = strlen
(vsitetop[i].bond[j].atom2), (!((size_t)(const void *)((atom1
) + 1) - (size_t)(const void *)(atom1) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].bond[j]
.atom2) + 1) - (size_t)(const void *)(vsitetop[i].bond[j].atom2
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom1, vsitetop
[i].bond[j].atom2) : (__builtin_constant_p (atom1) &&
((size_t)(const void *)((atom1) + 1) - (size_t)(const void *
)(atom1) == 1) && (__s1_len = strlen (atom1), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].bond[j].atom2) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom2) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom2) == 1) ? __builtin_strcmp
(atom1, vsitetop[i].bond[j].atom2) : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(vsitetop[i].bond[j].atom2); int __result = (((const unsigned
char *) (const char *) (atom1))[0] - __s2[0]); if (__s1_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (atom1))[1] - __s2[1]); if (__s1_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (atom1))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (atom1))[3] - __s2[3]); } } __result; }))) :
(__builtin_constant_p (vsitetop[i].bond[j].atom2) &&
((size_t)(const void *)((vsitetop[i].bond[j].atom2) + 1) - (
size_t)(const void *)(vsitetop[i].bond[j].atom2) == 1) &&
(__s2_len = strlen (vsitetop[i].bond[j].atom2), __s2_len <
4) ? (__builtin_constant_p (atom1) && ((size_t)(const
void *)((atom1) + 1) - (size_t)(const void *)(atom1) == 1) ?
__builtin_strcmp (atom1, vsitetop[i].bond[j].atom2) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (atom1); int __result = (((const unsigned char *) (const
char *) (vsitetop[i].bond[j].atom2))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (vsitetop[i].bond[j].atom2))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (vsitetop[i].bond
[j].atom2))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (vsitetop
[i].bond[j].atom2))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom1, vsitetop[i].bond[j].atom2)))); })))
386 {
387 j++;
388 }
389 if (j == vsitetop[i].nbonds)
390 {
391 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 391, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1, atom2, res);
392 }
393
394 return vsitetop[i].bond[j].value;
395}
396
397
398static real get_ddb_angle(t_vsitetop *vsitetop, int nvsitetop,
399 const char res[], const char atom1[],
400 const char atom2[], const char atom3[])
401{
402 int i, j;
403
404 i = 0;
405 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
406 {
407 i++;
408 }
409 if (i == nvsitetop)
410 {
411 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 411, "No vsite information for residue %s found in vsite database.\n", res);
412 }
413 j = 0;
414 while (j < vsitetop[i].nangles &&
415 ( strcmp(atom1, vsitetop[i].angle[j].atom1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom1) && __builtin_constant_p (vsitetop[i].angle[j
].atom1) && (__s1_len = strlen (atom1), __s2_len = strlen
(vsitetop[i].angle[j].atom1), (!((size_t)(const void *)((atom1
) + 1) - (size_t)(const void *)(atom1) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom1) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom1
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom1, vsitetop
[i].angle[j].atom1) : (__builtin_constant_p (atom1) &&
((size_t)(const void *)((atom1) + 1) - (size_t)(const void *
)(atom1) == 1) && (__s1_len = strlen (atom1), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom1)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom1
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom1) ==
1) ? __builtin_strcmp (atom1, vsitetop[i].angle[j].atom1) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom1); int __result
= (((const unsigned char *) (const char *) (atom1))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom1))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom1))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom1))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom1) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom1) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom1) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom1), __s2_len < 4) ? (__builtin_constant_p
(atom1) && ((size_t)(const void *)((atom1) + 1) - (size_t
)(const void *)(atom1) == 1) ? __builtin_strcmp (atom1, vsitetop
[i].angle[j].atom1) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom1); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom1))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom1))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom1))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom1))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom1, vsitetop[i].angle[j].atom1)))); }) ||
416 strcmp(atom2, vsitetop[i].angle[j].atom2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom2) && __builtin_constant_p (vsitetop[i].angle[j
].atom2) && (__s1_len = strlen (atom2), __s2_len = strlen
(vsitetop[i].angle[j].atom2), (!((size_t)(const void *)((atom2
) + 1) - (size_t)(const void *)(atom2) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom2) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom2
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom2, vsitetop
[i].angle[j].atom2) : (__builtin_constant_p (atom2) &&
((size_t)(const void *)((atom2) + 1) - (size_t)(const void *
)(atom2) == 1) && (__s1_len = strlen (atom2), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom2)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom2
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom2) ==
1) ? __builtin_strcmp (atom2, vsitetop[i].angle[j].atom2) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom2); int __result
= (((const unsigned char *) (const char *) (atom2))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom2))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom2))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom2))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom2) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom2) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom2) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom2), __s2_len < 4) ? (__builtin_constant_p
(atom2) && ((size_t)(const void *)((atom2) + 1) - (size_t
)(const void *)(atom2) == 1) ? __builtin_strcmp (atom2, vsitetop
[i].angle[j].atom2) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom2); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom2))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom2))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom2))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom2))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom2, vsitetop[i].angle[j].atom2)))); }) ||
417 strcmp(atom3, vsitetop[i].angle[j].atom3)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom3) && __builtin_constant_p (vsitetop[i].angle[j
].atom3) && (__s1_len = strlen (atom3), __s2_len = strlen
(vsitetop[i].angle[j].atom3), (!((size_t)(const void *)((atom3
) + 1) - (size_t)(const void *)(atom3) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom3) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom3
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom3, vsitetop
[i].angle[j].atom3) : (__builtin_constant_p (atom3) &&
((size_t)(const void *)((atom3) + 1) - (size_t)(const void *
)(atom3) == 1) && (__s1_len = strlen (atom3), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom3)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom3
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom3) ==
1) ? __builtin_strcmp (atom3, vsitetop[i].angle[j].atom3) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom3); int __result
= (((const unsigned char *) (const char *) (atom3))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom3))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom3))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom3))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom3) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom3) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom3) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom3), __s2_len < 4) ? (__builtin_constant_p
(atom3) && ((size_t)(const void *)((atom3) + 1) - (size_t
)(const void *)(atom3) == 1) ? __builtin_strcmp (atom3, vsitetop
[i].angle[j].atom3) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom3); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom3))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom3))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom3))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom3))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom3, vsitetop[i].angle[j].atom3)))); })) &&
418 ( strcmp(atom3, vsitetop[i].angle[j].atom1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom3) && __builtin_constant_p (vsitetop[i].angle[j
].atom1) && (__s1_len = strlen (atom3), __s2_len = strlen
(vsitetop[i].angle[j].atom1), (!((size_t)(const void *)((atom3
) + 1) - (size_t)(const void *)(atom3) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom1) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom1
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom3, vsitetop
[i].angle[j].atom1) : (__builtin_constant_p (atom3) &&
((size_t)(const void *)((atom3) + 1) - (size_t)(const void *
)(atom3) == 1) && (__s1_len = strlen (atom3), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom1)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom1
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom1) ==
1) ? __builtin_strcmp (atom3, vsitetop[i].angle[j].atom1) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom1); int __result
= (((const unsigned char *) (const char *) (atom3))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom3))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom3))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom3))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom1) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom1) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom1) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom1), __s2_len < 4) ? (__builtin_constant_p
(atom3) && ((size_t)(const void *)((atom3) + 1) - (size_t
)(const void *)(atom3) == 1) ? __builtin_strcmp (atom3, vsitetop
[i].angle[j].atom1) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom3); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom1))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom1))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom1))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom1))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom3, vsitetop[i].angle[j].atom1)))); }) ||
419 strcmp(atom2, vsitetop[i].angle[j].atom2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom2) && __builtin_constant_p (vsitetop[i].angle[j
].atom2) && (__s1_len = strlen (atom2), __s2_len = strlen
(vsitetop[i].angle[j].atom2), (!((size_t)(const void *)((atom2
) + 1) - (size_t)(const void *)(atom2) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom2) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom2
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom2, vsitetop
[i].angle[j].atom2) : (__builtin_constant_p (atom2) &&
((size_t)(const void *)((atom2) + 1) - (size_t)(const void *
)(atom2) == 1) && (__s1_len = strlen (atom2), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom2)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom2
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom2) ==
1) ? __builtin_strcmp (atom2, vsitetop[i].angle[j].atom2) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom2); int __result
= (((const unsigned char *) (const char *) (atom2))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom2))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom2))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom2))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom2) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom2) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom2) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom2), __s2_len < 4) ? (__builtin_constant_p
(atom2) && ((size_t)(const void *)((atom2) + 1) - (size_t
)(const void *)(atom2) == 1) ? __builtin_strcmp (atom2, vsitetop
[i].angle[j].atom2) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom2); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom2))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom2))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom2))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom2))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom2, vsitetop[i].angle[j].atom2)))); }) ||
420 strcmp(atom1, vsitetop[i].angle[j].atom3)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(atom1) && __builtin_constant_p (vsitetop[i].angle[j
].atom3) && (__s1_len = strlen (atom1), __s2_len = strlen
(vsitetop[i].angle[j].atom3), (!((size_t)(const void *)((atom1
) + 1) - (size_t)(const void *)(atom1) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((vsitetop[i].angle[j
].atom3) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom3
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (atom1, vsitetop
[i].angle[j].atom3) : (__builtin_constant_p (atom1) &&
((size_t)(const void *)((atom1) + 1) - (size_t)(const void *
)(atom1) == 1) && (__s1_len = strlen (atom1), __s1_len
< 4) ? (__builtin_constant_p (vsitetop[i].angle[j].atom3)
&& ((size_t)(const void *)((vsitetop[i].angle[j].atom3
) + 1) - (size_t)(const void *)(vsitetop[i].angle[j].atom3) ==
1) ? __builtin_strcmp (atom1, vsitetop[i].angle[j].atom3) : (
__extension__ ({ const unsigned char *__s2 = (const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom3); int __result
= (((const unsigned char *) (const char *) (atom1))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom1))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (atom1))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (atom1))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (vsitetop[i].
angle[j].atom3) && ((size_t)(const void *)((vsitetop[
i].angle[j].atom3) + 1) - (size_t)(const void *)(vsitetop[i].
angle[j].atom3) == 1) && (__s2_len = strlen (vsitetop
[i].angle[j].atom3), __s2_len < 4) ? (__builtin_constant_p
(atom1) && ((size_t)(const void *)((atom1) + 1) - (size_t
)(const void *)(atom1) == 1) ? __builtin_strcmp (atom1, vsitetop
[i].angle[j].atom3) : (- (__extension__ ({ const unsigned char
*__s2 = (const unsigned char *) (const char *) (atom1); int __result
= (((const unsigned char *) (const char *) (vsitetop[i].angle
[j].atom3))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
vsitetop[i].angle[j].atom3))[1] - __s2[1]); if (__s2_len >
1 && __result == 0) { __result = (((const unsigned char
*) (const char *) (vsitetop[i].angle[j].atom3))[2] - __s2[2]
); if (__s2_len > 2 && __result == 0) __result = (
((const unsigned char *) (const char *) (vsitetop[i].angle[j]
.atom3))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(atom1, vsitetop[i].angle[j].atom3)))); })))
421 {
422 j++;
423 }
424 if (j == vsitetop[i].nangles)
425 {
426 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 426, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1, atom2, atom3, res);
427 }
428
429 return vsitetop[i].angle[j].value;
430}
431
432
433static void count_bonds(int atom, t_params *psb, char ***atomname,
434 int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
435 int *nrheavies, int heavies[])
436{
437 int i, heavy, other, nrb, nrH, nrhv;
438
439 /* find heavy atom bound to this hydrogen */
440 heavy = NOTSET-12345;
441 for (i = 0; (i < psb->nr) && (heavy == NOTSET-12345); i++)
442 {
443 if (psb->param[i].AIa[0] == atom)
444 {
445 heavy = psb->param[i].AJa[1];
446 }
447 else if (psb->param[i].AJa[1] == atom)
448 {
449 heavy = psb->param[i].AIa[0];
450 }
451 }
452 if (heavy == NOTSET-12345)
453 {
454 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 454, "unbound hydrogen atom %d", atom+1);
455 }
456 /* find all atoms bound to heavy atom */
457 other = NOTSET-12345;
458 nrb = 0;
459 nrH = 0;
460 nrhv = 0;
461 for (i = 0; i < psb->nr; i++)
462 {
463 if (psb->param[i].AIa[0] == heavy)
464 {
465 other = psb->param[i].AJa[1];
466 }
467 else if (psb->param[i].AJa[1] == heavy)
468 {
469 other = psb->param[i].AIa[0];
470 }
471 if (other != NOTSET-12345)
472 {
473 nrb++;
474 if (is_hydrogen(*(atomname[other])))
475 {
476 Hatoms[nrH] = other;
477 nrH++;
478 }
479 else
480 {
481 heavies[nrhv] = other;
482 nrhv++;
483 }
484 other = NOTSET-12345;
485 }
486 }
487 *Heavy = heavy;
488 *nrbonds = nrb;
489 *nrHatoms = nrH;
490 *nrheavies = nrhv;
491}
492
493static void print_bonds(FILE *fp, int o2n[],
494 int nrHatoms, int Hatoms[], int Heavy,
495 int nrheavies, int heavies[])
496{
497 int i;
498
499 fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
500 for (i = 0; i < nrHatoms; i++)
501 {
502 fprintf(fp, " %d", o2n[Hatoms[i]]+1);
503 }
504 fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
505 for (i = 0; i < nrheavies; i++)
506 {
507 fprintf(fp, " %d", o2n[heavies[i]]+1);
508 }
509 fprintf(fp, "\n");
510}
511
512static int get_atype(int atom, t_atoms *at, int nrtp, t_restp rtp[],
513 gmx_residuetype_t rt)
514{
515 int type;
516 gmx_bool bNterm;
517 int j;
518 t_restp *rtpp;
519
520 if (at->atom[atom].m)
521 {
522 type = at->atom[atom].type;
523 }
524 else
525 {
526 /* get type from rtp */
527 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
528 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
529 (at->atom[atom].resind == 0);
530 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
531 type = rtpp->atom[j].type;
532 }
533 return type;
534}
535
536static int vsite_nm2type(const char *name, gpp_atomtype_t atype)
537{
538 int tp;
539
540 tp = get_atomtype_type(name, atype);
541 if (tp == NOTSET-12345)
542 {
543 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 543, "Dummy mass type (%s) not found in atom type database",
544 name);
545 }
546
547 return tp;
548}
549
550static real get_amass(int atom, t_atoms *at, int nrtp, t_restp rtp[],
551 gmx_residuetype_t rt)
552{
553 real mass;
554 gmx_bool bNterm;
555 int j;
556 t_restp *rtpp;
557
558 if (at->atom[atom].m)
559 {
560 mass = at->atom[atom].m;
561 }
562 else
563 {
564 /* get mass from rtp */
565 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
566 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
567 (at->atom[atom].resind == 0);
568 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
569 mass = rtpp->atom[j].m;
570 }
571 return mass;
572}
573
574static void my_add_param(t_params *plist, int ai, int aj, real b)
575{
576 static real c[MAXFORCEPARAM12] =
577 { NOTSET-12345, NOTSET-12345, NOTSET-12345, NOTSET-12345, NOTSET-12345, NOTSET-12345 };
578
579 c[0] = b;
580 add_param(plist, ai, aj, c, NULL((void*)0));
581}
582
583static void add_vsites(t_params plist[], int vsite_type[],
584 int Heavy, int nrHatoms, int Hatoms[],
585 int nrheavies, int heavies[])
586{
587 int i, j, ftype, other, moreheavy, bb;
588 gmx_bool bSwapParity;
589
590 for (i = 0; i < nrHatoms; i++)
591 {
592 ftype = vsite_type[Hatoms[i]];
593 /* Errors in setting the vsite_type should really be caugth earlier,
594 * because here it's not possible to print any useful error message.
595 * But it's still better to print a message than to segfault.
596 */
597 if (ftype == NOTSET-12345)
598 {
599 gmx_incons("Undetected error in setting up virtual sites")_gmx_error("incons", "Undetected error in setting up virtual sites"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 599);
600 }
601 bSwapParity = (ftype < 0);
602 vsite_type[Hatoms[i]] = ftype = abs(ftype);
603 if (ftype == F_BONDS)
604 {
605 if ( (nrheavies != 1) && (nrHatoms != 1) )
606 {
607 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 607, "cannot make constraint in add_vsites for %d heavy "
608 "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
609 }
610 my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET-12345);
611 }
612 else
613 {
614 switch (ftype)
615 {
616 case F_VSITE3:
617 case F_VSITE3FD:
618 case F_VSITE3OUT:
619 if (nrheavies < 2)
620 {
621 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 621, "Not enough heavy atoms (%d) for %s (min 3)",
622 nrheavies+1,
623 interaction_function[vsite_type[Hatoms[i]]].name);
624 }
625 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
626 bSwapParity);
627 break;
628 case F_VSITE3FAD:
629 {
630 if (nrheavies > 1)
631 {
632 moreheavy = heavies[1];
633 }
634 else
635 {
636 /* find more heavy atoms */
637 other = moreheavy = NOTSET-12345;
638 for (j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET-12345); j++)
639 {
640 if (plist[F_BONDS].param[j].AIa[0] == heavies[0])
641 {
642 other = plist[F_BONDS].param[j].AJa[1];
643 }
644 else if (plist[F_BONDS].param[j].AJa[1] == heavies[0])
645 {
646 other = plist[F_BONDS].param[j].AIa[0];
647 }
648 if ( (other != NOTSET-12345) && (other != Heavy) )
649 {
650 moreheavy = other;
651 }
652 }
653 if (moreheavy == NOTSET-12345)
654 {
655 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 655, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
656 }
657 }
658 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
659 bSwapParity);
660 break;
661 }
662 case F_VSITE4FD:
663 case F_VSITE4FDN:
664 if (nrheavies < 3)
665 {
666 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 666, "Not enough heavy atoms (%d) for %s (min 4)",
667 nrheavies+1,
668 interaction_function[vsite_type[Hatoms[i]]].name);
669 }
670 add_vsite4_atoms(&plist[ftype],
671 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
672 break;
673
674 default:
675 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 675, "can't use add_vsites for interaction function %s",
676 interaction_function[vsite_type[Hatoms[i]]].name);
677 } /* switch ftype */
678 } /* else */
679 } /* for i */
680}
681
682#define ANGLE_6RING((3.14159265358979323846/180.0)*120) (DEG2RAD(3.14159265358979323846/180.0)*120)
683
684/* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
685/* get a^2 when a, b and alpha are given: */
686#define cosrule(b, c, alpha)( sqr(b) + sqr(c) - 2*b*c*cos(alpha) ) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
687/* get cos(alpha) when a, b and c are given: */
688#define acosrule(a, b, c)( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) ) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
689
690static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], t_params plist[],
691 int nrfound, int *ats, real bond_cc, real bond_ch,
692 real xcom, gmx_bool bDoZ)
693{
694 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
695 enum {
696 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
697 atCZ, atHZ, atNR
698 };
699
700 int i, nvsite;
701 real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
702 real xCG, yCG, xCE1, yCE1, xCE2, yCE2;
703 /* CG, CE1 and CE2 stay and each get a part of the total mass,
704 * so the c-o-m stays the same.
705 */
706
707 if (bDoZ)
708 {
709 if (atNR != nrfound)
710 {
711 gmx_incons("Generating vsites on 6-rings")_gmx_error("incons", "Generating vsites on 6-rings", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 711);
712 }
713 }
714
715 /* constraints between CG, CE1 and CE2: */
716 dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING)( sqr(bond_cc) + sqr(bond_cc) - 2*bond_cc*bond_cc*cos(((3.14159265358979323846
/180.0)*120)) ) );
717 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
718 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
719 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
720
721 /* rest will be vsite3 */
722 mtot = 0;
723 nvsite = 0;
724 for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
725 {
726 mtot += at->atom[ats[i]].m;
727 if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
728 {
729 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
730 (*vsite_type)[ats[i]] = F_VSITE3;
731 nvsite++;
732 }
733 }
734 /* Distribute mass so center-of-mass stays the same.
735 * The center-of-mass in the call is defined with x=0 at
736 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
737 */
738 xCG = -bond_cc+bond_cc*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
739 yCG = 0;
740 xCE1 = 0;
Value stored to 'xCE1' is never read
741 yCE1 = bond_cc*sin(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
742 xCE2 = 0;
743 yCE2 = -bond_cc*sin(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
744
745 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
746 mrest = mtot-mG;
747 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
748 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
749
750 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
751 tmp1 = dCGCE*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120)*0.5);
752 tmp2 = bond_cc*cos(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120)) + tmp1;
753 tmp1 *= 2;
754 a = b = -bond_ch / tmp1;
755 /* HE1 and HE2: */
756 add_vsite3_param(&plist[F_VSITE3],
757 ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
758 add_vsite3_param(&plist[F_VSITE3],
759 ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
760 /* CD1, CD2 and CZ: */
761 a = b = tmp2 / tmp1;
762 add_vsite3_param(&plist[F_VSITE3],
763 ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
764 add_vsite3_param(&plist[F_VSITE3],
765 ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
766 if (bDoZ)
767 {
768 add_vsite3_param(&plist[F_VSITE3],
769 ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
770 }
771 /* HD1, HD2 and HZ: */
772 a = b = ( bond_ch + tmp2 ) / tmp1;
773 add_vsite3_param(&plist[F_VSITE3],
774 ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
775 add_vsite3_param(&plist[F_VSITE3],
776 ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
777 if (bDoZ)
778 {
779 add_vsite3_param(&plist[F_VSITE3],
780 ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
781 }
782
783 return nvsite;
784}
785
786static int gen_vsites_phe(t_atoms *at, int *vsite_type[], t_params plist[],
787 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
788{
789 real bond_cc, bond_ch;
790 real xcom, mtot;
791 int i;
792 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
793 enum {
794 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
795 atCZ, atHZ, atNR
796 };
797 real x[atNR], y[atNR];
798 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
799 * (angle is always 120 degrees).
800 */
801 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "CE1");
802 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "HD1");
803
804 x[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
805 y[atCG] = 0;
806 x[atCD1] = -bond_cc;
807 y[atCD1] = bond_cc*sin(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
808 x[atHD1] = x[atCD1]+bond_ch*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
809 y[atHD1] = y[atCD1]+bond_ch*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120));
810 x[atCE1] = 0;
811 y[atCE1] = y[atCD1];
812 x[atHE1] = x[atCE1]-bond_ch*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
813 y[atHE1] = y[atCE1]+bond_ch*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120));
814 x[atCD2] = x[atCD1];
815 y[atCD2] = -y[atCD1];
816 x[atHD2] = x[atHD1];
817 y[atHD2] = -y[atHD1];
818 x[atCE2] = x[atCE1];
819 y[atCE2] = -y[atCE1];
820 x[atHE2] = x[atHE1];
821 y[atHE2] = -y[atHE1];
822 x[atCZ] = bond_cc*cos(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
823 y[atCZ] = 0;
824 x[atHZ] = x[atCZ]+bond_ch;
825 y[atHZ] = 0;
826
827 xcom = mtot = 0;
828 for (i = 0; i < atNR; i++)
829 {
830 xcom += x[i]*at->atom[ats[i]].m;
831 mtot += at->atom[ats[i]].m;
832 }
833 xcom /= mtot;
834
835 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE1);
836}
837
838static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
839 real xk, real yk, real *a, real *b)
840{
841 /* determine parameters by solving the equation system, since we know the
842 * virtual site coordinates here.
843 */
844 real dx_ij, dx_ik, dy_ij, dy_ik;
845 real b_ij, b_ik;
846
847 dx_ij = xj-xi;
848 dy_ij = yj-yi;
849 dx_ik = xk-xi;
850 dy_ik = yk-yi;
851 b_ij = sqrt(dx_ij*dx_ij+dy_ij*dy_ij);
852 b_ik = sqrt(dx_ik*dx_ik+dy_ik*dy_ik);
853
854 *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
855 *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
856}
857
858
859static int gen_vsites_trp(gpp_atomtype_t atype, rvec *newx[],
860 t_atom *newatom[], char ***newatomname[],
861 int *o2n[], int *newvsite_type[], int *newcgnr[],
862 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
863 t_atoms *at, int *vsite_type[], t_params plist[],
864 int nrfound, int *ats, int add_shift,
865 t_vsitetop *vsitetop, int nvsitetop)
866{
867#define NMASS 2
868 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
869 enum {
870 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
871 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
872 };
873 /* weights for determining the COM's of both rings (M1 and M2): */
874 real mw[NMASS][atNR] = {
875 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
876 0, 0, 0, 0, 0, 0 },
877 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
878 1, 1, 1, 1, 1, 1 }
879 };
880
881 real xi[atNR], yi[atNR];
882 real xcom[NMASS], ycom[NMASS], I, alpha;
883 real lineA, lineB, dist;
884 real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
885 real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
886 real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
887 real b_CG_CD1, b_CZ3_HZ3;
888 real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
889 real a_CB_CG_CD1, a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
890 real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
891 real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
892 real xM[NMASS];
893 int atM[NMASS], tpM, i, i0, j, nvsite;
894 real mwtot, mtot, mM[NMASS], dCBM1, dCBM2, dM1M2;
895 real a, b, c[MAXFORCEPARAM12];
896 rvec r_ij, r_ik, t1, t2;
897 char name[10];
898
899 if (atNR != nrfound)
900 {
901 gmx_incons("atom types in gen_vsites_trp")_gmx_error("incons", "atom types in gen_vsites_trp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 901);
902 }
903 /* Get geometry from database */
904 b_CD2_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE2");
905 b_NE1_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "CE2");
906 b_CG_CD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD1");
907 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD2");
908 b_CB_CG = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CB", "CG");
909 b_CE2_CZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE2", "CZ2");
910 b_CD2_CE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE3");
911 b_CE3_CZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "CZ3");
912 b_CZ2_CH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "CH2");
913
914 b_CD1_HD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD1", "HD1");
915 b_CZ2_HZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "HZ2");
916 b_NE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "HE1");
917 b_CH2_HH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CH2", "HH2");
918 b_CE3_HE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "HE3");
919 b_CZ3_HZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ3", "HZ3");
920
921 a_NE1_CE2_CD2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "NE1", "CE2", "CD2");
922 a_CE2_CD2_CG = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CG");
923 a_CB_CG_CD2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD2");
924 a_CD2_CG_CD1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CG", "CD1");
925 a_CB_CG_CD1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD1");
926
927 a_CE2_CD2_CE3 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CE3");
928 a_CD2_CE2_CZ2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE2", "CZ2");
929 a_CD2_CE3_CZ3 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "CZ3");
930 a_CE3_CZ3_HZ3 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE3", "CZ3", "HZ3");
931 a_CZ2_CH2_HH2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CZ2", "CH2", "HH2");
932 a_CE2_CZ2_HZ2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "HZ2");
933 a_CE2_CZ2_CH2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "CH2");
934 a_CG_CD1_HD1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CG", "CD1", "HD1");
935 a_HE1_NE1_CE2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "HE1", "NE1", "CE2");
936 a_CD2_CE3_HE3 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "HE3");
937
938 /* Calculate local coordinates.
939 * y-axis (x=0) is the bond CD2-CE2.
940 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
941 * intersects the middle of the bond.
942 */
943 xi[atCD2] = 0;
944 yi[atCD2] = -0.5*b_CD2_CE2;
945
946 xi[atCE2] = 0;
947 yi[atCE2] = 0.5*b_CD2_CE2;
948
949 xi[atNE1] = -b_NE1_CE2*sin(a_NE1_CE2_CD2);
950 yi[atNE1] = yi[atCE2]-b_NE1_CE2*cos(a_NE1_CE2_CD2);
951
952 xi[atCG] = -b_CG_CD2*sin(a_CE2_CD2_CG);
953 yi[atCG] = yi[atCD2]+b_CG_CD2*cos(a_CE2_CD2_CG);
954
955 alpha = a_CE2_CD2_CG + M_PI3.14159265358979323846 - a_CB_CG_CD2;
956 xi[atCB] = xi[atCG]-b_CB_CG*sin(alpha);
957 yi[atCB] = yi[atCG]+b_CB_CG*cos(alpha);
958
959 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI3.14159265358979323846;
960 xi[atCD1] = xi[atCG]-b_CG_CD1*sin(alpha);
961 yi[atCD1] = yi[atCG]+b_CG_CD1*cos(alpha);
962
963 xi[atCE3] = b_CD2_CE3*sin(a_CE2_CD2_CE3);
964 yi[atCE3] = yi[atCD2]+b_CD2_CE3*cos(a_CE2_CD2_CE3);
965
966 xi[atCZ2] = b_CE2_CZ2*sin(a_CD2_CE2_CZ2);
967 yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*cos(a_CD2_CE2_CZ2);
968
969 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI3.14159265358979323846;
970 xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*sin(alpha);
971 yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*cos(alpha);
972
973 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI3.14159265358979323846;
974 xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*sin(alpha);
975 yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*cos(alpha);
976
977 /* hydrogens */
978 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
979 xi[atHD1] = xi[atCD1]-b_CD1_HD1*sin(alpha);
980 yi[atHD1] = yi[atCD1]+b_CD1_HD1*cos(alpha);
981
982 alpha = a_NE1_CE2_CD2 + M_PI3.14159265358979323846 - a_HE1_NE1_CE2;
983 xi[atHE1] = xi[atNE1]-b_NE1_HE1*sin(alpha);
984 yi[atHE1] = yi[atNE1]-b_NE1_HE1*cos(alpha);
985
986 alpha = a_CE2_CD2_CE3 + M_PI3.14159265358979323846 - a_CD2_CE3_HE3;
987 xi[atHE3] = xi[atCE3]+b_CE3_HE3*sin(alpha);
988 yi[atHE3] = yi[atCE3]+b_CE3_HE3*cos(alpha);
989
990 alpha = a_CD2_CE2_CZ2 + M_PI3.14159265358979323846 - a_CE2_CZ2_HZ2;
991 xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*sin(alpha);
992 yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*cos(alpha);
993
994 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
995 xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*sin(alpha);
996 yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*cos(alpha);
997
998 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
999 xi[atHH2] = xi[atCH2]+b_CH2_HH2*sin(alpha);
1000 yi[atHH2] = yi[atCH2]-b_CH2_HH2*cos(alpha);
1001
1002 /* Determine coeff. for the line CB-CG */
1003 lineA = (yi[atCB]-yi[atCG])/(xi[atCB]-xi[atCG]);
1004 lineB = yi[atCG]-lineA*xi[atCG];
1005
1006 /* Calculate masses for each ring and put it on the dummy masses */
1007 for (j = 0; j < NMASS; j++)
1008 {
1009 mM[j] = xcom[j] = ycom[j] = 0;
1010 }
1011 for (i = 0; i < atNR; i++)
1012 {
1013 if (i != atCB)
1014 {
1015 for (j = 0; j < NMASS; j++)
1016 {
1017 mM[j] += mw[j][i] * at->atom[ats[i]].m;
1018 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
1019 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
1020 }
1021 }
1022 }
1023 for (j = 0; j < NMASS; j++)
1024 {
1025 xcom[j] /= mM[j];
1026 ycom[j] /= mM[j];
1027 }
1028
1029 /* get dummy mass type */
1030 tpM = vsite_nm2type("MW", atype);
1031 /* make space for 2 masses: shift all atoms starting with CB */
1032 i0 = ats[atCB];
1033 for (j = 0; j < NMASS; j++)
1034 {
1035 atM[j] = i0+*nadd+j;
1036 }
1037 if (debug)
1038 {
1039 fprintf(stderrstderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
1040 }
1041 *nadd += NMASS;
1042 for (j = i0; j < at->nr; j++)
1043 {
1044 (*o2n)[j] = j+*nadd;
1045 }
1046 srenew(*newx, at->nr+*nadd)(*newx) = save_realloc("*newx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1046, (*newx), (at->nr+*nadd), sizeof(*(*newx)));
1047 srenew(*newatom, at->nr+*nadd)(*newatom) = save_realloc("*newatom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1047, (*newatom), (at->nr+*nadd), sizeof(*(*newatom)));
1048 srenew(*newatomname, at->nr+*nadd)(*newatomname) = save_realloc("*newatomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1048, (*newatomname), (at->nr+*nadd), sizeof(*(*newatomname
)));
1049 srenew(*newvsite_type, at->nr+*nadd)(*newvsite_type) = save_realloc("*newvsite_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1049, (*newvsite_type), (at->nr+*nadd), sizeof(*(*newvsite_type
)));
1050 srenew(*newcgnr, at->nr+*nadd)(*newcgnr) = save_realloc("*newcgnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1050, (*newcgnr), (at->nr+*nadd), sizeof(*(*newcgnr)));
1051 for (j = 0; j < NMASS; j++)
1052 {
1053 (*newatomname)[at->nr+*nadd-1-j] = NULL((void*)0);
1054 }
1055
1056 /* Dummy masses will be placed at the center-of-mass in each ring. */
1057
1058 /* calc initial position for dummy masses in real (non-local) coordinates.
1059 * Cheat by using the routine to calculate virtual site parameters. It is
1060 * much easier when we have the coordinates expressed in terms of
1061 * CB, CG, CD2.
1062 */
1063 rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1064 rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1065 calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1066 xi[atCD2], yi[atCD2], &a, &b);
1067 svmul(a, r_ij, t1);
1068 svmul(b, r_ik, t2);
1069 rvec_add(t1, t2, t1);
1070 rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1071
1072 calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1073 xi[atCD2], yi[atCD2], &a, &b);
1074 svmul(a, r_ij, t1);
1075 svmul(b, r_ik, t2);
1076 rvec_add(t1, t2, t1);
1077 rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1078
1079 /* set parameters for the masses */
1080 for (j = 0; j < NMASS; j++)
1081 {
1082 sprintf(name, "MW%d", j+1);
1083 (*newatomname) [atM[j]] = put_symtab(symtab, name);
1084 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
1085 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
1086 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
1087 (*newatom) [atM[j]].ptype = eptAtom;
1088 (*newatom) [atM[j]].resind = at->atom[i0].resind;
1089 (*newatom) [atM[j]].elem[0] = 'M';
1090 (*newatom) [atM[j]].elem[1] = '\0';
1091 (*newvsite_type)[atM[j]] = NOTSET-12345;
1092 (*newcgnr) [atM[j]] = (*cgnr)[i0];
1093 }
1094 /* renumber cgnr: */
1095 for (i = i0; i < at->nr; i++)
1096 {
1097 (*cgnr)[i]++;
1098 }
1099
1100 /* constraints between CB, M1 and M2 */
1101 /* 'add_shift' says which atoms won't be renumbered afterwards */
1102 dCBM1 = sqrt( sqr(xcom[0]-xi[atCB]) + sqr(ycom[0]-yi[atCB]) );
1103 dM1M2 = sqrt( sqr(xcom[0]-xcom[1]) + sqr(ycom[0]-ycom[1]) );
1104 dCBM2 = sqrt( sqr(xcom[1]-xi[atCB]) + sqr(ycom[1]-yi[atCB]) );
1105 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[0], dCBM1);
1106 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[1], dCBM2);
1107 my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
1108
1109 /* rest will be vsite3 */
1110 nvsite = 0;
1111 for (i = 0; i < atNR; i++)
1112 {
1113 if (i != atCB)
1114 {
1115 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1116 (*vsite_type)[ats[i]] = F_VSITE3;
1117 nvsite++;
1118 }
1119 }
1120
1121 /* now define all vsites from M1, M2, CB, ie:
1122 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1123 for (i = 0; i < atNR; i++)
1124 {
1125 if ( (*vsite_type)[ats[i]] == F_VSITE3)
1126 {
1127 calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
1128 add_vsite3_param(&plist[F_VSITE3],
1129 ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
1130 }
1131 }
1132 return nvsite;
1133#undef NMASS
1134}
1135
1136
1137static int gen_vsites_tyr(gpp_atomtype_t atype, rvec *newx[],
1138 t_atom *newatom[], char ***newatomname[],
1139 int *o2n[], int *newvsite_type[], int *newcgnr[],
1140 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
1141 t_atoms *at, int *vsite_type[], t_params plist[],
1142 int nrfound, int *ats, int add_shift,
1143 t_vsitetop *vsitetop, int nvsitetop)
1144{
1145 int nvsite, i, i0, j, atM, tpM;
1146 real dCGCE, dCEOH, dCGM, tmp1, a, b;
1147 real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1148 real xcom, mtot;
1149 real vmass, vdist, mM;
1150 rvec r1;
1151 char name[10];
1152
1153 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1154 enum {
1155 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
1156 atCZ, atOH, atHH, atNR
1157 };
1158 real xi[atNR], yi[atNR];
1159 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1160 rest gets virtualized.
1161 Now we have two linked triangles with one improper keeping them flat */
1162 if (atNR != nrfound)
1163 {
1164 gmx_incons("Number of atom types in gen_vsites_tyr")_gmx_error("incons", "Number of atom types in gen_vsites_tyr"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1164);
1165 }
1166
1167 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1168 * for the ring part (angle is always 120 degrees).
1169 */
1170 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "CE1");
1171 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "HD1");
1172 bond_co = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CZ", "OH");
1173 bond_oh = get_ddb_bond(vsitetop, nvsitetop, "TYR", "OH", "HH");
1174 angle_coh = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, "TYR", "CZ", "OH", "HH");
1175
1176 xi[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
1177 yi[atCG] = 0;
1178 xi[atCD1] = -bond_cc;
1179 yi[atCD1] = bond_cc*sin(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
1180 xi[atHD1] = xi[atCD1]+bond_ch*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
1181 yi[atHD1] = yi[atCD1]+bond_ch*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120));
1182 xi[atCE1] = 0;
1183 yi[atCE1] = yi[atCD1];
1184 xi[atHE1] = xi[atCE1]-bond_ch*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120));
1185 yi[atHE1] = yi[atCE1]+bond_ch*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120));
1186 xi[atCD2] = xi[atCD1];
1187 yi[atCD2] = -yi[atCD1];
1188 xi[atHD2] = xi[atHD1];
1189 yi[atHD2] = -yi[atHD1];
1190 xi[atCE2] = xi[atCE1];
1191 yi[atCE2] = -yi[atCE1];
1192 xi[atHE2] = xi[atHE1];
1193 yi[atHE2] = -yi[atHE1];
1194 xi[atCZ] = bond_cc*cos(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120));
1195 yi[atCZ] = 0;
1196 xi[atOH] = xi[atCZ]+bond_co;
1197 yi[atOH] = 0;
1198
1199 xcom = mtot = 0;
1200 for (i = 0; i < atOH; i++)
1201 {
1202 xcom += xi[i]*at->atom[ats[i]].m;
1203 mtot += at->atom[ats[i]].m;
1204 }
1205 xcom /= mtot;
1206
1207 /* first do 6 ring as default,
1208 except CZ (we'll do that different) and HZ (we don't have that): */
1209 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE0);
1210
1211 /* then construct CZ from the 2nd triangle */
1212 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1213 a = b = 0.5 * bond_co / ( bond_co - bond_cc*cos(ANGLE_6RING((3.14159265358979323846/180.0)*120)) );
1214 add_vsite3_param(&plist[F_VSITE3],
1215 ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1216 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1217
1218 /* constraints between CE1, CE2 and OH */
1219 dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING)( sqr(bond_cc) + sqr(bond_cc) - 2*bond_cc*bond_cc*cos(((3.14159265358979323846
/180.0)*120)) ) );
1220 dCEOH = sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING)( sqr(bond_cc) + sqr(bond_co) - 2*bond_cc*bond_co*cos(((3.14159265358979323846
/180.0)*120)) ) );
1221 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1222 my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1223
1224 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1225 * we need to introduce a constraint to CG.
1226 * CG is much further away, so that will lead to instabilities in LINCS
1227 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1228 * the use of lincs_order=8 we introduce a dummy mass three times further
1229 * away from OH than HH. The mass is accordingly a third, with the remaining
1230 * 2/3 moved to OH. This shouldnt cause any problems since the forces will
1231 * apply to the HH constructed atom and not directly on the virtual mass.
1232 */
1233
1234 vdist = 2.0*bond_oh;
1235 mM = at->atom[ats[atHH]].m/2.0;
1236 at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
1237 at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1238 at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
1239
1240 /* get dummy mass type */
1241 tpM = vsite_nm2type("MW", atype);
1242 /* make space for 1 mass: shift HH only */
1243 i0 = ats[atHH];
1244 atM = i0+*nadd;
1245 if (debug)
1246 {
1247 fprintf(stderrstderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
1248 }
1249 (*nadd)++;
1250 for (j = i0; j < at->nr; j++)
1251 {
1252 (*o2n)[j] = j+*nadd;
1253 }
1254 srenew(*newx, at->nr+*nadd)(*newx) = save_realloc("*newx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1254, (*newx), (at->nr+*nadd), sizeof(*(*newx)));
1255 srenew(*newatom, at->nr+*nadd)(*newatom) = save_realloc("*newatom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1255, (*newatom), (at->nr+*nadd), sizeof(*(*newatom)));
1256 srenew(*newatomname, at->nr+*nadd)(*newatomname) = save_realloc("*newatomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1256, (*newatomname), (at->nr+*nadd), sizeof(*(*newatomname
)));
1257 srenew(*newvsite_type, at->nr+*nadd)(*newvsite_type) = save_realloc("*newvsite_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1257, (*newvsite_type), (at->nr+*nadd), sizeof(*(*newvsite_type
)));
1258 srenew(*newcgnr, at->nr+*nadd)(*newcgnr) = save_realloc("*newcgnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1258, (*newcgnr), (at->nr+*nadd), sizeof(*(*newcgnr)));
1259 (*newatomname)[at->nr+*nadd-1] = NULL((void*)0);
1260
1261 /* Calc the dummy mass initial position */
1262 rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1263 svmul(2.0, r1, r1);
1264 rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1265
1266 strcpy(name, "MW1");
1267 (*newatomname) [atM] = put_symtab(symtab, name);
1268 (*newatom) [atM].m = (*newatom)[atM].mB = mM;
1269 (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
1270 (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
1271 (*newatom) [atM].ptype = eptAtom;
1272 (*newatom) [atM].resind = at->atom[i0].resind;
1273 (*newatom) [atM].elem[0] = 'M';
1274 (*newatom) [atM].elem[1] = '\0';
1275 (*newvsite_type)[atM] = NOTSET-12345;
1276 (*newcgnr) [atM] = (*cgnr)[i0];
1277 /* renumber cgnr: */
1278 for (i = i0; i < at->nr; i++)
1279 {
1280 (*cgnr)[i]++;
1281 }
1282
1283 (*vsite_type)[ats[atHH]] = F_VSITE2;
1284 nvsite++;
1285 /* assume we also want the COH angle constrained: */
1286 tmp1 = bond_cc*cos(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120)) + dCGCE*sin(ANGLE_6RING((3.14159265358979323846/180.0)*120)*0.5) + bond_co;
1287 dCGM = sqrt( cosrule(tmp1, vdist, angle_coh)( sqr(tmp1) + sqr(vdist) - 2*tmp1*vdist*cos(angle_coh) ) );
1288 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
1289 my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
1290
1291 add_vsite2_param(&plist[F_VSITE2],
1292 ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
1293 return nvsite;
1294}
1295
1296static int gen_vsites_his(t_atoms *at, int *vsite_type[], t_params plist[],
1297 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
1298{
1299 int nvsite, i;
1300 real a, b, alpha, dCGCE1, dCGNE2;
1301 real sinalpha, cosalpha;
1302 real xcom, ycom, mtot;
1303 real mG, mrest, mCE1, mNE2;
1304 real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1305 real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1306 real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1307 real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1308 char resname[10];
1309
1310 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1311 enum {
1312 atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
1313 };
1314 real x[atNR], y[atNR];
1315
1316 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1317 rest gets virtualized */
1318 /* check number of atoms, 3 hydrogens may be missing: */
1319 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1320 * Don't understand the above logic. Shouldn't it be && rather than || ???
1321 */
1322 if ((nrfound < atNR-3) || (nrfound > atNR))
1323 {
1324 gmx_incons("Generating vsites for HIS")_gmx_error("incons", "Generating vsites for HIS", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1324);
1325 }
1326
1327 /* avoid warnings about uninitialized variables */
1328 b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
1329 a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
1330
1331 if (ats[atHD1] != NOTSET-12345)
1332 {
1333 if (ats[atHE2] != NOTSET-12345)
1334 {
1335 sprintf(resname, "HISH");
1336 }
1337 else
1338 {
1339 sprintf(resname, "HISA");
1340 }
1341 }
1342 else
1343 {
1344 sprintf(resname, "HISB");
1345 }
1346
1347 /* Get geometry from database */
1348 b_CG_ND1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "ND1");
1349 b_ND1_CE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "CE1");
1350 b_CE1_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "NE2");
1351 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "CD2");
1352 b_CD2_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "NE2");
1353 a_CG_ND1_CE1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "ND1", "CE1");
1354 a_CG_CD2_NE2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "CD2", "NE2");
1355 a_ND1_CE1_NE2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "ND1", "CE1", "NE2");
1356 a_CE1_NE2_CD2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "CD2");
1357
1358 if (ats[atHD1] != NOTSET-12345)
1359 {
1360 b_ND1_HD1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "HD1");
1361 a_CE1_ND1_HD1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "ND1", "HD1");
1362 }
1363 if (ats[atHE2] != NOTSET-12345)
1364 {
1365 b_NE2_HE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "NE2", "HE2");
1366 a_CE1_NE2_HE2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "HE2");
1367 }
1368 if (ats[atHD2] != NOTSET-12345)
1369 {
1370 b_CD2_HD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "HD2");
1371 a_NE2_CD2_HD2 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CD2", "HD2");
1372 }
1373 if (ats[atHE1] != NOTSET-12345)
1374 {
1375 b_CE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "HE1");
1376 a_NE2_CE1_HE1 = DEG2RAD(3.14159265358979323846/180.0)*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CE1", "HE1");
1377 }
1378
1379 /* constraints between CG, CE1 and NE1 */
1380 dCGCE1 = sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1)( sqr(b_CG_ND1) + sqr(b_ND1_CE1) - 2*b_CG_ND1*b_ND1_CE1*cos(a_CG_ND1_CE1
) ) );
1381 dCGNE2 = sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2)( sqr(b_CG_CD2) + sqr(b_CD2_NE2) - 2*b_CG_CD2*b_CD2_NE2*cos(a_CG_CD2_NE2
) ) );
1382
1383 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1384 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1385 /* we already have a constraint CE1-NE2, so we don't add it again */
1386
1387 /* calculate the positions in a local frame of reference.
1388 * The x-axis is the line from CG that makes a right angle
1389 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1390 */
1391 /* First calculate the x-axis intersection with y-axis (=yCE1).
1392 * Get cos(angle CG-CE1-NE2) :
1393 */
1394 cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2)( (sqr(dCGCE1)+sqr(b_CE1_NE2)-sqr(dCGNE2))/(2*dCGCE1*b_CE1_NE2
) );
1395 x[atCE1] = 0;
1396 y[atCE1] = cosalpha*dCGCE1;
1397 x[atNE2] = 0;
1398 y[atNE2] = y[atCE1]-b_CE1_NE2;
1399 sinalpha = sqrt(1-cosalpha*cosalpha);
1400 x[atCG] = -sinalpha*dCGCE1;
1401 y[atCG] = 0;
1402 x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1403 y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1404
1405 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1406
1407 x[atND1] = -b_ND1_CE1*sin(a_ND1_CE1_NE2);
1408 y[atND1] = y[atCE1]-b_ND1_CE1*cos(a_ND1_CE1_NE2);
1409
1410 x[atCD2] = -b_CD2_NE2*sin(a_CE1_NE2_CD2);
1411 y[atCD2] = y[atNE2]+b_CD2_NE2*cos(a_CE1_NE2_CD2);
1412
1413 /* And finally the hydrogen positions */
1414 if (ats[atHE1] != NOTSET-12345)
1415 {
1416 x[atHE1] = x[atCE1] + b_CE1_HE1*sin(a_NE2_CE1_HE1);
1417 y[atHE1] = y[atCE1] - b_CE1_HE1*cos(a_NE2_CE1_HE1);
1418 }
1419 /* HD2 - first get (ccw) angle from (positive) y-axis */
1420 if (ats[atHD2] != NOTSET-12345)
1421 {
1422 alpha = a_CE1_NE2_CD2 + M_PI3.14159265358979323846 - a_NE2_CD2_HD2;
1423 x[atHD2] = x[atCD2] - b_CD2_HD2*sin(alpha);
1424 y[atHD2] = y[atCD2] + b_CD2_HD2*cos(alpha);
1425 }
1426 if (ats[atHD1] != NOTSET-12345)
1427 {
1428 /* HD1 - first get (cw) angle from (positive) y-axis */
1429 alpha = a_ND1_CE1_NE2 + M_PI3.14159265358979323846 - a_CE1_ND1_HD1;
1430 x[atHD1] = x[atND1] - b_ND1_HD1*sin(alpha);
1431 y[atHD1] = y[atND1] - b_ND1_HD1*cos(alpha);
1432 }
1433 if (ats[atHE2] != NOTSET-12345)
1434 {
1435 x[atHE2] = x[atNE2] + b_NE2_HE2*sin(a_CE1_NE2_HE2);
1436 y[atHE2] = y[atNE2] + b_NE2_HE2*cos(a_CE1_NE2_HE2);
1437 }
1438 /* Have all coordinates now */
1439
1440 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1441 * set the rest to vsite3
1442 */
1443 mtot = xcom = ycom = 0;
1444 nvsite = 0;
1445 for (i = 0; i < atNR; i++)
1446 {
1447 if (ats[i] != NOTSET-12345)
1448 {
1449 mtot += at->atom[ats[i]].m;
1450 xcom += x[i]*at->atom[ats[i]].m;
1451 ycom += y[i]*at->atom[ats[i]].m;
1452 if (i != atCG && i != atCE1 && i != atNE2)
1453 {
1454 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1455 (*vsite_type)[ats[i]] = F_VSITE3;
1456 nvsite++;
1457 }
1458 }
1459 }
1460 if (nvsite+3 != nrfound)
1461 {
1462 gmx_incons("Generating vsites for HIS")_gmx_error("incons", "Generating vsites for HIS", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1462);
1463 }
1464
1465 xcom /= mtot;
1466 ycom /= mtot;
1467
1468 /* distribute mass so that com stays the same */
1469 mG = xcom*mtot/x[atCG];
1470 mrest = mtot-mG;
1471 mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1472 mNE2 = mrest-mCE1;
1473
1474 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1475 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1476 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1477
1478 /* HE1 */
1479 if (ats[atHE1] != NOTSET-12345)
1480 {
1481 calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1482 x[atCG], y[atCG], &a, &b);
1483 add_vsite3_param(&plist[F_VSITE3],
1484 ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1485 }
1486 /* HE2 */
1487 if (ats[atHE2] != NOTSET-12345)
1488 {
1489 calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1490 x[atCG], y[atCG], &a, &b);
1491 add_vsite3_param(&plist[F_VSITE3],
1492 ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1493 }
1494
1495 /* ND1 */
1496 calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1497 x[atCG], y[atCG], &a, &b);
1498 add_vsite3_param(&plist[F_VSITE3],
1499 ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1500
1501 /* CD2 */
1502 calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1503 x[atCG], y[atCG], &a, &b);
1504 add_vsite3_param(&plist[F_VSITE3],
1505 ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1506
1507 /* HD1 */
1508 if (ats[atHD1] != NOTSET-12345)
1509 {
1510 calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1511 x[atCG], y[atCG], &a, &b);
1512 add_vsite3_param(&plist[F_VSITE3],
1513 ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1514 }
1515 /* HD2 */
1516 if (ats[atHD2] != NOTSET-12345)
1517 {
1518 calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1519 x[atCG], y[atCG], &a, &b);
1520 add_vsite3_param(&plist[F_VSITE3],
1521 ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1522 }
1523 return nvsite;
1524}
1525
1526static gmx_bool is_vsite(int vsite_type)
1527{
1528 if (vsite_type == NOTSET-12345)
1529 {
1530 return FALSE0;
1531 }
1532 switch (abs(vsite_type) )
1533 {
1534 case F_VSITE3:
1535 case F_VSITE3FD:
1536 case F_VSITE3OUT:
1537 case F_VSITE3FAD:
1538 case F_VSITE4FD:
1539 case F_VSITE4FDN:
1540 return TRUE1;
1541 default:
1542 return FALSE0;
1543 }
1544}
1545
1546static char atomnamesuffix[] = "1234";
1547
1548void do_vsites(int nrtp, t_restp rtp[], gpp_atomtype_t atype,
1549 t_atoms *at, t_symtab *symtab, rvec *x[],
1550 t_params plist[], int *vsite_type[], int *cgnr[],
1551 real mHmult, gmx_bool bVsiteAromatics,
1552 const char *ffdir)
1553{
1554#define MAXATOMSPERRESIDUE16 16
1555 int i, j, k, m, i0, ni0, whatres, resind, add_shift, ftype, nvsite, nadd;
1556 int ai, aj, ak, al;
1557 int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1558 int Hatoms[4], heavies[4], bb;
1559 gmx_bool bWARNING, bAddVsiteParam, bFirstWater;
1560 matrix tmpmat;
1561 gmx_bool *bResProcessed;
1562 real mHtot, mtot, fact, fact2;
1563 rvec rpar, rperp, temp;
1564 char name[10], tpname[32], nexttpname[32], *ch;
1565 rvec *newx;
1566 int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE16];
1567 t_atom *newatom;
1568 t_params *params;
1569 char ***newatomname;
1570 char *resnm = NULL((void*)0);
1571 int ndb, f;
1572 char **db;
1573 int nvsiteconf, nvsitetop, cmplength;
1574 gmx_bool isN, planarN, bFound;
1575 gmx_residuetype_t rt;
1576
1577 t_vsiteconf *vsiteconflist;
1578 /* pointer to a list of CH3/NH3/NH2 configuration entries.
1579 * See comments in read_vsite_database. It isnt beautiful,
1580 * but it had to be fixed, and I dont even want to try to
1581 * maintain this part of the code...
1582 */
1583 t_vsitetop *vsitetop;
1584 /* Pointer to a list of geometry (bond/angle) entries for
1585 * residues like PHE, TRP, TYR, HIS, etc., where we need
1586 * to know the geometry to construct vsite aromatics.
1587 * Note that equilibrium geometry isnt necessarily the same
1588 * as the individual bond and angle values given in the
1589 * force field (rings can be strained).
1590 */
1591
1592 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1593 PHE, TRP, TYR and HIS to a construction of virtual sites */
1594 enum {
1595 resPHE, resTRP, resTYR, resHIS, resNR
1596 };
1597 const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1598 /* Amber03 alternative names for termini */
1599 const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1600 const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1601 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1602 gmx_bool bPartial[resNR] = { FALSE0, FALSE0, FALSE0, TRUE1 };
1603 /* the atnms for every residue MUST correspond to the enums in the
1604 gen_vsites_* (one for each residue) routines! */
1605 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1606 const char *atnms[resNR][MAXATOMSPERRESIDUE16+1] = {
1607 { "CG", /* PHE */
1608 "CD1", "HD1", "CD2", "HD2",
1609 "CE1", "HE1", "CE2", "HE2",
1610 "CZ", "HZ", NULL((void*)0) },
1611 { "CB", /* TRP */
1612 "CG",
1613 "CD1", "HD1", "CD2",
1614 "NE1", "HE1", "CE2", "CE3", "HE3",
1615 "CZ2", "HZ2", "CZ3", "HZ3",
1616 "CH2", "HH2", NULL((void*)0) },
1617 { "CG", /* TYR */
1618 "CD1", "HD1", "CD2", "HD2",
1619 "CE1", "HE1", "CE2", "HE2",
1620 "CZ", "OH", "HH", NULL((void*)0) },
1621 { "CG", /* HIS */
1622 "ND1", "HD1", "CD2", "HD2",
1623 "CE1", "HE1", "NE2", "HE2", NULL((void*)0) }
1624 };
1625
1626 if (debug)
1627 {
1628 printf("Searching for atoms to make virtual sites ...\n");
1629 fprintf(debug, "# # # VSITES # # #\n");
1630 }
1631
1632 ndb = fflib_search_file_end(ffdir, ".vsd", FALSE0, &db);
1633 nvsiteconf = 0;
1634 vsiteconflist = NULL((void*)0);
1635 nvsitetop = 0;
1636 vsitetop = NULL((void*)0);
1637 for (f = 0; f < ndb; f++)
1638 {
1639 read_vsite_database(db[f], &vsiteconflist, &nvsiteconf, &vsitetop, &nvsitetop);
1640 sfree(db[f])save_free("db[f]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1640, (db[f]));
1641 }
1642 sfree(db)save_free("db", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1642, (db));
1643
1644 bFirstWater = TRUE1;
1645 nvsite = 0;
1646 nadd = 0;
1647 /* we need a marker for which atoms should *not* be renumbered afterwards */
1648 add_shift = 10*at->nr;
1649 /* make arrays where masses can be inserted into */
1650 snew(newx, at->nr)(newx) = save_calloc("newx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1650, (at->nr), sizeof(*(newx)));
1651 snew(newatom, at->nr)(newatom) = save_calloc("newatom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1651, (at->nr), sizeof(*(newatom)));
1652 snew(newatomname, at->nr)(newatomname) = save_calloc("newatomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1652, (at->nr), sizeof(*(newatomname)));
1653 snew(newvsite_type, at->nr)(newvsite_type) = save_calloc("newvsite_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1653, (at->nr), sizeof(*(newvsite_type)));
1654 snew(newcgnr, at->nr)(newcgnr) = save_calloc("newcgnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1654, (at->nr), sizeof(*(newcgnr)));
1655 /* make index array to tell where the atoms go to when masses are inserted */
1656 snew(o2n, at->nr)(o2n) = save_calloc("o2n", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1656, (at->nr), sizeof(*(o2n)));
1657 for (i = 0; i < at->nr; i++)
1658 {
1659 o2n[i] = i;
1660 }
1661 /* make index to tell which residues were already processed */
1662 snew(bResProcessed, at->nres)(bResProcessed) = save_calloc("bResProcessed", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1662, (at->nres), sizeof(*(bResProcessed)));
1663
1664 gmx_residuetype_init(&rt);
1665
1666 /* generate vsite constructions */
1667 /* loop over all atoms */
1668 resind = -1;
1669 for (i = 0; (i < at->nr); i++)
1670 {
1671 if (at->atom[i].resind != resind)
1672 {
1673 resind = at->atom[i].resind;
1674 resnm = *(at->resinfo[resind].name);
1675 }
1676 /* first check for aromatics to virtualize */
1677 /* don't waste our effort on DNA, water etc. */
1678 /* Only do the vsite aromatic stuff when we reach the
1679 * CA atom, since there might be an X2/X3 group on the
1680 * N-terminus that must be treated first.
1681 */
1682 if (bVsiteAromatics &&
1683 !strcmp(*(at->atomname[i]), "CA")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*(at->atomname[i])) && __builtin_constant_p ("CA"
) && (__s1_len = strlen (*(at->atomname[i])), __s2_len
= strlen ("CA"), (!((size_t)(const void *)((*(at->atomname
[i])) + 1) - (size_t)(const void *)(*(at->atomname[i])) ==
1) || __s1_len >= 4) && (!((size_t)(const void *)
(("CA") + 1) - (size_t)(const void *)("CA") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (*(at->atomname[i]), "CA") :
(__builtin_constant_p (*(at->atomname[i])) && ((size_t
)(const void *)((*(at->atomname[i])) + 1) - (size_t)(const
void *)(*(at->atomname[i])) == 1) && (__s1_len = strlen
(*(at->atomname[i])), __s1_len < 4) ? (__builtin_constant_p
("CA") && ((size_t)(const void *)(("CA") + 1) - (size_t
)(const void *)("CA") == 1) ? __builtin_strcmp (*(at->atomname
[i]), "CA") : (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) ("CA"); int __result = (
((const unsigned char *) (const char *) (*(at->atomname[i]
)))[0] - __s2[0]); if (__s1_len > 0 && __result ==
0) { __result = (((const unsigned char *) (const char *) (*(
at->atomname[i])))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (*(at->atomname[i])))[2] - __s2[2]); if (__s1_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (*(at->atomname[i])))[3] - __s2[3]
); } } __result; }))) : (__builtin_constant_p ("CA") &&
((size_t)(const void *)(("CA") + 1) - (size_t)(const void *)
("CA") == 1) && (__s2_len = strlen ("CA"), __s2_len <
4) ? (__builtin_constant_p (*(at->atomname[i])) &&
((size_t)(const void *)((*(at->atomname[i])) + 1) - (size_t
)(const void *)(*(at->atomname[i])) == 1) ? __builtin_strcmp
(*(at->atomname[i]), "CA") : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (*(at->
atomname[i])); int __result = (((const unsigned char *) (const
char *) ("CA"))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("CA"))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("CA"))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) ("CA"))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(*(at->atomname[i]), "CA")))); }) &&
1684 !bResProcessed[resind] &&
1685 gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name)) )
1686 {
1687 /* mark this residue */
1688 bResProcessed[resind] = TRUE1;
1689 /* find out if this residue needs converting */
1690 whatres = NOTSET-12345;
1691 for (j = 0; j < resNR && whatres == NOTSET-12345; j++)
1692 {
1693
1694 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1695
1696 bFound = ((gmx_strncasecmp(resnm, resnms[j], cmplength) == 0) ||
1697 (gmx_strncasecmp(resnm, resnmsN[j], cmplength) == 0) ||
1698 (gmx_strncasecmp(resnm, resnmsC[j], cmplength) == 0));
1699
1700 if (bFound)
1701 {
1702 whatres = j;
1703 /* get atoms we will be needing for the conversion */
1704 nrfound = 0;
1705 for (k = 0; atnms[j][k]; k++)
1706 {
1707 ats[k] = NOTSET-12345;
1708 for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET-12345; m++)
1709 {
1710 if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1711 {
1712 ats[k] = m;
1713 nrfound++;
1714 }
1715 }
1716 }
1717
1718 /* now k is number of atom names in atnms[j] */
1719 if (j == resHIS)
1720 {
1721 needed = k-3;
1722 }
1723 else
1724 {
1725 needed = k;
1726 }
1727 if (nrfound < needed)
1728 {
1729 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1729, "not enough atoms found (%d, need %d) in "
1730 "residue %s %d while\n "
1731 "generating aromatics virtual site construction",
1732 nrfound, needed, resnm, at->resinfo[resind].nr);
1733 }
1734 /* Advance overall atom counter */
1735 i++;
1736 }
1737 }
1738 /* the enums for every residue MUST correspond to atnms[residue] */
1739 switch (whatres)
1740 {
1741 case resPHE:
1742 if (debug)
1743 {
1744 fprintf(stderrstderr, "PHE at %d\n", o2n[ats[0]]+1);
1745 }
1746 nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1747 break;
1748 case resTRP:
1749 if (debug)
1750 {
1751 fprintf(stderrstderr, "TRP at %d\n", o2n[ats[0]]+1);
1752 }
1753 nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1754 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1755 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1756 break;
1757 case resTYR:
1758 if (debug)
1759 {
1760 fprintf(stderrstderr, "TYR at %d\n", o2n[ats[0]]+1);
1761 }
1762 nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1763 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1764 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1765 break;
1766 case resHIS:
1767 if (debug)
1768 {
1769 fprintf(stderrstderr, "HIS at %d\n", o2n[ats[0]]+1);
1770 }
1771 nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1772 break;
1773 case NOTSET-12345:
1774 /* this means this residue won't be processed */
1775 break;
1776 default:
1777 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1777, "DEATH HORROR in do_vsites (%s:%d)",
1778 __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c", __LINE__1778);
1779 } /* switch whatres */
1780 /* skip back to beginning of residue */
1781 while (i > 0 && at->atom[i-1].resind == resind)
1782 {
1783 i--;
1784 }
1785 } /* if bVsiteAromatics & is protein */
1786
1787 /* now process the rest of the hydrogens */
1788 /* only process hydrogen atoms which are not already set */
1789 if ( ((*vsite_type)[i] == NOTSET-12345) && is_hydrogen(*(at->atomname[i])))
1790 {
1791 /* find heavy atom, count #bonds from it and #H atoms bound to it
1792 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1793 count_bonds(i, &plist[F_BONDS], at->atomname,
1794 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1795 /* get Heavy atom type */
1796 tpHeavy = get_atype(Heavy, at, nrtp, rtp, rt);
1797 strcpy(tpname, get_atomtype_name(tpHeavy, atype));
1798
1799 bWARNING = FALSE0;
1800 bAddVsiteParam = TRUE1;
1801 /* nested if's which check nrHatoms, nrbonds and atomname */
1802 if (nrHatoms == 1)
1803 {
1804 switch (nrbonds)
1805 {
1806 case 2: /* -O-H */
1807 (*vsite_type)[i] = F_BONDS;
1808 break;
1809 case 3: /* =CH-, -NH- or =NH+- */
1810 (*vsite_type)[i] = F_VSITE3FD;
1811 break;
1812 case 4: /* --CH- (tert) */
1813 /* The old type 4FD had stability issues, so
1814 * all new constructs should use 4FDN
1815 */
1816 (*vsite_type)[i] = F_VSITE4FDN;
1817
1818 /* Check parity of heavy atoms from coordinates */
1819 ai = Heavy;
1820 aj = heavies[0];
1821 ak = heavies[1];
1822 al = heavies[2];
1823 rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1824 rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1825 rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1826
1827 if (det(tmpmat) > 0)
1828 {
1829 /* swap parity */
1830 heavies[1] = aj;
1831 heavies[0] = ak;
1832 }
1833
1834 break;
1835 default: /* nrbonds != 2, 3 or 4 */
1836 bWARNING = TRUE1;
1837 }
1838
1839 }
1840 else if ( /*(nrHatoms == 2) && (nrbonds == 2) && REMOVED this test
1841 DvdS 19-01-04 */
1842 (gmx_strncasecmp(*at->atomname[Heavy], "OW", 2) == 0) )
1843 {
1844 bAddVsiteParam = FALSE0; /* this is water: skip these hydrogens */
1845 if (bFirstWater)
1846 {
1847 bFirstWater = FALSE0;
1848 if (debug)
1849 {
1850 fprintf(debug,
1851 "Not converting hydrogens in water to virtual sites\n");
1852 }
1853 }
1854 }
1855 else if ( (nrHatoms == 2) && (nrbonds == 4) )
1856 {
1857 /* -CH2- , -NH2+- */
1858 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1859 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1860 }
1861 else
1862 {
1863 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1864 * If it is a nitrogen, first check if it is planar.
1865 */
1866 isN = planarN = FALSE0;
1867 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1868 {
1869 isN = TRUE1;
1870 j = nitrogen_is_planar(vsiteconflist, nvsiteconf, tpname);
1871 if (j < 0)
1872 {
1873 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1873, "No vsite database NH2 entry for type %s\n", tpname);
1874 }
1875 planarN = (j == 1);
1876 }
1877 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
1878 {
1879 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1880 (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1881 (*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
1882 }
1883 else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1884 ( isN && !planarN ) ) ||
1885 ( (nrHatoms == 3) && (nrbonds == 4) ) )
1886 {
1887 /* CH3, NH3 or non-planar NH2 group */
1888 int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1889 gmx_bool Hat_SwapParity[3] = { FALSE0, TRUE1, FALSE0 };
1890
1891 if (debug)
1892 {
1893 fprintf(stderrstderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
1894 }
1895 bAddVsiteParam = FALSE0; /* we'll do this ourselves! */
1896 /* -NH2 (umbrella), -NH3+ or -CH3 */
1897 (*vsite_type)[Heavy] = F_VSITE3;
1898 for (j = 0; j < nrHatoms; j++)
1899 {
1900 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1901 }
1902 /* get dummy mass type from first char of heavy atom type (N or C) */
1903
1904 strcpy(nexttpname, get_atomtype_name(get_atype(heavies[0], at, nrtp, rtp, rt), atype));
1905 ch = get_dummymass_name(vsiteconflist, nvsiteconf, tpname, nexttpname);
1906
1907 if (ch == NULL((void*)0))
1908 {
1909 if (ndb > 0)
1910 {
1911 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1911, "Can't find dummy mass for type %s bonded to type %s in the virtual site database (.vsd files). Add it to the database!\n", tpname, nexttpname);
1912 }
1913 else
1914 {
1915 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1915, "A dummy mass for type %s bonded to type %s is required, but no virtual site database (.vsd) files where found.\n", tpname, nexttpname);
1916 }
1917 }
1918 else
1919 {
1920 strcpy(name, ch);
1921 }
1922
1923 tpM = vsite_nm2type(name, atype);
1924 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1925#define NMASS 2
1926 i0 = Heavy;
1927 ni0 = i0+nadd;
1928 if (debug)
1929 {
1930 fprintf(stderrstderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
1931 }
1932 nadd += NMASS;
1933 for (j = i0; j < at->nr; j++)
1934 {
1935 o2n[j] = j+nadd;
1936 }
1937
1938 srenew(newx, at->nr+nadd)(newx) = save_realloc("newx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1938, (newx), (at->nr+nadd), sizeof(*(newx)));
1939 srenew(newatom, at->nr+nadd)(newatom) = save_realloc("newatom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1939, (newatom), (at->nr+nadd), sizeof(*(newatom)));
1940 srenew(newatomname, at->nr+nadd)(newatomname) = save_realloc("newatomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1940, (newatomname), (at->nr+nadd), sizeof(*(newatomname
)));
1941 srenew(newvsite_type, at->nr+nadd)(newvsite_type) = save_realloc("newvsite_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1941, (newvsite_type), (at->nr+nadd), sizeof(*(newvsite_type
)));
1942 srenew(newcgnr, at->nr+nadd)(newcgnr) = save_realloc("newcgnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 1942, (newcgnr), (at->nr+nadd), sizeof(*(newcgnr)));
1943
1944 for (j = 0; j < NMASS; j++)
1945 {
1946 newatomname[at->nr+nadd-1-j] = NULL((void*)0);
1947 }
1948
1949 /* calculate starting position for the masses */
1950 mHtot = 0;
1951 /* get atom masses, and set Heavy and Hatoms mass to zero */
1952 for (j = 0; j < nrHatoms; j++)
1953 {
1954 mHtot += get_amass(Hatoms[j], at, nrtp, rtp, rt);
1955 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1956 }
1957 mtot = mHtot + get_amass(Heavy, at, nrtp, rtp, rt);
1958 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1959 if (mHmult != 1.0)
1960 {
1961 mHtot *= mHmult;
1962 }
1963 fact2 = mHtot/mtot;
1964 fact = sqrt(fact2);
1965 /* generate vectors parallel and perpendicular to rotational axis:
1966 * rpar = Heavy -> Hcom
1967 * rperp = Hcom -> H1 */
1968 clear_rvec(rpar);
1969 for (j = 0; j < nrHatoms; j++)
1970 {
1971 rvec_inc(rpar, (*x)[Hatoms[j]]);
1972 }
1973 svmul(1.0/nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
1974 rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
1975 rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
1976 rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
1977 /* calc mass positions */
1978 svmul(fact2, rpar, temp);
1979 for (j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1980 {
1981 rvec_add((*x)[Heavy], temp, newx[ni0+j]);
1982 }
1983 svmul(fact, rperp, temp);
1984 rvec_inc(newx[ni0 ], temp);
1985 rvec_dec(newx[ni0+1], temp);
1986 /* set atom parameters for the masses */
1987 for (j = 0; (j < NMASS); j++)
1988 {
1989 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1990 name[0] = 'M';
1991 for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
1992 {
1993 name[k+1] = (*at->atomname[Heavy])[k];
1994 }
1995 name[k+1] = atomnamesuffix[j];
1996 name[k+2] = '\0';
1997 newatomname[ni0+j] = put_symtab(symtab, name);
1998 newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
1999 newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
2000 newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
2001 newatom[ni0+j].ptype = eptAtom;
2002 newatom[ni0+j].resind = at->atom[i0].resind;
2003 newatom[ni0+j].elem[0] = 'M';
2004 newatom[ni0+j].elem[1] = '\0';
2005 newvsite_type[ni0+j] = NOTSET-12345;
2006 newcgnr[ni0+j] = (*cgnr)[i0];
2007 }
2008 /* add constraints between dummy masses and to heavies[0] */
2009 /* 'add_shift' says which atoms won't be renumbered afterwards */
2010 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0, NOTSET-12345);
2011 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0+1, NOTSET-12345);
2012 my_add_param(&(plist[F_CONSTRNC]), add_shift+ni0, add_shift+ni0+1, NOTSET-12345);
2013
2014 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
2015 /* note that vsite_type cannot be NOTSET, because we just set it */
2016 add_vsite3_atoms (&plist[(*vsite_type)[Heavy]],
2017 Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
2018 FALSE0);
2019 for (j = 0; j < nrHatoms; j++)
2020 {
2021 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
2022 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
2023 Hat_SwapParity[j]);
2024 }
2025#undef NMASS
2026 }
2027 else
2028 {
2029 bWARNING = TRUE1;
2030 }
2031
2032 }
2033 if (bWARNING)
2034 {
2035 fprintf(stderrstderr,
2036 "Warning: cannot convert atom %d %s (bound to a heavy atom "
2037 "%s with \n"
2038 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2039 i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2040 }
2041 if (bAddVsiteParam)
2042 {
2043 /* add vsite parameters to topology,
2044 also get rid of negative vsite_types */
2045 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
2046 nrheavies, heavies);
2047 /* transfer mass of virtual site to Heavy atom */
2048 for (j = 0; j < nrHatoms; j++)
2049 {
2050 if (is_vsite((*vsite_type)[Hatoms[j]]))
2051 {
2052 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
2053 at->atom[Heavy].mB = at->atom[Heavy].m;
2054 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2055 }
2056 }
2057 }
2058 nvsite += nrHatoms;
2059 if (debug)
2060 {
2061 fprintf(debug, "atom %d: ", o2n[i]+1);
2062 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2063 }
2064 } /* if vsite NOTSET & is hydrogen */
2065
2066 } /* for i < at->nr */
2067
2068 gmx_residuetype_destroy(rt);
2069
2070 if (debug)
2071 {
2072 fprintf(debug, "Before inserting new atoms:\n");
2073 for (i = 0; i < at->nr; i++)
2074 {
2075 fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i+1, o2n[i]+1,
2076 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2077 at->resinfo[at->atom[i].resind].nr,
2078 at->resinfo[at->atom[i].resind].name ?
2079 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2080 (*cgnr)[i],
2081 ((*vsite_type)[i] == NOTSET-12345) ?
2082 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2083 }
2084 fprintf(debug, "new atoms to be inserted:\n");
2085 for (i = 0; i < at->nr+nadd; i++)
2086 {
2087 if (newatomname[i])
2088 {
2089 fprintf(debug, "%4d %4s %4d %6d %-10s\n", i+1,
2090 newatomname[i] ? *(newatomname[i]) : "(NULL)",
2091 newatom[i].resind, newcgnr[i],
2092 (newvsite_type[i] == NOTSET-12345) ?
2093 "NOTSET" : interaction_function[newvsite_type[i]].name);
2094 }
2095 }
2096 }
2097
2098 /* add all original atoms to the new arrays, using o2n index array */
2099 for (i = 0; i < at->nr; i++)
2100 {
2101 newatomname [o2n[i]] = at->atomname [i];
2102 newatom [o2n[i]] = at->atom [i];
2103 newvsite_type[o2n[i]] = (*vsite_type)[i];
2104 newcgnr [o2n[i]] = (*cgnr) [i];
2105 copy_rvec((*x)[i], newx[o2n[i]]);
2106 }
2107 /* throw away old atoms */
2108 sfree(at->atom)save_free("at->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2108, (at->atom));
2109 sfree(at->atomname)save_free("at->atomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2109, (at->atomname));
2110 sfree(*vsite_type)save_free("*vsite_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2110, (*vsite_type));
2111 sfree(*cgnr)save_free("*cgnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2111, (*cgnr));
2112 sfree(*x)save_free("*x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2112, (*x));
2113 /* put in the new ones */
2114 at->nr += nadd;
2115 at->atom = newatom;
2116 at->atomname = newatomname;
2117 *vsite_type = newvsite_type;
2118 *cgnr = newcgnr;
2119 *x = newx;
2120 if (at->nr > add_shift)
2121 {
2122 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2122, "Added impossible amount of dummy masses "
2123 "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
2124 }
2125
2126 if (debug)
2127 {
2128 fprintf(debug, "After inserting new atoms:\n");
2129 for (i = 0; i < at->nr; i++)
2130 {
2131 fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i+1,
2132 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2133 at->resinfo[at->atom[i].resind].nr,
2134 at->resinfo[at->atom[i].resind].name ?
2135 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2136 (*cgnr)[i],
2137 ((*vsite_type)[i] == NOTSET-12345) ?
2138 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2139 }
2140 }
2141
2142 /* now renumber all the interactions because of the added atoms */
2143 for (ftype = 0; ftype < F_NRE; ftype++)
2144 {
2145 params = &(plist[ftype]);
2146 if (debug)
2147 {
2148 fprintf(debug, "Renumbering %d %s\n", params->nr,
2149 interaction_function[ftype].longname);
2150 }
2151 for (i = 0; i < params->nr; i++)
2152 {
2153 for (j = 0; j < NRAL(ftype)(interaction_function[(ftype)].nratoms); j++)
2154 {
2155 if (params->param[i].a[j] >= add_shift)
2156 {
2157 if (debug)
2158 {
2159 fprintf(debug, " [%u -> %u]", params->param[i].a[j],
2160 params->param[i].a[j]-add_shift);
2161 }
2162 params->param[i].a[j] = params->param[i].a[j]-add_shift;
2163 }
2164 else
2165 {
2166 if (debug)
2167 {
2168 fprintf(debug, " [%u -> %d]", params->param[i].a[j],
2169 o2n[params->param[i].a[j]]);
2170 }
2171 params->param[i].a[j] = o2n[params->param[i].a[j]];
2172 }
2173 }
2174 if (debug)
2175 {
2176 fprintf(debug, "\n");
2177 }
2178 }
2179 }
2180 /* now check if atoms in the added constraints are in increasing order */
2181 params = &(plist[F_CONSTRNC]);
2182 for (i = 0; i < params->nr; i++)
2183 {
2184 if (params->param[i].AIa[0] > params->param[i].AJa[1])
2185 {
2186 j = params->param[i].AJa[1];
2187 params->param[i].AJa[1] = params->param[i].AIa[0];
2188 params->param[i].AIa[0] = j;
2189 }
2190 }
2191
2192 /* clean up */
2193 sfree(o2n)save_free("o2n", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2193, (o2n));
2194
2195 /* tell the user what we did */
2196 fprintf(stderrstderr, "Marked %d virtual sites\n", nvsite);
2197 fprintf(stderrstderr, "Added %d dummy masses\n", nadd);
2198 fprintf(stderrstderr, "Added %d new constraints\n", plist[F_CONSTRNC].nr);
2199}
2200
2201void do_h_mass(t_params *psb, int vsite_type[], t_atoms *at, real mHmult,
2202 gmx_bool bDeuterate)
2203{
2204 int i, j, a;
2205
2206 /* loop over all atoms */
2207 for (i = 0; i < at->nr; i++)
2208 {
2209 /* adjust masses if i is hydrogen and not a virtual site */
2210 if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
2211 {
2212 /* find bonded heavy atom */
2213 a = NOTSET-12345;
2214 for (j = 0; (j < psb->nr) && (a == NOTSET-12345); j++)
2215 {
2216 /* if other atom is not a virtual site, it is the one we want */
2217 if ( (psb->param[j].AIa[0] == i) &&
2218 !is_vsite(vsite_type[psb->param[j].AJa[1]]) )
2219 {
2220 a = psb->param[j].AJa[1];
2221 }
2222 else if ( (psb->param[j].AJa[1] == i) &&
2223 !is_vsite(vsite_type[psb->param[j].AIa[0]]) )
2224 {
2225 a = psb->param[j].AIa[0];
2226 }
2227 }
2228 if (a == NOTSET-12345)
2229 {
2230 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/gen_vsite.c"
, 2230, "Unbound hydrogen atom (%d) found while adjusting mass",
2231 i+1);
2232 }
2233
2234 /* adjust mass of i (hydrogen) with mHmult
2235 and correct mass of a (bonded atom) with same amount */
2236 if (!bDeuterate)
2237 {
2238 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
2239 at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
2240 }
2241 at->atom[i].m *= mHmult;
2242 at->atom[i].mB *= mHmult;
2243 }
2244 }
2245}