/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c

Bug Summary

File:	gromacs/gmxpreprocess/toppush.c
Location:	line 2213, column 5
Description:	Value stored to 'ret' is never read

Annotated Source Code

1	/*
2	* This file is part of the GROMACS molecular simulation package.
3	*
4	* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5	* Copyright (c) 2001-2004, The GROMACS development team.
6	* Copyright (c) 2013,2014, by the GROMACS development team, led by
7	* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8	* and including many others, as listed in the AUTHORS file in the
9	* top-level source directory and at http://www.gromacs.org.
10	*
11	* GROMACS is free software; you can redistribute it and/or
12	* modify it under the terms of the GNU Lesser General Public License
13	* as published by the Free Software Foundation; either version 2.1
14	* of the License, or (at your option) any later version.
15	*
16	* GROMACS is distributed in the hope that it will be useful,
17	* but WITHOUT ANY WARRANTY; without even the implied warranty of
18	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19	* Lesser General Public License for more details.
20	*
21	* You should have received a copy of the GNU Lesser General Public
22	* License along with GROMACS; if not, see
23	* http://www.gnu.org/licenses, or write to the Free Software Foundation,
24	* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25	*
26	* If you want to redistribute modifications to GROMACS, please
27	* consider that scientific software is very special. Version
28	* control is crucial - bugs must be traceable. We will be happy to
29	* consider code for inclusion in the official distribution, but
30	* derived work must not be called official GROMACS. Details are found
31	* in the README & COPYING files - if they are missing, get the
32	* official version at http://www.gromacs.org.
33	*
34	* To help us fund GROMACS development, we humbly ask that you cite
35	* the research papers on the package. Check out http://www.gromacs.org.
36	*/
37	#ifdef HAVE_CONFIG_H1
38	#include <config.h>
39	#endif
40
41	#include <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, nrnr)(plist->param) = save_calloc("plist->param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/toppush.c" , 72, (nrnr), 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, &param, FALSE0, &param_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, &param, TRUE1, &param_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, &param, 0, FALSE0, bGenPairs);
1848	bFoundB = default_nb_params(ftype, bondtype, at, &param, 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, &param, 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], &param);
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], &param);
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;
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, &param, 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], &param);
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);
	Value stored to 'ret' is never read
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], &param);
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", &param, -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], &param);
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	}