2 * This file is part of the GROMACS molecular simulation package.
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,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/topology/forcefieldparameters.h"
45 #include "gromacs/topology/ifunc.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/utility/stringstream.h"
49 #include "gromacs/utility/textwriter.h"
50 #include "gromacs/utility/txtdump.h"
52 static void printHarmonicInteraction(gmx::TextWriter* writer,
53 const t_iparams& iparams,
57 writer->writeLineFormatted("%sA=%12.5e, %sA=%12.5e, %sB=%12.5e, %sB=%12.5e",
65 iparams.harmonic.krB);
68 void pr_iparams(FILE* fp, t_functype ftype, const t_iparams& iparams)
70 gmx::StringOutputStream stream;
72 gmx::TextWriter writer(&stream);
73 printInteractionParameters(&writer, ftype, iparams);
75 fputs(stream.toString().c_str(), fp);
78 void printInteractionParameters(gmx::TextWriter* writer, t_functype ftype, const t_iparams& iparams)
83 case F_G96ANGLES: printHarmonicInteraction(writer, iparams, "th", "ct"); break;
84 case F_CROSS_BOND_BONDS:
85 writer->writeLineFormatted("r1e=%15.8e, r2e=%15.8e, krr=%15.8e",
88 iparams.cross_bb.krr);
90 case F_CROSS_BOND_ANGLES:
91 writer->writeLineFormatted("r1e=%15.8e, r1e=%15.8e, r3e=%15.8e, krt=%15.8e",
95 iparams.cross_ba.krt);
98 writer->writeLineFormatted("klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e",
99 iparams.linangle.klinA,
101 iparams.linangle.klinB,
102 iparams.linangle.aB);
105 writer->writeLineFormatted(
106 "thetaA=%15.8e, kthetaA=%15.8e, r13A=%15.8e, kUBA=%15.8e, thetaB=%15.8e, "
107 "kthetaB=%15.8e, r13B=%15.8e, kUBB=%15.8e",
117 case F_QUARTIC_ANGLES:
118 writer->writeStringFormatted("theta=%15.8e", iparams.qangle.theta);
119 for (int i = 0; i < 5; i++)
121 writer->writeStringFormatted(", c%c=%15.8e", '0' + i, iparams.qangle.c[i]);
123 writer->ensureLineBreak();
126 writer->writeLineFormatted(
127 "a=%15.8e, b=%15.8e, c=%15.8e", iparams.bham.a, iparams.bham.b, iparams.bham.c);
131 case F_HARMONIC: printHarmonicInteraction(writer, iparams, "b0", "cb"); break;
132 case F_IDIHS: printHarmonicInteraction(writer, iparams, "xi", "cx"); break;
134 writer->writeLineFormatted(
135 "b0A=%15.8e, cbA=%15.8e, betaA=%15.8e, b0B=%15.8e, cbB=%15.8e, betaB=%15.8e",
141 iparams.morse.betaB);
144 writer->writeLineFormatted("b0=%15.8e, kb=%15.8e, kcub=%15.8e",
149 case F_CONNBONDS: writer->ensureEmptyLine(); break;
151 writer->writeLineFormatted("bm=%15.8e, kb=%15.8e", iparams.fene.bm, iparams.fene.kb);
154 writer->writeLineFormatted(
155 "lowA=%15.8e, up1A=%15.8e, up2A=%15.8e, kA=%15.8e, lowB=%15.8e, up1B=%15.8e, "
156 "up2B=%15.8e, kB=%15.8e,",
157 iparams.restraint.lowA,
158 iparams.restraint.up1A,
159 iparams.restraint.up2A,
160 iparams.restraint.kA,
161 iparams.restraint.lowB,
162 iparams.restraint.up1B,
163 iparams.restraint.up2B,
164 iparams.restraint.kB);
170 writer->writeLineFormatted(
171 "tab=%d, kA=%15.8e, kB=%15.8e", iparams.tab.table, iparams.tab.kA, iparams.tab.kB);
174 writer->writeLineFormatted("alpha=%15.8e", iparams.polarize.alpha);
177 writer->writeLineFormatted("alpha=%15.8e drcut=%15.8e khyp=%15.8e",
178 iparams.anharm_polarize.alpha,
179 iparams.anharm_polarize.drcut,
180 iparams.anharm_polarize.khyp);
183 writer->writeLineFormatted("a=%15.8e, alpha1=%15.8e, alpha2=%15.8e",
185 iparams.thole.alpha1,
186 iparams.thole.alpha2);
189 writer->writeLineFormatted(
190 "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f",
199 writer->writeLineFormatted("c6=%15.8e, c12=%15.8e", iparams.lj.c6, iparams.lj.c12);
202 writer->writeLineFormatted("c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e",
209 writer->writeLineFormatted("fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
217 writer->writeLineFormatted("qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
227 writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d",
235 writer->writeLineFormatted(
236 "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)",
237 iparams.disres.label,
242 iparams.disres.kfac);
245 writer->writeLineFormatted(
246 "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)",
248 iparams.orires.label,
249 iparams.orires.power,
252 iparams.orires.kfac);
255 writer->writeLineFormatted(
256 "phiA=%15.8e, dphiA=%15.8e, kfacA=%15.8e, phiB=%15.8e, dphiB=%15.8e, "
259 iparams.dihres.dphiA,
260 iparams.dihres.kfacA,
262 iparams.dihres.dphiB,
263 iparams.dihres.kfacB);
266 writer->writeLineFormatted(
267 "pos0A=(%15.8e,%15.8e,%15.8e), fcA=(%15.8e,%15.8e,%15.8e), "
268 "pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)",
269 iparams.posres.pos0A[XX],
270 iparams.posres.pos0A[YY],
271 iparams.posres.pos0A[ZZ],
272 iparams.posres.fcA[XX],
273 iparams.posres.fcA[YY],
274 iparams.posres.fcA[ZZ],
275 iparams.posres.pos0B[XX],
276 iparams.posres.pos0B[YY],
277 iparams.posres.pos0B[ZZ],
278 iparams.posres.fcB[XX],
279 iparams.posres.fcB[YY],
280 iparams.posres.fcB[ZZ]);
283 writer->writeLineFormatted(
284 "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e",
285 iparams.fbposres.pos0[XX],
286 iparams.fbposres.pos0[YY],
287 iparams.fbposres.pos0[ZZ],
288 iparams.fbposres.geom,
293 for (int i = 0; i < NR_RBDIHS; i++)
295 writer->writeStringFormatted(
296 "%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams.rbdihs.rbcA[i]);
298 writer->ensureLineBreak();
299 for (int i = 0; i < NR_RBDIHS; i++)
301 writer->writeStringFormatted(
302 "%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams.rbdihs.rbcB[i]);
304 writer->ensureLineBreak();
308 /* Use the OPLS -> Ryckaert-Bellemans formula backwards to get
309 * the OPLS potential constants back.
311 const real* rbcA = iparams.rbdihs.rbcA;
312 const real* rbcB = iparams.rbdihs.rbcB;
315 VA[3] = -0.25 * rbcA[4];
316 VA[2] = -0.5 * rbcA[3];
317 VA[1] = 4.0 * VA[3] - rbcA[2];
318 VA[0] = 3.0 * VA[2] - 2.0 * rbcA[1];
320 VB[3] = -0.25 * rbcB[4];
321 VB[2] = -0.5 * rbcB[3];
322 VB[1] = 4.0 * VB[3] - rbcB[2];
323 VB[0] = 3.0 * VB[2] - 2.0 * rbcB[1];
325 for (int i = 0; i < NR_FOURDIHS; i++)
327 writer->writeStringFormatted("%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
329 writer->ensureLineBreak();
330 for (int i = 0; i < NR_FOURDIHS; i++)
332 writer->writeStringFormatted("%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
334 writer->ensureLineBreak();
340 writer->writeLineFormatted("dA=%15.8e, dB=%15.8e", iparams.constr.dA, iparams.constr.dB);
343 writer->writeLineFormatted("doh=%15.8e, dhh=%15.8e", iparams.settle.doh, iparams.settle.dhh);
345 case F_VSITE1: writer->ensureEmptyLine(); break;
346 case F_VSITE2: writer->writeLineFormatted("a=%15.8e", iparams.vsite.a); break;
350 writer->writeLineFormatted("a=%15.8e, b=%15.8e", iparams.vsite.a, iparams.vsite.b);
355 writer->writeLineFormatted(
356 "a=%15.8e, b=%15.8e, c=%15.8e", iparams.vsite.a, iparams.vsite.b, iparams.vsite.c);
359 writer->writeLineFormatted("n=%2d, a=%15.8e", iparams.vsiten.n, iparams.vsiten.a);
361 case F_GB12_NOLONGERUSED:
362 case F_GB13_NOLONGERUSED:
363 case F_GB14_NOLONGERUSED:
364 // These could only be generated by grompp, not written in
365 // a .top file. Now that implicit solvent is not
366 // supported, they can't be generated, and the values are
367 // ignored if read from an old .tpr file. So there is
371 writer->writeLineFormatted("cmapA=%1d, cmapB=%1d", iparams.cmap.cmapA, iparams.cmap.cmapB);
373 case F_RESTRANGLES: printHarmonicInteraction(writer, iparams, "ktheta", "costheta0"); break;
375 writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e", iparams.pdihs.phiA, iparams.pdihs.cpA);
378 writer->writeLineFormatted("kphi=%15.8e", iparams.cbtdihs.cbtcA[0]);
379 for (int i = 1; i < NR_CBTDIHS; i++)
381 writer->writeStringFormatted(", cbtcA[%d]=%15.8e", i - 1, iparams.cbtdihs.cbtcA[i]);
383 writer->ensureLineBreak();
387 "unknown function type %d (%s) in %s line %d",
389 interaction_function[ftype].name,
396 static void printIlist(FILE* fp,
399 const t_functype* functype,
401 gmx_bool bShowNumbers,
402 gmx_bool bShowParameters,
403 const t_iparams* iparams)
405 indent = pr_title(fp, indent, title);
406 pr_indent(fp, indent);
407 fprintf(fp, "nr: %d\n", ilist.size());
410 pr_indent(fp, indent);
411 fprintf(fp, "iatoms:\n");
413 for (int i = 0; i < ilist.size();)
415 pr_indent(fp, indent + INDENT);
416 const int type = ilist.iatoms[i];
417 const int ftype = functype[type];
420 fprintf(fp, "%d type=%d ", j, type);
423 printf("(%s)", interaction_function[ftype].name);
424 for (int k = 0; k < interaction_function[ftype].nratoms; k++)
426 fprintf(fp, " %3d", ilist.iatoms[i + 1 + k]);
431 pr_iparams(fp, ftype, iparams[type]);
434 i += 1 + interaction_function[ftype].nratoms;
439 void pr_ilist(FILE* fp,
442 const t_functype* functype,
443 const InteractionList& ilist,
444 gmx_bool bShowNumbers,
445 gmx_bool bShowParameters,
446 const t_iparams* iparams)
448 printIlist(fp, indent, title, functype, ilist, bShowNumbers, bShowParameters, iparams);
451 void pr_idef(FILE* fp, int indent, const char* title, const t_idef* idef, gmx_bool bShowNumbers, gmx_bool bShowParameters)
453 if (available(fp, idef, indent, title))
455 indent = pr_title(fp, indent, title);
456 pr_indent(fp, indent);
457 fprintf(fp, "atnr=%d\n", idef->atnr);
458 pr_indent(fp, indent);
459 fprintf(fp, "ntypes=%d\n", idef->ntypes);
460 for (int i = 0; i < idef->ntypes; i++)
462 pr_indent(fp, indent + INDENT);
465 bShowNumbers ? i : -1,
466 interaction_function[idef->functype[i]].name);
467 pr_iparams(fp, idef->functype[i], idef->iparams[i]);
469 pr_real(fp, indent, "fudgeQQ", idef->fudgeQQ);
471 for (int j = 0; (j < F_NRE); j++)
475 interaction_function[j].longname,
485 void init_idef(t_idef* idef)
489 idef->functype = nullptr;
490 idef->iparams = nullptr;
492 idef->iparams_posres = nullptr;
493 idef->iparams_fbposres = nullptr;
494 for (int f = 0; f < F_NRE; ++f)
496 idef->il[f].iatoms = nullptr;
497 idef->il[f].nalloc = 0;
502 InteractionDefinitions::InteractionDefinitions(const gmx_ffparams_t& ffparams) :
503 iparams(ffparams.iparams), functype(ffparams.functype), cmap_grid(ffparams.cmap_grid)
507 void InteractionDefinitions::clear()
509 /* Clear the counts */
510 for (auto& ilist : il)
514 iparams_posres.clear();
515 iparams_fbposres.clear();
518 void done_idef(t_idef* idef)
520 sfree(idef->functype);
521 sfree(idef->iparams);
522 sfree(idef->iparams_posres);
523 sfree(idef->iparams_fbposres);
524 for (int f = 0; f < F_NRE; ++f)
526 sfree(idef->il[f].iatoms);