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,2017, The GROMACS development team.
7 * Copyright (c) 2018,2019,2020, 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.
51 #include "gromacs/applied_forces/awh/read_params.h"
52 #include "gromacs/fileio/readinp.h"
53 #include "gromacs/fileio/warninp.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/calc_verletbuf.h"
60 #include "gromacs/mdrun/mdmodules.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/multipletimestepping.h"
64 #include "gromacs/mdtypes/pull_params.h"
65 #include "gromacs/options/options.h"
66 #include "gromacs/options/treesupport.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/selection/indexutil.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/topology/index.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/symtab.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/exceptions.h"
77 #include "gromacs/utility/fatalerror.h"
78 #include "gromacs/utility/filestream.h"
79 #include "gromacs/utility/gmxassert.h"
80 #include "gromacs/utility/ikeyvaluetreeerror.h"
81 #include "gromacs/utility/keyvaluetree.h"
82 #include "gromacs/utility/keyvaluetreebuilder.h"
83 #include "gromacs/utility/keyvaluetreemdpwriter.h"
84 #include "gromacs/utility/keyvaluetreetransform.h"
85 #include "gromacs/utility/mdmodulenotification.h"
86 #include "gromacs/utility/smalloc.h"
87 #include "gromacs/utility/strconvert.h"
88 #include "gromacs/utility/stringcompare.h"
89 #include "gromacs/utility/stringutil.h"
90 #include "gromacs/utility/textwriter.h"
95 /* Resource parameters
96 * Do not change any of these until you read the instruction
97 * in readinp.h. Some cpp's do not take spaces after the backslash
98 * (like the c-shell), which will give you a very weird compiler
102 struct gmx_inputrec_strings
104 char tcgrps
[STRLEN
], tau_t
[STRLEN
], ref_t
[STRLEN
], acc
[STRLEN
], accgrps
[STRLEN
], freeze
[STRLEN
],
105 frdim
[STRLEN
], energy
[STRLEN
], user1
[STRLEN
], user2
[STRLEN
], vcm
[STRLEN
],
106 x_compressed_groups
[STRLEN
], couple_moltype
[STRLEN
], orirefitgrp
[STRLEN
],
107 egptable
[STRLEN
], egpexcl
[STRLEN
], wall_atomtype
[STRLEN
], wall_density
[STRLEN
],
108 deform
[STRLEN
], QMMM
[STRLEN
], imd_grp
[STRLEN
];
109 char fep_lambda
[efptNR
][STRLEN
];
110 char lambda_weights
[STRLEN
];
111 std::vector
<std::string
> pullGroupNames
;
112 std::vector
<std::string
> rotateGroupNames
;
113 char anneal
[STRLEN
], anneal_npoints
[STRLEN
], anneal_time
[STRLEN
], anneal_temp
[STRLEN
];
116 static gmx_inputrec_strings
* inputrecStrings
= nullptr;
118 void init_inputrec_strings()
123 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
124 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
126 inputrecStrings
= new gmx_inputrec_strings();
129 void done_inputrec_strings()
131 delete inputrecStrings
;
132 inputrecStrings
= nullptr;
138 egrptpALL
, /* All particles have to be a member of a group. */
139 egrptpALL_GENREST
, /* A rest group with name is generated for particles *
140 * that are not part of any group. */
141 egrptpPART
, /* As egrptpALL_GENREST, but no name is generated *
142 * for the rest group. */
143 egrptpONE
/* Merge all selected groups into one group, *
144 * make a rest group for the remaining particles. */
147 static const char* constraints
[eshNR
+ 1] = { "none", "h-bonds", "all-bonds",
148 "h-angles", "all-angles", nullptr };
150 static const char* couple_lam
[ecouplamNR
+ 1] = { "vdw-q", "vdw", "q", "none", nullptr };
152 static void GetSimTemps(int ntemps
, t_simtemp
* simtemp
, double* temperature_lambdas
)
157 for (i
= 0; i
< ntemps
; i
++)
159 /* simple linear scaling -- allows more control */
160 if (simtemp
->eSimTempScale
== esimtempLINEAR
)
162 simtemp
->temperatures
[i
] =
164 + (simtemp
->simtemp_high
- simtemp
->simtemp_low
) * temperature_lambdas
[i
];
166 else if (simtemp
->eSimTempScale
167 == esimtempGEOMETRIC
) /* should give roughly equal acceptance for constant heat capacity . . . */
169 simtemp
->temperatures
[i
] = simtemp
->simtemp_low
170 * std::pow(simtemp
->simtemp_high
/ simtemp
->simtemp_low
,
171 static_cast<real
>((1.0 * i
) / (ntemps
- 1)));
173 else if (simtemp
->eSimTempScale
== esimtempEXPONENTIAL
)
175 simtemp
->temperatures
[i
] = simtemp
->simtemp_low
176 + (simtemp
->simtemp_high
- simtemp
->simtemp_low
)
177 * (std::expm1(temperature_lambdas
[i
]) / std::expm1(1.0));
182 sprintf(errorstr
, "eSimTempScale=%d not defined", simtemp
->eSimTempScale
);
183 gmx_fatal(FARGS
, "%s", errorstr
);
189 static void _low_check(bool b
, const char* s
, warninp_t wi
)
193 warning_error(wi
, s
);
197 static void check_nst(const char* desc_nst
, int nst
, const char* desc_p
, int* p
, warninp_t wi
)
201 if (*p
> 0 && *p
% nst
!= 0)
203 /* Round up to the next multiple of nst */
204 *p
= ((*p
) / nst
+ 1) * nst
;
205 sprintf(buf
, "%s should be a multiple of %s, changing %s to %d\n", desc_p
, desc_nst
, desc_p
, *p
);
210 static int lcd(int n1
, int n2
)
215 for (i
= 2; (i
<= n1
&& i
<= n2
); i
++)
217 if (n1
% i
== 0 && n2
% i
== 0)
226 //! Convert legacy mdp entries to modern ones.
227 static void process_interaction_modifier(int* eintmod
)
229 if (*eintmod
== eintmodPOTSHIFT_VERLET_UNSUPPORTED
)
231 *eintmod
= eintmodPOTSHIFT
;
235 void check_ir(const char* mdparin
,
236 const gmx::MdModulesNotifier
& mdModulesNotifier
,
240 /* Check internal consistency.
241 * NOTE: index groups are not set here yet, don't check things
242 * like temperature coupling group options here, but in triple_check
245 /* Strange macro: first one fills the err_buf, and then one can check
246 * the condition, which will print the message and increase the error
249 #define CHECK(b) _low_check(b, err_buf, wi)
250 char err_buf
[256], warn_buf
[STRLEN
];
253 t_lambda
* fep
= ir
->fepvals
;
254 t_expanded
* expand
= ir
->expandedvals
;
256 set_warning_line(wi
, mdparin
, -1);
258 /* We cannot check MTS requirements with an invalid MTS setup
259 * and we will already have generated errors with an invalid MTS setup.
261 if (gmx::haveValidMtsSetup(*ir
))
263 std::vector
<std::string
> errorMessages
= gmx::checkMtsRequirements(*ir
);
265 for (const auto& errorMessage
: errorMessages
)
267 warning_error(wi
, errorMessage
.c_str());
271 if (ir
->coulombtype
== eelRF_NEC_UNSUPPORTED
)
273 sprintf(warn_buf
, "%s electrostatics is no longer supported", eel_names
[eelRF_NEC_UNSUPPORTED
]);
274 warning_error(wi
, warn_buf
);
277 /* BASIC CUT-OFF STUFF */
278 if (ir
->rcoulomb
< 0)
280 warning_error(wi
, "rcoulomb should be >= 0");
284 warning_error(wi
, "rvdw should be >= 0");
286 if (ir
->rlist
< 0 && !(ir
->cutoff_scheme
== ecutsVERLET
&& ir
->verletbuf_tol
> 0))
288 warning_error(wi
, "rlist should be >= 0");
291 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
292 "neighbour-list update scheme for efficient buffering for improved energy "
293 "conservation, please use the Verlet cut-off scheme instead.)");
294 CHECK(ir
->nstlist
< 0);
296 process_interaction_modifier(&ir
->coulomb_modifier
);
297 process_interaction_modifier(&ir
->vdw_modifier
);
299 if (ir
->cutoff_scheme
== ecutsGROUP
)
302 "The group cutoff scheme has been removed since GROMACS 2020. "
303 "Please use the Verlet cutoff scheme.");
305 if (ir
->cutoff_scheme
== ecutsVERLET
)
309 /* Normal Verlet type neighbor-list, currently only limited feature support */
310 if (inputrec2nboundeddim(ir
) < 3)
312 warning_error(wi
, "With Verlet lists only full pbc or pbc=xy with walls is supported");
315 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
316 if (ir
->rcoulomb
!= ir
->rvdw
)
318 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
319 // for PME load balancing, we can support this exception.
320 bool bUsesPmeTwinRangeKernel
= (EEL_PME_EWALD(ir
->coulombtype
) && ir
->vdwtype
== evdwCUT
321 && ir
->rcoulomb
> ir
->rvdw
);
322 if (!bUsesPmeTwinRangeKernel
)
325 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
326 "rcoulomb>rvdw with PME electrostatics)");
330 if (ir
->vdwtype
== evdwSHIFT
|| ir
->vdwtype
== evdwSWITCH
)
332 if (ir
->vdw_modifier
== eintmodNONE
|| ir
->vdw_modifier
== eintmodPOTSHIFT
)
334 ir
->vdw_modifier
= (ir
->vdwtype
== evdwSHIFT
? eintmodFORCESWITCH
: eintmodPOTSWITCH
);
337 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
339 evdw_names
[ir
->vdwtype
], evdw_names
[evdwCUT
], eintmod_names
[ir
->vdw_modifier
]);
340 warning_note(wi
, warn_buf
);
342 ir
->vdwtype
= evdwCUT
;
346 sprintf(warn_buf
, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
347 evdw_names
[ir
->vdwtype
], eintmod_names
[ir
->vdw_modifier
]);
348 warning_error(wi
, warn_buf
);
352 if (!(ir
->vdwtype
== evdwCUT
|| ir
->vdwtype
== evdwPME
))
355 "With Verlet lists only cut-off and PME LJ interactions are supported");
357 if (!(ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
) || EEL_PME(ir
->coulombtype
)
358 || ir
->coulombtype
== eelEWALD
))
361 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
362 "electrostatics are supported");
364 if (!(ir
->coulomb_modifier
== eintmodNONE
|| ir
->coulomb_modifier
== eintmodPOTSHIFT
))
366 sprintf(warn_buf
, "coulomb_modifier=%s is not supported", eintmod_names
[ir
->coulomb_modifier
]);
367 warning_error(wi
, warn_buf
);
370 if (EEL_USER(ir
->coulombtype
))
372 sprintf(warn_buf
, "Coulomb type %s is not supported with the verlet scheme",
373 eel_names
[ir
->coulombtype
]);
374 warning_error(wi
, warn_buf
);
377 if (ir
->nstlist
<= 0)
379 warning_error(wi
, "With Verlet lists nstlist should be larger than 0");
382 if (ir
->nstlist
< 10)
385 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
386 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
390 rc_max
= std::max(ir
->rvdw
, ir
->rcoulomb
);
394 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
395 ir
->verletbuf_tol
= 0;
398 else if (ir
->verletbuf_tol
<= 0)
400 if (ir
->verletbuf_tol
== 0)
402 warning_error(wi
, "Can not have Verlet buffer tolerance of exactly 0");
405 if (ir
->rlist
< rc_max
)
408 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
411 if (ir
->rlist
== rc_max
&& ir
->nstlist
> 1)
415 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
416 "buffer. The cluster pair list does have a buffering effect, but choosing "
417 "a larger rlist might be necessary for good energy conservation.");
422 if (ir
->rlist
> rc_max
)
425 "You have set rlist larger than the interaction cut-off, but you also "
426 "have verlet-buffer-tolerance > 0. Will set rlist using "
427 "verlet-buffer-tolerance.");
430 if (ir
->nstlist
== 1)
432 /* No buffer required */
437 if (EI_DYNAMICS(ir
->eI
))
439 if (inputrec2nboundeddim(ir
) < 3)
442 "The box volume is required for calculating rlist from the "
443 "energy drift with verlet-buffer-tolerance > 0. You are "
444 "using at least one unbounded dimension, so no volume can be "
445 "computed. Either use a finite box, or set rlist yourself "
446 "together with verlet-buffer-tolerance = -1.");
448 /* Set rlist temporarily so we can continue processing */
453 /* Set the buffer to 5% of the cut-off */
454 ir
->rlist
= (1.0 + verlet_buffer_ratio_nodynamics
) * rc_max
;
460 /* GENERAL INTEGRATOR STUFF */
463 if (ir
->etc
!= etcNO
)
465 if (EI_RANDOM(ir
->eI
))
468 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
469 "implicitly. See the documentation for more information on which "
470 "parameters affect temperature for %s.",
471 etcoupl_names
[ir
->etc
], ei_names
[ir
->eI
], ei_names
[ir
->eI
]);
476 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
478 etcoupl_names
[ir
->etc
], ei_names
[ir
->eI
]);
480 warning_note(wi
, warn_buf
);
484 if (ir
->eI
== eiVVAK
)
487 "Integrator method %s is implemented primarily for validation purposes; for "
488 "molecular dynamics, you should probably be using %s or %s",
489 ei_names
[eiVVAK
], ei_names
[eiMD
], ei_names
[eiVV
]);
490 warning_note(wi
, warn_buf
);
492 if (!EI_DYNAMICS(ir
->eI
))
494 if (ir
->epc
!= epcNO
)
497 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
498 epcoupl_names
[ir
->epc
], ei_names
[ir
->eI
]);
499 warning_note(wi
, warn_buf
);
503 if (EI_DYNAMICS(ir
->eI
))
505 if (ir
->nstcalcenergy
< 0)
507 ir
->nstcalcenergy
= ir_optimal_nstcalcenergy(ir
);
508 if (ir
->nstenergy
!= 0 && ir
->nstenergy
< ir
->nstcalcenergy
)
510 /* nstcalcenergy larger than nstener does not make sense.
511 * We ideally want nstcalcenergy=nstener.
515 ir
->nstcalcenergy
= lcd(ir
->nstenergy
, ir
->nstlist
);
519 ir
->nstcalcenergy
= ir
->nstenergy
;
523 else if ((ir
->nstenergy
> 0 && ir
->nstcalcenergy
> ir
->nstenergy
)
524 || (ir
->efep
!= efepNO
&& ir
->fepvals
->nstdhdl
> 0
525 && (ir
->nstcalcenergy
> ir
->fepvals
->nstdhdl
)))
528 const char* nsten
= "nstenergy";
529 const char* nstdh
= "nstdhdl";
530 const char* min_name
= nsten
;
531 int min_nst
= ir
->nstenergy
;
533 /* find the smallest of ( nstenergy, nstdhdl ) */
534 if (ir
->efep
!= efepNO
&& ir
->fepvals
->nstdhdl
> 0
535 && (ir
->nstenergy
== 0 || ir
->fepvals
->nstdhdl
< ir
->nstenergy
))
537 min_nst
= ir
->fepvals
->nstdhdl
;
540 /* If the user sets nstenergy small, we should respect that */
541 sprintf(warn_buf
, "Setting nstcalcenergy (%d) equal to %s (%d)", ir
->nstcalcenergy
,
543 warning_note(wi
, warn_buf
);
544 ir
->nstcalcenergy
= min_nst
;
547 if (ir
->epc
!= epcNO
)
549 if (ir
->nstpcouple
< 0)
551 ir
->nstpcouple
= ir_optimal_nstpcouple(ir
);
553 if (ir
->useMts
&& ir
->nstpcouple
% ir
->mtsLevels
.back().stepFactor
!= 0)
556 "With multiple time stepping, nstpcouple should be a mutiple of "
561 if (ir
->nstcalcenergy
> 0)
563 if (ir
->efep
!= efepNO
)
565 /* nstdhdl should be a multiple of nstcalcenergy */
566 check_nst("nstcalcenergy", ir
->nstcalcenergy
, "nstdhdl", &ir
->fepvals
->nstdhdl
, wi
);
570 /* nstexpanded should be a multiple of nstcalcenergy */
571 check_nst("nstcalcenergy", ir
->nstcalcenergy
, "nstexpanded",
572 &ir
->expandedvals
->nstexpanded
, wi
);
574 /* for storing exact averages nstenergy should be
575 * a multiple of nstcalcenergy
577 check_nst("nstcalcenergy", ir
->nstcalcenergy
, "nstenergy", &ir
->nstenergy
, wi
);
580 // Inquire all MdModules, if their parameters match with the energy
581 // calculation frequency
582 gmx::EnergyCalculationFrequencyErrors
energyCalculationFrequencyErrors(ir
->nstcalcenergy
);
583 mdModulesNotifier
.preProcessingNotifications_
.notify(&energyCalculationFrequencyErrors
);
585 // Emit all errors from the energy calculation frequency checks
586 for (const std::string
& energyFrequencyErrorMessage
:
587 energyCalculationFrequencyErrors
.errorMessages())
589 warning_error(wi
, energyFrequencyErrorMessage
);
593 if (ir
->nsteps
== 0 && !ir
->bContinuation
)
596 "For a correct single-point energy evaluation with nsteps = 0, use "
597 "continuation = yes to avoid constraining the input coordinates.");
601 if ((EI_SD(ir
->eI
) || ir
->eI
== eiBD
) && ir
->bContinuation
&& ir
->ld_seed
!= -1)
604 "You are doing a continuation with SD or BD, make sure that ld_seed is "
605 "different from the previous run (using ld_seed=-1 will ensure this)");
611 sprintf(err_buf
, "TPI only works with pbc = %s", c_pbcTypeNames
[PbcType::Xyz
].c_str());
612 CHECK(ir
->pbcType
!= PbcType::Xyz
);
613 sprintf(err_buf
, "with TPI nstlist should be larger than zero");
614 CHECK(ir
->nstlist
<= 0);
615 sprintf(err_buf
, "TPI does not work with full electrostatics other than PME");
616 CHECK(EEL_FULL(ir
->coulombtype
) && !EEL_PME(ir
->coulombtype
));
620 if ((opts
->nshake
> 0) && (opts
->bMorse
))
622 sprintf(warn_buf
, "Using morse bond-potentials while constraining bonds is useless");
623 warning(wi
, warn_buf
);
626 if ((EI_SD(ir
->eI
) || ir
->eI
== eiBD
) && ir
->bContinuation
&& ir
->ld_seed
!= -1)
629 "You are doing a continuation with SD or BD, make sure that ld_seed is "
630 "different from the previous run (using ld_seed=-1 will ensure this)");
632 /* verify simulated tempering options */
636 bool bAllTempZero
= TRUE
;
637 for (i
= 0; i
< fep
->n_lambda
; i
++)
639 sprintf(err_buf
, "Entry %d for %s must be between 0 and 1, instead is %g", i
,
640 efpt_names
[efptTEMPERATURE
], fep
->all_lambda
[efptTEMPERATURE
][i
]);
641 CHECK((fep
->all_lambda
[efptTEMPERATURE
][i
] < 0) || (fep
->all_lambda
[efptTEMPERATURE
][i
] > 1));
642 if (fep
->all_lambda
[efptTEMPERATURE
][i
] > 0)
644 bAllTempZero
= FALSE
;
647 sprintf(err_buf
, "if simulated tempering is on, temperature-lambdas may not be all zero");
648 CHECK(bAllTempZero
== TRUE
);
650 sprintf(err_buf
, "Simulated tempering is currently only compatible with md-vv");
651 CHECK(ir
->eI
!= eiVV
);
653 /* check compatability of the temperature coupling with simulated tempering */
655 if (ir
->etc
== etcNOSEHOOVER
)
658 "Nose-Hoover based temperature control such as [%s] my not be "
659 "entirelyconsistent with simulated tempering",
660 etcoupl_names
[ir
->etc
]);
661 warning_note(wi
, warn_buf
);
664 /* check that the temperatures make sense */
667 "Higher simulated tempering temperature (%g) must be >= than the simulated "
668 "tempering lower temperature (%g)",
669 ir
->simtempvals
->simtemp_high
, ir
->simtempvals
->simtemp_low
);
670 CHECK(ir
->simtempvals
->simtemp_high
<= ir
->simtempvals
->simtemp_low
);
672 sprintf(err_buf
, "Higher simulated tempering temperature (%g) must be >= zero",
673 ir
->simtempvals
->simtemp_high
);
674 CHECK(ir
->simtempvals
->simtemp_high
<= 0);
676 sprintf(err_buf
, "Lower simulated tempering temperature (%g) must be >= zero",
677 ir
->simtempvals
->simtemp_low
);
678 CHECK(ir
->simtempvals
->simtemp_low
<= 0);
681 /* verify free energy options */
683 if (ir
->efep
!= efepNO
)
686 sprintf(err_buf
, "The soft-core power is %d and can only be 1 or 2", fep
->sc_power
);
687 CHECK(fep
->sc_alpha
!= 0 && fep
->sc_power
!= 1 && fep
->sc_power
!= 2);
690 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
692 static_cast<int>(fep
->sc_r_power
));
693 CHECK(fep
->sc_alpha
!= 0 && fep
->sc_r_power
!= 6.0);
696 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
699 CHECK(fep
->delta_lambda
> 0 && ((fep
->init_fep_state
> 0) || (fep
->init_lambda
> 0)));
701 sprintf(err_buf
, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
703 CHECK(fep
->delta_lambda
> 0 && (ir
->efep
== efepEXPANDED
));
705 sprintf(err_buf
, "Can only use expanded ensemble with md-vv (for now)");
706 CHECK(!(EI_VV(ir
->eI
)) && (ir
->efep
== efepEXPANDED
));
708 sprintf(err_buf
, "Free-energy not implemented for Ewald");
709 CHECK(ir
->coulombtype
== eelEWALD
);
711 /* check validty of lambda inputs */
712 if (fep
->n_lambda
== 0)
714 /* Clear output in case of no states:*/
715 sprintf(err_buf
, "init-lambda-state set to %d: no lambda states are defined.",
716 fep
->init_fep_state
);
717 CHECK((fep
->init_fep_state
>= 0) && (fep
->n_lambda
== 0));
721 sprintf(err_buf
, "initial thermodynamic state %d does not exist, only goes to %d",
722 fep
->init_fep_state
, fep
->n_lambda
- 1);
723 CHECK((fep
->init_fep_state
>= fep
->n_lambda
));
727 "Lambda state must be set, either with init-lambda-state or with init-lambda");
728 CHECK((fep
->init_fep_state
< 0) && (fep
->init_lambda
< 0));
731 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
732 "init-lambda-state or with init-lambda, but not both",
733 fep
->init_lambda
, fep
->init_fep_state
);
734 CHECK((fep
->init_fep_state
>= 0) && (fep
->init_lambda
>= 0));
737 if ((fep
->init_lambda
>= 0) && (fep
->delta_lambda
== 0))
741 for (i
= 0; i
< efptNR
; i
++)
743 if (fep
->separate_dvdl
[i
])
748 if (n_lambda_terms
> 1)
751 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
752 "use init-lambda to set lambda state (except for slow growth). Use "
753 "init-lambda-state instead.");
754 warning(wi
, warn_buf
);
757 if (n_lambda_terms
< 2 && fep
->n_lambda
> 0)
760 "init-lambda is deprecated for setting lambda state (except for slow "
761 "growth). Use init-lambda-state instead.");
765 for (j
= 0; j
< efptNR
; j
++)
767 for (i
= 0; i
< fep
->n_lambda
; i
++)
769 sprintf(err_buf
, "Entry %d for %s must be between 0 and 1, instead is %g", i
,
770 efpt_names
[j
], fep
->all_lambda
[j
][i
]);
771 CHECK((fep
->all_lambda
[j
][i
] < 0) || (fep
->all_lambda
[j
][i
] > 1));
775 if ((fep
->sc_alpha
> 0) && (!fep
->bScCoul
))
777 for (i
= 0; i
< fep
->n_lambda
; i
++)
780 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
781 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
782 "crashes, and is not supported.",
783 i
, fep
->all_lambda
[efptVDW
][i
], fep
->all_lambda
[efptCOUL
][i
]);
784 CHECK((fep
->sc_alpha
> 0)
785 && (((fep
->all_lambda
[efptCOUL
][i
] > 0.0) && (fep
->all_lambda
[efptCOUL
][i
] < 1.0))
786 && ((fep
->all_lambda
[efptVDW
][i
] > 0.0) && (fep
->all_lambda
[efptVDW
][i
] < 1.0))));
790 if ((fep
->bScCoul
) && (EEL_PME(ir
->coulombtype
)))
792 real sigma
, lambda
, r_sc
;
795 /* Maximum estimate for A and B charges equal with lambda power 1 */
797 r_sc
= std::pow(lambda
* fep
->sc_alpha
* std::pow(sigma
/ ir
->rcoulomb
, fep
->sc_r_power
) + 1.0,
798 1.0 / fep
->sc_r_power
);
800 "With PME there is a minor soft core effect present at the cut-off, "
801 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
802 "energy conservation, but usually other effects dominate. With a common sigma "
803 "value of %g nm the fraction of the particle-particle potential at the cut-off "
804 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
805 fep
->sc_r_power
, sigma
, lambda
, r_sc
- 1.0, ir
->ewald_rtol
);
806 warning_note(wi
, warn_buf
);
809 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
810 be treated differently, but that's the next step */
812 for (i
= 0; i
< efptNR
; i
++)
814 for (j
= 0; j
< fep
->n_lambda
; j
++)
816 sprintf(err_buf
, "%s[%d] must be between 0 and 1", efpt_names
[i
], j
);
817 CHECK((fep
->all_lambda
[i
][j
] < 0) || (fep
->all_lambda
[i
][j
] > 1));
822 if ((ir
->bSimTemp
) || (ir
->efep
== efepEXPANDED
))
826 /* checking equilibration of weights inputs for validity */
829 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
831 expand
->equil_n_at_lam
, elmceq_names
[elmceqNUMATLAM
]);
832 CHECK((expand
->equil_n_at_lam
> 0) && (expand
->elmceq
!= elmceqNUMATLAM
));
835 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
837 expand
->equil_samples
, elmceq_names
[elmceqSAMPLES
]);
838 CHECK((expand
->equil_samples
> 0) && (expand
->elmceq
!= elmceqSAMPLES
));
841 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
842 expand
->equil_steps
, elmceq_names
[elmceqSTEPS
]);
843 CHECK((expand
->equil_steps
> 0) && (expand
->elmceq
!= elmceqSTEPS
));
846 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
847 expand
->equil_samples
, elmceq_names
[elmceqWLDELTA
]);
848 CHECK((expand
->equil_wl_delta
> 0) && (expand
->elmceq
!= elmceqWLDELTA
));
851 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
852 expand
->equil_ratio
, elmceq_names
[elmceqRATIO
]);
853 CHECK((expand
->equil_ratio
> 0) && (expand
->elmceq
!= elmceqRATIO
));
856 "weight-equil-number-all-lambda (%d) must be a positive integer if "
857 "lmc-weights-equil=%s",
858 expand
->equil_n_at_lam
, elmceq_names
[elmceqNUMATLAM
]);
859 CHECK((expand
->equil_n_at_lam
<= 0) && (expand
->elmceq
== elmceqNUMATLAM
));
862 "weight-equil-number-samples (%d) must be a positive integer if "
863 "lmc-weights-equil=%s",
864 expand
->equil_samples
, elmceq_names
[elmceqSAMPLES
]);
865 CHECK((expand
->equil_samples
<= 0) && (expand
->elmceq
== elmceqSAMPLES
));
868 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
869 expand
->equil_steps
, elmceq_names
[elmceqSTEPS
]);
870 CHECK((expand
->equil_steps
<= 0) && (expand
->elmceq
== elmceqSTEPS
));
872 sprintf(err_buf
, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
873 expand
->equil_wl_delta
, elmceq_names
[elmceqWLDELTA
]);
874 CHECK((expand
->equil_wl_delta
<= 0) && (expand
->elmceq
== elmceqWLDELTA
));
876 sprintf(err_buf
, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
877 expand
->equil_ratio
, elmceq_names
[elmceqRATIO
]);
878 CHECK((expand
->equil_ratio
<= 0) && (expand
->elmceq
== elmceqRATIO
));
880 sprintf(err_buf
, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
881 elmceq_names
[elmceqWLDELTA
], elamstats_names
[elamstatsWL
], elamstats_names
[elamstatsWWL
]);
882 CHECK((expand
->elmceq
== elmceqWLDELTA
) && (!EWL(expand
->elamstats
)));
884 sprintf(err_buf
, "lmc-repeats (%d) must be greater than 0", expand
->lmc_repeats
);
885 CHECK((expand
->lmc_repeats
<= 0));
886 sprintf(err_buf
, "minimum-var-min (%d) must be greater than 0", expand
->minvarmin
);
887 CHECK((expand
->minvarmin
<= 0));
888 sprintf(err_buf
, "weight-c-range (%d) must be greater or equal to 0", expand
->c_range
);
889 CHECK((expand
->c_range
< 0));
891 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
893 fep
->init_fep_state
, expand
->lmc_forced_nstart
);
894 CHECK((fep
->init_fep_state
!= 0) && (expand
->lmc_forced_nstart
> 0)
895 && (expand
->elmcmove
!= elmcmoveNO
));
896 sprintf(err_buf
, "lmc-forced-nstart (%d) must not be negative", expand
->lmc_forced_nstart
);
897 CHECK((expand
->lmc_forced_nstart
< 0));
898 sprintf(err_buf
, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
899 fep
->init_fep_state
);
900 CHECK((fep
->init_fep_state
< 0) || (fep
->init_fep_state
>= fep
->n_lambda
));
902 sprintf(err_buf
, "init-wl-delta (%f) must be greater than or equal to 0", expand
->init_wl_delta
);
903 CHECK((expand
->init_wl_delta
< 0));
904 sprintf(err_buf
, "wl-ratio (%f) must be between 0 and 1", expand
->wl_ratio
);
905 CHECK((expand
->wl_ratio
<= 0) || (expand
->wl_ratio
>= 1));
906 sprintf(err_buf
, "wl-scale (%f) must be between 0 and 1", expand
->wl_scale
);
907 CHECK((expand
->wl_scale
<= 0) || (expand
->wl_scale
>= 1));
909 /* if there is no temperature control, we need to specify an MC temperature */
910 if (!integratorHasReferenceTemperature(ir
) && (expand
->elmcmove
!= elmcmoveNO
)
911 && (expand
->mc_temp
<= 0.0))
914 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
915 "to a positive number");
916 warning_error(wi
, err_buf
);
918 if (expand
->nstTij
> 0)
920 sprintf(err_buf
, "nstlog must be non-zero");
921 CHECK(ir
->nstlog
== 0);
922 // Avoid modulus by zero in the case that already triggered an error exit.
926 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
927 expand
->nstTij
, ir
->nstlog
);
928 CHECK((expand
->nstTij
% ir
->nstlog
) != 0);
934 sprintf(err_buf
, "walls only work with pbc=%s", c_pbcTypeNames
[PbcType::XY
].c_str());
935 CHECK(ir
->nwall
&& ir
->pbcType
!= PbcType::XY
);
938 if (ir
->pbcType
!= PbcType::Xyz
&& ir
->nwall
!= 2)
940 if (ir
->pbcType
== PbcType::No
)
942 if (ir
->epc
!= epcNO
)
944 warning(wi
, "Turning off pressure coupling for vacuum system");
950 sprintf(err_buf
, "Can not have pressure coupling with pbc=%s",
951 c_pbcTypeNames
[ir
->pbcType
].c_str());
952 CHECK(ir
->epc
!= epcNO
);
954 sprintf(err_buf
, "Can not have Ewald with pbc=%s", c_pbcTypeNames
[ir
->pbcType
].c_str());
955 CHECK(EEL_FULL(ir
->coulombtype
));
957 sprintf(err_buf
, "Can not have dispersion correction with pbc=%s",
958 c_pbcTypeNames
[ir
->pbcType
].c_str());
959 CHECK(ir
->eDispCorr
!= edispcNO
);
962 if (ir
->rlist
== 0.0)
965 "can only have neighborlist cut-off zero (=infinite)\n"
966 "with coulombtype = %s or coulombtype = %s\n"
967 "without periodic boundary conditions (pbc = %s) and\n"
968 "rcoulomb and rvdw set to zero",
969 eel_names
[eelCUT
], eel_names
[eelUSER
], c_pbcTypeNames
[PbcType::No
].c_str());
970 CHECK(((ir
->coulombtype
!= eelCUT
) && (ir
->coulombtype
!= eelUSER
))
971 || (ir
->pbcType
!= PbcType::No
) || (ir
->rcoulomb
!= 0.0) || (ir
->rvdw
!= 0.0));
976 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
977 "nstype=simple and only one MPI rank");
982 if (ir
->nstcomm
== 0)
984 // TODO Change this behaviour. There should be exactly one way
985 // to turn off an algorithm.
986 ir
->comm_mode
= ecmNO
;
988 if (ir
->comm_mode
!= ecmNO
)
992 // TODO Such input was once valid. Now that we've been
993 // helpful for a few years, we should reject such input,
994 // lest we have to support every historical decision
997 "If you want to remove the rotation around the center of mass, you should set "
998 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
999 "its absolute value");
1000 ir
->nstcomm
= abs(ir
->nstcomm
);
1003 if (ir
->nstcalcenergy
> 0 && ir
->nstcomm
< ir
->nstcalcenergy
)
1006 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1007 "nstcomm to nstcalcenergy");
1008 ir
->nstcomm
= ir
->nstcalcenergy
;
1011 if (ir
->comm_mode
== ecmANGULAR
)
1014 "Can not remove the rotation around the center of mass with periodic "
1016 CHECK(ir
->bPeriodicMols
);
1017 if (ir
->pbcType
!= PbcType::No
)
1020 "Removing the rotation around the center of mass in a periodic system, "
1021 "this can lead to artifacts. Only use this on a single (cluster of) "
1022 "molecules. This cluster should not cross periodic boundaries.");
1027 if (EI_STATE_VELOCITY(ir
->eI
) && !EI_SD(ir
->eI
) && ir
->pbcType
== PbcType::No
&& ir
->comm_mode
!= ecmANGULAR
)
1030 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1031 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1034 warning_note(wi
, warn_buf
);
1037 /* TEMPERATURE COUPLING */
1038 if (ir
->etc
== etcYES
)
1040 ir
->etc
= etcBERENDSEN
;
1042 "Old option for temperature coupling given: "
1043 "changing \"yes\" to \"Berendsen\"\n");
1046 if ((ir
->etc
== etcNOSEHOOVER
) || (ir
->epc
== epcMTTK
))
1048 if (ir
->opts
.nhchainlength
< 1)
1051 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1053 ir
->opts
.nhchainlength
);
1054 ir
->opts
.nhchainlength
= 1;
1055 warning(wi
, warn_buf
);
1058 if (ir
->etc
== etcNOSEHOOVER
&& !EI_VV(ir
->eI
) && ir
->opts
.nhchainlength
> 1)
1062 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1063 ir
->opts
.nhchainlength
= 1;
1068 ir
->opts
.nhchainlength
= 0;
1071 if (ir
->eI
== eiVVAK
)
1074 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1077 CHECK((ir
->nsttcouple
!= 1) || (ir
->nstpcouple
!= 1));
1080 if (ETC_ANDERSEN(ir
->etc
))
1082 sprintf(err_buf
, "%s temperature control not supported for integrator %s.",
1083 etcoupl_names
[ir
->etc
], ei_names
[ir
->eI
]);
1084 CHECK(!(EI_VV(ir
->eI
)));
1086 if (ir
->nstcomm
> 0 && (ir
->etc
== etcANDERSEN
))
1089 "Center of mass removal not necessary for %s. All velocities of coupled "
1090 "groups are rerandomized periodically, so flying ice cube errors will not "
1092 etcoupl_names
[ir
->etc
]);
1093 warning_note(wi
, warn_buf
);
1097 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1098 "randomized every time step",
1099 ir
->nstcomm
, etcoupl_names
[ir
->etc
]);
1100 CHECK(ir
->nstcomm
> 1 && (ir
->etc
== etcANDERSEN
));
1103 if (ir
->etc
== etcBERENDSEN
)
1106 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1107 "might want to consider using the %s thermostat.",
1108 ETCOUPLTYPE(ir
->etc
), ETCOUPLTYPE(etcVRESCALE
));
1109 warning_note(wi
, warn_buf
);
1112 if ((ir
->etc
== etcNOSEHOOVER
|| ETC_ANDERSEN(ir
->etc
)) && ir
->epc
== epcBERENDSEN
)
1115 "Using Berendsen pressure coupling invalidates the "
1116 "true ensemble for the thermostat");
1117 warning(wi
, warn_buf
);
1120 /* PRESSURE COUPLING */
1121 if (ir
->epc
== epcISOTROPIC
)
1123 ir
->epc
= epcBERENDSEN
;
1125 "Old option for pressure coupling given: "
1126 "changing \"Isotropic\" to \"Berendsen\"\n");
1129 if (ir
->epc
!= epcNO
)
1131 dt_pcoupl
= ir
->nstpcouple
* ir
->delta_t
;
1133 sprintf(err_buf
, "tau-p must be > 0 instead of %g\n", ir
->tau_p
);
1134 CHECK(ir
->tau_p
<= 0);
1136 if (ir
->tau_p
/ dt_pcoupl
< pcouple_min_integration_steps(ir
->epc
) - 10 * GMX_REAL_EPS
)
1139 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1140 "times larger than nstpcouple*dt (%g)",
1141 EPCOUPLTYPE(ir
->epc
), ir
->tau_p
, pcouple_min_integration_steps(ir
->epc
), dt_pcoupl
);
1142 warning(wi
, warn_buf
);
1146 "compressibility must be > 0 when using pressure"
1148 EPCOUPLTYPE(ir
->epc
));
1149 CHECK(ir
->compress
[XX
][XX
] < 0 || ir
->compress
[YY
][YY
] < 0 || ir
->compress
[ZZ
][ZZ
] < 0
1150 || (trace(ir
->compress
) == 0 && ir
->compress
[YY
][XX
] <= 0 && ir
->compress
[ZZ
][XX
] <= 0
1151 && ir
->compress
[ZZ
][YY
] <= 0));
1153 if (epcPARRINELLORAHMAN
== ir
->epc
&& opts
->bGenVel
)
1156 "You are generating velocities so I am assuming you "
1157 "are equilibrating a system. You are using "
1158 "%s pressure coupling, but this can be "
1159 "unstable for equilibration. If your system crashes, try "
1160 "equilibrating first with Berendsen pressure coupling. If "
1161 "you are not equilibrating the system, you can probably "
1162 "ignore this warning.",
1163 epcoupl_names
[ir
->epc
]);
1164 warning(wi
, warn_buf
);
1170 if (ir
->epc
== epcMTTK
)
1172 warning_error(wi
, "MTTK pressure coupling requires a Velocity-verlet integrator");
1176 /* ELECTROSTATICS */
1177 /* More checks are in triple check (grompp.c) */
1179 if (ir
->coulombtype
== eelSWITCH
)
1182 "coulombtype = %s is only for testing purposes and can lead to serious "
1183 "artifacts, advice: use coulombtype = %s",
1184 eel_names
[ir
->coulombtype
], eel_names
[eelRF_ZERO
]);
1185 warning(wi
, warn_buf
);
1188 if (EEL_RF(ir
->coulombtype
) && ir
->epsilon_rf
== 1 && ir
->epsilon_r
!= 1)
1191 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1192 "format and exchanging epsilon-r and epsilon-rf",
1194 warning(wi
, warn_buf
);
1195 ir
->epsilon_rf
= ir
->epsilon_r
;
1196 ir
->epsilon_r
= 1.0;
1199 if (ir
->epsilon_r
== 0)
1202 "It is pointless to use long-range electrostatics with infinite relative "
1204 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1206 CHECK(EEL_FULL(ir
->coulombtype
));
1209 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1211 sprintf(err_buf
, "epsilon-r must be >= 0 instead of %g\n", ir
->epsilon_r
);
1212 CHECK(ir
->epsilon_r
< 0);
1215 if (EEL_RF(ir
->coulombtype
))
1217 /* reaction field (at the cut-off) */
1219 if (ir
->coulombtype
== eelRF_ZERO
&& ir
->epsilon_rf
!= 0)
1222 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1223 eel_names
[ir
->coulombtype
]);
1224 warning(wi
, warn_buf
);
1225 ir
->epsilon_rf
= 0.0;
1228 sprintf(err_buf
, "epsilon-rf must be >= epsilon-r");
1229 CHECK((ir
->epsilon_rf
< ir
->epsilon_r
&& ir
->epsilon_rf
!= 0) || (ir
->epsilon_r
== 0));
1230 if (ir
->epsilon_rf
== ir
->epsilon_r
)
1232 sprintf(warn_buf
, "Using epsilon-rf = epsilon-r with %s does not make sense",
1233 eel_names
[ir
->coulombtype
]);
1234 warning(wi
, warn_buf
);
1237 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1238 * means the interaction is zero outside rcoulomb, but it helps to
1239 * provide accurate energy conservation.
1241 if (ir_coulomb_might_be_zero_at_cutoff(ir
))
1243 if (ir_coulomb_switched(ir
))
1246 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1247 "potential modifier options!",
1248 eel_names
[ir
->coulombtype
]);
1249 CHECK(ir
->rcoulomb_switch
>= ir
->rcoulomb
);
1253 if (ir
->coulombtype
== eelSWITCH
|| ir
->coulombtype
== eelSHIFT
)
1256 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1257 "secondary coulomb-modifier.");
1258 CHECK(ir
->coulomb_modifier
!= eintmodNONE
);
1260 if (ir
->vdwtype
== evdwSWITCH
|| ir
->vdwtype
== evdwSHIFT
)
1263 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1264 "secondary vdw-modifier.");
1265 CHECK(ir
->vdw_modifier
!= eintmodNONE
);
1268 if (ir
->coulombtype
== eelSWITCH
|| ir
->coulombtype
== eelSHIFT
|| ir
->vdwtype
== evdwSWITCH
1269 || ir
->vdwtype
== evdwSHIFT
)
1272 "The switch/shift interaction settings are just for compatibility; you will get "
1274 "performance from applying potential modifiers to your interactions!\n");
1275 warning_note(wi
, warn_buf
);
1278 if (ir
->coulombtype
== eelPMESWITCH
|| ir
->coulomb_modifier
== eintmodPOTSWITCH
)
1280 if (ir
->rcoulomb_switch
/ ir
->rcoulomb
< 0.9499)
1282 real percentage
= 100 * (ir
->rcoulomb
- ir
->rcoulomb_switch
) / ir
->rcoulomb
;
1284 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1285 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1286 "will be good regardless, since ewald_rtol = %g.",
1287 percentage
, ir
->rcoulomb_switch
, ir
->rcoulomb
, ir
->ewald_rtol
);
1288 warning(wi
, warn_buf
);
1292 if (ir
->vdwtype
== evdwSWITCH
|| ir
->vdw_modifier
== eintmodPOTSWITCH
)
1294 if (ir
->rvdw_switch
== 0)
1297 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1298 "potential. This suggests it was not set in the mdp, which can lead to large "
1299 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1300 "switching range.");
1301 warning(wi
, warn_buf
);
1305 if (EEL_FULL(ir
->coulombtype
))
1307 if (ir
->coulombtype
== eelPMESWITCH
|| ir
->coulombtype
== eelPMEUSER
1308 || ir
->coulombtype
== eelPMEUSERSWITCH
)
1310 sprintf(err_buf
, "With coulombtype = %s, rcoulomb must be <= rlist",
1311 eel_names
[ir
->coulombtype
]);
1312 CHECK(ir
->rcoulomb
> ir
->rlist
);
1316 if (EEL_PME(ir
->coulombtype
) || EVDW_PME(ir
->vdwtype
))
1318 // TODO: Move these checks into the ewald module with the options class
1320 int orderMax
= (ir
->coulombtype
== eelP3M_AD
? 8 : 12);
1322 if (ir
->pme_order
< orderMin
|| ir
->pme_order
> orderMax
)
1324 sprintf(warn_buf
, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1325 eel_names
[ir
->coulombtype
], orderMin
, orderMax
);
1326 warning_error(wi
, warn_buf
);
1330 if (ir
->nwall
== 2 && EEL_FULL(ir
->coulombtype
))
1332 if (ir
->ewald_geometry
== eewg3D
)
1334 sprintf(warn_buf
, "With pbc=%s you should use ewald-geometry=%s",
1335 c_pbcTypeNames
[ir
->pbcType
].c_str(), eewg_names
[eewg3DC
]);
1336 warning(wi
, warn_buf
);
1338 /* This check avoids extra pbc coding for exclusion corrections */
1339 sprintf(err_buf
, "wall-ewald-zfac should be >= 2");
1340 CHECK(ir
->wall_ewald_zfac
< 2);
1342 if ((ir
->ewald_geometry
== eewg3DC
) && (ir
->pbcType
!= PbcType::XY
) && EEL_FULL(ir
->coulombtype
))
1344 sprintf(warn_buf
, "With %s and ewald_geometry = %s you should use pbc = %s",
1345 eel_names
[ir
->coulombtype
], eewg_names
[eewg3DC
], c_pbcTypeNames
[PbcType::XY
].c_str());
1346 warning(wi
, warn_buf
);
1348 if ((ir
->epsilon_surface
!= 0) && EEL_FULL(ir
->coulombtype
))
1350 sprintf(err_buf
, "Cannot have periodic molecules with epsilon_surface > 0");
1351 CHECK(ir
->bPeriodicMols
);
1352 sprintf(warn_buf
, "With epsilon_surface > 0 all molecules should be neutral.");
1353 warning_note(wi
, warn_buf
);
1355 "With epsilon_surface > 0 you can only use domain decomposition "
1356 "when there are only small molecules with all bonds constrained (mdrun will check "
1358 warning_note(wi
, warn_buf
);
1361 if (ir_vdw_switched(ir
))
1363 sprintf(err_buf
, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1364 CHECK(ir
->rvdw_switch
>= ir
->rvdw
);
1366 if (ir
->rvdw_switch
< 0.5 * ir
->rvdw
)
1369 "You are applying a switch function to vdw forces or potentials from %g to %g "
1370 "nm, which is more than half the interaction range, whereas switch functions "
1371 "are intended to act only close to the cut-off.",
1372 ir
->rvdw_switch
, ir
->rvdw
);
1373 warning_note(wi
, warn_buf
);
1377 if (ir
->vdwtype
== evdwPME
)
1379 if (!(ir
->vdw_modifier
== eintmodNONE
|| ir
->vdw_modifier
== eintmodPOTSHIFT
))
1381 sprintf(err_buf
, "With vdwtype = %s, the only supported modifiers are %s and %s",
1382 evdw_names
[ir
->vdwtype
], eintmod_names
[eintmodPOTSHIFT
], eintmod_names
[eintmodNONE
]);
1383 warning_error(wi
, err_buf
);
1387 if (ir
->vdwtype
== evdwUSER
&& ir
->eDispCorr
!= edispcNO
)
1390 "You have selected user tables with dispersion correction, the dispersion "
1391 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1392 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1393 "really want dispersion correction to -C6/r^6.");
1396 if (ir
->eI
== eiLBFGS
&& (ir
->coulombtype
== eelCUT
|| ir
->vdwtype
== evdwCUT
) && ir
->rvdw
!= 0)
1398 warning(wi
, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1401 if (ir
->eI
== eiLBFGS
&& ir
->nbfgscorr
<= 0)
1403 warning(wi
, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1406 /* IMPLICIT SOLVENT */
1407 if (ir
->coulombtype
== eelGB_NOTUSED
)
1409 sprintf(warn_buf
, "Invalid option %s for coulombtype", eel_names
[ir
->coulombtype
]);
1410 warning_error(wi
, warn_buf
);
1415 warning_error(wi
, "The QMMM integration you are trying to use is no longer supported");
1420 gmx_fatal(FARGS
, "AdResS simulations are no longer supported");
1424 /* interpret a number of doubles from a string and put them in an array,
1425 after allocating space for them.
1426 str = the input string
1427 n = the (pre-allocated) number of doubles read
1428 r = the output array of doubles. */
1429 static void parse_n_real(char* str
, int* n
, real
** r
, warninp_t wi
)
1431 auto values
= gmx::splitString(str
);
1435 for (int i
= 0; i
< *n
; i
++)
1439 (*r
)[i
] = gmx::fromString
<real
>(values
[i
]);
1441 catch (gmx::GromacsException
&)
1443 warning_error(wi
, "Invalid value " + values
[i
]
1444 + " in string in mdp file. Expected a real number.");
1450 static void do_fep_params(t_inputrec
* ir
, char fep_lambda
[][STRLEN
], char weights
[STRLEN
], warninp_t wi
)
1453 int i
, j
, max_n_lambda
, nweights
, nfep
[efptNR
];
1454 t_lambda
* fep
= ir
->fepvals
;
1455 t_expanded
* expand
= ir
->expandedvals
;
1456 real
** count_fep_lambdas
;
1457 bool bOneLambda
= TRUE
;
1459 snew(count_fep_lambdas
, efptNR
);
1461 /* FEP input processing */
1462 /* first, identify the number of lambda values for each type.
1463 All that are nonzero must have the same number */
1465 for (i
= 0; i
< efptNR
; i
++)
1467 parse_n_real(fep_lambda
[i
], &(nfep
[i
]), &(count_fep_lambdas
[i
]), wi
);
1470 /* now, determine the number of components. All must be either zero, or equal. */
1473 for (i
= 0; i
< efptNR
; i
++)
1475 if (nfep
[i
] > max_n_lambda
)
1477 max_n_lambda
= nfep
[i
]; /* here's a nonzero one. All of them
1478 must have the same number if its not zero.*/
1483 for (i
= 0; i
< efptNR
; i
++)
1487 ir
->fepvals
->separate_dvdl
[i
] = FALSE
;
1489 else if (nfep
[i
] == max_n_lambda
)
1491 if (i
!= efptTEMPERATURE
) /* we treat this differently -- not really a reason to compute
1492 the derivative with respect to the temperature currently */
1494 ir
->fepvals
->separate_dvdl
[i
] = TRUE
;
1500 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1502 nfep
[i
], efpt_names
[i
], max_n_lambda
);
1505 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1506 ir
->fepvals
->separate_dvdl
[efptTEMPERATURE
] = FALSE
;
1508 /* the number of lambdas is the number we've read in, which is either zero
1509 or the same for all */
1510 fep
->n_lambda
= max_n_lambda
;
1512 /* allocate space for the array of lambda values */
1513 snew(fep
->all_lambda
, efptNR
);
1514 /* if init_lambda is defined, we need to set lambda */
1515 if ((fep
->init_lambda
> 0) && (fep
->n_lambda
== 0))
1517 ir
->fepvals
->separate_dvdl
[efptFEP
] = TRUE
;
1519 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1520 for (i
= 0; i
< efptNR
; i
++)
1522 snew(fep
->all_lambda
[i
], fep
->n_lambda
);
1523 if (nfep
[i
] > 0) /* if it's zero, then the count_fep_lambda arrays
1526 for (j
= 0; j
< fep
->n_lambda
; j
++)
1528 fep
->all_lambda
[i
][j
] = static_cast<double>(count_fep_lambdas
[i
][j
]);
1530 sfree(count_fep_lambdas
[i
]);
1533 sfree(count_fep_lambdas
);
1535 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1536 internal bookkeeping -- for now, init_lambda */
1538 if ((nfep
[efptFEP
] == 0) && (fep
->init_lambda
>= 0))
1540 for (i
= 0; i
< fep
->n_lambda
; i
++)
1542 fep
->all_lambda
[efptFEP
][i
] = fep
->init_lambda
;
1546 /* check to see if only a single component lambda is defined, and soft core is defined.
1547 In this case, turn on coulomb soft core */
1549 if (max_n_lambda
== 0)
1555 for (i
= 0; i
< efptNR
; i
++)
1557 if ((nfep
[i
] != 0) && (i
!= efptFEP
))
1563 if ((bOneLambda
) && (fep
->sc_alpha
> 0))
1565 fep
->bScCoul
= TRUE
;
1568 /* Fill in the others with the efptFEP if they are not explicitly
1569 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1570 they are all zero. */
1572 for (i
= 0; i
< efptNR
; i
++)
1574 if ((nfep
[i
] == 0) && (i
!= efptFEP
))
1576 for (j
= 0; j
< fep
->n_lambda
; j
++)
1578 fep
->all_lambda
[i
][j
] = fep
->all_lambda
[efptFEP
][j
];
1584 /* now read in the weights */
1585 parse_n_real(weights
, &nweights
, &(expand
->init_lambda_weights
), wi
);
1588 snew(expand
->init_lambda_weights
, fep
->n_lambda
); /* initialize to zero */
1590 else if (nweights
!= fep
->n_lambda
)
1592 gmx_fatal(FARGS
, "Number of weights (%d) is not equal to number of lambda values (%d)",
1593 nweights
, fep
->n_lambda
);
1595 if ((expand
->nstexpanded
< 0) && (ir
->efep
!= efepNO
))
1597 expand
->nstexpanded
= fep
->nstdhdl
;
1598 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1603 static void do_simtemp_params(t_inputrec
* ir
)
1606 snew(ir
->simtempvals
->temperatures
, ir
->fepvals
->n_lambda
);
1607 GetSimTemps(ir
->fepvals
->n_lambda
, ir
->simtempvals
, ir
->fepvals
->all_lambda
[efptTEMPERATURE
]);
1610 template<typename T
>
1611 void convertInts(warninp_t wi
, gmx::ArrayRef
<const std::string
> inputs
, const char* name
, T
* outputs
)
1614 for (const auto& input
: inputs
)
1618 outputs
[i
] = gmx::fromStdString
<T
>(input
);
1620 catch (gmx::GromacsException
&)
1622 auto message
= gmx::formatString(
1623 "Invalid value for mdp option %s. %s should only consist of integers separated "
1626 warning_error(wi
, message
);
1632 static void convertReals(warninp_t wi
, gmx::ArrayRef
<const std::string
> inputs
, const char* name
, real
* outputs
)
1635 for (const auto& input
: inputs
)
1639 outputs
[i
] = gmx::fromString
<real
>(input
);
1641 catch (gmx::GromacsException
&)
1643 auto message
= gmx::formatString(
1644 "Invalid value for mdp option %s. %s should only consist of real numbers "
1645 "separated by spaces.",
1647 warning_error(wi
, message
);
1653 static void convertRvecs(warninp_t wi
, gmx::ArrayRef
<const std::string
> inputs
, const char* name
, rvec
* outputs
)
1656 for (const auto& input
: inputs
)
1660 outputs
[i
][d
] = gmx::fromString
<real
>(input
);
1662 catch (gmx::GromacsException
&)
1664 auto message
= gmx::formatString(
1665 "Invalid value for mdp option %s. %s should only consist of real numbers "
1666 "separated by spaces.",
1668 warning_error(wi
, message
);
1679 static void do_wall_params(t_inputrec
* ir
, char* wall_atomtype
, char* wall_density
, t_gromppopts
* opts
, warninp_t wi
)
1681 opts
->wall_atomtype
[0] = nullptr;
1682 opts
->wall_atomtype
[1] = nullptr;
1684 ir
->wall_atomtype
[0] = -1;
1685 ir
->wall_atomtype
[1] = -1;
1686 ir
->wall_density
[0] = 0;
1687 ir
->wall_density
[1] = 0;
1691 auto wallAtomTypes
= gmx::splitString(wall_atomtype
);
1692 if (wallAtomTypes
.size() != size_t(ir
->nwall
))
1694 gmx_fatal(FARGS
, "Expected %d elements for wall_atomtype, found %zu", ir
->nwall
,
1695 wallAtomTypes
.size());
1697 GMX_RELEASE_ASSERT(ir
->nwall
< 3, "Invalid number of walls");
1698 for (int i
= 0; i
< ir
->nwall
; i
++)
1700 opts
->wall_atomtype
[i
] = gmx_strdup(wallAtomTypes
[i
].c_str());
1703 if (ir
->wall_type
== ewt93
|| ir
->wall_type
== ewt104
)
1705 auto wallDensity
= gmx::splitString(wall_density
);
1706 if (wallDensity
.size() != size_t(ir
->nwall
))
1708 gmx_fatal(FARGS
, "Expected %d elements for wall-density, found %zu", ir
->nwall
,
1709 wallDensity
.size());
1711 convertReals(wi
, wallDensity
, "wall-density", ir
->wall_density
);
1712 for (int i
= 0; i
< ir
->nwall
; i
++)
1714 if (ir
->wall_density
[i
] <= 0)
1716 gmx_fatal(FARGS
, "wall-density[%d] = %f\n", i
, ir
->wall_density
[i
]);
1723 static void add_wall_energrps(SimulationGroups
* groups
, int nwall
, t_symtab
* symtab
)
1727 AtomGroupIndices
* grps
= &(groups
->groups
[SimulationAtomGroupType::EnergyOutput
]);
1728 for (int i
= 0; i
< nwall
; i
++)
1730 groups
->groupNames
.emplace_back(put_symtab(symtab
, gmx::formatString("wall%d", i
).c_str()));
1731 grps
->emplace_back(groups
->groupNames
.size() - 1);
1736 static void read_expandedparams(std::vector
<t_inpfile
>* inp
, t_expanded
* expand
, warninp_t wi
)
1738 /* read expanded ensemble parameters */
1739 printStringNewline(inp
, "expanded ensemble variables");
1740 expand
->nstexpanded
= get_eint(inp
, "nstexpanded", -1, wi
);
1741 expand
->elamstats
= get_eeenum(inp
, "lmc-stats", elamstats_names
, wi
);
1742 expand
->elmcmove
= get_eeenum(inp
, "lmc-move", elmcmove_names
, wi
);
1743 expand
->elmceq
= get_eeenum(inp
, "lmc-weights-equil", elmceq_names
, wi
);
1744 expand
->equil_n_at_lam
= get_eint(inp
, "weight-equil-number-all-lambda", -1, wi
);
1745 expand
->equil_samples
= get_eint(inp
, "weight-equil-number-samples", -1, wi
);
1746 expand
->equil_steps
= get_eint(inp
, "weight-equil-number-steps", -1, wi
);
1747 expand
->equil_wl_delta
= get_ereal(inp
, "weight-equil-wl-delta", -1, wi
);
1748 expand
->equil_ratio
= get_ereal(inp
, "weight-equil-count-ratio", -1, wi
);
1749 printStringNewline(inp
, "Seed for Monte Carlo in lambda space");
1750 expand
->lmc_seed
= get_eint(inp
, "lmc-seed", -1, wi
);
1751 expand
->mc_temp
= get_ereal(inp
, "mc-temperature", -1, wi
);
1752 expand
->lmc_repeats
= get_eint(inp
, "lmc-repeats", 1, wi
);
1753 expand
->gibbsdeltalam
= get_eint(inp
, "lmc-gibbsdelta", -1, wi
);
1754 expand
->lmc_forced_nstart
= get_eint(inp
, "lmc-forced-nstart", 0, wi
);
1755 expand
->bSymmetrizedTMatrix
=
1756 (get_eeenum(inp
, "symmetrized-transition-matrix", yesno_names
, wi
) != 0);
1757 expand
->nstTij
= get_eint(inp
, "nst-transition-matrix", -1, wi
);
1758 expand
->minvarmin
= get_eint(inp
, "mininum-var-min", 100, wi
); /*default is reasonable */
1759 expand
->c_range
= get_eint(inp
, "weight-c-range", 0, wi
); /* default is just C=0 */
1760 expand
->wl_scale
= get_ereal(inp
, "wl-scale", 0.8, wi
);
1761 expand
->wl_ratio
= get_ereal(inp
, "wl-ratio", 0.8, wi
);
1762 expand
->init_wl_delta
= get_ereal(inp
, "init-wl-delta", 1.0, wi
);
1763 expand
->bWLoneovert
= (get_eeenum(inp
, "wl-oneovert", yesno_names
, wi
) != 0);
1766 /*! \brief Return whether an end state with the given coupling-lambda
1767 * value describes fully-interacting VDW.
1769 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1770 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1772 static bool couple_lambda_has_vdw_on(int couple_lambda_value
)
1774 return (couple_lambda_value
== ecouplamVDW
|| couple_lambda_value
== ecouplamVDWQ
);
1780 class MdpErrorHandler
: public gmx::IKeyValueTreeErrorHandler
1783 explicit MdpErrorHandler(warninp_t wi
) : wi_(wi
), mapping_(nullptr) {}
1785 void setBackMapping(const gmx::IKeyValueTreeBackMapping
& mapping
) { mapping_
= &mapping
; }
1787 bool onError(gmx::UserInputError
* ex
, const gmx::KeyValueTreePath
& context
) override
1790 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context
).c_str()));
1791 std::string message
= gmx::formatExceptionMessageToString(*ex
);
1792 warning_error(wi_
, message
.c_str());
1797 std::string
getOptionName(const gmx::KeyValueTreePath
& context
)
1799 if (mapping_
!= nullptr)
1801 gmx::KeyValueTreePath path
= mapping_
->originalPath(context
);
1802 GMX_ASSERT(path
.size() == 1, "Inconsistent mapping back to mdp options");
1805 GMX_ASSERT(context
.size() == 1, "Inconsistent context for mdp option parsing");
1810 const gmx::IKeyValueTreeBackMapping
* mapping_
;
1815 void get_ir(const char* mdparin
,
1816 const char* mdparout
,
1817 gmx::MDModules
* mdModules
,
1820 WriteMdpHeader writeMdpHeader
,
1824 double dumdub
[2][6];
1826 char warn_buf
[STRLEN
];
1827 t_lambda
* fep
= ir
->fepvals
;
1828 t_expanded
* expand
= ir
->expandedvals
;
1830 const char* no_names
[] = { "no", nullptr };
1832 init_inputrec_strings();
1833 gmx::TextInputFile
stream(mdparin
);
1834 std::vector
<t_inpfile
> inp
= read_inpfile(&stream
, mdparin
, wi
);
1836 snew(dumstr
[0], STRLEN
);
1837 snew(dumstr
[1], STRLEN
);
1839 /* ignore the following deprecated commands */
1840 replace_inp_entry(inp
, "title", nullptr);
1841 replace_inp_entry(inp
, "cpp", nullptr);
1842 replace_inp_entry(inp
, "domain-decomposition", nullptr);
1843 replace_inp_entry(inp
, "andersen-seed", nullptr);
1844 replace_inp_entry(inp
, "dihre", nullptr);
1845 replace_inp_entry(inp
, "dihre-fc", nullptr);
1846 replace_inp_entry(inp
, "dihre-tau", nullptr);
1847 replace_inp_entry(inp
, "nstdihreout", nullptr);
1848 replace_inp_entry(inp
, "nstcheckpoint", nullptr);
1849 replace_inp_entry(inp
, "optimize-fft", nullptr);
1850 replace_inp_entry(inp
, "adress_type", nullptr);
1851 replace_inp_entry(inp
, "adress_const_wf", nullptr);
1852 replace_inp_entry(inp
, "adress_ex_width", nullptr);
1853 replace_inp_entry(inp
, "adress_hy_width", nullptr);
1854 replace_inp_entry(inp
, "adress_ex_forcecap", nullptr);
1855 replace_inp_entry(inp
, "adress_interface_correction", nullptr);
1856 replace_inp_entry(inp
, "adress_site", nullptr);
1857 replace_inp_entry(inp
, "adress_reference_coords", nullptr);
1858 replace_inp_entry(inp
, "adress_tf_grp_names", nullptr);
1859 replace_inp_entry(inp
, "adress_cg_grp_names", nullptr);
1860 replace_inp_entry(inp
, "adress_do_hybridpairs", nullptr);
1861 replace_inp_entry(inp
, "rlistlong", nullptr);
1862 replace_inp_entry(inp
, "nstcalclr", nullptr);
1863 replace_inp_entry(inp
, "pull-print-com2", nullptr);
1864 replace_inp_entry(inp
, "gb-algorithm", nullptr);
1865 replace_inp_entry(inp
, "nstgbradii", nullptr);
1866 replace_inp_entry(inp
, "rgbradii", nullptr);
1867 replace_inp_entry(inp
, "gb-epsilon-solvent", nullptr);
1868 replace_inp_entry(inp
, "gb-saltconc", nullptr);
1869 replace_inp_entry(inp
, "gb-obc-alpha", nullptr);
1870 replace_inp_entry(inp
, "gb-obc-beta", nullptr);
1871 replace_inp_entry(inp
, "gb-obc-gamma", nullptr);
1872 replace_inp_entry(inp
, "gb-dielectric-offset", nullptr);
1873 replace_inp_entry(inp
, "sa-algorithm", nullptr);
1874 replace_inp_entry(inp
, "sa-surface-tension", nullptr);
1875 replace_inp_entry(inp
, "ns-type", nullptr);
1877 /* replace the following commands with the clearer new versions*/
1878 replace_inp_entry(inp
, "unconstrained-start", "continuation");
1879 replace_inp_entry(inp
, "foreign-lambda", "fep-lambdas");
1880 replace_inp_entry(inp
, "verlet-buffer-drift", "verlet-buffer-tolerance");
1881 replace_inp_entry(inp
, "nstxtcout", "nstxout-compressed");
1882 replace_inp_entry(inp
, "xtc-grps", "compressed-x-grps");
1883 replace_inp_entry(inp
, "xtc-precision", "compressed-x-precision");
1884 replace_inp_entry(inp
, "pull-print-com1", "pull-print-com");
1886 printStringNewline(&inp
, "VARIOUS PREPROCESSING OPTIONS");
1887 printStringNoNewline(&inp
, "Preprocessor information: use cpp syntax.");
1888 printStringNoNewline(&inp
, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1889 setStringEntry(&inp
, "include", opts
->include
, nullptr);
1890 printStringNoNewline(
1891 &inp
, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1892 setStringEntry(&inp
, "define", opts
->define
, nullptr);
1894 printStringNewline(&inp
, "RUN CONTROL PARAMETERS");
1895 ir
->eI
= get_eeenum(&inp
, "integrator", ei_names
, wi
);
1896 printStringNoNewline(&inp
, "Start time and timestep in ps");
1897 ir
->init_t
= get_ereal(&inp
, "tinit", 0.0, wi
);
1898 ir
->delta_t
= get_ereal(&inp
, "dt", 0.001, wi
);
1899 ir
->nsteps
= get_eint64(&inp
, "nsteps", 0, wi
);
1900 printStringNoNewline(&inp
, "For exact run continuation or redoing part of a run");
1901 ir
->init_step
= get_eint64(&inp
, "init-step", 0, wi
);
1902 printStringNoNewline(
1903 &inp
, "Part index is updated automatically on checkpointing (keeps files separate)");
1904 ir
->simulation_part
= get_eint(&inp
, "simulation-part", 1, wi
);
1905 printStringNoNewline(&inp
, "Multiple time-stepping");
1906 ir
->useMts
= (get_eeenum(&inp
, "mts", yesno_names
, wi
) != 0);
1909 gmx::GromppMtsOpts
& mtsOpts
= opts
->mtsOpts
;
1910 mtsOpts
.numLevels
= get_eint(&inp
, "mts-levels", 2, wi
);
1911 ir
->mtsLevels
.resize(2);
1912 mtsOpts
.level2Forces
= setStringEntry(&inp
, "mts-level2-forces",
1913 "longrange-nonbonded nonbonded pair dihedral");
1914 mtsOpts
.level2Factor
= get_eint(&inp
, "mts-level2-factor", 2, wi
);
1916 // We clear after reading without dynamics to not force the user to remove MTS mdp options
1917 if (!EI_DYNAMICS(ir
->eI
))
1922 printStringNoNewline(&inp
, "mode for center of mass motion removal");
1923 ir
->comm_mode
= get_eeenum(&inp
, "comm-mode", ecm_names
, wi
);
1924 printStringNoNewline(&inp
, "number of steps for center of mass motion removal");
1925 ir
->nstcomm
= get_eint(&inp
, "nstcomm", 100, wi
);
1926 printStringNoNewline(&inp
, "group(s) for center of mass motion removal");
1927 setStringEntry(&inp
, "comm-grps", inputrecStrings
->vcm
, nullptr);
1929 printStringNewline(&inp
, "LANGEVIN DYNAMICS OPTIONS");
1930 printStringNoNewline(&inp
, "Friction coefficient (amu/ps) and random seed");
1931 ir
->bd_fric
= get_ereal(&inp
, "bd-fric", 0.0, wi
);
1932 ir
->ld_seed
= get_eint64(&inp
, "ld-seed", -1, wi
);
1935 printStringNewline(&inp
, "ENERGY MINIMIZATION OPTIONS");
1936 printStringNoNewline(&inp
, "Force tolerance and initial step-size");
1937 ir
->em_tol
= get_ereal(&inp
, "emtol", 10.0, wi
);
1938 ir
->em_stepsize
= get_ereal(&inp
, "emstep", 0.01, wi
);
1939 printStringNoNewline(&inp
, "Max number of iterations in relax-shells");
1940 ir
->niter
= get_eint(&inp
, "niter", 20, wi
);
1941 printStringNoNewline(&inp
, "Step size (ps^2) for minimization of flexible constraints");
1942 ir
->fc_stepsize
= get_ereal(&inp
, "fcstep", 0, wi
);
1943 printStringNoNewline(&inp
, "Frequency of steepest descents steps when doing CG");
1944 ir
->nstcgsteep
= get_eint(&inp
, "nstcgsteep", 1000, wi
);
1945 ir
->nbfgscorr
= get_eint(&inp
, "nbfgscorr", 10, wi
);
1947 printStringNewline(&inp
, "TEST PARTICLE INSERTION OPTIONS");
1948 ir
->rtpi
= get_ereal(&inp
, "rtpi", 0.05, wi
);
1950 /* Output options */
1951 printStringNewline(&inp
, "OUTPUT CONTROL OPTIONS");
1952 printStringNoNewline(&inp
, "Output frequency for coords (x), velocities (v) and forces (f)");
1953 ir
->nstxout
= get_eint(&inp
, "nstxout", 0, wi
);
1954 ir
->nstvout
= get_eint(&inp
, "nstvout", 0, wi
);
1955 ir
->nstfout
= get_eint(&inp
, "nstfout", 0, wi
);
1956 printStringNoNewline(&inp
, "Output frequency for energies to log file and energy file");
1957 ir
->nstlog
= get_eint(&inp
, "nstlog", 1000, wi
);
1958 ir
->nstcalcenergy
= get_eint(&inp
, "nstcalcenergy", 100, wi
);
1959 ir
->nstenergy
= get_eint(&inp
, "nstenergy", 1000, wi
);
1960 printStringNoNewline(&inp
, "Output frequency and precision for .xtc file");
1961 ir
->nstxout_compressed
= get_eint(&inp
, "nstxout-compressed", 0, wi
);
1962 ir
->x_compression_precision
= get_ereal(&inp
, "compressed-x-precision", 1000.0, wi
);
1963 printStringNoNewline(&inp
, "This selects the subset of atoms for the compressed");
1964 printStringNoNewline(&inp
, "trajectory file. You can select multiple groups. By");
1965 printStringNoNewline(&inp
, "default, all atoms will be written.");
1966 setStringEntry(&inp
, "compressed-x-grps", inputrecStrings
->x_compressed_groups
, nullptr);
1967 printStringNoNewline(&inp
, "Selection of energy groups");
1968 setStringEntry(&inp
, "energygrps", inputrecStrings
->energy
, nullptr);
1970 /* Neighbor searching */
1971 printStringNewline(&inp
, "NEIGHBORSEARCHING PARAMETERS");
1972 printStringNoNewline(&inp
, "cut-off scheme (Verlet: particle based cut-offs)");
1973 ir
->cutoff_scheme
= get_eeenum(&inp
, "cutoff-scheme", ecutscheme_names
, wi
);
1974 printStringNoNewline(&inp
, "nblist update frequency");
1975 ir
->nstlist
= get_eint(&inp
, "nstlist", 10, wi
);
1976 printStringNoNewline(&inp
, "Periodic boundary conditions: xyz, no, xy");
1977 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
1978 std::vector
<const char*> pbcTypesNamesChar
;
1979 for (const auto& pbcTypeName
: c_pbcTypeNames
)
1981 pbcTypesNamesChar
.push_back(pbcTypeName
.c_str());
1983 ir
->pbcType
= static_cast<PbcType
>(get_eeenum(&inp
, "pbc", pbcTypesNamesChar
.data(), wi
));
1984 ir
->bPeriodicMols
= get_eeenum(&inp
, "periodic-molecules", yesno_names
, wi
) != 0;
1985 printStringNoNewline(&inp
,
1986 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1987 printStringNoNewline(&inp
, "a value of -1 means: use rlist");
1988 ir
->verletbuf_tol
= get_ereal(&inp
, "verlet-buffer-tolerance", 0.005, wi
);
1989 printStringNoNewline(&inp
, "nblist cut-off");
1990 ir
->rlist
= get_ereal(&inp
, "rlist", 1.0, wi
);
1991 printStringNoNewline(&inp
, "long-range cut-off for switched potentials");
1993 /* Electrostatics */
1994 printStringNewline(&inp
, "OPTIONS FOR ELECTROSTATICS AND VDW");
1995 printStringNoNewline(&inp
, "Method for doing electrostatics");
1996 ir
->coulombtype
= get_eeenum(&inp
, "coulombtype", eel_names
, wi
);
1997 ir
->coulomb_modifier
= get_eeenum(&inp
, "coulomb-modifier", eintmod_names
, wi
);
1998 printStringNoNewline(&inp
, "cut-off lengths");
1999 ir
->rcoulomb_switch
= get_ereal(&inp
, "rcoulomb-switch", 0.0, wi
);
2000 ir
->rcoulomb
= get_ereal(&inp
, "rcoulomb", 1.0, wi
);
2001 printStringNoNewline(&inp
,
2002 "Relative dielectric constant for the medium and the reaction field");
2003 ir
->epsilon_r
= get_ereal(&inp
, "epsilon-r", 1.0, wi
);
2004 ir
->epsilon_rf
= get_ereal(&inp
, "epsilon-rf", 0.0, wi
);
2005 printStringNoNewline(&inp
, "Method for doing Van der Waals");
2006 ir
->vdwtype
= get_eeenum(&inp
, "vdw-type", evdw_names
, wi
);
2007 ir
->vdw_modifier
= get_eeenum(&inp
, "vdw-modifier", eintmod_names
, wi
);
2008 printStringNoNewline(&inp
, "cut-off lengths");
2009 ir
->rvdw_switch
= get_ereal(&inp
, "rvdw-switch", 0.0, wi
);
2010 ir
->rvdw
= get_ereal(&inp
, "rvdw", 1.0, wi
);
2011 printStringNoNewline(&inp
, "Apply long range dispersion corrections for Energy and Pressure");
2012 ir
->eDispCorr
= get_eeenum(&inp
, "DispCorr", edispc_names
, wi
);
2013 printStringNoNewline(&inp
, "Extension of the potential lookup tables beyond the cut-off");
2014 ir
->tabext
= get_ereal(&inp
, "table-extension", 1.0, wi
);
2015 printStringNoNewline(&inp
, "Separate tables between energy group pairs");
2016 setStringEntry(&inp
, "energygrp-table", inputrecStrings
->egptable
, nullptr);
2017 printStringNoNewline(&inp
, "Spacing for the PME/PPPM FFT grid");
2018 ir
->fourier_spacing
= get_ereal(&inp
, "fourierspacing", 0.12, wi
);
2019 printStringNoNewline(&inp
, "FFT grid size, when a value is 0 fourierspacing will be used");
2020 ir
->nkx
= get_eint(&inp
, "fourier-nx", 0, wi
);
2021 ir
->nky
= get_eint(&inp
, "fourier-ny", 0, wi
);
2022 ir
->nkz
= get_eint(&inp
, "fourier-nz", 0, wi
);
2023 printStringNoNewline(&inp
, "EWALD/PME/PPPM parameters");
2024 ir
->pme_order
= get_eint(&inp
, "pme-order", 4, wi
);
2025 ir
->ewald_rtol
= get_ereal(&inp
, "ewald-rtol", 0.00001, wi
);
2026 ir
->ewald_rtol_lj
= get_ereal(&inp
, "ewald-rtol-lj", 0.001, wi
);
2027 ir
->ljpme_combination_rule
= get_eeenum(&inp
, "lj-pme-comb-rule", eljpme_names
, wi
);
2028 ir
->ewald_geometry
= get_eeenum(&inp
, "ewald-geometry", eewg_names
, wi
);
2029 ir
->epsilon_surface
= get_ereal(&inp
, "epsilon-surface", 0.0, wi
);
2031 /* Implicit solvation is no longer supported, but we need grompp
2032 to be able to refuse old .mdp files that would have built a tpr
2033 to run it. Thus, only "no" is accepted. */
2034 ir
->implicit_solvent
= (get_eeenum(&inp
, "implicit-solvent", no_names
, wi
) != 0);
2036 /* Coupling stuff */
2037 printStringNewline(&inp
, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2038 printStringNoNewline(&inp
, "Temperature coupling");
2039 ir
->etc
= get_eeenum(&inp
, "tcoupl", etcoupl_names
, wi
);
2040 ir
->nsttcouple
= get_eint(&inp
, "nsttcouple", -1, wi
);
2041 ir
->opts
.nhchainlength
= get_eint(&inp
, "nh-chain-length", 10, wi
);
2042 ir
->bPrintNHChains
= (get_eeenum(&inp
, "print-nose-hoover-chain-variables", yesno_names
, wi
) != 0);
2043 printStringNoNewline(&inp
, "Groups to couple separately");
2044 setStringEntry(&inp
, "tc-grps", inputrecStrings
->tcgrps
, nullptr);
2045 printStringNoNewline(&inp
, "Time constant (ps) and reference temperature (K)");
2046 setStringEntry(&inp
, "tau-t", inputrecStrings
->tau_t
, nullptr);
2047 setStringEntry(&inp
, "ref-t", inputrecStrings
->ref_t
, nullptr);
2048 printStringNoNewline(&inp
, "pressure coupling");
2049 ir
->epc
= get_eeenum(&inp
, "pcoupl", epcoupl_names
, wi
);
2050 ir
->epct
= get_eeenum(&inp
, "pcoupltype", epcoupltype_names
, wi
);
2051 ir
->nstpcouple
= get_eint(&inp
, "nstpcouple", -1, wi
);
2052 printStringNoNewline(&inp
, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2053 ir
->tau_p
= get_ereal(&inp
, "tau-p", 1.0, wi
);
2054 setStringEntry(&inp
, "compressibility", dumstr
[0], nullptr);
2055 setStringEntry(&inp
, "ref-p", dumstr
[1], nullptr);
2056 printStringNoNewline(&inp
, "Scaling of reference coordinates, No, All or COM");
2057 ir
->refcoord_scaling
= get_eeenum(&inp
, "refcoord-scaling", erefscaling_names
, wi
);
2060 printStringNewline(&inp
, "OPTIONS FOR QMMM calculations");
2061 ir
->bQMMM
= (get_eeenum(&inp
, "QMMM", yesno_names
, wi
) != 0);
2062 printStringNoNewline(&inp
, "Groups treated with MiMiC");
2063 setStringEntry(&inp
, "QMMM-grps", inputrecStrings
->QMMM
, nullptr);
2065 /* Simulated annealing */
2066 printStringNewline(&inp
, "SIMULATED ANNEALING");
2067 printStringNoNewline(&inp
, "Type of annealing for each temperature group (no/single/periodic)");
2068 setStringEntry(&inp
, "annealing", inputrecStrings
->anneal
, nullptr);
2069 printStringNoNewline(&inp
,
2070 "Number of time points to use for specifying annealing in each group");
2071 setStringEntry(&inp
, "annealing-npoints", inputrecStrings
->anneal_npoints
, nullptr);
2072 printStringNoNewline(&inp
, "List of times at the annealing points for each group");
2073 setStringEntry(&inp
, "annealing-time", inputrecStrings
->anneal_time
, nullptr);
2074 printStringNoNewline(&inp
, "Temp. at each annealing point, for each group.");
2075 setStringEntry(&inp
, "annealing-temp", inputrecStrings
->anneal_temp
, nullptr);
2078 printStringNewline(&inp
, "GENERATE VELOCITIES FOR STARTUP RUN");
2079 opts
->bGenVel
= (get_eeenum(&inp
, "gen-vel", yesno_names
, wi
) != 0);
2080 opts
->tempi
= get_ereal(&inp
, "gen-temp", 300.0, wi
);
2081 opts
->seed
= get_eint(&inp
, "gen-seed", -1, wi
);
2084 printStringNewline(&inp
, "OPTIONS FOR BONDS");
2085 opts
->nshake
= get_eeenum(&inp
, "constraints", constraints
, wi
);
2086 printStringNoNewline(&inp
, "Type of constraint algorithm");
2087 ir
->eConstrAlg
= get_eeenum(&inp
, "constraint-algorithm", econstr_names
, wi
);
2088 printStringNoNewline(&inp
, "Do not constrain the start configuration");
2089 ir
->bContinuation
= (get_eeenum(&inp
, "continuation", yesno_names
, wi
) != 0);
2090 printStringNoNewline(&inp
,
2091 "Use successive overrelaxation to reduce the number of shake iterations");
2092 ir
->bShakeSOR
= (get_eeenum(&inp
, "Shake-SOR", yesno_names
, wi
) != 0);
2093 printStringNoNewline(&inp
, "Relative tolerance of shake");
2094 ir
->shake_tol
= get_ereal(&inp
, "shake-tol", 0.0001, wi
);
2095 printStringNoNewline(&inp
, "Highest order in the expansion of the constraint coupling matrix");
2096 ir
->nProjOrder
= get_eint(&inp
, "lincs-order", 4, wi
);
2097 printStringNoNewline(&inp
, "Number of iterations in the final step of LINCS. 1 is fine for");
2098 printStringNoNewline(&inp
, "normal simulations, but use 2 to conserve energy in NVE runs.");
2099 printStringNoNewline(&inp
, "For energy minimization with constraints it should be 4 to 8.");
2100 ir
->nLincsIter
= get_eint(&inp
, "lincs-iter", 1, wi
);
2101 printStringNoNewline(&inp
, "Lincs will write a warning to the stderr if in one step a bond");
2102 printStringNoNewline(&inp
, "rotates over more degrees than");
2103 ir
->LincsWarnAngle
= get_ereal(&inp
, "lincs-warnangle", 30.0, wi
);
2104 printStringNoNewline(&inp
, "Convert harmonic bonds to morse potentials");
2105 opts
->bMorse
= (get_eeenum(&inp
, "morse", yesno_names
, wi
) != 0);
2107 /* Energy group exclusions */
2108 printStringNewline(&inp
, "ENERGY GROUP EXCLUSIONS");
2109 printStringNoNewline(
2110 &inp
, "Pairs of energy groups for which all non-bonded interactions are excluded");
2111 setStringEntry(&inp
, "energygrp-excl", inputrecStrings
->egpexcl
, nullptr);
2114 printStringNewline(&inp
, "WALLS");
2115 printStringNoNewline(
2116 &inp
, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2117 ir
->nwall
= get_eint(&inp
, "nwall", 0, wi
);
2118 ir
->wall_type
= get_eeenum(&inp
, "wall-type", ewt_names
, wi
);
2119 ir
->wall_r_linpot
= get_ereal(&inp
, "wall-r-linpot", -1, wi
);
2120 setStringEntry(&inp
, "wall-atomtype", inputrecStrings
->wall_atomtype
, nullptr);
2121 setStringEntry(&inp
, "wall-density", inputrecStrings
->wall_density
, nullptr);
2122 ir
->wall_ewald_zfac
= get_ereal(&inp
, "wall-ewald-zfac", 3, wi
);
2125 printStringNewline(&inp
, "COM PULLING");
2126 ir
->bPull
= (get_eeenum(&inp
, "pull", yesno_names
, wi
) != 0);
2129 ir
->pull
= std::make_unique
<pull_params_t
>();
2130 inputrecStrings
->pullGroupNames
= read_pullparams(&inp
, ir
->pull
.get(), wi
);
2134 for (int c
= 0; c
< ir
->pull
->ncoord
; c
++)
2136 if (ir
->pull
->coord
[c
].eType
== epullCONSTRAINT
)
2139 "Constraint COM pulling is not supported in combination with "
2140 "multiple time stepping");
2148 NOTE: needs COM pulling or free energy input */
2149 printStringNewline(&inp
, "AWH biasing");
2150 ir
->bDoAwh
= (get_eeenum(&inp
, "awh", yesno_names
, wi
) != 0);
2153 ir
->awhParams
= gmx::readAwhParams(&inp
, wi
);
2156 /* Enforced rotation */
2157 printStringNewline(&inp
, "ENFORCED ROTATION");
2158 printStringNoNewline(&inp
, "Enforced rotation: No or Yes");
2159 ir
->bRot
= (get_eeenum(&inp
, "rotation", yesno_names
, wi
) != 0);
2163 inputrecStrings
->rotateGroupNames
= read_rotparams(&inp
, ir
->rot
, wi
);
2166 /* Interactive MD */
2168 printStringNewline(&inp
, "Group to display and/or manipulate in interactive MD session");
2169 setStringEntry(&inp
, "IMD-group", inputrecStrings
->imd_grp
, nullptr);
2170 if (inputrecStrings
->imd_grp
[0] != '\0')
2177 printStringNewline(&inp
, "NMR refinement stuff");
2178 printStringNoNewline(&inp
, "Distance restraints type: No, Simple or Ensemble");
2179 ir
->eDisre
= get_eeenum(&inp
, "disre", edisre_names
, wi
);
2180 printStringNoNewline(
2181 &inp
, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2182 ir
->eDisreWeighting
= get_eeenum(&inp
, "disre-weighting", edisreweighting_names
, wi
);
2183 printStringNoNewline(&inp
, "Use sqrt of the time averaged times the instantaneous violation");
2184 ir
->bDisreMixed
= (get_eeenum(&inp
, "disre-mixed", yesno_names
, wi
) != 0);
2185 ir
->dr_fc
= get_ereal(&inp
, "disre-fc", 1000.0, wi
);
2186 ir
->dr_tau
= get_ereal(&inp
, "disre-tau", 0.0, wi
);
2187 printStringNoNewline(&inp
, "Output frequency for pair distances to energy file");
2188 ir
->nstdisreout
= get_eint(&inp
, "nstdisreout", 100, wi
);
2189 printStringNoNewline(&inp
, "Orientation restraints: No or Yes");
2190 opts
->bOrire
= (get_eeenum(&inp
, "orire", yesno_names
, wi
) != 0);
2191 printStringNoNewline(&inp
, "Orientation restraints force constant and tau for time averaging");
2192 ir
->orires_fc
= get_ereal(&inp
, "orire-fc", 0.0, wi
);
2193 ir
->orires_tau
= get_ereal(&inp
, "orire-tau", 0.0, wi
);
2194 setStringEntry(&inp
, "orire-fitgrp", inputrecStrings
->orirefitgrp
, nullptr);
2195 printStringNoNewline(&inp
, "Output frequency for trace(SD) and S to energy file");
2196 ir
->nstorireout
= get_eint(&inp
, "nstorireout", 100, wi
);
2198 /* free energy variables */
2199 printStringNewline(&inp
, "Free energy variables");
2200 ir
->efep
= get_eeenum(&inp
, "free-energy", efep_names
, wi
);
2201 setStringEntry(&inp
, "couple-moltype", inputrecStrings
->couple_moltype
, nullptr);
2202 opts
->couple_lam0
= get_eeenum(&inp
, "couple-lambda0", couple_lam
, wi
);
2203 opts
->couple_lam1
= get_eeenum(&inp
, "couple-lambda1", couple_lam
, wi
);
2204 opts
->bCoupleIntra
= (get_eeenum(&inp
, "couple-intramol", yesno_names
, wi
) != 0);
2206 fep
->init_lambda
= get_ereal(&inp
, "init-lambda", -1, wi
); /* start with -1 so
2208 it was not entered */
2209 fep
->init_fep_state
= get_eint(&inp
, "init-lambda-state", -1, wi
);
2210 fep
->delta_lambda
= get_ereal(&inp
, "delta-lambda", 0.0, wi
);
2211 fep
->nstdhdl
= get_eint(&inp
, "nstdhdl", 50, wi
);
2212 setStringEntry(&inp
, "fep-lambdas", inputrecStrings
->fep_lambda
[efptFEP
], nullptr);
2213 setStringEntry(&inp
, "mass-lambdas", inputrecStrings
->fep_lambda
[efptMASS
], nullptr);
2214 setStringEntry(&inp
, "coul-lambdas", inputrecStrings
->fep_lambda
[efptCOUL
], nullptr);
2215 setStringEntry(&inp
, "vdw-lambdas", inputrecStrings
->fep_lambda
[efptVDW
], nullptr);
2216 setStringEntry(&inp
, "bonded-lambdas", inputrecStrings
->fep_lambda
[efptBONDED
], nullptr);
2217 setStringEntry(&inp
, "restraint-lambdas", inputrecStrings
->fep_lambda
[efptRESTRAINT
], nullptr);
2218 setStringEntry(&inp
, "temperature-lambdas", inputrecStrings
->fep_lambda
[efptTEMPERATURE
], nullptr);
2219 fep
->lambda_neighbors
= get_eint(&inp
, "calc-lambda-neighbors", 1, wi
);
2220 setStringEntry(&inp
, "init-lambda-weights", inputrecStrings
->lambda_weights
, nullptr);
2221 fep
->edHdLPrintEnergy
= get_eeenum(&inp
, "dhdl-print-energy", edHdLPrintEnergy_names
, wi
);
2222 fep
->sc_alpha
= get_ereal(&inp
, "sc-alpha", 0.0, wi
);
2223 fep
->sc_power
= get_eint(&inp
, "sc-power", 1, wi
);
2224 fep
->sc_r_power
= get_ereal(&inp
, "sc-r-power", 6.0, wi
);
2225 fep
->sc_sigma
= get_ereal(&inp
, "sc-sigma", 0.3, wi
);
2226 fep
->bScCoul
= (get_eeenum(&inp
, "sc-coul", yesno_names
, wi
) != 0);
2227 fep
->dh_hist_size
= get_eint(&inp
, "dh_hist_size", 0, wi
);
2228 fep
->dh_hist_spacing
= get_ereal(&inp
, "dh_hist_spacing", 0.1, wi
);
2229 fep
->separate_dhdl_file
= get_eeenum(&inp
, "separate-dhdl-file", separate_dhdl_file_names
, wi
);
2230 fep
->dhdl_derivatives
= get_eeenum(&inp
, "dhdl-derivatives", dhdl_derivatives_names
, wi
);
2231 fep
->dh_hist_size
= get_eint(&inp
, "dh_hist_size", 0, wi
);
2232 fep
->dh_hist_spacing
= get_ereal(&inp
, "dh_hist_spacing", 0.1, wi
);
2234 /* Non-equilibrium MD stuff */
2235 printStringNewline(&inp
, "Non-equilibrium MD stuff");
2236 setStringEntry(&inp
, "acc-grps", inputrecStrings
->accgrps
, nullptr);
2237 setStringEntry(&inp
, "accelerate", inputrecStrings
->acc
, nullptr);
2238 setStringEntry(&inp
, "freezegrps", inputrecStrings
->freeze
, nullptr);
2239 setStringEntry(&inp
, "freezedim", inputrecStrings
->frdim
, nullptr);
2240 ir
->cos_accel
= get_ereal(&inp
, "cos-acceleration", 0, wi
);
2241 setStringEntry(&inp
, "deform", inputrecStrings
->deform
, nullptr);
2243 /* simulated tempering variables */
2244 printStringNewline(&inp
, "simulated tempering variables");
2245 ir
->bSimTemp
= (get_eeenum(&inp
, "simulated-tempering", yesno_names
, wi
) != 0);
2246 ir
->simtempvals
->eSimTempScale
= get_eeenum(&inp
, "simulated-tempering-scaling", esimtemp_names
, wi
);
2247 ir
->simtempvals
->simtemp_low
= get_ereal(&inp
, "sim-temp-low", 300.0, wi
);
2248 ir
->simtempvals
->simtemp_high
= get_ereal(&inp
, "sim-temp-high", 300.0, wi
);
2250 /* expanded ensemble variables */
2251 if (ir
->efep
== efepEXPANDED
|| ir
->bSimTemp
)
2253 read_expandedparams(&inp
, expand
, wi
);
2256 /* Electric fields */
2258 gmx::KeyValueTreeObject convertedValues
= flatKeyValueTreeFromInpFile(inp
);
2259 gmx::KeyValueTreeTransformer transform
;
2260 transform
.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive
);
2261 mdModules
->initMdpTransform(transform
.rules());
2262 for (const auto& path
: transform
.mappedPaths())
2264 GMX_ASSERT(path
.size() == 1, "Inconsistent mapping back to mdp options");
2265 mark_einp_set(inp
, path
[0].c_str());
2267 MdpErrorHandler
errorHandler(wi
);
2268 auto result
= transform
.transform(convertedValues
, &errorHandler
);
2269 ir
->params
= new gmx::KeyValueTreeObject(result
.object());
2270 mdModules
->adjustInputrecBasedOnModules(ir
);
2271 errorHandler
.setBackMapping(result
.backMapping());
2272 mdModules
->assignOptionsToModules(*ir
->params
, &errorHandler
);
2275 /* Ion/water position swapping ("computational electrophysiology") */
2276 printStringNewline(&inp
,
2277 "Ion/water position swapping for computational electrophysiology setups");
2278 printStringNoNewline(&inp
, "Swap positions along direction: no, X, Y, Z");
2279 ir
->eSwapCoords
= get_eeenum(&inp
, "swapcoords", eSwapTypes_names
, wi
);
2280 if (ir
->eSwapCoords
!= eswapNO
)
2287 printStringNoNewline(&inp
, "Swap attempt frequency");
2288 ir
->swap
->nstswap
= get_eint(&inp
, "swap-frequency", 1, wi
);
2289 printStringNoNewline(&inp
, "Number of ion types to be controlled");
2290 nIonTypes
= get_eint(&inp
, "iontypes", 1, wi
);
2293 warning_error(wi
, "You need to provide at least one ion type for position exchanges.");
2295 ir
->swap
->ngrp
= nIonTypes
+ eSwapFixedGrpNR
;
2296 snew(ir
->swap
->grp
, ir
->swap
->ngrp
);
2297 for (i
= 0; i
< ir
->swap
->ngrp
; i
++)
2299 snew(ir
->swap
->grp
[i
].molname
, STRLEN
);
2301 printStringNoNewline(&inp
,
2302 "Two index groups that contain the compartment-partitioning atoms");
2303 setStringEntry(&inp
, "split-group0", ir
->swap
->grp
[eGrpSplit0
].molname
, nullptr);
2304 setStringEntry(&inp
, "split-group1", ir
->swap
->grp
[eGrpSplit1
].molname
, nullptr);
2305 printStringNoNewline(&inp
,
2306 "Use center of mass of split groups (yes/no), otherwise center of "
2307 "geometry is used");
2308 ir
->swap
->massw_split
[0] = (get_eeenum(&inp
, "massw-split0", yesno_names
, wi
) != 0);
2309 ir
->swap
->massw_split
[1] = (get_eeenum(&inp
, "massw-split1", yesno_names
, wi
) != 0);
2311 printStringNoNewline(&inp
, "Name of solvent molecules");
2312 setStringEntry(&inp
, "solvent-group", ir
->swap
->grp
[eGrpSolvent
].molname
, nullptr);
2314 printStringNoNewline(&inp
,
2315 "Split cylinder: radius, upper and lower extension (nm) (this will "
2316 "define the channels)");
2317 printStringNoNewline(&inp
,
2318 "Note that the split cylinder settings do not have an influence on "
2319 "the swapping protocol,");
2320 printStringNoNewline(
2322 "however, if correctly defined, the permeation events are recorded per channel");
2323 ir
->swap
->cyl0r
= get_ereal(&inp
, "cyl0-r", 2.0, wi
);
2324 ir
->swap
->cyl0u
= get_ereal(&inp
, "cyl0-up", 1.0, wi
);
2325 ir
->swap
->cyl0l
= get_ereal(&inp
, "cyl0-down", 1.0, wi
);
2326 ir
->swap
->cyl1r
= get_ereal(&inp
, "cyl1-r", 2.0, wi
);
2327 ir
->swap
->cyl1u
= get_ereal(&inp
, "cyl1-up", 1.0, wi
);
2328 ir
->swap
->cyl1l
= get_ereal(&inp
, "cyl1-down", 1.0, wi
);
2330 printStringNoNewline(
2332 "Average the number of ions per compartment over these many swap attempt steps");
2333 ir
->swap
->nAverage
= get_eint(&inp
, "coupl-steps", 10, wi
);
2335 printStringNoNewline(
2336 &inp
, "Names of the ion types that can be exchanged with solvent molecules,");
2337 printStringNoNewline(
2338 &inp
, "and the requested number of ions of this type in compartments A and B");
2339 printStringNoNewline(&inp
, "-1 means fix the numbers as found in step 0");
2340 for (i
= 0; i
< nIonTypes
; i
++)
2342 int ig
= eSwapFixedGrpNR
+ i
;
2344 sprintf(buf
, "iontype%d-name", i
);
2345 setStringEntry(&inp
, buf
, ir
->swap
->grp
[ig
].molname
, nullptr);
2346 sprintf(buf
, "iontype%d-in-A", i
);
2347 ir
->swap
->grp
[ig
].nmolReq
[0] = get_eint(&inp
, buf
, -1, wi
);
2348 sprintf(buf
, "iontype%d-in-B", i
);
2349 ir
->swap
->grp
[ig
].nmolReq
[1] = get_eint(&inp
, buf
, -1, wi
);
2352 printStringNoNewline(
2354 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2355 printStringNoNewline(
2357 "at maximum distance (= bulk concentration) to the split group layers. However,");
2358 printStringNoNewline(&inp
,
2359 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2360 "layer from the middle at 0.0");
2361 printStringNoNewline(&inp
,
2362 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2363 ir
->swap
->bulkOffset
[0] = get_ereal(&inp
, "bulk-offsetA", 0.0, wi
);
2364 ir
->swap
->bulkOffset
[1] = get_ereal(&inp
, "bulk-offsetB", 0.0, wi
);
2365 if (!(ir
->swap
->bulkOffset
[0] > -1.0 && ir
->swap
->bulkOffset
[0] < 1.0)
2366 || !(ir
->swap
->bulkOffset
[1] > -1.0 && ir
->swap
->bulkOffset
[1] < 1.0))
2368 warning_error(wi
, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2371 printStringNoNewline(
2372 &inp
, "Start to swap ions if threshold difference to requested count is reached");
2373 ir
->swap
->threshold
= get_ereal(&inp
, "threshold", 1.0, wi
);
2376 /* AdResS is no longer supported, but we need grompp to be able to
2377 refuse to process old .mdp files that used it. */
2378 ir
->bAdress
= (get_eeenum(&inp
, "adress", no_names
, wi
) != 0);
2380 /* User defined thingies */
2381 printStringNewline(&inp
, "User defined thingies");
2382 setStringEntry(&inp
, "user1-grps", inputrecStrings
->user1
, nullptr);
2383 setStringEntry(&inp
, "user2-grps", inputrecStrings
->user2
, nullptr);
2384 ir
->userint1
= get_eint(&inp
, "userint1", 0, wi
);
2385 ir
->userint2
= get_eint(&inp
, "userint2", 0, wi
);
2386 ir
->userint3
= get_eint(&inp
, "userint3", 0, wi
);
2387 ir
->userint4
= get_eint(&inp
, "userint4", 0, wi
);
2388 ir
->userreal1
= get_ereal(&inp
, "userreal1", 0, wi
);
2389 ir
->userreal2
= get_ereal(&inp
, "userreal2", 0, wi
);
2390 ir
->userreal3
= get_ereal(&inp
, "userreal3", 0, wi
);
2391 ir
->userreal4
= get_ereal(&inp
, "userreal4", 0, wi
);
2395 gmx::TextOutputFile
stream(mdparout
);
2396 write_inpfile(&stream
, mdparout
, &inp
, FALSE
, writeMdpHeader
, wi
);
2398 // Transform module data into a flat key-value tree for output.
2399 gmx::KeyValueTreeBuilder builder
;
2400 gmx::KeyValueTreeObjectBuilder builderObject
= builder
.rootObject();
2401 mdModules
->buildMdpOutput(&builderObject
);
2403 gmx::TextWriter
writer(&stream
);
2404 writeKeyValueTreeAsMdp(&writer
, builder
.build());
2409 /* Process options if necessary */
2410 for (m
= 0; m
< 2; m
++)
2412 for (i
= 0; i
< 2 * DIM
; i
++)
2421 if (sscanf(dumstr
[m
], "%lf", &(dumdub
[m
][XX
])) != 1)
2425 "Pressure coupling incorrect number of values (I need exactly 1)");
2427 dumdub
[m
][YY
] = dumdub
[m
][ZZ
] = dumdub
[m
][XX
];
2429 case epctSEMIISOTROPIC
:
2430 case epctSURFACETENSION
:
2431 if (sscanf(dumstr
[m
], "%lf%lf", &(dumdub
[m
][XX
]), &(dumdub
[m
][ZZ
])) != 2)
2435 "Pressure coupling incorrect number of values (I need exactly 2)");
2437 dumdub
[m
][YY
] = dumdub
[m
][XX
];
2439 case epctANISOTROPIC
:
2440 if (sscanf(dumstr
[m
], "%lf%lf%lf%lf%lf%lf", &(dumdub
[m
][XX
]), &(dumdub
[m
][YY
]),
2441 &(dumdub
[m
][ZZ
]), &(dumdub
[m
][3]), &(dumdub
[m
][4]), &(dumdub
[m
][5]))
2446 "Pressure coupling incorrect number of values (I need exactly 6)");
2450 gmx_fatal(FARGS
, "Pressure coupling type %s not implemented yet",
2451 epcoupltype_names
[ir
->epct
]);
2455 clear_mat(ir
->ref_p
);
2456 clear_mat(ir
->compress
);
2457 for (i
= 0; i
< DIM
; i
++)
2459 ir
->ref_p
[i
][i
] = dumdub
[1][i
];
2460 ir
->compress
[i
][i
] = dumdub
[0][i
];
2462 if (ir
->epct
== epctANISOTROPIC
)
2464 ir
->ref_p
[XX
][YY
] = dumdub
[1][3];
2465 ir
->ref_p
[XX
][ZZ
] = dumdub
[1][4];
2466 ir
->ref_p
[YY
][ZZ
] = dumdub
[1][5];
2467 if (ir
->ref_p
[XX
][YY
] != 0 && ir
->ref_p
[XX
][ZZ
] != 0 && ir
->ref_p
[YY
][ZZ
] != 0)
2470 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2471 "apply a threefold shear stress?\n");
2473 ir
->compress
[XX
][YY
] = dumdub
[0][3];
2474 ir
->compress
[XX
][ZZ
] = dumdub
[0][4];
2475 ir
->compress
[YY
][ZZ
] = dumdub
[0][5];
2476 for (i
= 0; i
< DIM
; i
++)
2478 for (m
= 0; m
< i
; m
++)
2480 ir
->ref_p
[i
][m
] = ir
->ref_p
[m
][i
];
2481 ir
->compress
[i
][m
] = ir
->compress
[m
][i
];
2486 if (ir
->comm_mode
== ecmNO
)
2491 opts
->couple_moltype
= nullptr;
2492 if (strlen(inputrecStrings
->couple_moltype
) > 0)
2494 if (ir
->efep
!= efepNO
)
2496 opts
->couple_moltype
= gmx_strdup(inputrecStrings
->couple_moltype
);
2497 if (opts
->couple_lam0
== opts
->couple_lam1
)
2499 warning(wi
, "The lambda=0 and lambda=1 states for coupling are identical");
2501 if (ir
->eI
== eiMD
&& (opts
->couple_lam0
== ecouplamNONE
|| opts
->couple_lam1
== ecouplamNONE
))
2505 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2512 "Free energy is turned off, so we will not decouple the molecule listed "
2516 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2517 if (ir
->efep
!= efepNO
)
2519 if (fep
->delta_lambda
!= 0)
2521 ir
->efep
= efepSLOWGROWTH
;
2525 if (fep
->edHdLPrintEnergy
== edHdLPrintEnergyYES
)
2527 fep
->edHdLPrintEnergy
= edHdLPrintEnergyTOTAL
;
2529 "Old option for dhdl-print-energy given: "
2530 "changing \"yes\" to \"total\"\n");
2533 if (ir
->bSimTemp
&& (fep
->edHdLPrintEnergy
== edHdLPrintEnergyNO
))
2535 /* always print out the energy to dhdl if we are doing
2536 expanded ensemble, since we need the total energy for
2537 analysis if the temperature is changing. In some
2538 conditions one may only want the potential energy, so
2539 we will allow that if the appropriate mdp setting has
2540 been enabled. Otherwise, total it is:
2542 fep
->edHdLPrintEnergy
= edHdLPrintEnergyTOTAL
;
2545 if ((ir
->efep
!= efepNO
) || ir
->bSimTemp
)
2547 ir
->bExpanded
= FALSE
;
2548 if ((ir
->efep
== efepEXPANDED
) || ir
->bSimTemp
)
2550 ir
->bExpanded
= TRUE
;
2552 do_fep_params(ir
, inputrecStrings
->fep_lambda
, inputrecStrings
->lambda_weights
, wi
);
2553 if (ir
->bSimTemp
) /* done after fep params */
2555 do_simtemp_params(ir
);
2558 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2559 * setup and not on the old way of specifying the free-energy setup,
2560 * we should check for using soft-core when not needed, since that
2561 * can complicate the sampling significantly.
2562 * Note that we only check for the automated coupling setup.
2563 * If the (advanced) user does FEP through manual topology changes,
2564 * this check will not be triggered.
2566 if (ir
->efep
!= efepNO
&& ir
->fepvals
->n_lambda
== 0 && ir
->fepvals
->sc_alpha
!= 0
2567 && (couple_lambda_has_vdw_on(opts
->couple_lam0
) && couple_lambda_has_vdw_on(opts
->couple_lam1
)))
2570 "You are using soft-core interactions while the Van der Waals interactions are "
2571 "not decoupled (note that the sc-coul option is only active when using lambda "
2572 "states). Although this will not lead to errors, you will need much more "
2573 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2578 ir
->fepvals
->n_lambda
= 0;
2581 /* WALL PARAMETERS */
2583 do_wall_params(ir
, inputrecStrings
->wall_atomtype
, inputrecStrings
->wall_density
, opts
, wi
);
2585 /* ORIENTATION RESTRAINT PARAMETERS */
2587 if (opts
->bOrire
&& gmx::splitString(inputrecStrings
->orirefitgrp
).size() != 1)
2589 warning_error(wi
, "ERROR: Need one orientation restraint fit group\n");
2592 /* DEFORMATION PARAMETERS */
2594 clear_mat(ir
->deform
);
2595 for (i
= 0; i
< 6; i
++)
2600 double gmx_unused canary
;
2601 int ndeform
= sscanf(inputrecStrings
->deform
, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub
[0][0]),
2602 &(dumdub
[0][1]), &(dumdub
[0][2]), &(dumdub
[0][3]), &(dumdub
[0][4]),
2603 &(dumdub
[0][5]), &canary
);
2605 if (strlen(inputrecStrings
->deform
) > 0 && ndeform
!= 6)
2609 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2610 inputrecStrings
->deform
)
2613 for (i
= 0; i
< 3; i
++)
2615 ir
->deform
[i
][i
] = dumdub
[0][i
];
2617 ir
->deform
[YY
][XX
] = dumdub
[0][3];
2618 ir
->deform
[ZZ
][XX
] = dumdub
[0][4];
2619 ir
->deform
[ZZ
][YY
] = dumdub
[0][5];
2620 if (ir
->epc
!= epcNO
)
2622 for (i
= 0; i
< 3; i
++)
2624 for (j
= 0; j
<= i
; j
++)
2626 if (ir
->deform
[i
][j
] != 0 && ir
->compress
[i
][j
] != 0)
2628 warning_error(wi
, "A box element has deform set and compressibility > 0");
2632 for (i
= 0; i
< 3; i
++)
2634 for (j
= 0; j
< i
; j
++)
2636 if (ir
->deform
[i
][j
] != 0)
2638 for (m
= j
; m
< DIM
; m
++)
2640 if (ir
->compress
[m
][j
] != 0)
2643 "An off-diagonal box element has deform set while "
2644 "compressibility > 0 for the same component of another box "
2645 "vector, this might lead to spurious periodicity effects.");
2646 warning(wi
, warn_buf
);
2654 /* Ion/water position swapping checks */
2655 if (ir
->eSwapCoords
!= eswapNO
)
2657 if (ir
->swap
->nstswap
< 1)
2659 warning_error(wi
, "swap_frequency must be 1 or larger when ion swapping is requested");
2661 if (ir
->swap
->nAverage
< 1)
2663 warning_error(wi
, "coupl_steps must be 1 or larger.\n");
2665 if (ir
->swap
->threshold
< 1.0)
2667 warning_error(wi
, "Ion count threshold must be at least 1.\n");
2671 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2674 std::vector
<std::string
> errorMessages
;
2675 ir
->mtsLevels
= gmx::setupMtsLevels(opts
->mtsOpts
, &errorMessages
);
2677 for (const auto& errorMessage
: errorMessages
)
2679 warning_error(wi
, errorMessage
.c_str());
2685 gmx::checkAwhParams(ir
->awhParams
, ir
, wi
);
2692 /* We would like gn to be const as well, but C doesn't allow this */
2693 /* TODO this is utility functionality (search for the index of a
2694 string in a collection), so should be refactored and located more
2696 int search_string(const char* s
, int ng
, char* gn
[])
2700 for (i
= 0; (i
< ng
); i
++)
2702 if (gmx_strcasecmp(s
, gn
[i
]) == 0)
2709 "Group %s referenced in the .mdp file was not found in the index file.\n"
2710 "Group names must match either [moleculetype] names or custom index group\n"
2711 "names, in which case you must supply an index file to the '-n' option\n"
2716 static void do_numbering(int natoms
,
2717 SimulationGroups
* groups
,
2718 gmx::ArrayRef
<std::string
> groupsFromMdpFile
,
2721 SimulationAtomGroupType gtype
,
2727 unsigned short* cbuf
;
2728 AtomGroupIndices
* grps
= &(groups
->groups
[gtype
]);
2729 int j
, gid
, aj
, ognr
, ntot
= 0;
2731 char warn_buf
[STRLEN
];
2733 title
= shortName(gtype
);
2736 /* Mark all id's as not set */
2737 for (int i
= 0; (i
< natoms
); i
++)
2742 for (int i
= 0; i
!= groupsFromMdpFile
.ssize(); ++i
)
2744 /* Lookup the group name in the block structure */
2745 gid
= search_string(groupsFromMdpFile
[i
].c_str(), block
->nr
, gnames
);
2746 if ((grptp
!= egrptpONE
) || (i
== 0))
2748 grps
->emplace_back(gid
);
2751 /* Now go over the atoms in the group */
2752 for (j
= block
->index
[gid
]; (j
< block
->index
[gid
+ 1]); j
++)
2757 /* Range checking */
2758 if ((aj
< 0) || (aj
>= natoms
))
2760 gmx_fatal(FARGS
, "Invalid atom number %d in indexfile", aj
+ 1);
2762 /* Lookup up the old group number */
2766 gmx_fatal(FARGS
, "Atom %d in multiple %s groups (%d and %d)", aj
+ 1, title
,
2771 /* Store the group number in buffer */
2772 if (grptp
== egrptpONE
)
2785 /* Now check whether we have done all atoms */
2788 if (grptp
== egrptpALL
)
2790 gmx_fatal(FARGS
, "%d atoms are not part of any of the %s groups", natoms
- ntot
, title
);
2792 else if (grptp
== egrptpPART
)
2794 sprintf(warn_buf
, "%d atoms are not part of any of the %s groups", natoms
- ntot
, title
);
2795 warning_note(wi
, warn_buf
);
2797 /* Assign all atoms currently unassigned to a rest group */
2798 for (j
= 0; (j
< natoms
); j
++)
2800 if (cbuf
[j
] == NOGID
)
2802 cbuf
[j
] = grps
->size();
2805 if (grptp
!= egrptpPART
)
2809 fprintf(stderr
, "Making dummy/rest group for %s containing %d elements\n", title
,
2812 /* Add group name "rest" */
2813 grps
->emplace_back(restnm
);
2815 /* Assign the rest name to all atoms not currently assigned to a group */
2816 for (j
= 0; (j
< natoms
); j
++)
2818 if (cbuf
[j
] == NOGID
)
2820 // group size was not updated before this here, so need to use -1.
2821 cbuf
[j
] = grps
->size() - 1;
2827 if (grps
->size() == 1 && (ntot
== 0 || ntot
== natoms
))
2829 /* All atoms are part of one (or no) group, no index required */
2830 groups
->groupNumbers
[gtype
].clear();
2834 for (int j
= 0; (j
< natoms
); j
++)
2836 groups
->groupNumbers
[gtype
].emplace_back(cbuf
[j
]);
2843 static void calc_nrdf(const gmx_mtop_t
* mtop
, t_inputrec
* ir
, char** gnames
)
2846 pull_params_t
* pull
;
2847 int natoms
, imin
, jmin
;
2848 int * nrdf2
, *na_vcm
, na_tot
;
2849 double * nrdf_tc
, *nrdf_vcm
, nrdf_uc
, *nrdf_vcm_sub
;
2854 * First calc 3xnr-atoms for each group
2855 * then subtract half a degree of freedom for each constraint
2857 * Only atoms and nuclei contribute to the degrees of freedom...
2862 const SimulationGroups
& groups
= mtop
->groups
;
2863 natoms
= mtop
->natoms
;
2865 /* Allocate one more for a possible rest group */
2866 /* We need to sum degrees of freedom into doubles,
2867 * since floats give too low nrdf's above 3 million atoms.
2869 snew(nrdf_tc
, groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
].size() + 1);
2870 snew(nrdf_vcm
, groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].size() + 1);
2871 snew(dof_vcm
, groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].size() + 1);
2872 snew(na_vcm
, groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].size() + 1);
2873 snew(nrdf_vcm_sub
, groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].size() + 1);
2875 for (gmx::index i
= 0; i
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
]); i
++)
2879 for (gmx::index i
= 0;
2880 i
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
]) + 1; i
++)
2883 clear_ivec(dof_vcm
[i
]);
2885 nrdf_vcm_sub
[i
] = 0;
2887 snew(nrdf2
, natoms
);
2888 for (const AtomProxy atomP
: AtomRange(*mtop
))
2890 const t_atom
& local
= atomP
.atom();
2891 int i
= atomP
.globalAtomNumber();
2893 if (local
.ptype
== eptAtom
|| local
.ptype
== eptNucleus
)
2895 int g
= getGroupType(groups
, SimulationAtomGroupType::Freeze
, i
);
2896 for (int d
= 0; d
< DIM
; d
++)
2898 if (opts
->nFreeze
[g
][d
] == 0)
2900 /* Add one DOF for particle i (counted as 2*1) */
2902 /* VCM group i has dim d as a DOF */
2903 dof_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, i
)][d
] =
2907 nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, i
)] +=
2909 nrdf_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, i
)] +=
2915 for (const gmx_molblock_t
& molb
: mtop
->molblock
)
2917 const gmx_moltype_t
& molt
= mtop
->moltype
[molb
.type
];
2918 const t_atom
* atom
= molt
.atoms
.atom
;
2919 for (int mol
= 0; mol
< molb
.nmol
; mol
++)
2921 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
2923 gmx::ArrayRef
<const int> ia
= molt
.ilist
[ftype
].iatoms
;
2924 for (int i
= 0; i
< molt
.ilist
[ftype
].size();)
2926 /* Subtract degrees of freedom for the constraints,
2927 * if the particles still have degrees of freedom left.
2928 * If one of the particles is a vsite or a shell, then all
2929 * constraint motion will go there, but since they do not
2930 * contribute to the constraints the degrees of freedom do not
2933 int ai
= as
+ ia
[i
+ 1];
2934 int aj
= as
+ ia
[i
+ 2];
2935 if (((atom
[ia
[i
+ 1]].ptype
== eptNucleus
) || (atom
[ia
[i
+ 1]].ptype
== eptAtom
))
2936 && ((atom
[ia
[i
+ 2]].ptype
== eptNucleus
) || (atom
[ia
[i
+ 2]].ptype
== eptAtom
)))
2954 imin
= std::min(imin
, nrdf2
[ai
]);
2955 jmin
= std::min(jmin
, nrdf2
[aj
]);
2958 nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
)] -=
2960 nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, aj
)] -=
2962 nrdf_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, ai
)] -=
2964 nrdf_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, aj
)] -=
2967 i
+= interaction_function
[ftype
].nratoms
+ 1;
2970 gmx::ArrayRef
<const int> ia
= molt
.ilist
[F_SETTLE
].iatoms
;
2971 for (int i
= 0; i
< molt
.ilist
[F_SETTLE
].size();)
2973 /* Subtract 1 dof from every atom in the SETTLE */
2974 for (int j
= 0; j
< 3; j
++)
2976 int ai
= as
+ ia
[i
+ 1 + j
];
2977 imin
= std::min(2, nrdf2
[ai
]);
2979 nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
)] -=
2981 nrdf_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, ai
)] -=
2986 as
+= molt
.atoms
.nr
;
2992 /* Correct nrdf for the COM constraints.
2993 * We correct using the TC and VCM group of the first atom
2994 * in the reference and pull group. If atoms in one pull group
2995 * belong to different TC or VCM groups it is anyhow difficult
2996 * to determine the optimal nrdf assignment.
2998 pull
= ir
->pull
.get();
3000 for (int i
= 0; i
< pull
->ncoord
; i
++)
3002 if (pull
->coord
[i
].eType
!= epullCONSTRAINT
)
3009 for (int j
= 0; j
< 2; j
++)
3011 const t_pull_group
* pgrp
;
3013 pgrp
= &pull
->group
[pull
->coord
[i
].group
[j
]];
3015 if (!pgrp
->ind
.empty())
3017 /* Subtract 1/2 dof from each group */
3018 int ai
= pgrp
->ind
[0];
3019 nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
)] -=
3021 nrdf_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, ai
)] -=
3023 if (nrdf_tc
[getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
)] < 0)
3026 "Center of mass pulling constraints caused the number of degrees "
3027 "of freedom for temperature coupling group %s to be negative",
3028 gnames
[groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
][getGroupType(
3029 groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
)]]);
3034 /* We need to subtract the whole DOF from group j=1 */
3041 if (ir
->nstcomm
!= 0)
3043 GMX_RELEASE_ASSERT(!groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].empty(),
3044 "Expect at least one group when removing COM motion");
3046 /* We remove COM motion up to dim ndof_com() */
3047 const int ndim_rm_vcm
= ndof_com(ir
);
3049 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3050 * the number of degrees of freedom in each vcm group when COM
3051 * translation is removed and 6 when rotation is removed as well.
3052 * Note that we do not and should not include the rest group here.
3054 for (gmx::index j
= 0;
3055 j
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
]); j
++)
3057 switch (ir
->comm_mode
)
3060 case ecmLINEAR_ACCELERATION_CORRECTION
:
3061 nrdf_vcm_sub
[j
] = 0;
3062 for (int d
= 0; d
< ndim_rm_vcm
; d
++)
3070 case ecmANGULAR
: nrdf_vcm_sub
[j
] = 6; break;
3071 default: gmx_incons("Checking comm_mode");
3075 for (gmx::index i
= 0;
3076 i
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
]); i
++)
3078 /* Count the number of atoms of TC group i for every VCM group */
3079 for (gmx::index j
= 0;
3080 j
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
]) + 1; j
++)
3085 for (int ai
= 0; ai
< natoms
; ai
++)
3087 if (getGroupType(groups
, SimulationAtomGroupType::TemperatureCoupling
, ai
) == i
)
3089 na_vcm
[getGroupType(groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, ai
)]++;
3093 /* Correct for VCM removal according to the fraction of each VCM
3094 * group present in this TC group.
3096 nrdf_uc
= nrdf_tc
[i
];
3098 for (gmx::index j
= 0;
3099 j
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
]) + 1; j
++)
3101 if (nrdf_vcm
[j
] > nrdf_vcm_sub
[j
])
3103 nrdf_tc
[i
] += nrdf_uc
* (static_cast<double>(na_vcm
[j
]) / static_cast<double>(na_tot
))
3104 * (nrdf_vcm
[j
] - nrdf_vcm_sub
[j
]) / nrdf_vcm
[j
];
3109 for (int i
= 0; (i
< gmx::ssize(groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
])); i
++)
3111 opts
->nrdf
[i
] = nrdf_tc
[i
];
3112 if (opts
->nrdf
[i
] < 0)
3116 fprintf(stderr
, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3117 gnames
[groups
.groups
[SimulationAtomGroupType::TemperatureCoupling
][i
]], opts
->nrdf
[i
]);
3125 sfree(nrdf_vcm_sub
);
3128 static bool do_egp_flag(t_inputrec
* ir
, SimulationGroups
* groups
, const char* option
, const char* val
, int flag
)
3130 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3131 * But since this is much larger than STRLEN, such a line can not be parsed.
3132 * The real maximum is the number of names that fit in a string: STRLEN/2.
3134 #define EGP_MAX (STRLEN / 2)
3138 auto names
= gmx::splitString(val
);
3139 if (names
.size() % 2 != 0)
3141 gmx_fatal(FARGS
, "The number of groups for %s is odd", option
);
3143 nr
= groups
->groups
[SimulationAtomGroupType::EnergyOutput
].size();
3145 for (size_t i
= 0; i
< names
.size() / 2; i
++)
3147 // TODO this needs to be replaced by a solution using std::find_if
3151 names
[2 * i
].c_str(),
3152 *(groups
->groupNames
[groups
->groups
[SimulationAtomGroupType::EnergyOutput
][j
]])))
3158 gmx_fatal(FARGS
, "%s in %s is not an energy group\n", names
[2 * i
].c_str(), option
);
3163 names
[2 * i
+ 1].c_str(),
3164 *(groups
->groupNames
[groups
->groups
[SimulationAtomGroupType::EnergyOutput
][k
]])))
3170 gmx_fatal(FARGS
, "%s in %s is not an energy group\n", names
[2 * i
+ 1].c_str(), option
);
3172 if ((j
< nr
) && (k
< nr
))
3174 ir
->opts
.egp_flags
[nr
* j
+ k
] |= flag
;
3175 ir
->opts
.egp_flags
[nr
* k
+ j
] |= flag
;
3184 static void make_swap_groups(t_swapcoords
* swap
, t_blocka
* grps
, char** gnames
)
3186 int ig
= -1, i
= 0, gind
;
3190 /* Just a quick check here, more thorough checks are in mdrun */
3191 if (strcmp(swap
->grp
[eGrpSplit0
].molname
, swap
->grp
[eGrpSplit1
].molname
) == 0)
3193 gmx_fatal(FARGS
, "The split groups can not both be '%s'.", swap
->grp
[eGrpSplit0
].molname
);
3196 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3197 for (ig
= 0; ig
< swap
->ngrp
; ig
++)
3199 swapg
= &swap
->grp
[ig
];
3200 gind
= search_string(swap
->grp
[ig
].molname
, grps
->nr
, gnames
);
3201 swapg
->nat
= grps
->index
[gind
+ 1] - grps
->index
[gind
];
3205 fprintf(stderr
, "%s group '%s' contains %d atoms.\n",
3206 ig
< 3 ? eSwapFixedGrp_names
[ig
] : "Swap", swap
->grp
[ig
].molname
, swapg
->nat
);
3207 snew(swapg
->ind
, swapg
->nat
);
3208 for (i
= 0; i
< swapg
->nat
; i
++)
3210 swapg
->ind
[i
] = grps
->a
[grps
->index
[gind
] + i
];
3215 gmx_fatal(FARGS
, "Swap group %s does not contain any atoms.", swap
->grp
[ig
].molname
);
3221 static void make_IMD_group(t_IMD
* IMDgroup
, char* IMDgname
, t_blocka
* grps
, char** gnames
)
3226 ig
= search_string(IMDgname
, grps
->nr
, gnames
);
3227 IMDgroup
->nat
= grps
->index
[ig
+ 1] - grps
->index
[ig
];
3229 if (IMDgroup
->nat
> 0)
3232 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3234 IMDgname
, IMDgroup
->nat
);
3235 snew(IMDgroup
->ind
, IMDgroup
->nat
);
3236 for (i
= 0; i
< IMDgroup
->nat
; i
++)
3238 IMDgroup
->ind
[i
] = grps
->a
[grps
->index
[ig
] + i
];
3243 /* Checks whether atoms are both part of a COM removal group and frozen.
3244 * If a fully frozen atom is part of a COM removal group, it is removed
3245 * from the COM removal group. A note is issued if such atoms are present.
3246 * A warning is issued for atom with one or two dimensions frozen that
3247 * are part of a COM removal group (mdrun would need to compute COM mass
3248 * per dimension to handle this correctly).
3249 * Also issues a warning when non-frozen atoms are not part of a COM
3250 * removal group while COM removal is active.
3252 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups
* groups
,
3254 const t_grpopts
& opts
,
3257 const int vcmRestGroup
=
3258 std::max(int(groups
->groups
[SimulationAtomGroupType::MassCenterVelocityRemoval
].size()), 1);
3260 int numFullyFrozenVcmAtoms
= 0;
3261 int numPartiallyFrozenVcmAtoms
= 0;
3262 int numNonVcmAtoms
= 0;
3263 for (int a
= 0; a
< numAtoms
; a
++)
3265 const int freezeGroup
= getGroupType(*groups
, SimulationAtomGroupType::Freeze
, a
);
3266 int numFrozenDims
= 0;
3267 for (int d
= 0; d
< DIM
; d
++)
3269 numFrozenDims
+= opts
.nFreeze
[freezeGroup
][d
];
3272 const int vcmGroup
= getGroupType(*groups
, SimulationAtomGroupType::MassCenterVelocityRemoval
, a
);
3273 if (vcmGroup
< vcmRestGroup
)
3275 if (numFrozenDims
== DIM
)
3277 /* Do not remove COM motion for this fully frozen atom */
3278 if (groups
->groupNumbers
[SimulationAtomGroupType::MassCenterVelocityRemoval
].empty())
3280 groups
->groupNumbers
[SimulationAtomGroupType::MassCenterVelocityRemoval
].resize(
3283 groups
->groupNumbers
[SimulationAtomGroupType::MassCenterVelocityRemoval
][a
] = vcmRestGroup
;
3284 numFullyFrozenVcmAtoms
++;
3286 else if (numFrozenDims
> 0)
3288 numPartiallyFrozenVcmAtoms
++;
3291 else if (numFrozenDims
< DIM
)
3297 if (numFullyFrozenVcmAtoms
> 0)
3299 std::string warningText
= gmx::formatString(
3300 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3301 "removing these atoms from the COMM removal group(s)",
3302 numFullyFrozenVcmAtoms
);
3303 warning_note(wi
, warningText
.c_str());
3305 if (numPartiallyFrozenVcmAtoms
> 0 && numPartiallyFrozenVcmAtoms
< numAtoms
)
3307 std::string warningText
= gmx::formatString(
3308 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3309 "removal group(s), due to limitations in the code these still contribute to the "
3310 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3312 numPartiallyFrozenVcmAtoms
, DIM
);
3313 warning(wi
, warningText
.c_str());
3315 if (numNonVcmAtoms
> 0)
3317 std::string warningText
= gmx::formatString(
3318 "%d atoms are not part of any center of mass motion removal group.\n"
3319 "This may lead to artifacts.\n"
3320 "In most cases one should use one group for the whole system.",
3322 warning(wi
, warningText
.c_str());
3326 void do_index(const char* mdparin
,
3330 const gmx::MdModulesNotifier
& notifier
,
3334 t_blocka
* defaultIndexGroups
;
3342 int i
, j
, k
, restnm
;
3343 bool bExcl
, bTable
, bAnneal
;
3344 char warn_buf
[STRLEN
];
3348 fprintf(stderr
, "processing index file...\n");
3352 snew(defaultIndexGroups
, 1);
3353 snew(defaultIndexGroups
->index
, 1);
3355 atoms_all
= gmx_mtop_global_atoms(mtop
);
3356 analyse(&atoms_all
, defaultIndexGroups
, &gnames
, FALSE
, TRUE
);
3357 done_atom(&atoms_all
);
3361 defaultIndexGroups
= init_index(ndx
, &gnames
);
3364 SimulationGroups
* groups
= &mtop
->groups
;
3365 natoms
= mtop
->natoms
;
3366 symtab
= &mtop
->symtab
;
3368 for (int i
= 0; (i
< defaultIndexGroups
->nr
); i
++)
3370 groups
->groupNames
.emplace_back(put_symtab(symtab
, gnames
[i
]));
3372 groups
->groupNames
.emplace_back(put_symtab(symtab
, "rest"));
3373 restnm
= groups
->groupNames
.size() - 1;
3374 GMX_RELEASE_ASSERT(restnm
== defaultIndexGroups
->nr
, "Size of allocations must match");
3375 srenew(gnames
, defaultIndexGroups
->nr
+ 1);
3376 gnames
[restnm
] = *(groups
->groupNames
.back());
3378 set_warning_line(wi
, mdparin
, -1);
3380 auto temperatureCouplingTauValues
= gmx::splitString(inputrecStrings
->tau_t
);
3381 auto temperatureCouplingReferenceValues
= gmx::splitString(inputrecStrings
->ref_t
);
3382 auto temperatureCouplingGroupNames
= gmx::splitString(inputrecStrings
->tcgrps
);
3383 if (temperatureCouplingTauValues
.size() != temperatureCouplingGroupNames
.size()
3384 || temperatureCouplingReferenceValues
.size() != temperatureCouplingGroupNames
.size())
3387 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3389 temperatureCouplingGroupNames
.size(), temperatureCouplingReferenceValues
.size(),
3390 temperatureCouplingTauValues
.size());
3393 const bool useReferenceTemperature
= integratorHasReferenceTemperature(ir
);
3394 do_numbering(natoms
, groups
, temperatureCouplingGroupNames
, defaultIndexGroups
, gnames
,
3395 SimulationAtomGroupType::TemperatureCoupling
, restnm
,
3396 useReferenceTemperature
? egrptpALL
: egrptpALL_GENREST
, bVerbose
, wi
);
3397 nr
= groups
->groups
[SimulationAtomGroupType::TemperatureCoupling
].size();
3399 snew(ir
->opts
.nrdf
, nr
);
3400 snew(ir
->opts
.tau_t
, nr
);
3401 snew(ir
->opts
.ref_t
, nr
);
3402 if (ir
->eI
== eiBD
&& ir
->bd_fric
== 0)
3404 fprintf(stderr
, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3407 if (useReferenceTemperature
)
3409 if (size_t(nr
) != temperatureCouplingReferenceValues
.size())
3411 gmx_fatal(FARGS
, "Not enough ref-t and tau-t values!");
3415 convertReals(wi
, temperatureCouplingTauValues
, "tau-t", ir
->opts
.tau_t
);
3416 for (i
= 0; (i
< nr
); i
++)
3418 if ((ir
->eI
== eiBD
) && ir
->opts
.tau_t
[i
] <= 0)
3420 sprintf(warn_buf
, "With integrator %s tau-t should be larger than 0", ei_names
[ir
->eI
]);
3421 warning_error(wi
, warn_buf
);
3424 if (ir
->etc
!= etcVRESCALE
&& ir
->opts
.tau_t
[i
] == 0)
3428 "tau-t = -1 is the value to signal that a group should not have "
3429 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3432 if (ir
->opts
.tau_t
[i
] >= 0)
3434 tau_min
= std::min(tau_min
, ir
->opts
.tau_t
[i
]);
3437 if (ir
->etc
!= etcNO
&& ir
->nsttcouple
== -1)
3439 ir
->nsttcouple
= ir_optimal_nsttcouple(ir
);
3444 if ((ir
->etc
== etcNOSEHOOVER
) && (ir
->epc
== epcBERENDSEN
))
3447 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3448 "md-vv; use either vrescale temperature with berendsen pressure or "
3449 "Nose-Hoover temperature with MTTK pressure");
3451 if (ir
->epc
== epcMTTK
)
3453 if (ir
->etc
!= etcNOSEHOOVER
)
3456 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3461 if (ir
->nstpcouple
!= ir
->nsttcouple
)
3463 int mincouple
= std::min(ir
->nstpcouple
, ir
->nsttcouple
);
3464 ir
->nstpcouple
= ir
->nsttcouple
= mincouple
;
3466 "for current Trotter decomposition methods with vv, nsttcouple and "
3467 "nstpcouple must be equal. Both have been reset to "
3468 "min(nsttcouple,nstpcouple) = %d",
3470 warning_note(wi
, warn_buf
);
3475 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3476 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3478 if (ETC_ANDERSEN(ir
->etc
))
3480 if (ir
->nsttcouple
!= 1)
3484 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3485 "need for larger nsttcouple > 1, since no global parameters are computed. "
3486 "nsttcouple has been reset to 1");
3487 warning_note(wi
, warn_buf
);
3490 nstcmin
= tcouple_min_integration_steps(ir
->etc
);
3493 if (tau_min
/ (ir
->delta_t
* ir
->nsttcouple
) < nstcmin
- 10 * GMX_REAL_EPS
)
3496 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3497 "least %d times larger than nsttcouple*dt (%g)",
3498 ETCOUPLTYPE(ir
->etc
), tau_min
, nstcmin
, ir
->nsttcouple
* ir
->delta_t
);
3499 warning(wi
, warn_buf
);
3502 convertReals(wi
, temperatureCouplingReferenceValues
, "ref-t", ir
->opts
.ref_t
);
3503 for (i
= 0; (i
< nr
); i
++)
3505 if (ir
->opts
.ref_t
[i
] < 0)
3507 gmx_fatal(FARGS
, "ref-t for group %d negative", i
);
3510 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3511 if we are in this conditional) if mc_temp is negative */
3512 if (ir
->expandedvals
->mc_temp
< 0)
3514 ir
->expandedvals
->mc_temp
= ir
->opts
.ref_t
[0]; /*for now, set to the first reft */
3518 /* Simulated annealing for each group. There are nr groups */
3519 auto simulatedAnnealingGroupNames
= gmx::splitString(inputrecStrings
->anneal
);
3520 if (simulatedAnnealingGroupNames
.size() == 1
3521 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames
[0], "N", 1))
3523 simulatedAnnealingGroupNames
.resize(0);
3525 if (!simulatedAnnealingGroupNames
.empty() && gmx::ssize(simulatedAnnealingGroupNames
) != nr
)
3527 gmx_fatal(FARGS
, "Wrong number of annealing values: %zu (for %d groups)\n",
3528 simulatedAnnealingGroupNames
.size(), nr
);
3532 snew(ir
->opts
.annealing
, nr
);
3533 snew(ir
->opts
.anneal_npoints
, nr
);
3534 snew(ir
->opts
.anneal_time
, nr
);
3535 snew(ir
->opts
.anneal_temp
, nr
);
3536 for (i
= 0; i
< nr
; i
++)
3538 ir
->opts
.annealing
[i
] = eannNO
;
3539 ir
->opts
.anneal_npoints
[i
] = 0;
3540 ir
->opts
.anneal_time
[i
] = nullptr;
3541 ir
->opts
.anneal_temp
[i
] = nullptr;
3543 if (!simulatedAnnealingGroupNames
.empty())
3546 for (i
= 0; i
< nr
; i
++)
3548 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames
[i
], "N", 1))
3550 ir
->opts
.annealing
[i
] = eannNO
;
3552 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames
[i
], "S", 1))
3554 ir
->opts
.annealing
[i
] = eannSINGLE
;
3557 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames
[i
], "P", 1))
3559 ir
->opts
.annealing
[i
] = eannPERIODIC
;
3565 /* Read the other fields too */
3566 auto simulatedAnnealingPoints
= gmx::splitString(inputrecStrings
->anneal_npoints
);
3567 if (simulatedAnnealingPoints
.size() != simulatedAnnealingGroupNames
.size())
3569 gmx_fatal(FARGS
, "Found %zu annealing-npoints values for %zu groups\n",
3570 simulatedAnnealingPoints
.size(), simulatedAnnealingGroupNames
.size());
3572 convertInts(wi
, simulatedAnnealingPoints
, "annealing points", ir
->opts
.anneal_npoints
);
3573 size_t numSimulatedAnnealingFields
= 0;
3574 for (i
= 0; i
< nr
; i
++)
3576 if (ir
->opts
.anneal_npoints
[i
] == 1)
3580 "Please specify at least a start and an end point for annealing\n");
3582 snew(ir
->opts
.anneal_time
[i
], ir
->opts
.anneal_npoints
[i
]);
3583 snew(ir
->opts
.anneal_temp
[i
], ir
->opts
.anneal_npoints
[i
]);
3584 numSimulatedAnnealingFields
+= ir
->opts
.anneal_npoints
[i
];
3587 auto simulatedAnnealingTimes
= gmx::splitString(inputrecStrings
->anneal_time
);
3589 if (simulatedAnnealingTimes
.size() != numSimulatedAnnealingFields
)
3591 gmx_fatal(FARGS
, "Found %zu annealing-time values, wanted %zu\n",
3592 simulatedAnnealingTimes
.size(), numSimulatedAnnealingFields
);
3594 auto simulatedAnnealingTemperatures
= gmx::splitString(inputrecStrings
->anneal_temp
);
3595 if (simulatedAnnealingTemperatures
.size() != numSimulatedAnnealingFields
)
3597 gmx_fatal(FARGS
, "Found %zu annealing-temp values, wanted %zu\n",
3598 simulatedAnnealingTemperatures
.size(), numSimulatedAnnealingFields
);
3601 std::vector
<real
> allSimulatedAnnealingTimes(numSimulatedAnnealingFields
);
3602 std::vector
<real
> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields
);
3603 convertReals(wi
, simulatedAnnealingTimes
, "anneal-time",
3604 allSimulatedAnnealingTimes
.data());
3605 convertReals(wi
, simulatedAnnealingTemperatures
, "anneal-temp",
3606 allSimulatedAnnealingTemperatures
.data());
3607 for (i
= 0, k
= 0; i
< nr
; i
++)
3609 for (j
= 0; j
< ir
->opts
.anneal_npoints
[i
]; j
++)
3611 ir
->opts
.anneal_time
[i
][j
] = allSimulatedAnnealingTimes
[k
];
3612 ir
->opts
.anneal_temp
[i
][j
] = allSimulatedAnnealingTemperatures
[k
];
3615 if (ir
->opts
.anneal_time
[i
][0] > (ir
->init_t
+ GMX_REAL_EPS
))
3617 gmx_fatal(FARGS
, "First time point for annealing > init_t.\n");
3623 if (ir
->opts
.anneal_time
[i
][j
] < ir
->opts
.anneal_time
[i
][j
- 1])
3626 "Annealing timepoints out of order: t=%f comes after "
3628 ir
->opts
.anneal_time
[i
][j
], ir
->opts
.anneal_time
[i
][j
- 1]);
3631 if (ir
->opts
.anneal_temp
[i
][j
] < 0)
3633 gmx_fatal(FARGS
, "Found negative temperature in annealing: %f\n",
3634 ir
->opts
.anneal_temp
[i
][j
]);
3639 /* Print out some summary information, to make sure we got it right */
3640 for (i
= 0; i
< nr
; i
++)
3642 if (ir
->opts
.annealing
[i
] != eannNO
)
3644 j
= groups
->groups
[SimulationAtomGroupType::TemperatureCoupling
][i
];
3645 fprintf(stderr
, "Simulated annealing for group %s: %s, %d timepoints\n",
3646 *(groups
->groupNames
[j
]), eann_names
[ir
->opts
.annealing
[i
]],
3647 ir
->opts
.anneal_npoints
[i
]);
3648 fprintf(stderr
, "Time (ps) Temperature (K)\n");
3649 /* All terms except the last one */
3650 for (j
= 0; j
< (ir
->opts
.anneal_npoints
[i
] - 1); j
++)
3652 fprintf(stderr
, "%9.1f %5.1f\n", ir
->opts
.anneal_time
[i
][j
],
3653 ir
->opts
.anneal_temp
[i
][j
]);
3656 /* Finally the last one */
3657 j
= ir
->opts
.anneal_npoints
[i
] - 1;
3658 if (ir
->opts
.annealing
[i
] == eannSINGLE
)
3660 fprintf(stderr
, "%9.1f- %5.1f\n", ir
->opts
.anneal_time
[i
][j
],
3661 ir
->opts
.anneal_temp
[i
][j
]);
3665 fprintf(stderr
, "%9.1f %5.1f\n", ir
->opts
.anneal_time
[i
][j
],
3666 ir
->opts
.anneal_temp
[i
][j
]);
3667 if (std::fabs(ir
->opts
.anneal_temp
[i
][j
] - ir
->opts
.anneal_temp
[i
][0]) > GMX_REAL_EPS
)
3670 "There is a temperature jump when your annealing "
3682 process_pull_groups(ir
->pull
->group
, inputrecStrings
->pullGroupNames
, defaultIndexGroups
, gnames
);
3684 checkPullCoords(ir
->pull
->group
, ir
->pull
->coord
);
3689 make_rotation_groups(ir
->rot
, inputrecStrings
->rotateGroupNames
, defaultIndexGroups
, gnames
);
3692 if (ir
->eSwapCoords
!= eswapNO
)
3694 make_swap_groups(ir
->swap
, defaultIndexGroups
, gnames
);
3697 /* Make indices for IMD session */
3700 make_IMD_group(ir
->imd
, inputrecStrings
->imd_grp
, defaultIndexGroups
, gnames
);
3703 gmx::IndexGroupsAndNames
defaultIndexGroupsAndNames(
3704 *defaultIndexGroups
, gmx::arrayRefFromArray(gnames
, defaultIndexGroups
->nr
));
3705 notifier
.preProcessingNotifications_
.notify(defaultIndexGroupsAndNames
);
3707 auto accelerations
= gmx::splitString(inputrecStrings
->acc
);
3708 auto accelerationGroupNames
= gmx::splitString(inputrecStrings
->accgrps
);
3709 if (accelerationGroupNames
.size() * DIM
!= accelerations
.size())
3711 gmx_fatal(FARGS
, "Invalid Acceleration input: %zu groups and %zu acc. values",
3712 accelerationGroupNames
.size(), accelerations
.size());
3714 do_numbering(natoms
, groups
, accelerationGroupNames
, defaultIndexGroups
, gnames
,
3715 SimulationAtomGroupType::Acceleration
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3716 nr
= groups
->groups
[SimulationAtomGroupType::Acceleration
].size();
3717 snew(ir
->opts
.acc
, nr
);
3718 ir
->opts
.ngacc
= nr
;
3720 convertRvecs(wi
, accelerations
, "anneal-time", ir
->opts
.acc
);
3722 auto freezeDims
= gmx::splitString(inputrecStrings
->frdim
);
3723 auto freezeGroupNames
= gmx::splitString(inputrecStrings
->freeze
);
3724 if (freezeDims
.size() != DIM
* freezeGroupNames
.size())
3726 gmx_fatal(FARGS
, "Invalid Freezing input: %zu groups and %zu freeze values",
3727 freezeGroupNames
.size(), freezeDims
.size());
3729 do_numbering(natoms
, groups
, freezeGroupNames
, defaultIndexGroups
, gnames
,
3730 SimulationAtomGroupType::Freeze
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3731 nr
= groups
->groups
[SimulationAtomGroupType::Freeze
].size();
3732 ir
->opts
.ngfrz
= nr
;
3733 snew(ir
->opts
.nFreeze
, nr
);
3734 for (i
= k
= 0; (size_t(i
) < freezeGroupNames
.size()); i
++)
3736 for (j
= 0; (j
< DIM
); j
++, k
++)
3738 ir
->opts
.nFreeze
[i
][j
] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims
[k
], "Y", 1));
3739 if (!ir
->opts
.nFreeze
[i
][j
])
3741 if (!gmx::equalCaseInsensitive(freezeDims
[k
], "N", 1))
3744 "Please use Y(ES) or N(O) for freezedim only "
3746 freezeDims
[k
].c_str());
3747 warning(wi
, warn_buf
);
3752 for (; (i
< nr
); i
++)
3754 for (j
= 0; (j
< DIM
); j
++)
3756 ir
->opts
.nFreeze
[i
][j
] = 0;
3760 auto energyGroupNames
= gmx::splitString(inputrecStrings
->energy
);
3761 do_numbering(natoms
, groups
, energyGroupNames
, defaultIndexGroups
, gnames
,
3762 SimulationAtomGroupType::EnergyOutput
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3763 add_wall_energrps(groups
, ir
->nwall
, symtab
);
3764 ir
->opts
.ngener
= groups
->groups
[SimulationAtomGroupType::EnergyOutput
].size();
3765 auto vcmGroupNames
= gmx::splitString(inputrecStrings
->vcm
);
3766 do_numbering(natoms
, groups
, vcmGroupNames
, defaultIndexGroups
, gnames
,
3767 SimulationAtomGroupType::MassCenterVelocityRemoval
, restnm
,
3768 vcmGroupNames
.empty() ? egrptpALL_GENREST
: egrptpPART
, bVerbose
, wi
);
3770 if (ir
->comm_mode
!= ecmNO
)
3772 checkAndUpdateVcmFreezeGroupConsistency(groups
, natoms
, ir
->opts
, wi
);
3775 /* Now we have filled the freeze struct, so we can calculate NRDF */
3776 calc_nrdf(mtop
, ir
, gnames
);
3778 auto user1GroupNames
= gmx::splitString(inputrecStrings
->user1
);
3779 do_numbering(natoms
, groups
, user1GroupNames
, defaultIndexGroups
, gnames
,
3780 SimulationAtomGroupType::User1
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3781 auto user2GroupNames
= gmx::splitString(inputrecStrings
->user2
);
3782 do_numbering(natoms
, groups
, user2GroupNames
, defaultIndexGroups
, gnames
,
3783 SimulationAtomGroupType::User2
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3784 auto compressedXGroupNames
= gmx::splitString(inputrecStrings
->x_compressed_groups
);
3785 do_numbering(natoms
, groups
, compressedXGroupNames
, defaultIndexGroups
, gnames
,
3786 SimulationAtomGroupType::CompressedPositionOutput
, restnm
, egrptpONE
, bVerbose
, wi
);
3787 auto orirefFitGroupNames
= gmx::splitString(inputrecStrings
->orirefitgrp
);
3788 do_numbering(natoms
, groups
, orirefFitGroupNames
, defaultIndexGroups
, gnames
,
3789 SimulationAtomGroupType::OrientationRestraintsFit
, restnm
, egrptpALL_GENREST
,
3792 /* MiMiC QMMM input processing */
3793 auto qmGroupNames
= gmx::splitString(inputrecStrings
->QMMM
);
3794 if (qmGroupNames
.size() > 1)
3796 gmx_fatal(FARGS
, "Currently, having more than one QM group in MiMiC is not supported");
3798 /* group rest, if any, is always MM! */
3799 do_numbering(natoms
, groups
, qmGroupNames
, defaultIndexGroups
, gnames
,
3800 SimulationAtomGroupType::QuantumMechanics
, restnm
, egrptpALL_GENREST
, bVerbose
, wi
);
3801 ir
->opts
.ngQM
= qmGroupNames
.size();
3803 /* end of MiMiC QMMM input */
3807 for (auto group
: gmx::keysOf(groups
->groups
))
3809 fprintf(stderr
, "%-16s has %zu element(s):", shortName(group
), groups
->groups
[group
].size());
3810 for (const auto& entry
: groups
->groups
[group
])
3812 fprintf(stderr
, " %s", *(groups
->groupNames
[entry
]));
3814 fprintf(stderr
, "\n");
3818 nr
= groups
->groups
[SimulationAtomGroupType::EnergyOutput
].size();
3819 snew(ir
->opts
.egp_flags
, nr
* nr
);
3821 bExcl
= do_egp_flag(ir
, groups
, "energygrp-excl", inputrecStrings
->egpexcl
, EGP_EXCL
);
3822 if (bExcl
&& ir
->cutoff_scheme
== ecutsVERLET
)
3824 warning_error(wi
, "Energy group exclusions are currently not supported");
3826 if (bExcl
&& EEL_FULL(ir
->coulombtype
))
3828 warning(wi
, "Can not exclude the lattice Coulomb energy between energy groups");
3831 bTable
= do_egp_flag(ir
, groups
, "energygrp-table", inputrecStrings
->egptable
, EGP_TABLE
);
3832 if (bTable
&& !(ir
->vdwtype
== evdwUSER
) && !(ir
->coulombtype
== eelUSER
)
3833 && !(ir
->coulombtype
== eelPMEUSER
) && !(ir
->coulombtype
== eelPMEUSERSWITCH
))
3836 "Can only have energy group pair tables in combination with user tables for VdW "
3840 /* final check before going out of scope if simulated tempering variables
3841 * need to be set to default values.
3843 if ((ir
->expandedvals
->nstexpanded
< 0) && ir
->bSimTemp
)
3845 ir
->expandedvals
->nstexpanded
= 2 * static_cast<int>(ir
->opts
.tau_t
[0] / ir
->delta_t
);
3846 warning(wi
, gmx::formatString(
3847 "the value for nstexpanded was not specified for "
3848 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3849 "by default, but it is recommended to set it to an explicit value!",
3850 ir
->expandedvals
->nstexpanded
));
3852 for (i
= 0; (i
< defaultIndexGroups
->nr
); i
++)
3857 done_blocka(defaultIndexGroups
);
3858 sfree(defaultIndexGroups
);
3862 static void check_disre(const gmx_mtop_t
* mtop
)
3864 if (gmx_mtop_ftype_count(mtop
, F_DISRES
) > 0)
3866 const gmx_ffparams_t
& ffparams
= mtop
->ffparams
;
3869 for (int i
= 0; i
< ffparams
.numTypes(); i
++)
3871 int ftype
= ffparams
.functype
[i
];
3872 if (ftype
== F_DISRES
)
3874 int label
= ffparams
.iparams
[i
].disres
.label
;
3875 if (label
== old_label
)
3877 fprintf(stderr
, "Distance restraint index %d occurs twice\n", label
);
3886 "Found %d double distance restraint indices,\n"
3887 "probably the parameters for multiple pairs in one restraint "
3888 "are not identical\n",
3894 static bool absolute_reference(const t_inputrec
* ir
, const gmx_mtop_t
* sys
, const bool posres_only
, ivec AbsRef
)
3897 gmx_mtop_ilistloop_t iloop
;
3899 const t_iparams
* pr
;
3906 for (d
= 0; d
< DIM
; d
++)
3908 AbsRef
[d
] = (d
< ndof_com(ir
) ? 0 : 1);
3910 /* Check for freeze groups */
3911 for (g
= 0; g
< ir
->opts
.ngfrz
; g
++)
3913 for (d
= 0; d
< DIM
; d
++)
3915 if (ir
->opts
.nFreeze
[g
][d
] != 0)
3923 /* Check for position restraints */
3924 iloop
= gmx_mtop_ilistloop_init(sys
);
3925 while (const InteractionLists
* ilist
= gmx_mtop_ilistloop_next(iloop
, &nmol
))
3927 if (nmol
> 0 && (AbsRef
[XX
] == 0 || AbsRef
[YY
] == 0 || AbsRef
[ZZ
] == 0))
3929 for (i
= 0; i
< (*ilist
)[F_POSRES
].size(); i
+= 2)
3931 pr
= &sys
->ffparams
.iparams
[(*ilist
)[F_POSRES
].iatoms
[i
]];
3932 for (d
= 0; d
< DIM
; d
++)
3934 if (pr
->posres
.fcA
[d
] != 0)
3940 for (i
= 0; i
< (*ilist
)[F_FBPOSRES
].size(); i
+= 2)
3942 /* Check for flat-bottom posres */
3943 pr
= &sys
->ffparams
.iparams
[(*ilist
)[F_FBPOSRES
].iatoms
[i
]];
3944 if (pr
->fbposres
.k
!= 0)
3946 switch (pr
->fbposres
.geom
)
3948 case efbposresSPHERE
: AbsRef
[XX
] = AbsRef
[YY
] = AbsRef
[ZZ
] = 1; break;
3949 case efbposresCYLINDERX
: AbsRef
[YY
] = AbsRef
[ZZ
] = 1; break;
3950 case efbposresCYLINDERY
: AbsRef
[XX
] = AbsRef
[ZZ
] = 1; break;
3951 case efbposresCYLINDER
:
3952 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3953 case efbposresCYLINDERZ
: AbsRef
[XX
] = AbsRef
[YY
] = 1; break;
3954 case efbposresX
: /* d=XX */
3955 case efbposresY
: /* d=YY */
3956 case efbposresZ
: /* d=ZZ */
3957 d
= pr
->fbposres
.geom
- efbposresX
;
3962 " Invalid geometry for flat-bottom position restraint.\n"
3963 "Expected nr between 1 and %d. Found %d\n",
3964 efbposresNR
- 1, pr
->fbposres
.geom
);
3971 return (AbsRef
[XX
] != 0 && AbsRef
[YY
] != 0 && AbsRef
[ZZ
] != 0);
3974 static void check_combination_rule_differences(const gmx_mtop_t
* mtop
,
3976 bool* bC6ParametersWorkWithGeometricRules
,
3977 bool* bC6ParametersWorkWithLBRules
,
3978 bool* bLBRulesPossible
)
3980 int ntypes
, tpi
, tpj
;
3983 double c6i
, c6j
, c12i
, c12j
;
3984 double c6
, c6_geometric
, c6_LB
;
3985 double sigmai
, sigmaj
, epsi
, epsj
;
3986 bool bCanDoLBRules
, bCanDoGeometricRules
;
3989 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3990 * force-field floating point parameters.
3993 ptr
= getenv("GMX_LJCOMB_TOL");
3997 double gmx_unused canary
;
3999 if (sscanf(ptr
, "%lf%lf", &dbl
, &canary
) != 1)
4002 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr
);
4007 *bC6ParametersWorkWithLBRules
= TRUE
;
4008 *bC6ParametersWorkWithGeometricRules
= TRUE
;
4009 bCanDoLBRules
= TRUE
;
4010 ntypes
= mtop
->ffparams
.atnr
;
4011 snew(typecount
, ntypes
);
4012 gmx_mtop_count_atomtypes(mtop
, state
, typecount
);
4013 *bLBRulesPossible
= TRUE
;
4014 for (tpi
= 0; tpi
< ntypes
; ++tpi
)
4016 c6i
= mtop
->ffparams
.iparams
[(ntypes
+ 1) * tpi
].lj
.c6
;
4017 c12i
= mtop
->ffparams
.iparams
[(ntypes
+ 1) * tpi
].lj
.c12
;
4018 for (tpj
= tpi
; tpj
< ntypes
; ++tpj
)
4020 c6j
= mtop
->ffparams
.iparams
[(ntypes
+ 1) * tpj
].lj
.c6
;
4021 c12j
= mtop
->ffparams
.iparams
[(ntypes
+ 1) * tpj
].lj
.c12
;
4022 c6
= mtop
->ffparams
.iparams
[ntypes
* tpi
+ tpj
].lj
.c6
;
4023 c6_geometric
= std::sqrt(c6i
* c6j
);
4024 if (!gmx_numzero(c6_geometric
))
4026 if (!gmx_numzero(c12i
) && !gmx_numzero(c12j
))
4028 sigmai
= gmx::sixthroot(c12i
/ c6i
);
4029 sigmaj
= gmx::sixthroot(c12j
/ c6j
);
4030 epsi
= c6i
* c6i
/ (4.0 * c12i
);
4031 epsj
= c6j
* c6j
/ (4.0 * c12j
);
4032 c6_LB
= 4.0 * std::sqrt(epsi
* epsj
) * gmx::power6(0.5 * (sigmai
+ sigmaj
));
4036 *bLBRulesPossible
= FALSE
;
4037 c6_LB
= c6_geometric
;
4039 bCanDoLBRules
= gmx_within_tol(c6_LB
, c6
, tol
);
4044 *bC6ParametersWorkWithLBRules
= FALSE
;
4047 bCanDoGeometricRules
= gmx_within_tol(c6_geometric
, c6
, tol
);
4049 if (!bCanDoGeometricRules
)
4051 *bC6ParametersWorkWithGeometricRules
= FALSE
;
4058 static void check_combination_rules(const t_inputrec
* ir
, const gmx_mtop_t
* mtop
, warninp_t wi
)
4060 bool bLBRulesPossible
, bC6ParametersWorkWithGeometricRules
, bC6ParametersWorkWithLBRules
;
4062 check_combination_rule_differences(mtop
, 0, &bC6ParametersWorkWithGeometricRules
,
4063 &bC6ParametersWorkWithLBRules
, &bLBRulesPossible
);
4064 if (ir
->ljpme_combination_rule
== eljpmeLB
)
4066 if (!bC6ParametersWorkWithLBRules
|| !bLBRulesPossible
)
4069 "You are using arithmetic-geometric combination rules "
4070 "in LJ-PME, but your non-bonded C6 parameters do not "
4071 "follow these rules.");
4076 if (!bC6ParametersWorkWithGeometricRules
)
4078 if (ir
->eDispCorr
!= edispcNO
)
4081 "You are using geometric combination rules in "
4082 "LJ-PME, but your non-bonded C6 parameters do "
4083 "not follow these rules. "
4084 "This will introduce very small errors in the forces and energies in "
4085 "your simulations. Dispersion correction will correct total energy "
4086 "and/or pressure for isotropic systems, but not forces or surface "
4092 "You are using geometric combination rules in "
4093 "LJ-PME, but your non-bonded C6 parameters do "
4094 "not follow these rules. "
4095 "This will introduce very small errors in the forces and energies in "
4096 "your simulations. If your system is homogeneous, consider using "
4097 "dispersion correction "
4098 "for the total energy and pressure.");
4104 void triple_check(const char* mdparin
, t_inputrec
* ir
, gmx_mtop_t
* sys
, warninp_t wi
)
4106 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4107 GMX_ASSERT(gmx::checkMtsRequirements(*ir
).empty(), "All MTS requirements should be met here");
4109 char err_buf
[STRLEN
];
4114 gmx_mtop_atomloop_block_t aloopb
;
4116 char warn_buf
[STRLEN
];
4118 set_warning_line(wi
, mdparin
, -1);
4120 if (absolute_reference(ir
, sys
, false, AbsRef
))
4123 "Removing center of mass motion in the presence of position restraints might "
4124 "cause artifacts. When you are using position restraints to equilibrate a "
4125 "macro-molecule, the artifacts are usually negligible.");
4128 if (ir
->cutoff_scheme
== ecutsVERLET
&& ir
->verletbuf_tol
> 0 && ir
->nstlist
> 1
4129 && ((EI_MD(ir
->eI
) || EI_SD(ir
->eI
)) && (ir
->etc
== etcVRESCALE
|| ir
->etc
== etcBERENDSEN
)))
4131 /* Check if a too small Verlet buffer might potentially
4132 * cause more drift than the thermostat can couple off.
4134 /* Temperature error fraction for warning and suggestion */
4135 const real T_error_warn
= 0.002;
4136 const real T_error_suggest
= 0.001;
4137 /* For safety: 2 DOF per atom (typical with constraints) */
4138 const real nrdf_at
= 2;
4139 real T
, tau
, max_T_error
;
4144 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
4146 T
= std::max(T
, ir
->opts
.ref_t
[i
]);
4147 tau
= std::max(tau
, ir
->opts
.tau_t
[i
]);
4151 /* This is a worst case estimate of the temperature error,
4152 * assuming perfect buffer estimation and no cancelation
4153 * of errors. The factor 0.5 is because energy distributes
4154 * equally over Ekin and Epot.
4156 max_T_error
= 0.5 * tau
* ir
->verletbuf_tol
/ (nrdf_at
* BOLTZ
* T
);
4157 if (max_T_error
> T_error_warn
)
4160 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4161 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4162 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4163 "%.0e or decrease tau_t.",
4164 ir
->verletbuf_tol
, T
, tau
, 100 * max_T_error
, 100 * T_error_suggest
,
4165 ir
->verletbuf_tol
* T_error_suggest
/ max_T_error
);
4166 warning(wi
, warn_buf
);
4171 if (ETC_ANDERSEN(ir
->etc
))
4175 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
4178 "all tau_t must currently be equal using Andersen temperature control, "
4179 "violated for group %d",
4181 CHECK(ir
->opts
.tau_t
[0] != ir
->opts
.tau_t
[i
]);
4183 "all tau_t must be positive using Andersen temperature control, "
4185 i
, ir
->opts
.tau_t
[i
]);
4186 CHECK(ir
->opts
.tau_t
[i
] < 0);
4189 if (ir
->etc
== etcANDERSENMASSIVE
&& ir
->comm_mode
!= ecmNO
)
4191 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
4193 int nsteps
= gmx::roundToInt(ir
->opts
.tau_t
[i
] / ir
->delta_t
);
4195 "tau_t/delta_t for group %d for temperature control method %s must be a "
4196 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4197 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4199 i
, etcoupl_names
[ir
->etc
], ir
->nstcomm
, ir
->opts
.tau_t
[i
], nsteps
);
4200 CHECK(nsteps
% ir
->nstcomm
!= 0);
4205 if (EI_DYNAMICS(ir
->eI
) && !EI_SD(ir
->eI
) && ir
->eI
!= eiBD
&& ir
->comm_mode
== ecmNO
4206 && !(absolute_reference(ir
, sys
, FALSE
, AbsRef
) || ir
->nsteps
<= 10) && !ETC_ANDERSEN(ir
->etc
))
4209 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4210 "rounding errors can lead to build up of kinetic energy of the center of mass");
4213 if (ir
->epc
== epcPARRINELLORAHMAN
&& ir
->etc
== etcNOSEHOOVER
)
4216 for (int g
= 0; g
< ir
->opts
.ngtc
; g
++)
4218 tau_t_max
= std::max(tau_t_max
, ir
->opts
.tau_t
[g
]);
4220 if (ir
->tau_p
< 1.9 * tau_t_max
)
4222 std::string message
= gmx::formatString(
4223 "With %s T-coupling and %s p-coupling, "
4224 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4225 etcoupl_names
[ir
->etc
], epcoupl_names
[ir
->epc
], "tau-p", ir
->tau_p
, "tau-t",
4227 warning(wi
, message
.c_str());
4231 /* Check for pressure coupling with absolute position restraints */
4232 if (ir
->epc
!= epcNO
&& ir
->refcoord_scaling
== erscNO
)
4234 absolute_reference(ir
, sys
, TRUE
, AbsRef
);
4236 for (m
= 0; m
< DIM
; m
++)
4238 if (AbsRef
[m
] && norm2(ir
->compress
[m
]) > 0)
4241 "You are using pressure coupling with absolute position restraints, "
4242 "this will give artifacts. Use the refcoord_scaling option.");
4250 aloopb
= gmx_mtop_atomloop_block_init(sys
);
4252 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
4254 if (atom
->q
!= 0 || atom
->qB
!= 0)
4262 if (EEL_FULL(ir
->coulombtype
))
4265 "You are using full electrostatics treatment %s for a system without charges.\n"
4266 "This costs a lot of performance for just processing zeros, consider using %s "
4268 EELTYPE(ir
->coulombtype
), EELTYPE(eelCUT
));
4269 warning(wi
, err_buf
);
4274 if (ir
->coulombtype
== eelCUT
&& ir
->rcoulomb
> 0)
4277 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4278 "You might want to consider using %s electrostatics.\n",
4280 warning_note(wi
, err_buf
);
4284 /* Check if combination rules used in LJ-PME are the same as in the force field */
4285 if (EVDW_PME(ir
->vdwtype
))
4287 check_combination_rules(ir
, sys
, wi
);
4290 /* Generalized reaction field */
4291 if (ir
->coulombtype
== eelGRF_NOTUSED
)
4294 "Generalized reaction-field electrostatics is no longer supported. "
4295 "You can use normal reaction-field instead and compute the reaction-field "
4296 "constant by hand.");
4300 for (int i
= 0; (i
< gmx::ssize(sys
->groups
.groups
[SimulationAtomGroupType::Acceleration
])); i
++)
4302 for (m
= 0; (m
< DIM
); m
++)
4304 if (fabs(ir
->opts
.acc
[i
][m
]) > 1e-6)
4313 snew(mgrp
, sys
->groups
.groups
[SimulationAtomGroupType::Acceleration
].size());
4314 for (const AtomProxy atomP
: AtomRange(*sys
))
4316 const t_atom
& local
= atomP
.atom();
4317 int i
= atomP
.globalAtomNumber();
4318 mgrp
[getGroupType(sys
->groups
, SimulationAtomGroupType::Acceleration
, i
)] += local
.m
;
4321 for (i
= 0; (i
< gmx::ssize(sys
->groups
.groups
[SimulationAtomGroupType::Acceleration
])); i
++)
4323 for (m
= 0; (m
< DIM
); m
++)
4325 acc
[m
] += ir
->opts
.acc
[i
][m
] * mgrp
[i
];
4329 for (m
= 0; (m
< DIM
); m
++)
4331 if (fabs(acc
[m
]) > 1e-6)
4333 const char* dim
[DIM
] = { "X", "Y", "Z" };
4334 fprintf(stderr
, "Net Acceleration in %s direction, will %s be corrected\n", dim
[m
],
4335 ir
->nstcomm
!= 0 ? "" : "not");
4336 if (ir
->nstcomm
!= 0 && m
< ndof_com(ir
))
4340 (i
< gmx::ssize(sys
->groups
.groups
[SimulationAtomGroupType::Acceleration
])); i
++)
4342 ir
->opts
.acc
[i
][m
] -= acc
[m
];
4350 if (ir
->efep
!= efepNO
&& ir
->fepvals
->sc_alpha
!= 0
4351 && !gmx_within_tol(sys
->ffparams
.reppow
, 12.0, 10 * GMX_DOUBLE_EPS
))
4353 gmx_fatal(FARGS
, "Soft-core interactions are only supported with VdW repulsion power 12");
4361 for (i
= 0; i
< ir
->pull
->ncoord
&& !bWarned
; i
++)
4363 if (ir
->pull
->coord
[i
].group
[0] == 0 || ir
->pull
->coord
[i
].group
[1] == 0)
4365 absolute_reference(ir
, sys
, FALSE
, AbsRef
);
4366 for (m
= 0; m
< DIM
; m
++)
4368 if (ir
->pull
->coord
[i
].dim
[m
] && !AbsRef
[m
])
4371 "You are using an absolute reference for pulling, but the rest of "
4372 "the system does not have an absolute reference. This will lead to "
4381 for (i
= 0; i
< 3; i
++)
4383 for (m
= 0; m
<= i
; m
++)
4385 if ((ir
->epc
!= epcNO
&& ir
->compress
[i
][m
] != 0) || ir
->deform
[i
][m
] != 0)
4387 for (c
= 0; c
< ir
->pull
->ncoord
; c
++)
4389 if (ir
->pull
->coord
[c
].eGeom
== epullgDIRPBC
&& ir
->pull
->coord
[c
].vec
[m
] != 0)
4392 "Can not have dynamic box while using pull geometry '%s' "
4394 EPULLGEOM(ir
->pull
->coord
[c
].eGeom
), 'x' + m
);
4405 void double_check(t_inputrec
* ir
, matrix box
, bool bHasNormalConstraints
, bool bHasAnyConstraints
, warninp_t wi
)
4407 char warn_buf
[STRLEN
];
4410 ptr
= check_box(ir
->pbcType
, box
);
4413 warning_error(wi
, ptr
);
4416 if (bHasNormalConstraints
&& ir
->eConstrAlg
== econtSHAKE
)
4418 if (ir
->shake_tol
<= 0.0)
4420 sprintf(warn_buf
, "ERROR: shake-tol must be > 0 instead of %g\n", ir
->shake_tol
);
4421 warning_error(wi
, warn_buf
);
4425 if ((ir
->eConstrAlg
== econtLINCS
) && bHasNormalConstraints
)
4427 /* If we have Lincs constraints: */
4428 if (ir
->eI
== eiMD
&& ir
->etc
== etcNO
&& ir
->eConstrAlg
== econtLINCS
&& ir
->nLincsIter
== 1)
4431 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4432 warning_note(wi
, warn_buf
);
4435 if ((ir
->eI
== eiCG
|| ir
->eI
== eiLBFGS
) && (ir
->nProjOrder
< 8))
4438 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4440 warning_note(wi
, warn_buf
);
4442 if (ir
->epc
== epcMTTK
)
4444 warning_error(wi
, "MTTK not compatible with lincs -- use shake instead.");
4448 if (bHasAnyConstraints
&& ir
->epc
== epcMTTK
)
4450 warning_error(wi
, "Constraints are not implemented with MTTK pressure control.");
4453 if (ir
->LincsWarnAngle
> 90.0)
4455 sprintf(warn_buf
, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4456 warning(wi
, warn_buf
);
4457 ir
->LincsWarnAngle
= 90.0;
4460 if (ir
->pbcType
!= PbcType::No
)
4462 if (ir
->nstlist
== 0)
4465 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4466 "atoms might cause the simulation to crash.");
4468 if (gmx::square(ir
->rlist
) >= max_cutoff2(ir
->pbcType
, box
))
4471 "ERROR: The cut-off length is longer than half the shortest box vector or "
4472 "longer than the smallest box diagonal element. Increase the box size or "
4473 "decrease rlist.\n");
4474 warning_error(wi
, warn_buf
);