File: | gromacs/gmxpreprocess/gen_vsite.c |
Location: | line 742, column 5 |
Description: | Value stored to 'xCE2' is never read |
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 | |
67 | typedef 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 | */ |
84 | typedef 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 | |
102 | enum { |
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 | |
107 | typedef char t_dirname[STRLEN4096]; |
108 | |
109 | static 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 | |
121 | static 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 | |
143 | static 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 | |
322 | static 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 | |
345 | static 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 | |
367 | static 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 | |
398 | static 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 | |
433 | static 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 | |
493 | static 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 | |
512 | static 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 | |
536 | static 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 | |
550 | static 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 | |
574 | static 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 | |
583 | static 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 | |
690 | static 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; |
741 | yCE1 = bond_cc*sin(0.5*ANGLE_6RING((3.14159265358979323846/180.0)*120)); |
742 | xCE2 = 0; |
Value stored to 'xCE2' is never read | |
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 | |
786 | static 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 | |
838 | static 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 | |
859 | static 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 | |
1137 | static 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 | |
1296 | static 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 | |
1526 | static 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 | |
1546 | static char atomnamesuffix[] = "1234"; |
1547 | |
1548 | void 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 | |
2201 | void 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 | } |