File: | gromacs/gmxpreprocess/toppush.c |
Location: | line 2110, column 5 |
Description: | Value stored to 'nr' 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 <assert.h> |
42 | #include <ctype.h> |
43 | #include <math.h> |
44 | #include <stdlib.h> |
45 | |
46 | #include "macros.h" |
47 | #include "names.h" |
48 | #include "toputil.h" |
49 | #include "toppush.h" |
50 | #include "topdirs.h" |
51 | #include "readir.h" |
52 | #include "symtab.h" |
53 | #include "warninp.h" |
54 | #include "gpp_atomtype.h" |
55 | #include "gpp_bond_atomtype.h" |
56 | |
57 | #include "gromacs/utility/cstringutil.h" |
58 | #include "gromacs/utility/fatalerror.h" |
59 | #include "gromacs/utility/smalloc.h" |
60 | |
61 | void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype, |
62 | warninp_t wi) |
63 | { |
64 | int i, j, k = -1, nf; |
65 | int nr, nrfp; |
66 | real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2; |
67 | char errbuf[256]; |
68 | |
69 | /* Lean mean shortcuts */ |
70 | nr = get_atomtype_ntypes(atype); |
71 | nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); |
72 | snew(plist->param, nr*nr)(plist->param) = save_calloc("plist->param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 72, (nr*nr), sizeof(*(plist->param))); |
73 | plist->nr = nr*nr; |
74 | |
75 | /* Fill the matrix with force parameters */ |
76 | switch (ftype) |
77 | { |
78 | case F_LJ: |
79 | switch (comb) |
80 | { |
81 | case eCOMB_GEOMETRIC: |
82 | /* Gromos rules */ |
83 | for (i = k = 0; (i < nr); i++) |
84 | { |
85 | for (j = 0; (j < nr); j++, k++) |
86 | { |
87 | for (nf = 0; (nf < nrfp); nf++) |
88 | { |
89 | ci = get_atomtype_nbparam(i, nf, atype); |
90 | cj = get_atomtype_nbparam(j, nf, atype); |
91 | c = sqrt(ci * cj); |
92 | plist->param[k].c[nf] = c; |
93 | } |
94 | } |
95 | } |
96 | break; |
97 | |
98 | case eCOMB_ARITHMETIC: |
99 | /* c0 and c1 are epsilon and sigma */ |
100 | for (i = k = 0; (i < nr); i++) |
101 | { |
102 | for (j = 0; (j < nr); j++, k++) |
103 | { |
104 | ci0 = get_atomtype_nbparam(i, 0, atype); |
105 | cj0 = get_atomtype_nbparam(j, 0, atype); |
106 | ci1 = get_atomtype_nbparam(i, 1, atype); |
107 | cj1 = get_atomtype_nbparam(j, 1, atype); |
108 | plist->param[k].c[0] = (ci0+cj0)*0.5; |
109 | plist->param[k].c[1] = sqrt(ci1*cj1); |
110 | } |
111 | } |
112 | |
113 | break; |
114 | case eCOMB_GEOM_SIG_EPS: |
115 | /* c0 and c1 are epsilon and sigma */ |
116 | for (i = k = 0; (i < nr); i++) |
117 | { |
118 | for (j = 0; (j < nr); j++, k++) |
119 | { |
120 | ci0 = get_atomtype_nbparam(i, 0, atype); |
121 | cj0 = get_atomtype_nbparam(j, 0, atype); |
122 | ci1 = get_atomtype_nbparam(i, 1, atype); |
123 | cj1 = get_atomtype_nbparam(j, 1, atype); |
124 | plist->param[k].c[0] = sqrt(ci0*cj0); |
125 | plist->param[k].c[1] = sqrt(ci1*cj1); |
126 | } |
127 | } |
128 | |
129 | break; |
130 | default: |
131 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 131, "No such combination rule %d", comb); |
132 | } |
133 | if (plist->nr != k) |
134 | { |
135 | gmx_incons("Topology processing, generate nb parameters")_gmx_error("incons", "Topology processing, generate nb parameters" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 135); |
136 | } |
137 | break; |
138 | |
139 | case F_BHAM: |
140 | /* Buckingham rules */ |
141 | for (i = k = 0; (i < nr); i++) |
142 | { |
143 | for (j = 0; (j < nr); j++, k++) |
144 | { |
145 | ci0 = get_atomtype_nbparam(i, 0, atype); |
146 | cj0 = get_atomtype_nbparam(j, 0, atype); |
147 | ci2 = get_atomtype_nbparam(i, 2, atype); |
148 | cj2 = get_atomtype_nbparam(j, 2, atype); |
149 | bi = get_atomtype_nbparam(i, 1, atype); |
150 | bj = get_atomtype_nbparam(j, 1, atype); |
151 | plist->param[k].c[0] = sqrt(ci0 * cj0); |
152 | if ((bi == 0) || (bj == 0)) |
153 | { |
154 | plist->param[k].c[1] = 0; |
155 | } |
156 | else |
157 | { |
158 | plist->param[k].c[1] = 2.0/(1/bi+1/bj); |
159 | } |
160 | plist->param[k].c[2] = sqrt(ci2 * cj2); |
161 | } |
162 | } |
163 | |
164 | break; |
165 | default: |
166 | sprintf(errbuf, "Invalid nonbonded type %s", |
167 | interaction_function[ftype].longname); |
168 | warning_error(wi, errbuf); |
169 | } |
170 | } |
171 | |
172 | static void realloc_nb_params(gpp_atomtype_t at, |
173 | t_nbparam ***nbparam, t_nbparam ***pair) |
174 | { |
175 | /* Add space in the non-bonded parameters matrix */ |
176 | int atnr = get_atomtype_ntypes(at); |
177 | srenew(*nbparam, atnr)(*nbparam) = save_realloc("*nbparam", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 177, (*nbparam), (atnr), sizeof(*(*nbparam))); |
178 | snew((*nbparam)[atnr-1], atnr)((*nbparam)[atnr-1]) = save_calloc("(*nbparam)[atnr-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 178, (atnr), sizeof(*((*nbparam)[atnr-1]))); |
179 | if (pair) |
180 | { |
181 | srenew(*pair, atnr)(*pair) = save_realloc("*pair", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 181, (*pair), (atnr), sizeof(*(*pair))); |
182 | snew((*pair)[atnr-1], atnr)((*pair)[atnr-1]) = save_calloc("(*pair)[atnr-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 182, (atnr), sizeof(*((*pair)[atnr-1]))); |
183 | } |
184 | } |
185 | |
186 | static void copy_B_from_A(int ftype, double *c) |
187 | { |
188 | int nrfpA, nrfpB, i; |
189 | |
190 | nrfpA = NRFPA(ftype)(interaction_function[(ftype)].nrfpA); |
191 | nrfpB = NRFPB(ftype)(interaction_function[(ftype)].nrfpB); |
192 | |
193 | /* Copy the B parameters from the first nrfpB A parameters */ |
194 | for (i = 0; (i < nrfpB); i++) |
195 | { |
196 | c[nrfpA+i] = c[i]; |
197 | } |
198 | } |
199 | |
200 | void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat, |
201 | char *line, int nb_funct, |
202 | t_nbparam ***nbparam, t_nbparam ***pair, |
203 | warninp_t wi) |
204 | { |
205 | typedef struct { |
206 | const char *entry; |
207 | int ptype; |
208 | } t_xlate; |
209 | t_xlate xl[eptNR] = { |
210 | { "A", eptAtom }, |
211 | { "N", eptNucleus }, |
212 | { "S", eptShell }, |
213 | { "B", eptBond }, |
214 | { "V", eptVSite }, |
215 | }; |
216 | |
217 | int nr, i, nfields, j, pt, nfp0 = -1; |
218 | int batype_nr, nread; |
219 | char type[STRLEN4096], btype[STRLEN4096], ptype[STRLEN4096]; |
220 | double m, q; |
221 | double c[MAXFORCEPARAM12]; |
222 | double radius, vol, surftens, gb_radius, S_hct; |
223 | char tmpfield[12][100]; /* Max 12 fields of width 100 */ |
224 | char errbuf[256]; |
225 | t_atom *atom; |
226 | t_param *param; |
227 | int atomnr; |
228 | gmx_bool have_atomic_number; |
229 | gmx_bool have_bonded_type; |
230 | |
231 | snew(atom, 1)(atom) = save_calloc("atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 231, (1), sizeof(*(atom))); |
232 | snew(param, 1)(param) = save_calloc("param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 232, (1), sizeof(*(param))); |
233 | |
234 | /* First assign input line to temporary array */ |
235 | nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s", |
236 | tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5], |
237 | tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]); |
238 | |
239 | /* Comments on optional fields in the atomtypes section: |
240 | * |
241 | * The force field format is getting a bit old. For OPLS-AA we needed |
242 | * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff |
243 | * we also needed the atomic numbers. |
244 | * To avoid making all old or user-generated force fields unusable we |
245 | * have introduced both these quantities as optional columns, and do some |
246 | * acrobatics to check whether they are present or not. |
247 | * This will all look much nicer when we switch to XML... sigh. |
248 | * |
249 | * Field 0 (mandatory) is the nonbonded type name. (string) |
250 | * Field 1 (optional) is the bonded type (string) |
251 | * Field 2 (optional) is the atomic number (int) |
252 | * Field 3 (mandatory) is the mass (numerical) |
253 | * Field 4 (mandatory) is the charge (numerical) |
254 | * Field 5 (mandatory) is the particle type (single character) |
255 | * This is followed by a number of nonbonded parameters. |
256 | * |
257 | * The safest way to identify the format is the particle type field. |
258 | * |
259 | * So, here is what we do: |
260 | * |
261 | * A. Read in the first six fields as strings |
262 | * B. If field 3 (starting from 0) is a single char, we have neither |
263 | * bonded_type or atomic numbers. |
264 | * C. If field 5 is a single char we have both. |
265 | * D. If field 4 is a single char we check field 1. If this begins with |
266 | * an alphabetical character we have bonded types, otherwise atomic numbers. |
267 | */ |
268 | |
269 | if (nfields < 6) |
270 | { |
271 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 271); |
272 | return; |
273 | } |
274 | |
275 | if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0])((*__ctype_b_loc ())[(int) ((tmpfield[5][0]))] & (unsigned short int) _ISalpha) ) |
276 | { |
277 | have_bonded_type = TRUE1; |
278 | have_atomic_number = TRUE1; |
279 | } |
280 | else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0])((*__ctype_b_loc ())[(int) ((tmpfield[3][0]))] & (unsigned short int) _ISalpha) ) |
281 | { |
282 | have_bonded_type = FALSE0; |
283 | have_atomic_number = FALSE0; |
284 | } |
285 | else |
286 | { |
287 | have_bonded_type = ( isalpha(tmpfield[1][0])((*__ctype_b_loc ())[(int) ((tmpfield[1][0]))] & (unsigned short int) _ISalpha) != 0 ); |
288 | have_atomic_number = !have_bonded_type; |
289 | } |
290 | |
291 | /* optional fields */ |
292 | surftens = -1; |
293 | vol = -1; |
294 | radius = -1; |
295 | gb_radius = -1; |
296 | atomnr = -1; |
297 | S_hct = -1; |
298 | |
299 | switch (nb_funct) |
300 | { |
301 | |
302 | case F_LJ: |
303 | nfp0 = 2; |
304 | |
305 | if (have_atomic_number) |
306 | { |
307 | if (have_bonded_type) |
308 | { |
309 | nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf", |
310 | type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], |
311 | &radius, &vol, &surftens, &gb_radius); |
312 | if (nread < 8) |
313 | { |
314 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 314); |
315 | return; |
316 | } |
317 | } |
318 | else |
319 | { |
320 | /* have_atomic_number && !have_bonded_type */ |
321 | nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf", |
322 | type, &atomnr, &m, &q, ptype, &c[0], &c[1], |
323 | &radius, &vol, &surftens, &gb_radius); |
324 | if (nread < 7) |
325 | { |
326 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 326); |
327 | return; |
328 | } |
329 | } |
330 | } |
331 | else |
332 | { |
333 | if (have_bonded_type) |
334 | { |
335 | /* !have_atomic_number && have_bonded_type */ |
336 | nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf", |
337 | type, btype, &m, &q, ptype, &c[0], &c[1], |
338 | &radius, &vol, &surftens, &gb_radius); |
339 | if (nread < 7) |
340 | { |
341 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 341); |
342 | return; |
343 | } |
344 | } |
345 | else |
346 | { |
347 | /* !have_atomic_number && !have_bonded_type */ |
348 | nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf", |
349 | type, &m, &q, ptype, &c[0], &c[1], |
350 | &radius, &vol, &surftens, &gb_radius); |
351 | if (nread < 6) |
352 | { |
353 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 353); |
354 | return; |
355 | } |
356 | } |
357 | } |
358 | |
359 | if (!have_bonded_type) |
360 | { |
361 | strcpy(btype, type); |
362 | } |
363 | |
364 | if (!have_atomic_number) |
365 | { |
366 | atomnr = -1; |
367 | } |
368 | |
369 | break; |
370 | |
371 | case F_BHAM: |
372 | nfp0 = 3; |
373 | |
374 | if (have_atomic_number) |
375 | { |
376 | if (have_bonded_type) |
377 | { |
378 | nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", |
379 | type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2], |
380 | &radius, &vol, &surftens, &gb_radius); |
381 | if (nread < 9) |
382 | { |
383 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 383); |
384 | return; |
385 | } |
386 | } |
387 | else |
388 | { |
389 | /* have_atomic_number && !have_bonded_type */ |
390 | nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", |
391 | type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2], |
392 | &radius, &vol, &surftens, &gb_radius); |
393 | if (nread < 8) |
394 | { |
395 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 395); |
396 | return; |
397 | } |
398 | } |
399 | } |
400 | else |
401 | { |
402 | if (have_bonded_type) |
403 | { |
404 | /* !have_atomic_number && have_bonded_type */ |
405 | nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", |
406 | type, btype, &m, &q, ptype, &c[0], &c[1], &c[2], |
407 | &radius, &vol, &surftens, &gb_radius); |
408 | if (nread < 8) |
409 | { |
410 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 410); |
411 | return; |
412 | } |
413 | } |
414 | else |
415 | { |
416 | /* !have_atomic_number && !have_bonded_type */ |
417 | nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", |
418 | type, &m, &q, ptype, &c[0], &c[1], &c[2], |
419 | &radius, &vol, &surftens, &gb_radius); |
420 | if (nread < 7) |
421 | { |
422 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 422); |
423 | return; |
424 | } |
425 | } |
426 | } |
427 | |
428 | if (!have_bonded_type) |
429 | { |
430 | strcpy(btype, type); |
431 | } |
432 | |
433 | if (!have_atomic_number) |
434 | { |
435 | atomnr = -1; |
436 | } |
437 | |
438 | break; |
439 | |
440 | default: |
441 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 441, "Invalid function type %d in push_at %s %d", nb_funct, |
442 | __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c", __LINE__442); |
443 | } |
444 | for (j = nfp0; (j < MAXFORCEPARAM12); j++) |
445 | { |
446 | c[j] = 0.0; |
447 | } |
448 | |
449 | if (strlen(type) == 1 && isdigit(type[0])((*__ctype_b_loc ())[(int) ((type[0]))] & (unsigned short int) _ISdigit)) |
450 | { |
451 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 451, "Atom type names can't be single digits."); |
452 | } |
453 | |
454 | if (strlen(btype) == 1 && isdigit(btype[0])((*__ctype_b_loc ())[(int) ((btype[0]))] & (unsigned short int) _ISdigit)) |
455 | { |
456 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 456, "Bond atom type names can't be single digits."); |
457 | } |
458 | |
459 | /* Hack to read old topologies */ |
460 | if (gmx_strcasecmp(ptype, "D") == 0) |
461 | { |
462 | sprintf(ptype, "V"); |
463 | } |
464 | for (j = 0; (j < eptNR); j++) |
465 | { |
466 | if (gmx_strcasecmp(ptype, xl[j].entry) == 0) |
467 | { |
468 | break; |
469 | } |
470 | } |
471 | if (j == eptNR) |
472 | { |
473 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 473, "Invalid particle type %s on line %s", |
474 | ptype, line); |
475 | } |
476 | pt = xl[j].ptype; |
477 | if (debug) |
478 | { |
479 | fprintf(debug, "ptype: %s\n", ptype_str[pt]); |
480 | } |
481 | |
482 | atom->q = q; |
483 | atom->m = m; |
484 | atom->ptype = pt; |
485 | for (i = 0; (i < MAXFORCEPARAM12); i++) |
486 | { |
487 | param->c[i] = c[i]; |
488 | } |
489 | |
490 | if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET-12345) |
491 | { |
492 | add_bond_atomtype(bat, symtab, btype); |
493 | } |
494 | batype_nr = get_bond_atomtype_type(btype, bat); |
495 | |
496 | if ((nr = get_atomtype_type(type, at)) != NOTSET-12345) |
497 | { |
498 | sprintf(errbuf, "Overriding atomtype %s", type); |
499 | warning(wi, errbuf); |
500 | if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr, |
501 | radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET-12345) |
502 | { |
503 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 503, "Replacing atomtype %s failed", type); |
504 | } |
505 | } |
506 | else if ((nr = add_atomtype(at, symtab, atom, type, param, |
507 | batype_nr, radius, vol, |
508 | surftens, atomnr, gb_radius, S_hct)) == NOTSET-12345) |
509 | { |
510 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 510, "Adding atomtype %s failed", type); |
511 | } |
512 | else |
513 | { |
514 | /* Add space in the non-bonded parameters matrix */ |
515 | realloc_nb_params(at, nbparam, pair); |
516 | } |
517 | sfree(atom)save_free("atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 517, (atom)); |
518 | sfree(param)save_free("param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 518, (param)); |
519 | } |
520 | |
521 | static void push_bondtype(t_params * bt, |
522 | t_param * b, |
523 | int nral, |
524 | int ftype, |
525 | gmx_bool bAllowRepeat, |
526 | char * line, |
527 | warninp_t wi) |
528 | { |
529 | int i, j; |
530 | gmx_bool bTest, bFound, bCont, bId; |
531 | int nr = bt->nr; |
532 | int nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); |
533 | char errbuf[256]; |
534 | |
535 | /* If bAllowRepeat is TRUE, we allow multiple entries as long as they |
536 | are on directly _adjacent_ lines. |
537 | */ |
538 | |
539 | /* First check if our atomtypes are _identical_ (not reversed) to the previous |
540 | entry. If they are not identical we search for earlier duplicates. If they are |
541 | we can skip it, since we already searched for the first line |
542 | in this group. |
543 | */ |
544 | |
545 | bFound = FALSE0; |
546 | bCont = FALSE0; |
547 | |
548 | if (bAllowRepeat && nr > 1) |
549 | { |
550 | for (j = 0, bCont = TRUE1; (j < nral); j++) |
551 | { |
552 | bCont = bCont && (b->a[j] == bt->param[nr-2].a[j]); |
553 | } |
554 | } |
555 | |
556 | /* Search for earlier duplicates if this entry was not a continuation |
557 | from the previous line. |
558 | */ |
559 | if (!bCont) |
560 | { |
561 | bFound = FALSE0; |
562 | for (i = 0; (i < nr); i++) |
563 | { |
564 | bTest = TRUE1; |
565 | for (j = 0; (j < nral); j++) |
566 | { |
567 | bTest = (bTest && (b->a[j] == bt->param[i].a[j])); |
568 | } |
569 | if (!bTest) |
570 | { |
571 | bTest = TRUE1; |
572 | for (j = 0; (j < nral); j++) |
573 | { |
574 | bTest = (bTest && (b->a[nral-1-j] == bt->param[i].a[j])); |
575 | } |
576 | } |
577 | if (bTest) |
578 | { |
579 | if (!bFound) |
580 | { |
581 | bId = TRUE1; |
582 | for (j = 0; (j < nrfp); j++) |
583 | { |
584 | bId = bId && (bt->param[i].c[j] == b->c[j]); |
585 | } |
586 | if (!bId) |
587 | { |
588 | sprintf(errbuf, "Overriding %s parameters.%s", |
589 | interaction_function[ftype].longname, |
590 | (ftype == F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : ""); |
591 | warning(wi, errbuf); |
592 | fprintf(stderrstderr, " old:"); |
593 | for (j = 0; (j < nrfp); j++) |
594 | { |
595 | fprintf(stderrstderr, " %g", bt->param[i].c[j]); |
596 | } |
597 | fprintf(stderrstderr, " \n new: %s\n\n", line); |
598 | } |
599 | } |
600 | /* Overwrite it! */ |
601 | for (j = 0; (j < nrfp); j++) |
602 | { |
603 | bt->param[i].c[j] = b->c[j]; |
604 | } |
605 | bFound = TRUE1; |
606 | } |
607 | } |
608 | } |
609 | if (!bFound) |
610 | { |
611 | /* alloc */ |
612 | pr_alloc (2, bt); |
613 | |
614 | /* fill the arrays up and down */ |
615 | memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c)); |
616 | memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a)); |
617 | memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c)); |
618 | |
619 | /* The definitions of linear angles depend on the order of atoms, |
620 | * that means that for atoms i-j-k, with certain parameter a, the |
621 | * corresponding k-j-i angle will have parameter 1-a. |
622 | */ |
623 | if (ftype == F_LINEAR_ANGLES) |
624 | { |
625 | bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0]; |
626 | bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2]; |
627 | } |
628 | |
629 | for (j = 0; (j < nral); j++) |
630 | { |
631 | bt->param[bt->nr+1].a[j] = b->a[nral-1-j]; |
632 | } |
633 | |
634 | bt->nr += 2; |
635 | } |
636 | } |
637 | |
638 | void push_bt(directive d, t_params bt[], int nral, |
639 | gpp_atomtype_t at, |
640 | t_bond_atomtype bat, char *line, |
641 | warninp_t wi) |
642 | { |
643 | const char *formal[MAXATOMLIST6+1] = { |
644 | "%s", |
645 | "%s%s", |
646 | "%s%s%s", |
647 | "%s%s%s%s", |
648 | "%s%s%s%s%s", |
649 | "%s%s%s%s%s%s", |
650 | "%s%s%s%s%s%s%s" |
651 | }; |
652 | const char *formnl[MAXATOMLIST6+1] = { |
653 | "%*s", |
654 | "%*s%*s", |
655 | "%*s%*s%*s", |
656 | "%*s%*s%*s%*s", |
657 | "%*s%*s%*s%*s%*s", |
658 | "%*s%*s%*s%*s%*s%*s", |
659 | "%*s%*s%*s%*s%*s%*s%*s" |
660 | }; |
661 | const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; |
662 | int i, ft, ftype, nn, nrfp, nrfpA, nrfpB; |
663 | char f1[STRLEN4096]; |
664 | char alc[MAXATOMLIST6+1][20]; |
665 | /* One force parameter more, so we can check if we read too many */ |
666 | double c[MAXFORCEPARAM12+1]; |
667 | t_param p; |
668 | char errbuf[256]; |
669 | |
670 | if ((bat && at) || (!bat && !at)) |
671 | { |
672 | gmx_incons("You should pass either bat or at to push_bt")_gmx_error("incons", "You should pass either bat or at to push_bt" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 672); |
673 | } |
674 | |
675 | /* Make format string (nral ints+functype) */ |
676 | if ((nn = sscanf(line, formal[nral], |
677 | alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1) |
678 | { |
679 | sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral); |
680 | warning_error(wi, errbuf); |
681 | return; |
682 | } |
683 | |
684 | ft = strtol(alc[nral], NULL((void*)0), 10); |
685 | ftype = ifunc_index(d, ft); |
686 | nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); |
687 | nrfpA = interaction_function[ftype].nrfpA; |
688 | nrfpB = interaction_function[ftype].nrfpB; |
689 | strcpy(f1, formnl[nral]); |
690 | strcat(f1, formlf); |
691 | if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11], &c[12])) |
692 | != nrfp) |
693 | { |
694 | if (nn == nrfpA) |
695 | { |
696 | /* Copy the B-state from the A-state */ |
697 | copy_B_from_A(ftype, c); |
698 | } |
699 | else |
700 | { |
701 | if (nn < nrfpA) |
702 | { |
703 | warning_error(wi, "Not enough parameters"); |
704 | } |
705 | else if (nn > nrfpA && nn < nrfp) |
706 | { |
707 | warning_error(wi, "Too many parameters or not enough parameters for topology B"); |
708 | } |
709 | else if (nn > nrfp) |
710 | { |
711 | warning_error(wi, "Too many parameters"); |
712 | } |
713 | for (i = nn; (i < nrfp); i++) |
714 | { |
715 | c[i] = 0.0; |
716 | } |
717 | } |
718 | } |
719 | for (i = 0; (i < nral); i++) |
720 | { |
721 | if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET-12345)) |
722 | { |
723 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 723, "Unknown atomtype %s\n", alc[i]); |
724 | } |
725 | else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET-12345)) |
726 | { |
727 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 727, "Unknown bond_atomtype %s\n", alc[i]); |
728 | } |
729 | } |
730 | for (i = 0; (i < nrfp); i++) |
731 | { |
732 | p.c[i] = c[i]; |
733 | } |
734 | push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE0, line, wi); |
735 | } |
736 | |
737 | |
738 | void push_dihedraltype(directive d, t_params bt[], |
739 | t_bond_atomtype bat, char *line, |
740 | warninp_t wi) |
741 | { |
742 | const char *formal[MAXATOMLIST6+1] = { |
743 | "%s", |
744 | "%s%s", |
745 | "%s%s%s", |
746 | "%s%s%s%s", |
747 | "%s%s%s%s%s", |
748 | "%s%s%s%s%s%s", |
749 | "%s%s%s%s%s%s%s" |
750 | }; |
751 | const char *formnl[MAXATOMLIST6+1] = { |
752 | "%*s", |
753 | "%*s%*s", |
754 | "%*s%*s%*s", |
755 | "%*s%*s%*s%*s", |
756 | "%*s%*s%*s%*s%*s", |
757 | "%*s%*s%*s%*s%*s%*s", |
758 | "%*s%*s%*s%*s%*s%*s%*s" |
759 | }; |
760 | const char *formlf[MAXFORCEPARAM12] = { |
761 | "%lf", |
762 | "%lf%lf", |
763 | "%lf%lf%lf", |
764 | "%lf%lf%lf%lf", |
765 | "%lf%lf%lf%lf%lf", |
766 | "%lf%lf%lf%lf%lf%lf", |
767 | "%lf%lf%lf%lf%lf%lf%lf", |
768 | "%lf%lf%lf%lf%lf%lf%lf%lf", |
769 | "%lf%lf%lf%lf%lf%lf%lf%lf%lf", |
770 | "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", |
771 | "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", |
772 | "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", |
773 | }; |
774 | int i, ft, ftype, nn, nrfp, nrfpA, nrfpB, nral; |
775 | char f1[STRLEN4096]; |
776 | char alc[MAXATOMLIST6+1][20]; |
777 | double c[MAXFORCEPARAM12]; |
778 | t_param p; |
779 | gmx_bool bAllowRepeat; |
780 | char errbuf[256]; |
781 | |
782 | /* This routine accepts dihedraltypes defined from either 2 or 4 atoms. |
783 | * |
784 | * We first check for 2 atoms with the 3th column being an integer |
785 | * defining the type. If this isn't the case, we try it with 4 atoms |
786 | * and the 5th column defining the dihedral type. |
787 | */ |
788 | nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]); |
789 | if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0])((*__ctype_b_loc ())[(int) ((alc[2][0]))] & (unsigned short int) _ISdigit)) |
790 | { |
791 | nral = 2; |
792 | ft = strtol(alc[nral], NULL((void*)0), 10); |
793 | /* Move atom types around a bit and use 'X' for wildcard atoms |
794 | * to create a 4-atom dihedral definition with arbitrary atoms in |
795 | * position 1 and 4. |
796 | */ |
797 | if (alc[2][0] == '2') |
798 | { |
799 | /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */ |
800 | strcpy(alc[3], alc[1]); |
801 | sprintf(alc[2], "X"); |
802 | sprintf(alc[1], "X"); |
803 | /* alc[0] stays put */ |
804 | } |
805 | else |
806 | { |
807 | /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */ |
808 | sprintf(alc[3], "X"); |
809 | strcpy(alc[2], alc[1]); |
810 | strcpy(alc[1], alc[0]); |
811 | sprintf(alc[0], "X"); |
812 | } |
813 | } |
814 | else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0])((*__ctype_b_loc ())[(int) ((alc[4][0]))] & (unsigned short int) _ISdigit)) |
815 | { |
816 | nral = 4; |
817 | ft = strtol(alc[nral], NULL((void*)0), 10); |
818 | } |
819 | else |
820 | { |
821 | sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1); |
822 | warning_error(wi, errbuf); |
823 | return; |
824 | } |
825 | |
826 | if (ft == 9) |
827 | { |
828 | /* Previously, we have always overwritten parameters if e.g. a torsion |
829 | with the same atomtypes occurs on multiple lines. However, CHARMM and |
830 | some other force fields specify multiple dihedrals over some bonds, |
831 | including cosines with multiplicity 6 and somethimes even higher. |
832 | Thus, they cannot be represented with Ryckaert-Bellemans terms. |
833 | To add support for these force fields, Dihedral type 9 is identical to |
834 | normal proper dihedrals, but repeated entries are allowed. |
835 | */ |
836 | bAllowRepeat = TRUE1; |
837 | ft = 1; |
838 | } |
839 | else |
840 | { |
841 | bAllowRepeat = FALSE0; |
842 | } |
843 | |
844 | |
845 | ftype = ifunc_index(d, ft); |
846 | nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); |
847 | nrfpA = interaction_function[ftype].nrfpA; |
848 | nrfpB = interaction_function[ftype].nrfpB; |
849 | |
850 | strcpy(f1, formnl[nral]); |
851 | strcat(f1, formlf[nrfp-1]); |
852 | |
853 | /* Check number of parameters given */ |
854 | if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11])) |
855 | != nrfp) |
856 | { |
857 | if (nn == nrfpA) |
858 | { |
859 | /* Copy the B-state from the A-state */ |
860 | copy_B_from_A(ftype, c); |
861 | } |
862 | else |
863 | { |
864 | if (nn < nrfpA) |
865 | { |
866 | warning_error(wi, "Not enough parameters"); |
867 | } |
868 | else if (nn > nrfpA && nn < nrfp) |
869 | { |
870 | warning_error(wi, "Too many parameters or not enough parameters for topology B"); |
871 | } |
872 | else if (nn > nrfp) |
873 | { |
874 | warning_error(wi, "Too many parameters"); |
875 | } |
876 | for (i = nn; (i < nrfp); i++) |
877 | { |
878 | c[i] = 0.0; |
879 | } |
880 | } |
881 | } |
882 | |
883 | for (i = 0; (i < 4); i++) |
884 | { |
885 | if (!strcmp(alc[i], "X")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (alc[i]) && __builtin_constant_p ("X") && (__s1_len = strlen (alc[i]), __s2_len = strlen ("X"), (!((size_t)(const void *)((alc[i]) + 1) - (size_t)(const void *)(alc[i]) == 1) || __s1_len >= 4) && (!((size_t)(const void *)(("X" ) + 1) - (size_t)(const void *)("X") == 1) || __s2_len >= 4 )) ? __builtin_strcmp (alc[i], "X") : (__builtin_constant_p ( alc[i]) && ((size_t)(const void *)((alc[i]) + 1) - (size_t )(const void *)(alc[i]) == 1) && (__s1_len = strlen ( alc[i]), __s1_len < 4) ? (__builtin_constant_p ("X") && ((size_t)(const void *)(("X") + 1) - (size_t)(const void *)( "X") == 1) ? __builtin_strcmp (alc[i], "X") : (__extension__ ( { const unsigned char *__s2 = (const unsigned char *) (const char *) ("X"); int __result = (((const unsigned char *) (const char *) (alc[i]))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ( alc[i]))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ( alc[i]))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (alc [i]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("X") && ((size_t)(const void *)(("X") + 1) - (size_t )(const void *)("X") == 1) && (__s2_len = strlen ("X" ), __s2_len < 4) ? (__builtin_constant_p (alc[i]) && ((size_t)(const void *)((alc[i]) + 1) - (size_t)(const void * )(alc[i]) == 1) ? __builtin_strcmp (alc[i], "X") : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (alc[i]); int __result = (((const unsigned char *) ( const char *) ("X"))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ("X"))[1] - __s2[1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ("X"))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) ("X"))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (alc[i], "X")))); })) |
886 | { |
887 | p.a[i] = -1; |
888 | } |
889 | else |
890 | { |
891 | if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET-12345) |
892 | { |
893 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 893, "Unknown bond_atomtype %s", alc[i]); |
894 | } |
895 | } |
896 | } |
897 | for (i = 0; (i < nrfp); i++) |
898 | { |
899 | p.c[i] = c[i]; |
900 | } |
901 | /* Always use 4 atoms here, since we created two wildcard atoms |
902 | * if there wasn't of them 4 already. |
903 | */ |
904 | push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi); |
905 | } |
906 | |
907 | |
908 | void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype, |
909 | char *pline, int nb_funct, |
910 | warninp_t wi) |
911 | { |
912 | /* swap the atoms */ |
913 | const char *form2 = "%*s%*s%*s%lf%lf"; |
914 | const char *form3 = "%*s%*s%*s%lf%lf%lf"; |
915 | const char *form4 = "%*s%*s%*s%lf%lf%lf%lf"; |
916 | const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf"; |
917 | char a0[80], a1[80]; |
918 | int i, f, n, ftype, atnr, nrfp; |
919 | double c[4], dum; |
920 | real cr[4], sig6; |
921 | atom_id ai, aj; |
922 | t_nbparam *nbp; |
923 | gmx_bool bId; |
924 | char errbuf[256]; |
925 | |
926 | if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3) |
927 | { |
928 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 928); |
929 | return; |
930 | } |
931 | |
932 | ftype = ifunc_index(d, f); |
933 | |
934 | if (ftype != nb_funct) |
935 | { |
936 | sprintf(errbuf, "Trying to add %s while the default nonbond type is %s", |
937 | interaction_function[ftype].longname, |
938 | interaction_function[nb_funct].longname); |
939 | warning_error(wi, errbuf); |
940 | return; |
941 | } |
942 | |
943 | /* Get the force parameters */ |
944 | nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); |
945 | if (ftype == F_LJ14) |
946 | { |
947 | n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]); |
948 | if (n < 2) |
949 | { |
950 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 950); |
951 | return; |
952 | } |
953 | /* When the B topology parameters are not set, |
954 | * copy them from topology A |
955 | */ |
956 | assert(nrfp <= 4)((void) (0)); |
957 | for (i = n; i < nrfp; i++) |
958 | { |
959 | c[i] = c[i-2]; |
960 | } |
961 | } |
962 | else if (ftype == F_LJC14_Q) |
963 | { |
964 | n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum); |
965 | if (n != 4) |
966 | { |
967 | incorrect_n_param(wi)_incorrect_n_param(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 967); |
968 | return; |
969 | } |
970 | } |
971 | else if (nrfp == 2) |
972 | { |
973 | if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2) |
974 | { |
975 | incorrect_n_param(wi)_incorrect_n_param(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 975); |
976 | return; |
977 | } |
978 | } |
979 | else if (nrfp == 3) |
980 | { |
981 | if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3) |
982 | { |
983 | incorrect_n_param(wi)_incorrect_n_param(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 983); |
984 | return; |
985 | } |
986 | } |
987 | else |
988 | { |
989 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 989, "Number of force parameters for nonbonded interactions is %d" |
990 | " in file %s, line %d", nrfp, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c", __LINE__990); |
991 | } |
992 | for (i = 0; (i < nrfp); i++) |
993 | { |
994 | cr[i] = c[i]; |
995 | } |
996 | |
997 | /* Put the parameters in the matrix */ |
998 | if ((ai = get_atomtype_type (a0, atype)) == NOTSET-12345) |
999 | { |
1000 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1000, "Atomtype %s not found", a0); |
1001 | } |
1002 | if ((aj = get_atomtype_type (a1, atype)) == NOTSET-12345) |
1003 | { |
1004 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1004, "Atomtype %s not found", a1); |
1005 | } |
1006 | nbp = &(nbt[max(ai, aj)(((ai) > (aj)) ? (ai) : (aj) )][min(ai, aj)(((ai) < (aj)) ? (ai) : (aj) )]); |
1007 | |
1008 | if (nbp->bSet) |
1009 | { |
1010 | bId = TRUE1; |
1011 | for (i = 0; i < nrfp; i++) |
1012 | { |
1013 | bId = bId && (nbp->c[i] == cr[i]); |
1014 | } |
1015 | if (!bId) |
1016 | { |
1017 | sprintf(errbuf, "Overriding non-bonded parameters,"); |
1018 | warning(wi, errbuf); |
1019 | fprintf(stderrstderr, " old:"); |
1020 | for (i = 0; i < nrfp; i++) |
1021 | { |
1022 | fprintf(stderrstderr, " %g", nbp->c[i]); |
1023 | } |
1024 | fprintf(stderrstderr, " new\n%s\n", pline); |
1025 | } |
1026 | } |
1027 | nbp->bSet = TRUE1; |
1028 | for (i = 0; i < nrfp; i++) |
1029 | { |
1030 | nbp->c[i] = cr[i]; |
1031 | } |
1032 | } |
1033 | |
1034 | void |
1035 | push_gb_params (gpp_atomtype_t at, char *line, |
1036 | warninp_t wi) |
1037 | { |
1038 | int nfield; |
1039 | int atype; |
1040 | double radius, vol, surftens, gb_radius, S_hct; |
1041 | char atypename[STRLEN4096]; |
1042 | char errbuf[STRLEN4096]; |
1043 | |
1044 | if ( (nfield = sscanf(line, "%s%lf%lf%lf%lf%lf", atypename, &radius, &vol, &surftens, &gb_radius, &S_hct)) != 6) |
1045 | { |
1046 | sprintf(errbuf, "Too few gb parameters for type %s\n", atypename); |
1047 | warning(wi, errbuf); |
1048 | } |
1049 | |
1050 | /* Search for atomtype */ |
1051 | atype = get_atomtype_type(atypename, at); |
1052 | |
1053 | if (atype == NOTSET-12345) |
1054 | { |
1055 | printf("Couldn't find topology match for atomtype %s\n", atypename); |
1056 | abort(); |
1057 | } |
1058 | |
1059 | set_atomtype_gbparam(at, atype, radius, vol, surftens, gb_radius, S_hct); |
1060 | } |
1061 | |
1062 | void |
1063 | push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at, |
1064 | t_bond_atomtype bat, char *line, |
1065 | warninp_t wi) |
1066 | { |
1067 | const char *formal = "%s%s%s%s%s%s%s%s"; |
1068 | |
1069 | int i, j, ft, ftype, nn, nrfp, nrfpA, nrfpB; |
1070 | int start; |
1071 | int nxcmap, nycmap, ncmap, read_cmap, sl, nct; |
1072 | char s[20], alc[MAXATOMLIST6+2][20]; |
1073 | t_param p; |
1074 | gmx_bool bAllowRepeat; |
1075 | char errbuf[256]; |
1076 | |
1077 | /* Keep the compiler happy */ |
1078 | read_cmap = 0; |
1079 | start = 0; |
1080 | |
1081 | if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7])) != nral+3) |
1082 | { |
1083 | sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1); |
1084 | warning_error(wi, errbuf); |
1085 | return; |
1086 | } |
1087 | |
1088 | /* Compute an offset for each line where the cmap parameters start |
1089 | * ie. where the atom types and grid spacing information ends |
1090 | */ |
1091 | for (i = 0; i < nn; i++) |
1092 | { |
1093 | start += (int)strlen(alc[i]); |
1094 | } |
1095 | |
1096 | /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */ |
1097 | /* start is the position on the line where we start to read the actual cmap grid data from the itp file */ |
1098 | start = start + nn -1; |
1099 | |
1100 | ft = strtol(alc[nral], NULL((void*)0), 10); |
1101 | nxcmap = strtol(alc[nral+1], NULL((void*)0), 10); |
1102 | nycmap = strtol(alc[nral+2], NULL((void*)0), 10); |
1103 | |
1104 | /* Check for equal grid spacing in x and y dims */ |
1105 | if (nxcmap != nycmap) |
1106 | { |
1107 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1107, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap); |
1108 | } |
1109 | |
1110 | ncmap = nxcmap*nycmap; |
1111 | ftype = ifunc_index(d, ft); |
1112 | nrfpA = strtol(alc[6], NULL((void*)0), 10)*strtol(alc[6], NULL((void*)0), 10); |
1113 | nrfpB = strtol(alc[7], NULL((void*)0), 10)*strtol(alc[7], NULL((void*)0), 10); |
1114 | nrfp = nrfpA+nrfpB; |
1115 | |
1116 | /* Allocate memory for the CMAP grid */ |
1117 | bt->ncmap += nrfp; |
1118 | srenew(bt->cmap, bt->ncmap)(bt->cmap) = save_realloc("bt->cmap", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1118, (bt->cmap), (bt->ncmap), sizeof(*(bt->cmap)) ); |
1119 | |
1120 | /* Read in CMAP parameters */ |
1121 | sl = 0; |
1122 | for (i = 0; i < ncmap; i++) |
1123 | { |
1124 | while (isspace(*(line+start+sl))((*__ctype_b_loc ())[(int) ((*(line+start+sl)))] & (unsigned short int) _ISspace)) |
1125 | { |
1126 | sl++; |
1127 | } |
1128 | nn = sscanf(line+start+sl, " %s ", s); |
1129 | sl += strlen(s); |
1130 | bt->cmap[i+(bt->ncmap)-nrfp] = strtod(s, NULL((void*)0)); |
1131 | |
1132 | if (nn == 1) |
1133 | { |
1134 | read_cmap++; |
1135 | } |
1136 | else |
1137 | { |
1138 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1138, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]); |
1139 | } |
1140 | |
1141 | } |
1142 | |
1143 | /* Check do that we got the number of parameters we expected */ |
1144 | if (read_cmap == nrfpA) |
1145 | { |
1146 | for (i = 0; i < ncmap; i++) |
1147 | { |
1148 | bt->cmap[i+ncmap] = bt->cmap[i]; |
1149 | } |
1150 | } |
1151 | else |
1152 | { |
1153 | if (read_cmap < nrfpA) |
1154 | { |
1155 | warning_error(wi, "Not enough cmap parameters"); |
1156 | } |
1157 | else if (read_cmap > nrfpA && read_cmap < nrfp) |
1158 | { |
1159 | warning_error(wi, "Too many cmap parameters or not enough parameters for topology B"); |
1160 | } |
1161 | else if (read_cmap > nrfp) |
1162 | { |
1163 | warning_error(wi, "Too many cmap parameters"); |
1164 | } |
1165 | } |
1166 | |
1167 | |
1168 | /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids |
1169 | * so we can safely assign them each time |
1170 | */ |
1171 | bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */ |
1172 | bt->nc = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */ |
1173 | nct = (nral+1) * bt->nc; |
1174 | |
1175 | /* Allocate memory for the cmap_types information */ |
1176 | srenew(bt->cmap_types, nct)(bt->cmap_types) = save_realloc("bt->cmap_types", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1176, (bt->cmap_types), (nct), sizeof(*(bt->cmap_types ))); |
1177 | |
1178 | for (i = 0; (i < nral); i++) |
1179 | { |
1180 | if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET-12345)) |
1181 | { |
1182 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1182, "Unknown atomtype %s\n", alc[i]); |
1183 | } |
1184 | else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET-12345)) |
1185 | { |
1186 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1186, "Unknown bond_atomtype %s\n", alc[i]); |
1187 | } |
1188 | |
1189 | /* Assign a grid number to each cmap_type */ |
1190 | bt->cmap_types[bt->nct++] = get_bond_atomtype_type(alc[i], bat); |
1191 | } |
1192 | |
1193 | /* Assign a type number to this cmap */ |
1194 | bt->cmap_types[bt->nct++] = bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */ |
1195 | |
1196 | /* Check for the correct number of atoms (again) */ |
1197 | if (bt->nct != nct) |
1198 | { |
1199 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1199, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt->nc); |
1200 | } |
1201 | |
1202 | /* Is this correct?? */ |
1203 | for (i = 0; i < MAXFORCEPARAM12; i++) |
1204 | { |
1205 | p.c[i] = NOTSET-12345; |
1206 | } |
1207 | |
1208 | /* Push the bond to the bondlist */ |
1209 | push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE0, line, wi); |
1210 | } |
1211 | |
1212 | |
1213 | static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr, |
1214 | int atomicnumber, |
1215 | int type, char *ctype, int ptype, |
1216 | char *resnumberic, |
1217 | char *resname, char *name, real m0, real q0, |
1218 | int typeB, char *ctypeB, real mB, real qB) |
1219 | { |
1220 | int j, resind = 0, resnr; |
1221 | unsigned char ric; |
1222 | int nr = at->nr; |
1223 | |
1224 | if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1))) |
1225 | { |
1226 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1226, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr); |
1227 | } |
1228 | |
1229 | j = strlen(resnumberic) - 1; |
1230 | if (isdigit(resnumberic[j])((*__ctype_b_loc ())[(int) ((resnumberic[j]))] & (unsigned short int) _ISdigit)) |
1231 | { |
1232 | ric = ' '; |
1233 | } |
1234 | else |
1235 | { |
1236 | ric = resnumberic[j]; |
1237 | if (j == 0 || !isdigit(resnumberic[j-1])((*__ctype_b_loc ())[(int) ((resnumberic[j-1]))] & (unsigned short int) _ISdigit)) |
1238 | { |
1239 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1239, "Invalid residue number '%s' for atom %d", |
1240 | resnumberic, atomnr); |
1241 | } |
1242 | } |
1243 | resnr = strtol(resnumberic, NULL((void*)0), 10); |
1244 | |
1245 | if (nr > 0) |
1246 | { |
1247 | resind = at->atom[nr-1].resind; |
1248 | } |
1249 | if (nr == 0 || strcmp(resname, *at->resinfo[resind].name)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (resname) && __builtin_constant_p (*at->resinfo[resind ].name) && (__s1_len = strlen (resname), __s2_len = strlen (*at->resinfo[resind].name), (!((size_t)(const void *)((resname ) + 1) - (size_t)(const void *)(resname) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((*at->resinfo[resind ].name) + 1) - (size_t)(const void *)(*at->resinfo[resind] .name) == 1) || __s2_len >= 4)) ? __builtin_strcmp (resname , *at->resinfo[resind].name) : (__builtin_constant_p (resname ) && ((size_t)(const void *)((resname) + 1) - (size_t )(const void *)(resname) == 1) && (__s1_len = strlen ( resname), __s1_len < 4) ? (__builtin_constant_p (*at->resinfo [resind].name) && ((size_t)(const void *)((*at->resinfo [resind].name) + 1) - (size_t)(const void *)(*at->resinfo[ resind].name) == 1) ? __builtin_strcmp (resname, *at->resinfo [resind].name) : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (*at->resinfo[resind ].name); int __result = (((const unsigned char *) (const char *) (resname))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ( resname))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ( resname))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (resname ))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ( *at->resinfo[resind].name) && ((size_t)(const void *)((*at->resinfo[resind].name) + 1) - (size_t)(const void *)(*at->resinfo[resind].name) == 1) && (__s2_len = strlen (*at->resinfo[resind].name), __s2_len < 4) ? (__builtin_constant_p (resname) && ((size_t)(const void *)((resname) + 1) - (size_t)(const void *)(resname) == 1) ? __builtin_strcmp (resname , *at->resinfo[resind].name) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (resname ); int __result = (((const unsigned char *) (const char *) (* at->resinfo[resind].name))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (*at->resinfo[resind].name))[1] - __s2[ 1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (*at->resinfo[ resind].name))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (* at->resinfo[resind].name))[3] - __s2[3]); } } __result; }) ))) : __builtin_strcmp (resname, *at->resinfo[resind].name )))); }) != 0 || |
1250 | resnr != at->resinfo[resind].nr || |
1251 | ric != at->resinfo[resind].ic) |
1252 | { |
1253 | if (nr == 0) |
1254 | { |
1255 | resind = 0; |
1256 | } |
1257 | else |
1258 | { |
1259 | resind++; |
1260 | } |
1261 | at->nres = resind + 1; |
1262 | srenew(at->resinfo, at->nres)(at->resinfo) = save_realloc("at->resinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1262, (at->resinfo), (at->nres), sizeof(*(at->resinfo ))); |
1263 | at->resinfo[resind].name = put_symtab(symtab, resname); |
1264 | at->resinfo[resind].nr = resnr; |
1265 | at->resinfo[resind].ic = ric; |
1266 | } |
1267 | else |
1268 | { |
1269 | resind = at->atom[at->nr-1].resind; |
1270 | } |
1271 | |
1272 | /* New atom instance |
1273 | * get new space for arrays |
1274 | */ |
1275 | srenew(at->atom, nr+1)(at->atom) = save_realloc("at->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1275, (at->atom), (nr+1), sizeof(*(at->atom))); |
1276 | srenew(at->atomname, nr+1)(at->atomname) = save_realloc("at->atomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1276, (at->atomname), (nr+1), sizeof(*(at->atomname)) ); |
1277 | srenew(at->atomtype, nr+1)(at->atomtype) = save_realloc("at->atomtype", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1277, (at->atomtype), (nr+1), sizeof(*(at->atomtype)) ); |
1278 | srenew(at->atomtypeB, nr+1)(at->atomtypeB) = save_realloc("at->atomtypeB", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1278, (at->atomtypeB), (nr+1), sizeof(*(at->atomtypeB ))); |
1279 | |
1280 | /* fill the list */ |
1281 | at->atom[nr].type = type; |
1282 | at->atom[nr].ptype = ptype; |
1283 | at->atom[nr].q = q0; |
1284 | at->atom[nr].m = m0; |
1285 | at->atom[nr].typeB = typeB; |
1286 | at->atom[nr].qB = qB; |
1287 | at->atom[nr].mB = mB; |
1288 | |
1289 | at->atom[nr].resind = resind; |
1290 | at->atom[nr].atomnumber = atomicnumber; |
1291 | at->atomname[nr] = put_symtab(symtab, name); |
1292 | at->atomtype[nr] = put_symtab(symtab, ctype); |
1293 | at->atomtypeB[nr] = put_symtab(symtab, ctypeB); |
1294 | at->nr++; |
1295 | } |
1296 | |
1297 | void push_cg(t_block *block, int *lastindex, int index, int a) |
1298 | { |
1299 | if (debug) |
1300 | { |
1301 | fprintf (debug, "Index %d, Atom %d\n", index, a); |
1302 | } |
1303 | |
1304 | if (((block->nr) && (*lastindex != index)) || (!block->nr)) |
1305 | { |
1306 | /* add a new block */ |
1307 | block->nr++; |
1308 | srenew(block->index, block->nr+1)(block->index) = save_realloc("block->index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1308, (block->index), (block->nr+1), sizeof(*(block-> index))); |
1309 | } |
1310 | block->index[block->nr] = a + 1; |
1311 | *lastindex = index; |
1312 | } |
1313 | |
1314 | void push_atom(t_symtab *symtab, t_block *cgs, |
1315 | t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg, |
1316 | warninp_t wi) |
1317 | { |
1318 | int nr, ptype; |
1319 | int cgnumber, atomnr, type, typeB, nscan; |
1320 | char id[STRLEN4096], ctype[STRLEN4096], ctypeB[STRLEN4096], |
1321 | resnumberic[STRLEN4096], resname[STRLEN4096], name[STRLEN4096], check[STRLEN4096]; |
1322 | double m, q, mb, qb; |
1323 | real m0, q0, mB, qB; |
1324 | |
1325 | /* Make a shortcut for writing in this molecule */ |
1326 | nr = at->nr; |
1327 | |
1328 | /* Fixed parameters */ |
1329 | if (sscanf(line, "%s%s%s%s%s%d", |
1330 | id, ctype, resnumberic, resname, name, &cgnumber) != 6) |
1331 | { |
1332 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1332); |
1333 | return; |
1334 | } |
1335 | sscanf(id, "%d", &atomnr); |
1336 | if ((type = get_atomtype_type(ctype, atype)) == NOTSET-12345) |
1337 | { |
1338 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1338, "Atomtype %s not found", ctype); |
1339 | } |
1340 | ptype = get_atomtype_ptype(type, atype); |
1341 | |
1342 | /* Set default from type */ |
1343 | q0 = get_atomtype_qA(type, atype); |
1344 | m0 = get_atomtype_massA(type, atype); |
1345 | typeB = type; |
1346 | qB = q0; |
1347 | mB = m0; |
1348 | |
1349 | /* Optional parameters */ |
1350 | nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s", |
1351 | &q, &m, ctypeB, &qb, &mb, check); |
1352 | |
1353 | /* Nasty switch that falls thru all the way down! */ |
1354 | if (nscan > 0) |
1355 | { |
1356 | q0 = qB = q; |
1357 | if (nscan > 1) |
1358 | { |
1359 | m0 = mB = m; |
1360 | if (nscan > 2) |
1361 | { |
1362 | if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET-12345) |
1363 | { |
1364 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1364, "Atomtype %s not found", ctypeB); |
1365 | } |
1366 | qB = get_atomtype_qA(typeB, atype); |
1367 | mB = get_atomtype_massA(typeB, atype); |
1368 | if (nscan > 3) |
1369 | { |
1370 | qB = qb; |
1371 | if (nscan > 4) |
1372 | { |
1373 | mB = mb; |
1374 | if (nscan > 5) |
1375 | { |
1376 | warning_error(wi, "Too many parameters"); |
1377 | } |
1378 | } |
1379 | } |
1380 | } |
1381 | } |
1382 | } |
1383 | if (debug) |
1384 | { |
1385 | fprintf(debug, "mB=%g, qB=%g, typeB=%d\n", mB, qB, typeB); |
1386 | } |
1387 | |
1388 | push_cg(cgs, lastcg, cgnumber, nr); |
1389 | |
1390 | push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype), |
1391 | type, ctype, ptype, resnumberic, |
1392 | resname, name, m0, q0, typeB, |
1393 | typeB == type ? ctype : ctypeB, mB, qB); |
1394 | } |
1395 | |
1396 | void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line, |
1397 | warninp_t wi) |
1398 | { |
1399 | char type[STRLEN4096]; |
1400 | int nrexcl, i; |
1401 | t_molinfo *newmol; |
1402 | |
1403 | if ((sscanf(line, "%s%d", type, &nrexcl)) != 2) |
1404 | { |
1405 | warning_error(wi, "Expected a molecule type name and nrexcl"); |
1406 | } |
1407 | |
1408 | /* Test if this atomtype overwrites another */ |
1409 | i = 0; |
1410 | while (i < *nmol) |
1411 | { |
1412 | if (gmx_strcasecmp(*((*mol)[i].name), type) == 0) |
1413 | { |
1414 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1414, "moleculetype %s is redefined", type); |
1415 | } |
1416 | i++; |
1417 | } |
1418 | |
1419 | (*nmol)++; |
1420 | srenew(*mol, *nmol)(*mol) = save_realloc("*mol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1420, (*mol), (*nmol), sizeof(*(*mol))); |
1421 | newmol = &((*mol)[*nmol-1]); |
1422 | init_molinfo(newmol); |
1423 | |
1424 | /* Fill in the values */ |
1425 | newmol->name = put_symtab(symtab, type); |
1426 | newmol->nrexcl = nrexcl; |
1427 | newmol->excl_set = FALSE0; |
1428 | } |
1429 | |
1430 | static gmx_bool default_nb_params(int ftype, t_params bt[], t_atoms *at, |
1431 | t_param *p, int c_start, gmx_bool bB, gmx_bool bGenPairs) |
1432 | { |
1433 | int i, j, ti, tj, ntype; |
1434 | gmx_bool bFound; |
1435 | t_param *pi = NULL((void*)0); |
1436 | int nr = bt[ftype].nr; |
1437 | int nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); |
1438 | int nrfp = interaction_function[ftype].nrfpA; |
1439 | int nrfpB = interaction_function[ftype].nrfpB; |
1440 | |
1441 | if ((!bB && nrfp == 0) || (bB && nrfpB == 0)) |
1442 | { |
1443 | return TRUE1; |
1444 | } |
1445 | |
1446 | bFound = FALSE0; |
1447 | if (bGenPairs) |
1448 | { |
1449 | /* First test the generated-pair position to save |
1450 | * time when we have 1000*1000 entries for e.g. OPLS... |
1451 | */ |
1452 | ntype = sqrt(nr); |
1453 | if (bB) |
1454 | { |
1455 | ti = at->atom[p->a[0]].typeB; |
1456 | tj = at->atom[p->a[1]].typeB; |
1457 | } |
1458 | else |
1459 | { |
1460 | ti = at->atom[p->a[0]].type; |
1461 | tj = at->atom[p->a[1]].type; |
1462 | } |
1463 | pi = &(bt[ftype].param[ntype*ti+tj]); |
1464 | bFound = ((ti == pi->a[0]) && (tj == pi->a[1])); |
1465 | } |
1466 | |
1467 | /* Search explicitly if we didnt find it */ |
1468 | if (!bFound) |
1469 | { |
1470 | for (i = 0; ((i < nr) && !bFound); i++) |
1471 | { |
1472 | pi = &(bt[ftype].param[i]); |
1473 | if (bB) |
1474 | { |
1475 | for (j = 0; ((j < nral) && |
1476 | (at->atom[p->a[j]].typeB == pi->a[j])); j++) |
1477 | { |
1478 | ; |
1479 | } |
1480 | } |
1481 | else |
1482 | { |
1483 | for (j = 0; ((j < nral) && |
1484 | (at->atom[p->a[j]].type == pi->a[j])); j++) |
1485 | { |
1486 | ; |
1487 | } |
1488 | } |
1489 | bFound = (j == nral); |
1490 | } |
1491 | } |
1492 | |
1493 | if (bFound) |
1494 | { |
1495 | if (bB) |
1496 | { |
1497 | if (nrfp+nrfpB > MAXFORCEPARAM12) |
1498 | { |
1499 | gmx_incons("Too many force parameters")_gmx_error("incons", "Too many force parameters", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1499); |
1500 | } |
1501 | for (j = c_start; (j < nrfpB); j++) |
1502 | { |
1503 | p->c[nrfp+j] = pi->c[j]; |
1504 | } |
1505 | } |
1506 | else |
1507 | { |
1508 | for (j = c_start; (j < nrfp); j++) |
1509 | { |
1510 | p->c[j] = pi->c[j]; |
1511 | } |
1512 | } |
1513 | } |
1514 | else |
1515 | { |
1516 | for (j = c_start; (j < nrfp); j++) |
1517 | { |
1518 | p->c[j] = 0.0; |
1519 | } |
1520 | } |
1521 | return bFound; |
1522 | } |
1523 | |
1524 | static gmx_bool default_cmap_params(t_params bondtype[], |
1525 | t_atoms *at, gpp_atomtype_t atype, |
1526 | t_param *p, gmx_bool bB, |
1527 | int *cmap_type, int *nparam_def) |
1528 | { |
1529 | int i, j, nparam_found; |
1530 | int ct; |
1531 | gmx_bool bFound = FALSE0; |
1532 | |
1533 | nparam_found = 0; |
1534 | ct = 0; |
1535 | |
1536 | /* Match the current cmap angle against the list of cmap_types */ |
1537 | for (i = 0; i < bondtype->nct && !bFound; i += 6) |
1538 | { |
1539 | if (bB) |
1540 | { |
1541 | |
1542 | } |
1543 | else |
1544 | { |
1545 | if ( |
1546 | (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype->cmap_types[i]) && |
1547 | (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype->cmap_types[i+1]) && |
1548 | (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype->cmap_types[i+2]) && |
1549 | (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype->cmap_types[i+3]) && |
1550 | (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype->cmap_types[i+4])) |
1551 | { |
1552 | /* Found cmap torsion */ |
1553 | bFound = TRUE1; |
1554 | ct = bondtype->cmap_types[i+5]; |
1555 | nparam_found = 1; |
1556 | } |
1557 | } |
1558 | } |
1559 | |
1560 | /* If we did not find a matching type for this cmap torsion */ |
1561 | if (!bFound) |
1562 | { |
1563 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1563, "Unknown cmap torsion between atoms %d %d %d %d %d\n", |
1564 | p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1); |
1565 | } |
1566 | |
1567 | *nparam_def = nparam_found; |
1568 | *cmap_type = ct; |
1569 | |
1570 | return bFound; |
1571 | } |
1572 | |
1573 | static gmx_bool default_params(int ftype, t_params bt[], |
1574 | t_atoms *at, gpp_atomtype_t atype, |
1575 | t_param *p, gmx_bool bB, |
1576 | t_param **param_def, |
1577 | int *nparam_def) |
1578 | { |
1579 | int i, j, nparam_found; |
1580 | gmx_bool bFound, bSame; |
1581 | t_param *pi = NULL((void*)0); |
1582 | t_param *pj = NULL((void*)0); |
1583 | int nr = bt[ftype].nr; |
1584 | int nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); |
1585 | int nrfpA = interaction_function[ftype].nrfpA; |
1586 | int nrfpB = interaction_function[ftype].nrfpB; |
1587 | |
1588 | if ((!bB && nrfpA == 0) || (bB && nrfpB == 0)) |
1589 | { |
1590 | return TRUE1; |
1591 | } |
1592 | |
1593 | |
1594 | /* We allow wildcards now. The first type (with or without wildcards) that |
1595 | * fits is used, so you should probably put the wildcarded bondtypes |
1596 | * at the end of each section. |
1597 | */ |
1598 | bFound = FALSE0; |
1599 | nparam_found = 0; |
1600 | /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a |
1601 | * special case for this. Check for B state outside loop to speed it up. |
1602 | */ |
1603 | if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS) |
1604 | { |
1605 | if (bB) |
1606 | { |
1607 | for (i = 0; ((i < nr) && !bFound); i++) |
1608 | { |
1609 | pi = &(bt[ftype].param[i]); |
1610 | bFound = |
1611 | ( |
1612 | ((pi->AIa[0] == -1) || (get_atomtype_batype(at->atom[p->AIa[0]].typeB, atype) == pi->AIa[0])) && |
1613 | ((pi->AJa[1] == -1) || (get_atomtype_batype(at->atom[p->AJa[1]].typeB, atype) == pi->AJa[1])) && |
1614 | ((pi->AKa[2] == -1) || (get_atomtype_batype(at->atom[p->AKa[2]].typeB, atype) == pi->AKa[2])) && |
1615 | ((pi->ALa[3] == -1) || (get_atomtype_batype(at->atom[p->ALa[3]].typeB, atype) == pi->ALa[3])) |
1616 | ); |
1617 | } |
1618 | } |
1619 | else |
1620 | { |
1621 | /* State A */ |
1622 | for (i = 0; ((i < nr) && !bFound); i++) |
1623 | { |
1624 | pi = &(bt[ftype].param[i]); |
1625 | bFound = |
1626 | ( |
1627 | ((pi->AIa[0] == -1) || (get_atomtype_batype(at->atom[p->AIa[0]].type, atype) == pi->AIa[0])) && |
1628 | ((pi->AJa[1] == -1) || (get_atomtype_batype(at->atom[p->AJa[1]].type, atype) == pi->AJa[1])) && |
1629 | ((pi->AKa[2] == -1) || (get_atomtype_batype(at->atom[p->AKa[2]].type, atype) == pi->AKa[2])) && |
1630 | ((pi->ALa[3] == -1) || (get_atomtype_batype(at->atom[p->ALa[3]].type, atype) == pi->ALa[3])) |
1631 | ); |
1632 | } |
1633 | } |
1634 | /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm. |
1635 | * The rules in that case is that additional matches HAVE to be on adjacent lines! |
1636 | */ |
1637 | if (bFound == TRUE1) |
1638 | { |
1639 | nparam_found++; |
1640 | bSame = TRUE1; |
1641 | /* Continue from current i value */ |
1642 | for (j = i+1; j < nr && bSame; j += 2) |
1643 | { |
1644 | pj = &(bt[ftype].param[j]); |
1645 | bSame = (pi->AIa[0] == pj->AIa[0] && pi->AJa[1] == pj->AJa[1] && pi->AKa[2] == pj->AKa[2] && pi->ALa[3] == pj->ALa[3]); |
1646 | if (bSame) |
1647 | { |
1648 | nparam_found++; |
1649 | } |
1650 | /* nparam_found will be increased as long as the numbers match */ |
1651 | } |
1652 | } |
1653 | } |
1654 | else /* Not a dihedral */ |
1655 | { |
1656 | for (i = 0; ((i < nr) && !bFound); i++) |
1657 | { |
1658 | pi = &(bt[ftype].param[i]); |
1659 | if (bB) |
1660 | { |
1661 | for (j = 0; ((j < nral) && |
1662 | (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++) |
1663 | { |
1664 | ; |
1665 | } |
1666 | } |
1667 | else |
1668 | { |
1669 | for (j = 0; ((j < nral) && |
1670 | (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++) |
1671 | { |
1672 | ; |
1673 | } |
1674 | } |
1675 | bFound = (j == nral); |
1676 | } |
1677 | if (bFound) |
1678 | { |
1679 | nparam_found = 1; |
1680 | } |
1681 | } |
1682 | |
1683 | *param_def = pi; |
1684 | *nparam_def = nparam_found; |
1685 | |
1686 | return bFound; |
1687 | } |
1688 | |
1689 | |
1690 | |
1691 | void push_bond(directive d, t_params bondtype[], t_params bond[], |
1692 | t_atoms *at, gpp_atomtype_t atype, char *line, |
1693 | gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ, |
1694 | gmx_bool bZero, gmx_bool *bWarn_copy_A_B, |
1695 | warninp_t wi) |
1696 | { |
1697 | const char *aaformat[MAXATOMLIST6] = { |
1698 | "%d%d", |
1699 | "%d%d%d", |
1700 | "%d%d%d%d", |
1701 | "%d%d%d%d%d", |
1702 | "%d%d%d%d%d%d", |
1703 | "%d%d%d%d%d%d%d" |
1704 | }; |
1705 | const char *asformat[MAXATOMLIST6] = { |
1706 | "%*s%*s", |
1707 | "%*s%*s%*s", |
1708 | "%*s%*s%*s%*s", |
1709 | "%*s%*s%*s%*s%*s", |
1710 | "%*s%*s%*s%*s%*s%*s", |
1711 | "%*s%*s%*s%*s%*s%*s%*s" |
1712 | }; |
1713 | const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; |
1714 | int nr, i, j, nral, nral_fmt, nread, ftype; |
1715 | char format[STRLEN4096]; |
1716 | /* One force parameter more, so we can check if we read too many */ |
1717 | double cc[MAXFORCEPARAM12+1]; |
1718 | int aa[MAXATOMLIST6+1]; |
1719 | t_param param, paramB, *param_defA, *param_defB; |
1720 | gmx_bool bFoundA = FALSE0, bFoundB = FALSE0, bDef, bPert, bSwapParity = FALSE0; |
1721 | int nparam_defA, nparam_defB; |
1722 | char errbuf[256]; |
1723 | |
1724 | nparam_defA = nparam_defB = 0; |
1725 | |
1726 | ftype = ifunc_index(d, 1); |
1727 | nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); |
1728 | for (j = 0; j < MAXATOMLIST6; j++) |
1729 | { |
1730 | aa[j] = NOTSET-12345; |
1731 | } |
1732 | bDef = (NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)) > 0); |
1733 | |
1734 | if (ftype == F_SETTLE) |
1735 | { |
1736 | /* SETTLE acts on 3 atoms, but the topology format only specifies |
1737 | * the first atom (for historical reasons). |
1738 | */ |
1739 | nral_fmt = 1; |
1740 | } |
1741 | else |
1742 | { |
1743 | nral_fmt = nral; |
1744 | } |
1745 | |
1746 | nread = sscanf(line, aaformat[nral_fmt-1], |
1747 | &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]); |
1748 | |
1749 | if (ftype == F_SETTLE) |
1750 | { |
1751 | aa[3] = aa[1]; |
1752 | aa[1] = aa[0] + 1; |
1753 | aa[2] = aa[0] + 2; |
1754 | } |
1755 | |
1756 | if (nread < nral_fmt) |
1757 | { |
1758 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1758); |
1759 | return; |
1760 | } |
1761 | else if (nread > nral_fmt) |
1762 | { |
1763 | /* this is a hack to allow for virtual sites with swapped parity */ |
1764 | bSwapParity = (aa[nral] < 0); |
1765 | if (bSwapParity) |
1766 | { |
1767 | aa[nral] = -aa[nral]; |
1768 | } |
1769 | ftype = ifunc_index(d, aa[nral]); |
1770 | if (bSwapParity) |
1771 | { |
1772 | switch (ftype) |
1773 | { |
1774 | case F_VSITE3FAD: |
1775 | case F_VSITE3OUT: |
1776 | break; |
1777 | default: |
1778 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1778, "Negative function types only allowed for %s and %s", |
1779 | interaction_function[F_VSITE3FAD].longname, |
1780 | interaction_function[F_VSITE3OUT].longname); |
1781 | } |
1782 | } |
1783 | } |
1784 | |
1785 | |
1786 | /* Check for double atoms and atoms out of bounds */ |
1787 | for (i = 0; (i < nral); i++) |
1788 | { |
1789 | if (aa[i] < 1 || aa[i] > at->nr) |
1790 | { |
1791 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1791, "[ file %s, line %d ]:\n" |
1792 | "Atom index (%d) in %s out of bounds (1-%d).\n" |
1793 | "This probably means that you have inserted topology section \"%s\"\n" |
1794 | "in a part belonging to a different molecule than you intended to.\n" |
1795 | "In that case move the \"%s\" section to the right molecule.", |
1796 | get_warning_file(wi), get_warning_line(wi), |
1797 | aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d)); |
1798 | } |
1799 | for (j = i+1; (j < nral); j++) |
1800 | { |
1801 | if (aa[i] == aa[j]) |
1802 | { |
1803 | sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d)); |
1804 | warning(wi, errbuf); |
1805 | } |
1806 | } |
1807 | } |
1808 | |
1809 | /* default force parameters */ |
1810 | for (j = 0; (j < MAXATOMLIST6); j++) |
1811 | { |
1812 | param.a[j] = aa[j]-1; |
1813 | } |
1814 | for (j = 0; (j < MAXFORCEPARAM12); j++) |
1815 | { |
1816 | param.c[j] = 0.0; |
1817 | } |
1818 | |
1819 | /* Get force params for normal and free energy perturbation |
1820 | * studies, as determined by types! |
1821 | */ |
1822 | |
1823 | if (bBonded) |
1824 | { |
1825 | bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE0, ¶m_defA, &nparam_defA); |
1826 | if (bFoundA) |
1827 | { |
1828 | /* Copy the A-state and B-state default parameters. */ |
1829 | assert(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM)((void) (0)); |
1830 | for (j = 0; (j < NRFPA(ftype)(interaction_function[(ftype)].nrfpA)+NRFPB(ftype)(interaction_function[(ftype)].nrfpB)); j++) |
1831 | { |
1832 | param.c[j] = param_defA->c[j]; |
1833 | } |
1834 | } |
1835 | bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE1, ¶m_defB, &nparam_defB); |
1836 | if (bFoundB) |
1837 | { |
1838 | /* Copy only the B-state default parameters */ |
1839 | for (j = NRFPA(ftype)(interaction_function[(ftype)].nrfpA); (j < NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB))); j++) |
1840 | { |
1841 | param.c[j] = param_defB->c[j]; |
1842 | } |
1843 | } |
1844 | } |
1845 | else if (ftype == F_LJ14) |
1846 | { |
1847 | bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE0, bGenPairs); |
1848 | bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE1, bGenPairs); |
1849 | } |
1850 | else if (ftype == F_LJC14_Q) |
1851 | { |
1852 | param.c[0] = fudgeQQ; |
1853 | /* Fill in the A-state charges as default parameters */ |
1854 | param.c[1] = at->atom[param.a[0]].q; |
1855 | param.c[2] = at->atom[param.a[1]].q; |
1856 | /* The default LJ parameters are the standard 1-4 parameters */ |
1857 | bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE0, bGenPairs); |
1858 | bFoundB = TRUE1; |
1859 | } |
1860 | else if (ftype == F_LJC_PAIRS_NB) |
1861 | { |
1862 | /* Defaults are not supported here */ |
1863 | bFoundA = FALSE0; |
1864 | bFoundB = TRUE1; |
1865 | } |
1866 | else |
1867 | { |
1868 | gmx_incons("Unknown function type in push_bond")_gmx_error("incons", "Unknown function type in push_bond", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1868); |
1869 | } |
1870 | |
1871 | if (nread > nral_fmt) |
1872 | { |
1873 | /* Manually specified parameters - in this case we discard multiple torsion info! */ |
1874 | |
1875 | strcpy(format, asformat[nral_fmt-1]); |
1876 | strcat(format, ccformat); |
1877 | |
1878 | nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5], |
1879 | &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]); |
1880 | |
1881 | if ((nread == NRFPA(ftype)(interaction_function[(ftype)].nrfpA)) && (NRFPB(ftype)(interaction_function[(ftype)].nrfpB) != 0)) |
1882 | { |
1883 | /* We only have to issue a warning if these atoms are perturbed! */ |
1884 | bPert = FALSE0; |
1885 | for (j = 0; (j < nral); j++) |
1886 | { |
1887 | bPert = bPert || PERTURBED(at->atom[param.a[j]])(((at->atom[param.a[j]]).mB != (at->atom[param.a[j]]).m ) || ((at->atom[param.a[j]]).qB != (at->atom[param.a[j] ]).q) || ((at->atom[param.a[j]]).typeB != (at->atom[param .a[j]]).type)); |
1888 | } |
1889 | |
1890 | if (bPert && *bWarn_copy_A_B) |
1891 | { |
1892 | sprintf(errbuf, |
1893 | "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B"); |
1894 | warning(wi, errbuf); |
1895 | *bWarn_copy_A_B = FALSE0; |
1896 | } |
1897 | |
1898 | /* If only the A parameters were specified, copy them to the B state */ |
1899 | /* The B-state parameters correspond to the first nrfpB |
1900 | * A-state parameters. |
1901 | */ |
1902 | for (j = 0; (j < NRFPB(ftype)(interaction_function[(ftype)].nrfpB)); j++) |
1903 | { |
1904 | cc[nread++] = cc[j]; |
1905 | } |
1906 | } |
1907 | |
1908 | /* If nread was 0 or EOF, no parameters were read => use defaults. |
1909 | * If nread was nrfpA we copied above so nread=nrfp. |
1910 | * If nread was nrfp we are cool. |
1911 | * For F_LJC14_Q we allow supplying fudgeQQ only. |
1912 | * Anything else is an error! |
1913 | */ |
1914 | if ((nread != 0) && (nread != EOF(-1)) && (nread != NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB))) && |
1915 | !(ftype == F_LJC14_Q && nread == 1)) |
1916 | { |
1917 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 1917, "Incorrect number of parameters - found %d, expected %d or %d for %s.", |
1918 | nread, NRFPA(ftype)(interaction_function[(ftype)].nrfpA), NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)), |
1919 | interaction_function[ftype].longname); |
1920 | } |
1921 | |
1922 | for (j = 0; (j < nread); j++) |
1923 | { |
1924 | param.c[j] = cc[j]; |
1925 | } |
1926 | |
1927 | /* Check whether we have to use the defaults */ |
1928 | if (nread == NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB))) |
1929 | { |
1930 | bDef = FALSE0; |
1931 | } |
1932 | } |
1933 | else |
1934 | { |
1935 | nread = 0; |
1936 | } |
1937 | /* nread now holds the number of force parameters read! */ |
1938 | |
1939 | if (bDef) |
1940 | { |
1941 | /* Use defaults */ |
1942 | /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */ |
1943 | if (ftype == F_PDIHS) |
1944 | { |
1945 | if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB))) |
1946 | { |
1947 | sprintf(errbuf, |
1948 | "Cannot automatically perturb a torsion with multiple terms to different form.\n" |
1949 | "Please specify perturbed parameters manually for this torsion in your topology!"); |
1950 | warning_error(wi, errbuf); |
1951 | } |
1952 | } |
1953 | |
1954 | if (nread > 0 && nread < NRFPA(ftype)(interaction_function[(ftype)].nrfpA)) |
1955 | { |
1956 | /* Issue an error, do not use defaults */ |
1957 | sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype)(interaction_function[(ftype)].nrfpA)); |
1958 | warning_error(wi, errbuf); |
1959 | } |
1960 | |
1961 | if (nread == 0 || nread == EOF(-1)) |
1962 | { |
1963 | if (!bFoundA) |
1964 | { |
1965 | if (interaction_function[ftype].flags & IF_VSITE1<<1) |
1966 | { |
1967 | /* set them to NOTSET, will be calculated later */ |
1968 | for (j = 0; (j < MAXFORCEPARAM12); j++) |
1969 | { |
1970 | param.c[j] = NOTSET-12345; |
1971 | } |
1972 | |
1973 | if (bSwapParity) |
1974 | { |
1975 | param.C1c[1] = -1; /* flag to swap parity of vsite construction */ |
1976 | } |
1977 | } |
1978 | else |
1979 | { |
1980 | if (bZero) |
1981 | { |
1982 | fprintf(stderrstderr, "NOTE: No default %s types, using zeroes\n", |
1983 | interaction_function[ftype].longname); |
1984 | } |
1985 | else |
1986 | { |
1987 | sprintf(errbuf, "No default %s types", interaction_function[ftype].longname); |
1988 | warning_error(wi, errbuf); |
1989 | } |
1990 | } |
1991 | } |
1992 | else |
1993 | { |
1994 | if (bSwapParity) |
1995 | { |
1996 | switch (ftype) |
1997 | { |
1998 | case F_VSITE3FAD: |
1999 | param.C0c[0] = 360 - param.C0c[0]; |
2000 | break; |
2001 | case F_VSITE3OUT: |
2002 | param.C2c[2] = -param.C2c[2]; |
2003 | break; |
2004 | } |
2005 | } |
2006 | } |
2007 | if (!bFoundB) |
2008 | { |
2009 | /* We only have to issue a warning if these atoms are perturbed! */ |
2010 | bPert = FALSE0; |
2011 | for (j = 0; (j < nral); j++) |
2012 | { |
2013 | bPert = bPert || PERTURBED(at->atom[param.a[j]])(((at->atom[param.a[j]]).mB != (at->atom[param.a[j]]).m ) || ((at->atom[param.a[j]]).qB != (at->atom[param.a[j] ]).q) || ((at->atom[param.a[j]]).typeB != (at->atom[param .a[j]]).type)); |
2014 | } |
2015 | |
2016 | if (bPert) |
2017 | { |
2018 | sprintf(errbuf, "No default %s types for perturbed atoms, " |
2019 | "using normal values", interaction_function[ftype].longname); |
2020 | warning(wi, errbuf); |
2021 | } |
2022 | } |
2023 | } |
2024 | } |
2025 | |
2026 | if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) |
2027 | && param.c[5] != param.c[2]) |
2028 | { |
2029 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2029, "[ file %s, line %d ]:\n" |
2030 | " %s multiplicity can not be perturbed %f!=%f", |
2031 | get_warning_file(wi), get_warning_line(wi), |
2032 | interaction_function[ftype].longname, |
2033 | param.c[2], param.c[5]); |
2034 | } |
2035 | |
2036 | if (IS_TABULATED(ftype)(interaction_function[(ftype)].flags & 1<<6) && param.c[2] != param.c[0]) |
2037 | { |
2038 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2038, "[ file %s, line %d ]:\n" |
2039 | " %s table number can not be perturbed %d!=%d", |
2040 | get_warning_file(wi), get_warning_line(wi), |
2041 | interaction_function[ftype].longname, |
2042 | (int)(param.c[0]+0.5), (int)(param.c[2]+0.5)); |
2043 | } |
2044 | |
2045 | /* Dont add R-B dihedrals where all parameters are zero (no interaction) */ |
2046 | if (ftype == F_RBDIHS) |
2047 | { |
2048 | nr = 0; |
2049 | for (i = 0; i < NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[ (ftype)].nrfpB)); i++) |
2050 | { |
2051 | if (param.c[i] != 0) |
2052 | { |
2053 | nr++; |
2054 | } |
2055 | } |
2056 | if (nr == 0) |
2057 | { |
2058 | return; |
2059 | } |
2060 | } |
2061 | |
2062 | /* Put the values in the appropriate arrays */ |
2063 | add_param_to_list (&bond[ftype], ¶m); |
2064 | |
2065 | /* Push additional torsions from FF for ftype==9 if we have them. |
2066 | * We have already checked that the A/B states do not differ in this case, |
2067 | * so we do not have to double-check that again, or the vsite stuff. |
2068 | * In addition, those torsions cannot be automatically perturbed. |
2069 | */ |
2070 | if (bDef && ftype == F_PDIHS) |
2071 | { |
2072 | for (i = 1; i < nparam_defA; i++) |
2073 | { |
2074 | /* Advance pointer! */ |
2075 | param_defA += 2; |
2076 | for (j = 0; (j < NRFPA(ftype)(interaction_function[(ftype)].nrfpA)+NRFPB(ftype)(interaction_function[(ftype)].nrfpB)); j++) |
2077 | { |
2078 | param.c[j] = param_defA->c[j]; |
2079 | } |
2080 | /* And push the next term for this torsion */ |
2081 | add_param_to_list (&bond[ftype], ¶m); |
2082 | } |
2083 | } |
2084 | } |
2085 | |
2086 | void push_cmap(directive d, t_params bondtype[], t_params bond[], |
2087 | t_atoms *at, gpp_atomtype_t atype, char *line, |
2088 | warninp_t wi) |
2089 | { |
2090 | const char *aaformat[MAXATOMLIST6+1] = |
2091 | { |
2092 | "%d", |
2093 | "%d%d", |
2094 | "%d%d%d", |
2095 | "%d%d%d%d", |
2096 | "%d%d%d%d%d", |
2097 | "%d%d%d%d%d%d", |
2098 | "%d%d%d%d%d%d%d" |
2099 | }; |
2100 | |
2101 | int i, j, nr, ftype, nral, nread, ncmap_params; |
2102 | int cmap_type; |
2103 | int aa[MAXATOMLIST6+1]; |
2104 | char errbuf[256]; |
2105 | gmx_bool bFound; |
2106 | t_param param, paramB, *param_defA, *param_defB; |
2107 | |
2108 | ftype = ifunc_index(d, 1); |
2109 | nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); |
2110 | nr = bondtype[ftype].nr; |
Value stored to 'nr' is never read | |
2111 | ncmap_params = 0; |
2112 | |
2113 | nread = sscanf(line, aaformat[nral-1], |
2114 | &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]); |
2115 | |
2116 | if (nread < nral) |
2117 | { |
2118 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2118); |
2119 | return; |
2120 | } |
2121 | else if (nread == nral) |
2122 | { |
2123 | ftype = ifunc_index(d, 1); |
2124 | } |
2125 | |
2126 | /* Check for double atoms and atoms out of bounds */ |
2127 | for (i = 0; i < nral; i++) |
2128 | { |
2129 | if (aa[i] < 1 || aa[i] > at->nr) |
2130 | { |
2131 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2131, "[ file %s, line %d ]:\n" |
2132 | "Atom index (%d) in %s out of bounds (1-%d).\n" |
2133 | "This probably means that you have inserted topology section \"%s\"\n" |
2134 | "in a part belonging to a different molecule than you intended to.\n" |
2135 | "In that case move the \"%s\" section to the right molecule.", |
2136 | get_warning_file(wi), get_warning_line(wi), |
2137 | aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d)); |
2138 | } |
2139 | |
2140 | for (j = i+1; (j < nral); j++) |
2141 | { |
2142 | if (aa[i] == aa[j]) |
2143 | { |
2144 | sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d)); |
2145 | warning(wi, errbuf); |
2146 | } |
2147 | } |
2148 | } |
2149 | |
2150 | /* default force parameters */ |
2151 | for (j = 0; (j < MAXATOMLIST6); j++) |
2152 | { |
2153 | param.a[j] = aa[j]-1; |
2154 | } |
2155 | for (j = 0; (j < MAXFORCEPARAM12); j++) |
2156 | { |
2157 | param.c[j] = 0.0; |
2158 | } |
2159 | |
2160 | /* Get the cmap type for this cmap angle */ |
2161 | bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE0, &cmap_type, &ncmap_params); |
2162 | |
2163 | /* We want exactly one parameter (the cmap type in state A (currently no state B) back */ |
2164 | if (bFound && ncmap_params == 1) |
2165 | { |
2166 | /* Put the values in the appropriate arrays */ |
2167 | param.c[0] = cmap_type; |
2168 | add_param_to_list(&bond[ftype], ¶m); |
2169 | } |
2170 | else |
2171 | { |
2172 | /* This is essentially the same check as in default_cmap_params() done one more time */ |
2173 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2173, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n", |
2174 | param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1); |
2175 | } |
2176 | } |
2177 | |
2178 | |
2179 | |
2180 | void push_vsitesn(directive d, t_params bond[], |
2181 | t_atoms *at, char *line, |
2182 | warninp_t wi) |
2183 | { |
2184 | char *ptr; |
2185 | int type, ftype, j, n, ret, nj, a; |
2186 | int *atc = NULL((void*)0); |
2187 | double *weight = NULL((void*)0), weight_tot; |
2188 | t_param param; |
2189 | |
2190 | /* default force parameters */ |
2191 | for (j = 0; (j < MAXATOMLIST6); j++) |
2192 | { |
2193 | param.a[j] = NOTSET-12345; |
2194 | } |
2195 | for (j = 0; (j < MAXFORCEPARAM12); j++) |
2196 | { |
2197 | param.c[j] = 0.0; |
2198 | } |
2199 | |
2200 | ptr = line; |
2201 | ret = sscanf(ptr, "%d%n", &a, &n); |
2202 | ptr += n; |
2203 | if (ret == 0) |
2204 | { |
2205 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2205, "[ file %s, line %d ]:\n" |
2206 | " Expected an atom index in section \"%s\"", |
2207 | get_warning_file(wi), get_warning_line(wi), |
2208 | dir2str(d)); |
2209 | } |
2210 | |
2211 | param.a[0] = a - 1; |
2212 | |
2213 | ret = sscanf(ptr, "%d%n", &type, &n); |
2214 | ptr += n; |
2215 | ftype = ifunc_index(d, type); |
2216 | |
2217 | weight_tot = 0; |
2218 | nj = 0; |
2219 | do |
2220 | { |
2221 | ret = sscanf(ptr, "%d%n", &a, &n); |
2222 | ptr += n; |
2223 | if (ret > 0) |
2224 | { |
2225 | if (nj % 20 == 0) |
2226 | { |
2227 | srenew(atc, nj+20)(atc) = save_realloc("atc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2227, (atc), (nj+20), sizeof(*(atc))); |
2228 | srenew(weight, nj+20)(weight) = save_realloc("weight", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2228, (weight), (nj+20), sizeof(*(weight))); |
2229 | } |
2230 | atc[nj] = a - 1; |
2231 | switch (type) |
2232 | { |
2233 | case 1: |
2234 | weight[nj] = 1; |
2235 | break; |
2236 | case 2: |
2237 | /* Here we use the A-state mass as a parameter. |
2238 | * Note that the B-state mass has no influence. |
2239 | */ |
2240 | weight[nj] = at->atom[atc[nj]].m; |
2241 | break; |
2242 | case 3: |
2243 | weight[nj] = -1; |
2244 | ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n); |
2245 | ptr += n; |
2246 | if (weight[nj] < 0) |
2247 | { |
2248 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2248, "[ file %s, line %d ]:\n" |
2249 | " No weight or negative weight found for vsiten constructing atom %d (atom index %d)", |
2250 | get_warning_file(wi), get_warning_line(wi), |
2251 | nj+1, atc[nj]+1); |
2252 | } |
2253 | break; |
2254 | default: |
2255 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2255, "Unknown vsiten type %d", type); |
2256 | } |
2257 | weight_tot += weight[nj]; |
2258 | nj++; |
2259 | } |
2260 | } |
2261 | while (ret > 0); |
2262 | |
2263 | if (nj == 0) |
2264 | { |
2265 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2265, "[ file %s, line %d ]:\n" |
2266 | " Expected more than one atom index in section \"%s\"", |
2267 | get_warning_file(wi), get_warning_line(wi), |
2268 | dir2str(d)); |
2269 | } |
2270 | |
2271 | if (weight_tot == 0) |
2272 | { |
2273 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2273, "[ file %s, line %d ]:\n" |
2274 | " The total mass of the construting atoms is zero", |
2275 | get_warning_file(wi), get_warning_line(wi)); |
2276 | } |
2277 | |
2278 | for (j = 0; j < nj; j++) |
2279 | { |
2280 | param.a[1] = atc[j]; |
2281 | param.c[0] = nj; |
2282 | param.c[1] = weight[j]/weight_tot; |
2283 | /* Put the values in the appropriate arrays */ |
2284 | add_param_to_list (&bond[ftype], ¶m); |
2285 | } |
2286 | |
2287 | sfree(atc)save_free("atc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2287, (atc)); |
2288 | sfree(weight)save_free("weight", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2288, (weight)); |
2289 | } |
2290 | |
2291 | void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol, |
2292 | int *nrcopies, |
2293 | warninp_t wi) |
2294 | { |
2295 | int i, copies; |
2296 | char type[STRLEN4096]; |
2297 | |
2298 | *nrcopies = 0; |
2299 | if (sscanf(pline, "%s%d", type, &copies) != 2) |
2300 | { |
2301 | too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2301); |
2302 | return; |
2303 | } |
2304 | |
2305 | /* search moleculename */ |
2306 | for (i = 0; ((i < nrmols) && gmx_strcasecmp(type, *(mols[i].name))); i++) |
2307 | { |
2308 | ; |
2309 | } |
2310 | |
2311 | if (i < nrmols) |
2312 | { |
2313 | *nrcopies = copies; |
2314 | *whichmol = i; |
2315 | } |
2316 | else |
2317 | { |
2318 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2318, "No such moleculetype %s", type); |
2319 | } |
2320 | } |
2321 | |
2322 | void init_block2(t_block2 *b2, int natom) |
2323 | { |
2324 | int i; |
2325 | |
2326 | b2->nr = natom; |
2327 | snew(b2->nra, b2->nr)(b2->nra) = save_calloc("b2->nra", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2327, (b2->nr), sizeof(*(b2->nra))); |
2328 | snew(b2->a, b2->nr)(b2->a) = save_calloc("b2->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2328, (b2->nr), sizeof(*(b2->a))); |
2329 | for (i = 0; (i < b2->nr); i++) |
2330 | { |
2331 | b2->a[i] = NULL((void*)0); |
2332 | } |
2333 | } |
2334 | |
2335 | void done_block2(t_block2 *b2) |
2336 | { |
2337 | int i; |
2338 | |
2339 | if (b2->nr) |
2340 | { |
2341 | for (i = 0; (i < b2->nr); i++) |
2342 | { |
2343 | sfree(b2->a[i])save_free("b2->a[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2343, (b2->a[i])); |
2344 | } |
2345 | sfree(b2->a)save_free("b2->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2345, (b2->a)); |
2346 | sfree(b2->nra)save_free("b2->nra", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2346, (b2->nra)); |
2347 | b2->nr = 0; |
2348 | } |
2349 | } |
2350 | |
2351 | void push_excl(char *line, t_block2 *b2) |
2352 | { |
2353 | int i, j; |
2354 | int n; |
2355 | char base[STRLEN4096], format[STRLEN4096]; |
2356 | |
2357 | if (sscanf(line, "%d", &i) == 0) |
2358 | { |
2359 | return; |
2360 | } |
2361 | |
2362 | if ((1 <= i) && (i <= b2->nr)) |
2363 | { |
2364 | i--; |
2365 | } |
2366 | else |
2367 | { |
2368 | if (debug) |
2369 | { |
2370 | fprintf(debug, "Unbound atom %d\n", i-1); |
2371 | } |
2372 | return; |
2373 | } |
2374 | strcpy(base, "%*d"); |
2375 | do |
2376 | { |
2377 | strcpy(format, base); |
2378 | strcat(format, "%d"); |
2379 | n = sscanf(line, format, &j); |
2380 | if (n == 1) |
2381 | { |
2382 | if ((1 <= j) && (j <= b2->nr)) |
2383 | { |
2384 | j--; |
2385 | srenew(b2->a[i], ++(b2->nra[i]))(b2->a[i]) = save_realloc("b2->a[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2385, (b2->a[i]), (++(b2->nra[i])), sizeof(*(b2->a [i]))); |
2386 | b2->a[i][b2->nra[i]-1] = j; |
2387 | /* also add the reverse exclusion! */ |
2388 | srenew(b2->a[j], ++(b2->nra[j]))(b2->a[j]) = save_realloc("b2->a[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2388, (b2->a[j]), (++(b2->nra[j])), sizeof(*(b2->a [j]))); |
2389 | b2->a[j][b2->nra[j]-1] = i; |
2390 | strcat(base, "%*d"); |
2391 | } |
2392 | else |
2393 | { |
2394 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2394, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr); |
2395 | } |
2396 | } |
2397 | } |
2398 | while (n == 1); |
2399 | } |
2400 | |
2401 | void b_to_b2(t_blocka *b, t_block2 *b2) |
2402 | { |
2403 | int i; |
2404 | atom_id j, a; |
2405 | |
2406 | for (i = 0; (i < b->nr); i++) |
2407 | { |
2408 | for (j = b->index[i]; (j < b->index[i+1]); j++) |
2409 | { |
2410 | a = b->a[j]; |
2411 | srenew(b2->a[i], ++b2->nra[i])(b2->a[i]) = save_realloc("b2->a[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2411, (b2->a[i]), (++b2->nra[i]), sizeof(*(b2->a[i ]))); |
2412 | b2->a[i][b2->nra[i]-1] = a; |
2413 | } |
2414 | } |
2415 | } |
2416 | |
2417 | void b2_to_b(t_block2 *b2, t_blocka *b) |
2418 | { |
2419 | int i, nra; |
2420 | atom_id j; |
2421 | |
2422 | nra = 0; |
2423 | for (i = 0; (i < b2->nr); i++) |
2424 | { |
2425 | b->index[i] = nra; |
2426 | for (j = 0; (j < b2->nra[i]); j++) |
2427 | { |
2428 | b->a[nra+j] = b2->a[i][j]; |
2429 | } |
2430 | nra += b2->nra[i]; |
2431 | } |
2432 | /* terminate list */ |
2433 | b->index[i] = nra; |
2434 | } |
2435 | |
2436 | static int icomp(const void *v1, const void *v2) |
2437 | { |
2438 | return (*((atom_id *) v1))-(*((atom_id *) v2)); |
2439 | } |
2440 | |
2441 | void merge_excl(t_blocka *excl, t_block2 *b2) |
2442 | { |
2443 | int i, k; |
2444 | atom_id j; |
2445 | int nra; |
2446 | |
2447 | if (!b2->nr) |
2448 | { |
2449 | return; |
2450 | } |
2451 | else if (b2->nr != excl->nr) |
2452 | { |
2453 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2453, "DEATH HORROR: b2->nr = %d, while excl->nr = %d", |
2454 | b2->nr, excl->nr); |
2455 | } |
2456 | else if (debug) |
2457 | { |
2458 | fprintf(debug, "Entering merge_excl\n"); |
2459 | } |
2460 | |
2461 | /* First copy all entries from excl to b2 */ |
2462 | b_to_b2(excl, b2); |
2463 | |
2464 | /* Count and sort the exclusions */ |
2465 | nra = 0; |
2466 | for (i = 0; (i < b2->nr); i++) |
2467 | { |
2468 | if (b2->nra[i] > 0) |
2469 | { |
2470 | /* remove double entries */ |
2471 | qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp); |
2472 | k = 1; |
2473 | for (j = 1; (j < b2->nra[i]); j++) |
2474 | { |
2475 | if (b2->a[i][j] != b2->a[i][k-1]) |
2476 | { |
2477 | b2->a[i][k] = b2->a[i][j]; |
2478 | k++; |
2479 | } |
2480 | } |
2481 | b2->nra[i] = k; |
2482 | nra += b2->nra[i]; |
2483 | } |
2484 | } |
2485 | excl->nra = nra; |
2486 | srenew(excl->a, excl->nra)(excl->a) = save_realloc("excl->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2486, (excl->a), (excl->nra), sizeof(*(excl->a))); |
2487 | |
2488 | b2_to_b(b2, excl); |
2489 | } |
2490 | |
2491 | int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at, |
2492 | t_nbparam ***nbparam, t_nbparam ***pair) |
2493 | { |
2494 | t_atom atom; |
2495 | t_param param; |
2496 | int i, nr; |
2497 | |
2498 | /* Define an atom type with all parameters set to zero (no interactions) */ |
2499 | atom.q = 0.0; |
2500 | atom.m = 0.0; |
2501 | /* Type for decoupled atoms could be anything, |
2502 | * this should be changed automatically later when required. |
2503 | */ |
2504 | atom.ptype = eptAtom; |
2505 | for (i = 0; (i < MAXFORCEPARAM12); i++) |
2506 | { |
2507 | param.c[i] = 0.0; |
2508 | } |
2509 | |
2510 | nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0.0, 0.0, 0.0, 0, 0, 0); |
2511 | |
2512 | /* Add space in the non-bonded parameters matrix */ |
2513 | realloc_nb_params(at, nbparam, pair); |
2514 | |
2515 | return nr; |
2516 | } |
2517 | |
2518 | static void convert_pairs_to_pairsQ(t_params *plist, |
2519 | real fudgeQQ, t_atoms *atoms) |
2520 | { |
2521 | t_param *paramp1, *paramp2, *paramnew; |
2522 | int i, j, p1nr, p2nr, p2newnr; |
2523 | |
2524 | /* Add the pair list to the pairQ list */ |
2525 | p1nr = plist[F_LJ14].nr; |
2526 | p2nr = plist[F_LJC14_Q].nr; |
2527 | p2newnr = p1nr + p2nr; |
2528 | snew(paramnew, p2newnr)(paramnew) = save_calloc("paramnew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2528, (p2newnr), sizeof(*(paramnew))); |
2529 | |
2530 | paramp1 = plist[F_LJ14].param; |
2531 | paramp2 = plist[F_LJC14_Q].param; |
2532 | |
2533 | /* Fill in the new F_LJC14_Q array with the old one. NOTE: |
2534 | it may be possible to just ADD the converted F_LJ14 array |
2535 | to the old F_LJC14_Q array, but since we have to create |
2536 | a new sized memory structure, better just to deep copy it all. |
2537 | */ |
2538 | |
2539 | for (i = 0; i < p2nr; i++) |
2540 | { |
2541 | /* Copy over parameters */ |
2542 | for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */ |
2543 | { |
2544 | paramnew[i].c[j] = paramp2[i].c[j]; |
2545 | } |
2546 | |
2547 | /* copy over atoms */ |
2548 | for (j = 0; j < 2; j++) |
2549 | { |
2550 | paramnew[i].a[j] = paramp2[i].a[j]; |
2551 | } |
2552 | } |
2553 | |
2554 | for (i = p2nr; i < p2newnr; i++) |
2555 | { |
2556 | j = i-p2nr; |
2557 | |
2558 | /* Copy over parameters */ |
2559 | paramnew[i].c[0] = fudgeQQ; |
2560 | paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q; |
2561 | paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q; |
2562 | paramnew[i].c[3] = paramp1[j].c[0]; |
2563 | paramnew[i].c[4] = paramp1[j].c[1]; |
2564 | |
2565 | /* copy over atoms */ |
2566 | paramnew[i].a[0] = paramp1[j].a[0]; |
2567 | paramnew[i].a[1] = paramp1[j].a[1]; |
2568 | } |
2569 | |
2570 | /* free the old pairlists */ |
2571 | sfree(plist[F_LJC14_Q].param)save_free("plist[F_LJC14_Q].param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2571, (plist[F_LJC14_Q].param)); |
2572 | sfree(plist[F_LJ14].param)save_free("plist[F_LJ14].param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2572, (plist[F_LJ14].param)); |
2573 | |
2574 | /* now assign the new data to the F_LJC14_Q structure */ |
2575 | plist[F_LJC14_Q].param = paramnew; |
2576 | plist[F_LJC14_Q].nr = p2newnr; |
2577 | |
2578 | /* Empty the LJ14 pairlist */ |
2579 | plist[F_LJ14].nr = 0; |
2580 | plist[F_LJ14].param = NULL((void*)0); |
2581 | } |
2582 | |
2583 | static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp) |
2584 | { |
2585 | int n, ntype, i, j, k; |
2586 | t_atom *atom; |
2587 | t_blocka *excl; |
2588 | gmx_bool bExcl; |
2589 | t_param param; |
2590 | |
2591 | n = mol->atoms.nr; |
2592 | atom = mol->atoms.atom; |
2593 | |
2594 | ntype = sqrt(nbp->nr); |
2595 | |
2596 | for (i = 0; i < MAXATOMLIST6; i++) |
2597 | { |
2598 | param.a[i] = NOTSET-12345; |
2599 | } |
2600 | for (i = 0; i < MAXFORCEPARAM12; i++) |
2601 | { |
2602 | param.c[i] = NOTSET-12345; |
2603 | } |
2604 | |
2605 | /* Add a pair interaction for all non-excluded atom pairs */ |
2606 | excl = &mol->excls; |
2607 | for (i = 0; i < n; i++) |
2608 | { |
2609 | for (j = i+1; j < n; j++) |
2610 | { |
2611 | bExcl = FALSE0; |
2612 | for (k = excl->index[i]; k < excl->index[i+1]; k++) |
2613 | { |
2614 | if (excl->a[k] == j) |
2615 | { |
2616 | bExcl = TRUE1; |
2617 | } |
2618 | } |
2619 | if (!bExcl) |
2620 | { |
2621 | if (nb_funct != F_LJ) |
2622 | { |
2623 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2623, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones"); |
2624 | } |
2625 | param.a[0] = i; |
2626 | param.a[1] = j; |
2627 | param.c[0] = atom[i].q; |
2628 | param.c[1] = atom[j].q; |
2629 | param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0]; |
2630 | param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1]; |
2631 | add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m); |
2632 | } |
2633 | } |
2634 | } |
2635 | } |
2636 | |
2637 | static void set_excl_all(t_blocka *excl) |
2638 | { |
2639 | int nat, i, j, k; |
2640 | |
2641 | /* Get rid of the current exclusions and exclude all atom pairs */ |
2642 | nat = excl->nr; |
2643 | excl->nra = nat*nat; |
2644 | srenew(excl->a, excl->nra)(excl->a) = save_realloc("excl->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 2644, (excl->a), (excl->nra), sizeof(*(excl->a))); |
2645 | k = 0; |
2646 | for (i = 0; i < nat; i++) |
2647 | { |
2648 | excl->index[i] = k; |
2649 | for (j = 0; j < nat; j++) |
2650 | { |
2651 | excl->a[k++] = j; |
2652 | } |
2653 | } |
2654 | excl->index[nat] = k; |
2655 | } |
2656 | |
2657 | static void decouple_atoms(t_atoms *atoms, int atomtype_decouple, |
2658 | int couple_lam0, int couple_lam1) |
2659 | { |
2660 | int i; |
2661 | |
2662 | for (i = 0; i < atoms->nr; i++) |
2663 | { |
2664 | if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) |
2665 | { |
2666 | atoms->atom[i].q = 0.0; |
2667 | } |
2668 | if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) |
2669 | { |
2670 | atoms->atom[i].type = atomtype_decouple; |
2671 | } |
2672 | if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) |
2673 | { |
2674 | atoms->atom[i].qB = 0.0; |
2675 | } |
2676 | if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) |
2677 | { |
2678 | atoms->atom[i].typeB = atomtype_decouple; |
2679 | } |
2680 | } |
2681 | } |
2682 | |
2683 | void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ, |
2684 | int couple_lam0, int couple_lam1, |
2685 | gmx_bool bCoupleIntra, int nb_funct, t_params *nbp) |
2686 | { |
2687 | convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms); |
2688 | if (!bCoupleIntra) |
2689 | { |
2690 | generate_LJCpairsNB(mol, nb_funct, nbp); |
2691 | set_excl_all(&mol->excls); |
2692 | } |
2693 | decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1); |
2694 | } |