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,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
38 * \brief Defines the high-level constraint code.
40 * \author Berk Hess <hess@kth.se>
41 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_mdlib
55 #include "gromacs/domdec/domdec.h"
56 #include "gromacs/domdec/domdec_struct.h"
57 #include "gromacs/essentialdynamics/edsam.h"
58 #include "gromacs/fileio/confio.h"
59 #include "gromacs/fileio/gmxfio.h"
60 #include "gromacs/fileio/pdbio.h"
61 #include "gromacs/gmxlib/nrnb.h"
62 #include "gromacs/math/utilities.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/mdlib/gmx_omp_nthreads.h"
65 #include "gromacs/mdlib/lincs.h"
66 #include "gromacs/mdlib/mdrun.h"
67 #include "gromacs/mdlib/settle.h"
68 #include "gromacs/mdlib/shake.h"
69 #include "gromacs/mdtypes/commrec.h"
70 #include "gromacs/mdtypes/inputrec.h"
71 #include "gromacs/mdtypes/md_enums.h"
72 #include "gromacs/mdtypes/mdatom.h"
73 #include "gromacs/mdtypes/state.h"
74 #include "gromacs/pbcutil/pbc.h"
75 #include "gromacs/pulling/pull.h"
76 #include "gromacs/topology/block.h"
77 #include "gromacs/topology/ifunc.h"
78 #include "gromacs/topology/mtop_lookup.h"
79 #include "gromacs/topology/mtop_util.h"
80 #include "gromacs/utility/arrayref.h"
81 #include "gromacs/utility/exceptions.h"
82 #include "gromacs/utility/fatalerror.h"
83 #include "gromacs/utility/gmxassert.h"
84 #include "gromacs/utility/pleasecite.h"
85 #include "gromacs/utility/smalloc.h"
86 #include "gromacs/utility/txtdump.h"
91 /* \brief Impl class for Constraints
93 * \todo Members like md, idef are valid only for the lifetime of a
94 * domain, which would be good to make clearer in the structure of the
95 * code. It should not be possible to call apply() if setConstraints()
96 * has not been called. For example, this could be achieved if
97 * setConstraints returned a valid object with such a method. That
98 * still requires that we manage the lifetime of that object
99 * correctly, however. */
100 class Constraints::Impl
103 Impl(const gmx_mtop_t &mtop_p,
104 const t_inputrec &ir_p,
106 const t_mdatoms &md_p,
107 const t_commrec *cr_p,
108 const gmx_multisim_t &ms,
110 gmx_wallcycle *wcycle_p,
111 bool pbcHandlingRequired,
115 void setConstraints(const gmx_localtop_t &top,
116 const t_mdatoms &md);
117 bool apply(bool bLog,
130 ConstraintVariable econq);
131 //! The total number of constraints.
133 //! The number of flexible constraints.
135 //! A list of atoms to constraints for each moleculetype.
136 std::vector<t_blocka> at2con_mt;
137 //! The size of at2settle = number of moltypes
138 int n_at2settle_mt = 0;
139 //! A list of atoms to settles.
140 int **at2settle_mt = nullptr;
141 //! Whether any SETTLES cross charge-group boundaries.
142 bool bInterCGsettles = false;
144 Lincs *lincsd = nullptr; // TODO this should become a unique_ptr
146 shakedata *shaked = nullptr;
148 settledata *settled = nullptr;
149 //! The maximum number of warnings.
151 //! The number of warnings for LINCS.
152 int warncount_lincs = 0;
153 //! The number of warnings for SETTLE.
154 int warncount_settle = 0;
155 //! The essential dynamics data.
156 gmx_edsam * ed = nullptr;
158 //! Thread-local virial contribution.
159 tensor *vir_r_m_dr_th = {nullptr};
160 //! Did a settle error occur?
161 bool *bSettleErrorHasOccurred = nullptr;
163 //! Pointer to the global topology - only used for printing warnings.
164 const gmx_mtop_t &mtop;
165 //! Parameters for the interactions in this domain.
166 const t_idef *idef = nullptr;
167 //! Data about atoms in this domain.
169 //! Whether we need to do pbc for handling bonds.
170 bool pbcHandlingRequired_ = false;
174 //! Communication support.
175 const t_commrec *cr = nullptr;
176 //! Multi-sim support.
177 const gmx_multisim_t &ms;
178 /*!\brief Input options.
180 * \todo Replace with IMdpOptions */
181 const t_inputrec &ir;
182 //! Flop counting support.
183 t_nrnb *nrnb = nullptr;
184 //! Tracks wallcycle usage.
185 gmx_wallcycle *wcycle;
188 Constraints::~Constraints() = default;
190 int Constraints::numFlexibleConstraints() const
192 return impl_->nflexcon;
195 //! Clears constraint quantities for atoms in nonlocal region.
196 static void clear_constraint_quantity_nonlocal(gmx_domdec_t *dd, rvec *q)
198 int nonlocal_at_start, nonlocal_at_end, at;
200 dd_get_constraint_range(dd, &nonlocal_at_start, &nonlocal_at_end);
202 for (at = nonlocal_at_start; at < nonlocal_at_end; at++)
208 void too_many_constraint_warnings(int eConstrAlg, int warncount)
211 "Too many %s warnings (%d)\n"
212 "If you know what you are doing you can %s"
213 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
214 "but normally it is better to fix the problem",
215 (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE", warncount,
216 (eConstrAlg == econtLINCS) ?
217 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
220 //! Writes out coordinates.
221 static void write_constr_pdb(const char *fn, const char *title,
222 const gmx_mtop_t &mtop,
223 int start, int homenr, const t_commrec *cr,
224 const rvec x[], matrix box)
228 int dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
230 const char *anm, *resnm;
233 if (DOMAINDECOMP(cr))
236 dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
243 sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
247 sprintf(fname, "%s.pdb", fn);
250 out = gmx_fio_fopen(fname, "w");
252 fprintf(out, "TITLE %s\n", title);
253 gmx_write_pdb_box(out, -1, box);
255 for (i = start; i < start+homenr; i++)
259 if (i >= dd_numHomeAtoms(*dd) && i < dd_ac0)
263 ii = dd->globalAtomIndices[i];
269 mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, nullptr);
270 gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, anm, ' ', resnm, ' ', resnr, ' ',
271 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 0.0, "");
273 fprintf(out, "TER\n");
278 //! Writes out domain contents to help diagnose crashes.
279 static void dump_confs(FILE *log, int64_t step, const gmx_mtop_t &mtop,
280 int start, int homenr, const t_commrec *cr,
281 const rvec x[], rvec xprime[], matrix box)
283 char buf[STRLEN], buf2[22];
285 char *env = getenv("GMX_SUPPRESS_DUMP");
291 sprintf(buf, "step%sb", gmx_step_str(step, buf2));
292 write_constr_pdb(buf, "initial coordinates",
293 mtop, start, homenr, cr, x, box);
294 sprintf(buf, "step%sc", gmx_step_str(step, buf2));
295 write_constr_pdb(buf, "coordinates after constraining",
296 mtop, start, homenr, cr, xprime, box);
299 fprintf(log, "Wrote pdb files with previous and current coordinates\n");
301 fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
305 Constraints::apply(bool bLog,
318 ConstraintVariable econq)
320 return impl_->apply(bLog,
337 Constraints::Impl::apply(bool bLog,
350 ConstraintVariable econq)
356 real invdt, vir_fac = 0, t;
358 t_pbc pbc, *pbc_null;
362 wallcycle_start(wcycle, ewcCONSTR);
364 if (econq == ConstraintVariable::ForceDispl && !EI_ENERGY_MINIMIZATION(ir.eI))
366 gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
375 scaled_delta_t = step_scaling * ir.delta_t;
377 /* Prepare time step for use in constraint implementations, and
378 avoid generating inf when ir.delta_t = 0. */
385 invdt = 1.0/scaled_delta_t;
388 if (ir.efep != efepNO && EI_DYNAMICS(ir.eI))
390 /* Set the constraint lengths for the step at which this configuration
391 * is meant to be. The invmasses should not be changed.
393 lambda += delta_step*ir.fepvals->delta_lambda;
398 clear_mat(vir_r_m_dr);
400 const t_ilist *settle = &idef->il[F_SETTLE];
401 nsettle = settle->nr/(1+NRAL(F_SETTLE));
405 nth = gmx_omp_nthreads_get(emntSETTLE);
412 /* We do not need full pbc when constraints do not cross charge groups,
413 * i.e. when dd->constraint_comm==NULL.
414 * Note that PBC for constraints is different from PBC for bondeds.
415 * For constraints there is both forward and backward communication.
417 if (ir.ePBC != epbcNONE &&
418 (cr->dd || pbcHandlingRequired_) && !(cr->dd && cr->dd->constraint_comm == nullptr))
420 /* With pbc=screw the screw has been changed to a shift
421 * by the constraint coordinate communication routine,
422 * so that here we can use normal pbc.
424 pbc_null = set_pbc_dd(&pbc, ir.ePBC,
425 DOMAINDECOMP(cr) ? cr->dd->nc : nullptr,
433 /* Communicate the coordinates required for the non-local constraints
434 * for LINCS and/or SETTLE.
438 dd_move_x_constraints(cr->dd, box, x, xprime, econq == ConstraintVariable::Positions);
442 /* We need to initialize the non-local components of v.
443 * We never actually use these values, but we do increment them,
444 * so we should avoid uninitialized variables and overflows.
446 clear_constraint_quantity_nonlocal(cr->dd, v);
450 if (lincsd != nullptr)
452 bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms,
454 box, pbc_null, lambda, dvdlambda,
455 invdt, v, vir != nullptr, vir_r_m_dr,
457 maxwarn, &warncount_lincs);
458 if (!bOK && maxwarn < INT_MAX)
462 fprintf(log, "Constraint error in algorithm %s at step %s\n",
463 econstr_names[econtLINCS], gmx_step_str(step, buf));
469 if (shaked != nullptr)
471 bOK = constrain_shake(log, shaked,
473 *idef, ir, x, xprime, min_proj, nrnb,
475 invdt, v, vir != nullptr, vir_r_m_dr,
476 maxwarn < INT_MAX, econq);
478 if (!bOK && maxwarn < INT_MAX)
482 fprintf(log, "Constraint error in algorithm %s at step %s\n",
483 econstr_names[econtSHAKE], gmx_step_str(step, buf));
491 bool bSettleErrorHasOccurred0 = false;
495 case ConstraintVariable::Positions:
496 #pragma omp parallel for num_threads(nth) schedule(static)
497 for (th = 0; th < nth; th++)
503 clear_mat(vir_r_m_dr_th[th]);
510 invdt, v ? v[0] : nullptr,
512 th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
513 th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
515 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
517 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
520 inc_nrnb(nrnb, eNR_CONSTR_V, nsettle*3);
524 inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle*3);
527 case ConstraintVariable::Velocities:
528 case ConstraintVariable::Derivative:
529 case ConstraintVariable::Force:
530 case ConstraintVariable::ForceDispl:
531 #pragma omp parallel for num_threads(nth) schedule(static)
532 for (th = 0; th < nth; th++)
536 int calcvir_atom_end;
540 calcvir_atom_end = 0;
544 calcvir_atom_end = md.homenr;
549 clear_mat(vir_r_m_dr_th[th]);
552 int start_th = (nsettle* th )/nth;
553 int end_th = (nsettle*(th+1))/nth;
555 if (start_th >= 0 && end_th - start_th > 0)
557 settle_proj(settled, econq,
559 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
562 xprime, min_proj, calcvir_atom_end,
563 th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
566 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
568 /* This is an overestimate */
569 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
571 case ConstraintVariable::Deriv_FlexCon:
572 /* Nothing to do, since the are no flexible constraints in settles */
575 gmx_incons("Unknown constraint quantity for settle");
580 /* Reduce the virial contributions over the threads */
581 for (int th = 1; th < nth; th++)
583 m_add(vir_r_m_dr, vir_r_m_dr_th[th], vir_r_m_dr);
587 if (econq == ConstraintVariable::Positions)
589 for (int th = 1; th < nth; th++)
591 bSettleErrorHasOccurred0 = bSettleErrorHasOccurred0 || bSettleErrorHasOccurred[th];
594 if (bSettleErrorHasOccurred0)
598 "\nstep " "%" PRId64 ": One or more water molecules can not be settled.\n"
599 "Check for bad contacts and/or reduce the timestep if appropriate.\n",
603 fprintf(log, "%s", buf);
605 fprintf(stderr, "%s", buf);
607 if (warncount_settle > maxwarn)
609 too_many_constraint_warnings(-1, warncount_settle);
620 /* The normal uses of constrain() pass step_scaling = 1.0.
621 * The call to constrain() for SD1 that passes step_scaling =
622 * 0.5 also passes vir = NULL, so cannot reach this
623 * assertion. This assertion should remain until someone knows
624 * that this path works for their intended purpose, and then
625 * they can use scaled_delta_t instead of ir.delta_t
627 assert(gmx_within_tol(step_scaling, 1.0, GMX_REAL_EPS));
630 case ConstraintVariable::Positions:
631 vir_fac = 0.5/(ir.delta_t*ir.delta_t);
633 case ConstraintVariable::Velocities:
634 vir_fac = 0.5/ir.delta_t;
636 case ConstraintVariable::Force:
637 case ConstraintVariable::ForceDispl:
641 gmx_incons("Unsupported constraint quantity for virial");
646 vir_fac *= 2; /* only constraining over half the distance here */
648 for (int i = 0; i < DIM; i++)
650 for (int j = 0; j < DIM; j++)
652 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
659 dump_confs(log, step, mtop, start, homenr, cr, x, xprime, box);
662 if (econq == ConstraintVariable::Positions)
664 if (ir.bPull && pull_have_constraint(ir.pull_work))
666 if (EI_DYNAMICS(ir.eI))
668 t = ir.init_t + (step + delta_step)*ir.delta_t;
674 set_pbc(&pbc, ir.ePBC, box);
675 pull_constraint(ir.pull_work, &md, &pbc, cr, ir.delta_t, t, x, xprime, v, *vir);
677 if (ed && delta_step > 0)
679 /* apply the essential dynamics constraints here */
680 do_edsam(&ir, step, cr, xprime, v, box, ed);
683 wallcycle_stop(wcycle, ewcCONSTR);
685 if (v != nullptr && md.cFREEZE)
687 /* Set the velocities of frozen dimensions to zero */
689 int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
691 #pragma omp parallel for num_threads(numThreads) schedule(static)
692 for (int i = 0; i < md.homenr; i++)
694 int freezeGroup = md.cFREEZE[i];
696 for (int d = 0; d < DIM; d++)
698 if (ir.opts.nFreeze[freezeGroup][d])
709 ArrayRef<real> Constraints::rmsdData() const
713 return lincs_rmsdData(impl_->lincsd);
717 return EmptyArrayRef();
721 real Constraints::rmsd() const
725 return lincs_rmsd(impl_->lincsd);
733 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator)
735 if (haveDynamicsIntegrator)
737 return FlexibleConstraintTreatment::Include;
741 return FlexibleConstraintTreatment::Exclude;
745 t_blocka make_at2con(int numAtoms,
746 const t_ilist *ilists,
747 const t_iparams *iparams,
748 FlexibleConstraintTreatment flexibleConstraintTreatment)
750 GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams");
752 std::vector<int> count(numAtoms);
754 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
756 const t_ilist &ilist = ilists[ftype];
757 const int stride = 1 + NRAL(ftype);
758 for (int i = 0; i < ilist.nr; i += stride)
760 if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
761 !isConstraintFlexible(iparams, ilist.iatoms[i]))
763 for (int j = 1; j < 3; j++)
765 int a = ilist.iatoms[i + j];
773 at2con.nr = numAtoms;
774 at2con.nalloc_index = at2con.nr + 1;
775 snew(at2con.index, at2con.nalloc_index);
777 for (int a = 0; a < numAtoms; a++)
779 at2con.index[a + 1] = at2con.index[a] + count[a];
782 at2con.nra = at2con.index[at2con.nr];
783 at2con.nalloc_a = at2con.nra;
784 snew(at2con.a, at2con.nalloc_a);
786 /* The F_CONSTRNC constraints have constraint numbers
787 * that continue after the last F_CONSTR constraint.
789 int numConstraints = 0;
790 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
792 const t_ilist &ilist = ilists[ftype];
793 const int stride = 1 + NRAL(ftype);
794 for (int i = 0; i < ilist.nr; i += stride)
796 if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
797 !isConstraintFlexible(iparams, ilist.iatoms[i]))
799 for (int j = 1; j < 3; j++)
801 int a = ilist.iatoms[i + j];
802 at2con.a[at2con.index[a] + count[a]++] = numConstraints;
812 int countFlexibleConstraints(const t_ilist *ilist,
813 const t_iparams *iparams)
816 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
818 const int numIatomsPerConstraint = 3;
819 int ncon = ilist[ftype].nr / numIatomsPerConstraint;
820 t_iatom *ia = ilist[ftype].iatoms;
821 for (int con = 0; con < ncon; con++)
823 if (iparams[ia[0]].constr.dA == 0 &&
824 iparams[ia[0]].constr.dB == 0)
828 ia += numIatomsPerConstraint;
835 //! Returns the index of the settle to which each atom belongs.
836 static int *make_at2settle(int natoms, const t_ilist *ilist)
842 /* Set all to no settle */
843 for (a = 0; a < natoms; a++)
848 stride = 1 + NRAL(F_SETTLE);
850 for (s = 0; s < ilist->nr; s += stride)
852 at2s[ilist->iatoms[s+1]] = s/stride;
853 at2s[ilist->iatoms[s+2]] = s/stride;
854 at2s[ilist->iatoms[s+3]] = s/stride;
861 Constraints::Impl::setConstraints(const gmx_localtop_t &top,
868 /* With DD we might also need to call LINCS on a domain no constraints for
869 * communicating coordinates to other nodes that do have constraints.
871 if (ir.eConstrAlg == econtLINCS)
873 set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
875 if (ir.eConstrAlg == econtSHAKE)
879 // We are using the local topology, so there are only
880 // F_CONSTR constraints.
881 make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], &top.cgs, cr->dd);
885 make_shake_sblock_serial(shaked, idef, md);
892 settle_set_constraints(settled,
893 &idef->il[F_SETTLE], md);
896 /* Make a selection of the local atoms for essential dynamics */
899 dd_make_local_ed_indices(cr->dd, ed);
904 Constraints::setConstraints(const gmx_localtop_t &top,
907 impl_->setConstraints(top, md);
910 /*! \brief Makes a per-moleculetype container of mappings from atom
911 * indices to constraint indices.
913 * Note that flexible constraints are only enabled with a dynamical integrator. */
914 static std::vector<t_blocka>
915 makeAtomToConstraintMappings(const gmx_mtop_t &mtop,
916 FlexibleConstraintTreatment flexibleConstraintTreatment)
918 std::vector<t_blocka> mapping;
919 mapping.reserve(mtop.moltype.size());
920 for (const gmx_moltype_t &moltype : mtop.moltype)
922 mapping.push_back(make_at2con(moltype.atoms.nr,
924 mtop.ffparams.iparams,
925 flexibleConstraintTreatment));
930 Constraints::Constraints(const gmx_mtop_t &mtop,
931 const t_inputrec &ir,
935 const gmx_multisim_t &ms,
937 gmx_wallcycle *wcycle,
938 bool pbcHandlingRequired,
941 : impl_(new Impl(mtop,
955 Constraints::Impl::Impl(const gmx_mtop_t &mtop_p,
956 const t_inputrec &ir_p,
958 const t_mdatoms &md_p,
959 const t_commrec *cr_p,
960 const gmx_multisim_t &ms_p,
962 gmx_wallcycle *wcycle_p,
963 bool pbcHandlingRequired,
966 : ncon_tot(numConstraints),
969 pbcHandlingRequired_(pbcHandlingRequired),
977 if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
979 gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
983 if (numConstraints > 0)
985 at2con_mt = makeAtomToConstraintMappings(mtop,
986 flexibleConstraintTreatment(EI_DYNAMICS(ir.eI)));
988 for (const gmx_molblock_t &molblock : mtop.molblock)
990 int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist,
991 mtop.ffparams.iparams);
992 nflexcon += molblock.nmol*count;
999 fprintf(log, "There are %d flexible constraints\n",
1001 if (ir.fc_stepsize == 0)
1004 "WARNING: step size for flexible constraining = 0\n"
1005 " All flexible constraints will be rigid.\n"
1006 " Will try to keep all flexible constraints at their original length,\n"
1007 " but the lengths may exhibit some drift.\n\n");
1013 please_cite(log, "Hess2002");
1017 if (ir.eConstrAlg == econtLINCS)
1019 lincsd = init_lincs(log, mtop,
1020 nflexcon, at2con_mt,
1021 DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1022 ir.nLincsIter, ir.nProjOrder);
1025 if (ir.eConstrAlg == econtSHAKE)
1027 if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1029 gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1033 gmx_fatal(FARGS, "For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1035 please_cite(log, "Ryckaert77a");
1038 please_cite(log, "Barth95a");
1041 shaked = shake_init();
1047 please_cite(log, "Miyamoto92a");
1049 bInterCGsettles = inter_charge_group_settles(mtop);
1051 settled = settle_init(mtop);
1053 /* Make an atom to settle index for use in domain decomposition */
1054 n_at2settle_mt = mtop.moltype.size();
1055 snew(at2settle_mt, n_at2settle_mt);
1056 for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
1059 make_at2settle(mtop.moltype[mt].atoms.nr,
1060 &mtop.moltype[mt].ilist[F_SETTLE]);
1063 /* Allocate thread-local work arrays */
1064 int nthreads = gmx_omp_nthreads_get(emntSETTLE);
1065 if (nthreads > 1 && vir_r_m_dr_th == nullptr)
1067 snew(vir_r_m_dr_th, nthreads);
1068 snew(bSettleErrorHasOccurred, nthreads);
1073 char *env = getenv("GMX_MAXCONSTRWARN");
1077 sscanf(env, "%8d", &maxwarn);
1085 "Setting the maximum number of constraint warnings to %d\n",
1091 "Setting the maximum number of constraint warnings to %d\n",
1095 warncount_lincs = 0;
1096 warncount_settle = 0;
1099 Constraints::Impl::~Impl()
1104 void Constraints::saveEdsamPointer(gmx_edsam * ed)
1109 const ArrayRef<const t_blocka>
1110 Constraints::atom2constraints_moltype() const
1112 return impl_->at2con_mt;
1115 int *const* Constraints::atom2settle_moltype() const
1117 return impl_->at2settle_mt;
1121 bool inter_charge_group_constraints(const gmx_mtop_t &mtop)
1123 const gmx_moltype_t *molt;
1126 int *at2cg, cg, a, ftype, i;
1130 for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1132 molt = &mtop.moltype[mtop.molblock[mb].type];
1134 if (molt->ilist[F_CONSTR].nr > 0 ||
1135 molt->ilist[F_CONSTRNC].nr > 0 ||
1136 molt->ilist[F_SETTLE].nr > 0)
1139 snew(at2cg, molt->atoms.nr);
1140 for (cg = 0; cg < cgs->nr; cg++)
1142 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1148 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
1150 il = &molt->ilist[ftype];
1151 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype))
1153 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1167 bool inter_charge_group_settles(const gmx_mtop_t &mtop)
1169 const gmx_moltype_t *molt;
1172 int *at2cg, cg, a, ftype, i;
1176 for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
1178 molt = &mtop.moltype[mtop.molblock[mb].type];
1180 if (molt->ilist[F_SETTLE].nr > 0)
1183 snew(at2cg, molt->atoms.nr);
1184 for (cg = 0; cg < cgs->nr; cg++)
1186 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1192 for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
1194 il = &molt->ilist[ftype];
1195 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE))
1197 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1198 at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])