?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Split up compare.*
[gromacs.git] / src / gromacs / topology / topology.cpp
bloba94f36ee5a5e80d81b97adbc96bb1111a7009d43
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,2015,2016, 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 #include "gmxpre.h"
38
39 #include "topology.h"
40
41 #include <cstdio>
42
43 #include <algorithm>
44
45 #include "gromacs/math/vecdump.h"
46 #include "gromacs/topology/atoms.h"
47 #include "gromacs/topology/idef.h"
48 #include "gromacs/topology/ifunc.h"
49 #include "gromacs/topology/symtab.h"
50 #include "gromacs/utility/compare.h"
51 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/utility/stringutil.h"
53 #include "gromacs/utility/txtdump.h"
54
55 const char *gtypes[egcNR+1] = {
56 "T-Coupling", "Energy Mon.", "Acceleration", "Freeze",
57 "User1", "User2", "VCM", "Compressed X", "Or. Res. Fit", "QMMM", NULL
58 };
59
60 static void init_groups(gmx_groups_t *groups)
61 {
62 groups->ngrpname = 0;
63 groups->grpname = NULL;
64 for (int g = 0; g < egcNR; g++)
65 {
66 groups->grps[g].nm_ind = NULL;
67 groups->ngrpnr[g] = 0;
68 groups->grpnr[g] = NULL;
69 }
70
71 }
72
73 void init_mtop(gmx_mtop_t *mtop)
74 {
75 mtop->name = NULL;
76 mtop->nmoltype = 0;
77 mtop->moltype = NULL;
78 mtop->nmolblock = 0;
79 mtop->molblock = NULL;
80 mtop->maxres_renum = 0;
81 mtop->maxresnr = -1;
82 init_groups(&mtop->groups);
83 init_block(&mtop->mols);
84 open_symtab(&mtop->symtab);
85 }
86
87 void init_top(t_topology *top)
88 {
89 top->name = NULL;
90 init_atom(&(top->atoms));
91 init_atomtypes(&(top->atomtypes));
92 init_block(&top->cgs);
93 init_block(&top->mols);
94 init_blocka(&top->excls);
95 open_symtab(&top->symtab);
96 }
97
98
99 void done_moltype(gmx_moltype_t *molt)
100 {
101 done_atom(&molt->atoms);
102 done_block(&molt->cgs);
103 done_blocka(&molt->excls);
104
105 for (int f = 0; f < F_NRE; f++)
106 {
107 sfree(molt->ilist[f].iatoms);
108 molt->ilist[f].nalloc = 0;
109 }
110 }
111
112 void done_molblock(gmx_molblock_t *molb)
113 {
114 if (molb->nposres_xA > 0)
115 {
116 molb->nposres_xA = 0;
117 sfree(molb->posres_xA);
118 }
119 if (molb->nposres_xB > 0)
120 {
121 molb->nposres_xB = 0;
122 sfree(molb->posres_xB);
123 }
124 }
125
126 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
127 {
128 if (bDoneSymtab)
129 {
130 done_symtab(&mtop->symtab);
131 }
132
133 sfree(mtop->ffparams.functype);
134 sfree(mtop->ffparams.iparams);
135
136 for (int i = 0; i < mtop->nmoltype; i++)
137 {
138 done_moltype(&mtop->moltype[i]);
139 }
140 sfree(mtop->moltype);
141 for (int i = 0; i < mtop->nmolblock; i++)
142 {
143 done_molblock(&mtop->molblock[i]);
144 }
145 sfree(mtop->molblock);
146 done_block(&mtop->mols);
147 }
148
149 void done_top(t_topology *top)
150 {
151 sfree(top->idef.functype);
152 sfree(top->idef.iparams);
153 for (int f = 0; f < F_NRE; ++f)
154 {
155 sfree(top->idef.il[f].iatoms);
156 top->idef.il[f].iatoms = NULL;
157 top->idef.il[f].nalloc = 0;
158 }
159
160 done_atom(&(top->atoms));
161
162 /* For GB */
163 done_atomtypes(&(top->atomtypes));
164
165 done_symtab(&(top->symtab));
166 done_block(&(top->cgs));
167 done_block(&(top->mols));
168 done_blocka(&(top->excls));
169 }
170
171 static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
172 {
173 int i, j;
174
175 for (i = 0; (i < egcNR); i++)
176 {
177 fprintf(fp, "%s[%-12s] nr=%d, name=[", title, gtypes[i], grps[i].nr);
178 for (j = 0; (j < grps[i].nr); j++)
179 {
180 fprintf(fp, " %s", *(grpname[grps[i].nm_ind[j]]));
181 }
182 fprintf(fp, "]\n");
183 }
184 }
185
186 static void pr_groups(FILE *fp, int indent,
187 const gmx_groups_t *groups,
188 gmx_bool bShowNumbers)
189 {
190 int nat_max, i, g;
191
192 pr_grps(fp, "grp", groups->grps, groups->grpname);
193 pr_strings(fp, indent, "grpname", groups->grpname, groups->ngrpname, bShowNumbers);
194
195 pr_indent(fp, indent);
196 fprintf(fp, "groups ");
197 for (g = 0; g < egcNR; g++)
198 {
199 printf(" %5.5s", gtypes[g]);
200 }
201 printf("\n");
202
203 pr_indent(fp, indent);
204 fprintf(fp, "allocated ");
205 nat_max = 0;
206 for (g = 0; g < egcNR; g++)
207 {
208 printf(" %5d", groups->ngrpnr[g]);
209 nat_max = std::max(nat_max, groups->ngrpnr[g]);
210 }
211 printf("\n");
212
213 if (nat_max == 0)
214 {
215 pr_indent(fp, indent);
216 fprintf(fp, "groupnr[%5s] =", "*");
217 for (g = 0; g < egcNR; g++)
218 {
219 fprintf(fp, " %3d ", 0);
220 }
221 fprintf(fp, "\n");
222 }
223 else
224 {
225 for (i = 0; i < nat_max; i++)
226 {
227 pr_indent(fp, indent);
228 fprintf(fp, "groupnr[%5d] =", i);
229 for (g = 0; g < egcNR; g++)
230 {
231 fprintf(fp, " %3d ",
232 groups->grpnr[g] ? groups->grpnr[g][i] : 0);
233 }
234 fprintf(fp, "\n");
235 }
236 }
237 }
238
239 static void pr_moltype(FILE *fp, int indent, const char *title,
240 const gmx_moltype_t *molt, int n,
241 const gmx_ffparams_t *ffparams,
242 gmx_bool bShowNumbers)
243 {
244 int j;
245
246 indent = pr_title_n(fp, indent, title, n);
247 pr_indent(fp, indent);
248 fprintf(fp, "name=\"%s\"\n", *(molt->name));
249 pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
250 pr_block(fp, indent, "cgs", &molt->cgs, bShowNumbers);
251 pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers);
252 for (j = 0; (j < F_NRE); j++)
253 {
254 pr_ilist(fp, indent, interaction_function[j].longname,
255 ffparams->functype, &molt->ilist[j], bShowNumbers);
256 }
257 }
258
259 static void pr_molblock(FILE *fp, int indent, const char *title,
260 const gmx_molblock_t *molb, int n,
261 const gmx_moltype_t *molt)
262 {
263 indent = pr_title_n(fp, indent, title, n);
264 pr_indent(fp, indent);
265 fprintf(fp, "%-20s = %d \"%s\"\n",
266 "moltype", molb->type, *(molt[molb->type].name));
267 pr_int(fp, indent, "#molecules", molb->nmol);
268 pr_int(fp, indent, "#atoms_mol", molb->natoms_mol);
269 pr_int(fp, indent, "#posres_xA", molb->nposres_xA);
270 if (molb->nposres_xA > 0)
271 {
272 pr_rvecs(fp, indent, "posres_xA", molb->posres_xA, molb->nposres_xA);
273 }
274 pr_int(fp, indent, "#posres_xB", molb->nposres_xB);
275 if (molb->nposres_xB > 0)
276 {
277 pr_rvecs(fp, indent, "posres_xB", molb->posres_xB, molb->nposres_xB);
278 }
279 }
280
281 void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
282 gmx_bool bShowNumbers)
283 {
284 int mt, mb, j;
285
286 if (available(fp, mtop, indent, title))
287 {
288 indent = pr_title(fp, indent, title);
289 pr_indent(fp, indent);
290 fprintf(fp, "name=\"%s\"\n", *(mtop->name));
291 pr_int(fp, indent, "#atoms", mtop->natoms);
292 pr_int(fp, indent, "#molblock", mtop->nmolblock);
293 for (mb = 0; mb < mtop->nmolblock; mb++)
294 {
295 pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
296 }
297 pr_str(fp, indent, "bIntermolecularInteractions",
298 gmx::boolToString(mtop->bIntermolecularInteractions));
299 if (mtop->bIntermolecularInteractions)
300 {
301 for (j = 0; (j < F_NRE); j++)
302 {
303 pr_ilist(fp, indent, interaction_function[j].longname,
304 mtop->ffparams.functype,
305 &mtop->intermolecular_ilist[j], bShowNumbers);
306 }
307 }
308 pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
309 pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
310 for (mt = 0; mt < mtop->nmoltype; mt++)
311 {
312 pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt,
313 &mtop->ffparams, bShowNumbers);
314 }
315 pr_groups(fp, indent, &mtop->groups, bShowNumbers);
316 }
317 }
318
319 void pr_top(FILE *fp, int indent, const char *title, const t_topology *top, gmx_bool bShowNumbers)
320 {
321 if (available(fp, top, indent, title))
322 {
323 indent = pr_title(fp, indent, title);
324 pr_indent(fp, indent);
325 fprintf(fp, "name=\"%s\"\n", *(top->name));
326 pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
327 pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
328 pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
329 pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
330 pr_str(fp, indent, "bIntermolecularInteractions",
331 gmx::boolToString(top->bIntermolecularInteractions));
332 pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
333 pr_idef(fp, indent, "idef", &top->idef, bShowNumbers);
334 }
335 }
336
337 static void cmp_ilist(FILE *fp, int ftype, const t_ilist *il1, const t_ilist *il2)
338 {
339 int i;
340 char buf[256];
341
342 fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
343 sprintf(buf, "%s->nr", interaction_function[ftype].name);
344 cmp_int(fp, buf, -1, il1->nr, il2->nr);
345 sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
346 if (((il1->nr > 0) && (!il1->iatoms)) ||
347 ((il2->nr > 0) && (!il2->iatoms)) ||
348 ((il1->nr != il2->nr)))
349 {
350 fprintf(fp, "Comparing radically different topologies - %s is different\n",
351 buf);
352 }
353 else
354 {
355 for (i = 0; (i < il1->nr); i++)
356 {
357 cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
358 }
359 }
360 }
361
362 static void cmp_iparm(FILE *fp, const char *s, t_functype ft,
363 const t_iparams &ip1, const t_iparams &ip2, real ftol, real abstol)
364 {
365 int i;
366 gmx_bool bDiff;
367
368 bDiff = FALSE;
369 for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
370 {
371 bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
372 }
373 if (bDiff)
374 {
375 fprintf(fp, "%s1: ", s);
376 pr_iparams(fp, ft, &ip1);
377 fprintf(fp, "%s2: ", s);
378 pr_iparams(fp, ft, &ip2);
379 }
380 }
381
382 static void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft,
383 const t_iparams &ip1, real ftol, real abstol)
384 {
385 int nrfpA, nrfpB, p0, i;
386 gmx_bool bDiff;
387
388 /* Normally the first parameter is perturbable */
389 p0 = 0;
390 nrfpA = interaction_function[ft].nrfpA;
391 nrfpB = interaction_function[ft].nrfpB;
392 if (ft == F_PDIHS)
393 {
394 nrfpB = 2;
395 }
396 else if (interaction_function[ft].flags & IF_TABULATED)
397 {
398 /* For tabulated interactions only the second parameter is perturbable */
399 p0 = 1;
400 nrfpB = 1;
401 }
402 bDiff = FALSE;
403 for (i = 0; i < nrfpB && !bDiff; i++)
404 {
405 bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
406 }
407 if (bDiff)
408 {
409 fprintf(fp, "%s: ", s);
410 pr_iparams(fp, ft, &ip1);
411 }
412 }
413
414 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
415 {
416 cmp_int(fp, "cmap ngrid", -1, cmap1->ngrid, cmap2->ngrid);
417 cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
418 if (cmap1->ngrid == cmap2->ngrid &&
419 cmap1->grid_spacing == cmap2->grid_spacing)
420 {
421 int g;
422
423 for (g = 0; g < cmap1->ngrid; g++)
424 {
425 int i;
426
427 fprintf(fp, "comparing cmap %d\n", g);
428
429 for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
430 {
431 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
432 }
433 }
434 }
435 }
436
437 static void cmp_idef(FILE *fp, const t_idef *id1, const t_idef *id2, real ftol, real abstol)
438 {
439 int i;
440 char buf1[64], buf2[64];
441
442 fprintf(fp, "comparing idef\n");
443 if (id2)
444 {
445 cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
446 cmp_int(fp, "idef->atnr", -1, id1->atnr, id2->atnr);
447 for (i = 0; (i < std::min(id1->ntypes, id2->ntypes)); i++)
448 {
449 sprintf(buf1, "idef->functype[%d]", i);
450 sprintf(buf2, "idef->iparam[%d]", i);
451 cmp_int(fp, buf1, i, (int)id1->functype[i], (int)id2->functype[i]);
452 cmp_iparm(fp, buf2, id1->functype[i],
453 id1->iparams[i], id2->iparams[i], ftol, abstol);
454 }
455 cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
456 cmp_cmap(fp, &id1->cmap_grid, &id2->cmap_grid, ftol, abstol);
457 for (i = 0; (i < F_NRE); i++)
458 {
459 cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
460 }
461 }
462 else
463 {
464 for (i = 0; (i < id1->ntypes); i++)
465 {
466 cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
467 }
468 }
469 }
470
471 static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s)
472 {
473 char buf[32];
474
475 fprintf(fp, "comparing block %s\n", s);
476 sprintf(buf, "%s.nr", s);
477 cmp_int(fp, buf, -1, b1->nr, b2->nr);
478 }
479
480 static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s)
481 {
482 char buf[32];
483
484 fprintf(fp, "comparing blocka %s\n", s);
485 sprintf(buf, "%s.nr", s);
486 cmp_int(fp, buf, -1, b1->nr, b2->nr);
487 sprintf(buf, "%s.nra", s);
488 cmp_int(fp, buf, -1, b1->nra, b2->nra);
489 }
490
491 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol)
492 {
493 fprintf(fp, "comparing top\n");
494 if (t2)
495 {
496 cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
497 cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
498 cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
499 cmp_block(fp, &t1->mols, &t2->mols, "mols");
500 cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
501 cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
502 }
503 else
504 {
505 cmp_idef(fp, &(t1->idef), NULL, ftol, abstol);
506 cmp_atoms(fp, &(t1->atoms), NULL, ftol, abstol);
507 }
508 }
509
510 void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
511 int natoms0, int natoms1)
512 {
513 int i, j;
514 char buf[32];
515
516 fprintf(fp, "comparing groups\n");
517
518 for (i = 0; i < egcNR; i++)
519 {
520 sprintf(buf, "grps[%d].nr", i);
521 cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
522 if (g0->grps[i].nr == g1->grps[i].nr)
523 {
524 for (j = 0; j < g0->grps[i].nr; j++)
525 {
526 sprintf(buf, "grps[%d].name[%d]", i, j);
527 cmp_str(fp, buf, -1,
528 *g0->grpname[g0->grps[i].nm_ind[j]],
529 *g1->grpname[g1->grps[i].nm_ind[j]]);
530 }
531 }
532 cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
533 if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
534 (g0->grpnr[i] != NULL || g1->grpnr[i] != NULL))
535 {
536 for (j = 0; j < natoms0; j++)
537 {
538 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
539 }
540 }
541 }
542 /* We have compared the names in the groups lists,
543 * so we can skip the grpname list comparison.
544 */
545 }
The ultimate molecular dynamics simulation package