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, rfac=%15.8e",
185 iparams.thole.alpha1,
186 iparams.thole.alpha2,
190 writer->writeLineFormatted(
191 "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f",
200 writer->writeLineFormatted("c6=%15.8e, c12=%15.8e", iparams.lj.c6, iparams.lj.c12);
203 writer->writeLineFormatted("c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e",
210 writer->writeLineFormatted("fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
218 writer->writeLineFormatted("qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
228 writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d",
236 writer->writeLineFormatted(
237 "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)",
238 iparams.disres.label,
243 iparams.disres.kfac);
246 writer->writeLineFormatted(
247 "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)",
249 iparams.orires.label,
250 iparams.orires.power,
253 iparams.orires.kfac);
256 writer->writeLineFormatted(
257 "phiA=%15.8e, dphiA=%15.8e, kfacA=%15.8e, phiB=%15.8e, dphiB=%15.8e, "
260 iparams.dihres.dphiA,
261 iparams.dihres.kfacA,
263 iparams.dihres.dphiB,
264 iparams.dihres.kfacB);
267 writer->writeLineFormatted(
268 "pos0A=(%15.8e,%15.8e,%15.8e), fcA=(%15.8e,%15.8e,%15.8e), "
269 "pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)",
270 iparams.posres.pos0A[XX],
271 iparams.posres.pos0A[YY],
272 iparams.posres.pos0A[ZZ],
273 iparams.posres.fcA[XX],
274 iparams.posres.fcA[YY],
275 iparams.posres.fcA[ZZ],
276 iparams.posres.pos0B[XX],
277 iparams.posres.pos0B[YY],
278 iparams.posres.pos0B[ZZ],
279 iparams.posres.fcB[XX],
280 iparams.posres.fcB[YY],
281 iparams.posres.fcB[ZZ]);
284 writer->writeLineFormatted(
285 "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e",
286 iparams.fbposres.pos0[XX],
287 iparams.fbposres.pos0[YY],
288 iparams.fbposres.pos0[ZZ],
289 iparams.fbposres.geom,
294 for (int i = 0; i < NR_RBDIHS; i++)
296 writer->writeStringFormatted(
297 "%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams.rbdihs.rbcA[i]);
299 writer->ensureLineBreak();
300 for (int i = 0; i < NR_RBDIHS; i++)
302 writer->writeStringFormatted(
303 "%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams.rbdihs.rbcB[i]);
305 writer->ensureLineBreak();
309 /* Use the OPLS -> Ryckaert-Bellemans formula backwards to get
310 * the OPLS potential constants back.
312 const real* rbcA = iparams.rbdihs.rbcA;
313 const real* rbcB = iparams.rbdihs.rbcB;
316 VA[3] = -0.25 * rbcA[4];
317 VA[2] = -0.5 * rbcA[3];
318 VA[1] = 4.0 * VA[3] - rbcA[2];
319 VA[0] = 3.0 * VA[2] - 2.0 * rbcA[1];
321 VB[3] = -0.25 * rbcB[4];
322 VB[2] = -0.5 * rbcB[3];
323 VB[1] = 4.0 * VB[3] - rbcB[2];
324 VB[0] = 3.0 * VB[2] - 2.0 * rbcB[1];
326 for (int i = 0; i < NR_FOURDIHS; i++)
328 writer->writeStringFormatted("%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
330 writer->ensureLineBreak();
331 for (int i = 0; i < NR_FOURDIHS; i++)
333 writer->writeStringFormatted("%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
335 writer->ensureLineBreak();
341 writer->writeLineFormatted("dA=%15.8e, dB=%15.8e", iparams.constr.dA, iparams.constr.dB);
344 writer->writeLineFormatted("doh=%15.8e, dhh=%15.8e", iparams.settle.doh, iparams.settle.dhh);
346 case F_VSITE1: writer->ensureEmptyLine(); break;
347 case F_VSITE2: writer->writeLineFormatted("a=%15.8e", iparams.vsite.a); break;
351 writer->writeLineFormatted("a=%15.8e, b=%15.8e", iparams.vsite.a, iparams.vsite.b);
356 writer->writeLineFormatted(
357 "a=%15.8e, b=%15.8e, c=%15.8e", iparams.vsite.a, iparams.vsite.b, iparams.vsite.c);
360 writer->writeLineFormatted("n=%2d, a=%15.8e", iparams.vsiten.n, iparams.vsiten.a);
362 case F_GB12_NOLONGERUSED:
363 case F_GB13_NOLONGERUSED:
364 case F_GB14_NOLONGERUSED:
365 // These could only be generated by grompp, not written in
366 // a .top file. Now that implicit solvent is not
367 // supported, they can't be generated, and the values are
368 // ignored if read from an old .tpr file. So there is
372 writer->writeLineFormatted("cmapA=%1d, cmapB=%1d", iparams.cmap.cmapA, iparams.cmap.cmapB);
374 case F_RESTRANGLES: printHarmonicInteraction(writer, iparams, "ktheta", "costheta0"); break;
376 writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e", iparams.pdihs.phiA, iparams.pdihs.cpA);
379 writer->writeLineFormatted("kphi=%15.8e", iparams.cbtdihs.cbtcA[0]);
380 for (int i = 1; i < NR_CBTDIHS; i++)
382 writer->writeStringFormatted(", cbtcA[%d]=%15.8e", i - 1, iparams.cbtdihs.cbtcA[i]);
384 writer->ensureLineBreak();
388 "unknown function type %d (%s) in %s line %d",
390 interaction_function[ftype].name,
397 static void printIlist(FILE* fp,
400 const t_functype* functype,
402 gmx_bool bShowNumbers,
403 gmx_bool bShowParameters,
404 const t_iparams* iparams)
406 indent = pr_title(fp, indent, title);
407 pr_indent(fp, indent);
408 fprintf(fp, "nr: %d\n", ilist.size());
411 pr_indent(fp, indent);
412 fprintf(fp, "iatoms:\n");
414 for (int i = 0; i < ilist.size();)
416 pr_indent(fp, indent + INDENT);
417 const int type = ilist.iatoms[i];
418 const int ftype = functype[type];
421 fprintf(fp, "%d type=%d ", j, type);
424 printf("(%s)", interaction_function[ftype].name);
425 for (int k = 0; k < interaction_function[ftype].nratoms; k++)
427 fprintf(fp, " %3d", ilist.iatoms[i + 1 + k]);
432 pr_iparams(fp, ftype, iparams[type]);
435 i += 1 + interaction_function[ftype].nratoms;
440 void pr_ilist(FILE* fp,
443 const t_functype* functype,
444 const InteractionList& ilist,
445 gmx_bool bShowNumbers,
446 gmx_bool bShowParameters,
447 const t_iparams* iparams)
449 printIlist(fp, indent, title, functype, ilist, bShowNumbers, bShowParameters, iparams);
452 void pr_idef(FILE* fp, int indent, const char* title, const t_idef* idef, gmx_bool bShowNumbers, gmx_bool bShowParameters)
454 if (available(fp, idef, indent, title))
456 indent = pr_title(fp, indent, title);
457 pr_indent(fp, indent);
458 fprintf(fp, "atnr=%d\n", idef->atnr);
459 pr_indent(fp, indent);
460 fprintf(fp, "ntypes=%d\n", idef->ntypes);
461 for (int i = 0; i < idef->ntypes; i++)
463 pr_indent(fp, indent + INDENT);
466 bShowNumbers ? i : -1,
467 interaction_function[idef->functype[i]].name);
468 pr_iparams(fp, idef->functype[i], idef->iparams[i]);
470 pr_real(fp, indent, "fudgeQQ", idef->fudgeQQ);
472 for (int j = 0; (j < F_NRE); j++)
476 interaction_function[j].longname,
486 void init_idef(t_idef* idef)
490 idef->functype = nullptr;
491 idef->iparams = nullptr;
493 idef->iparams_posres = nullptr;
494 idef->iparams_fbposres = nullptr;
495 for (int f = 0; f < F_NRE; ++f)
497 idef->il[f].iatoms = nullptr;
498 idef->il[f].nalloc = 0;
503 InteractionDefinitions::InteractionDefinitions(const gmx_ffparams_t& ffparams) :
504 iparams(ffparams.iparams), functype(ffparams.functype), cmap_grid(ffparams.cmap_grid)
508 void InteractionDefinitions::clear()
510 /* Clear the counts */
511 for (auto& ilist : il)
515 iparams_posres.clear();
516 iparams_fbposres.clear();
519 void done_idef(t_idef* idef)
521 sfree(idef->functype);
522 sfree(idef->iparams);
523 sfree(idef->iparams_posres);
524 sfree(idef->iparams_fbposres);
525 for (int f = 0; f < F_NRE; ++f)
527 sfree(idef->il[f].iatoms);