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 by 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.
39 * \brief Defines the high-level constraint code.
41 * \author Berk Hess <hess@kth.se>
42 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 * \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/settle.h"
67 #include "gromacs/mdlib/shake.h"
68 #include "gromacs/mdtypes/commrec.h"
69 #include "gromacs/mdtypes/inputrec.h"
70 #include "gromacs/mdtypes/md_enums.h"
71 #include "gromacs/mdtypes/mdatom.h"
72 #include "gromacs/mdtypes/state.h"
73 #include "gromacs/pbcutil/pbc.h"
74 #include "gromacs/pulling/pull.h"
75 #include "gromacs/timing/wallcycle.h"
76 #include "gromacs/topology/ifunc.h"
77 #include "gromacs/topology/mtop_lookup.h"
78 #include "gromacs/topology/mtop_util.h"
79 #include "gromacs/utility/exceptions.h"
80 #include "gromacs/utility/fatalerror.h"
81 #include "gromacs/utility/gmxassert.h"
82 #include "gromacs/utility/listoflists.h"
83 #include "gromacs/utility/pleasecite.h"
84 #include "gromacs/utility/txtdump.h"
89 /* \brief Impl class for Constraints
91 * \todo Members like md, idef are valid only for the lifetime of a
92 * domain, which would be good to make clearer in the structure of the
93 * code. It should not be possible to call apply() if setConstraints()
94 * has not been called. For example, this could be achieved if
95 * setConstraints returned a valid object with such a method. That
96 * still requires that we manage the lifetime of that object
97 * correctly, however. */
98 class Constraints::Impl
101 Impl(const gmx_mtop_t
& mtop_p
,
102 const t_inputrec
& ir_p
,
105 const t_mdatoms
& md_p
,
106 const t_commrec
* cr_p
,
107 const gmx_multisim_t
* ms
,
109 gmx_wallcycle
* wcycle_p
,
110 bool pbcHandlingRequired
,
114 void setConstraints(gmx_localtop_t
* top
, const t_mdatoms
& md
);
115 bool apply(bool bLog
,
120 ArrayRefWithPadding
<RVec
> x
,
121 ArrayRefWithPadding
<RVec
> xprime
,
122 ArrayRef
<RVec
> min_proj
,
126 ArrayRefWithPadding
<RVec
> v
,
128 ConstraintVariable econq
);
129 //! The total number of constraints.
131 //! The number of flexible constraints.
133 //! A list of atoms to constraints for each moleculetype.
134 std::vector
<ListOfLists
<int>> at2con_mt
;
135 //! A list of atoms to settles for each moleculetype
136 std::vector
<std::vector
<int>> at2settle_mt
;
138 Lincs
* lincsd
= nullptr; // TODO this should become a unique_ptr
140 std::unique_ptr
<shakedata
> shaked
;
142 settledata
* settled
= nullptr;
143 //! The maximum number of warnings.
145 //! The number of warnings for LINCS.
146 int warncount_lincs
= 0;
147 //! The number of warnings for SETTLE.
148 int warncount_settle
= 0;
149 //! The essential dynamics data.
150 gmx_edsam
* ed
= nullptr;
152 //! Thread-local virial contribution.
153 tensor
* vir_r_m_dr_th
= { nullptr };
154 //! Did a settle error occur?
155 bool* bSettleErrorHasOccurred
= nullptr;
157 //! Pointer to the global topology - only used for printing warnings.
158 const gmx_mtop_t
& mtop
;
159 //! Parameters for the interactions in this domain.
160 const InteractionDefinitions
* idef
= nullptr;
161 //! Data about atoms in this domain.
163 //! Whether we need to do pbc for handling bonds.
164 bool pbcHandlingRequired_
= false;
168 //! Communication support.
169 const t_commrec
* cr
= nullptr;
170 //! Multi-sim support.
171 const gmx_multisim_t
* ms
= nullptr;
172 //! Pulling code object, if any.
173 pull_t
* pull_work
= nullptr;
174 /*!\brief Input options.
176 * \todo Replace with IMdpOptions */
177 const t_inputrec
& ir
;
178 //! Flop counting support.
179 t_nrnb
* nrnb
= nullptr;
180 //! Tracks wallcycle usage.
181 gmx_wallcycle
* wcycle
;
184 Constraints::~Constraints() = default;
186 int Constraints::numFlexibleConstraints() const
188 return impl_
->nflexcon
;
191 bool Constraints::havePerturbedConstraints() const
193 const gmx_ffparams_t
& ffparams
= impl_
->mtop
.ffparams
;
195 for (size_t i
= 0; i
< ffparams
.functype
.size(); i
++)
197 if ((ffparams
.functype
[i
] == F_CONSTR
|| ffparams
.functype
[i
] == F_CONSTRNC
)
198 && ffparams
.iparams
[i
].constr
.dA
!= ffparams
.iparams
[i
].constr
.dB
)
207 //! Clears constraint quantities for atoms in nonlocal region.
208 static void clear_constraint_quantity_nonlocal(gmx_domdec_t
* dd
, ArrayRef
<RVec
> q
)
210 int nonlocal_at_start
, nonlocal_at_end
, at
;
212 dd_get_constraint_range(dd
, &nonlocal_at_start
, &nonlocal_at_end
);
214 for (at
= nonlocal_at_start
; at
< nonlocal_at_end
; at
++)
220 void too_many_constraint_warnings(int eConstrAlg
, int warncount
)
224 "Too many %s warnings (%d)\n"
225 "If you know what you are doing you can %s"
226 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
227 "but normally it is better to fix the problem",
228 (eConstrAlg
== econtLINCS
) ? "LINCS" : "SETTLE", warncount
,
229 (eConstrAlg
== econtLINCS
) ? "adjust the lincs warning threshold in your mdp file\nor " : "\n");
232 //! Writes out coordinates.
233 static void write_constr_pdb(const char* fn
,
235 const gmx_mtop_t
& mtop
,
239 ArrayRef
<const RVec
> x
,
244 int dd_ac0
= 0, dd_ac1
= 0, i
, ii
, resnr
;
246 const char * anm
, *resnm
;
249 if (DOMAINDECOMP(cr
))
252 dd_get_constraint_range(dd
, &dd_ac0
, &dd_ac1
);
259 sprintf(fname
, "%s_n%d.pdb", fn
, cr
->sim_nodeid
);
263 sprintf(fname
, "%s.pdb", fn
);
266 out
= gmx_fio_fopen(fname
, "w");
268 fprintf(out
, "TITLE %s\n", title
);
269 gmx_write_pdb_box(out
, PbcType::Unset
, box
);
271 for (i
= start
; i
< start
+ homenr
; i
++)
275 if (i
>= dd_numHomeAtoms(*dd
) && i
< dd_ac0
)
279 ii
= dd
->globalAtomIndices
[i
];
285 mtopGetAtomAndResidueName(mtop
, ii
, &molb
, &anm
, &resnr
, &resnm
, nullptr);
286 gmx_fprintf_pdb_atomline(out
, epdbATOM
, ii
+ 1, anm
, ' ', resnm
, ' ', resnr
, ' ',
287 10 * x
[i
][XX
], 10 * x
[i
][YY
], 10 * x
[i
][ZZ
], 1.0, 0.0, "");
289 fprintf(out
, "TER\n");
294 //! Writes out domain contents to help diagnose crashes.
295 static void dump_confs(FILE* log
,
297 const gmx_mtop_t
& mtop
,
301 ArrayRef
<const RVec
> x
,
302 ArrayRef
<const RVec
> xprime
,
305 char buf
[STRLEN
], buf2
[22];
307 char* env
= getenv("GMX_SUPPRESS_DUMP");
313 sprintf(buf
, "step%sb", gmx_step_str(step
, buf2
));
314 write_constr_pdb(buf
, "initial coordinates", mtop
, start
, homenr
, cr
, x
, box
);
315 sprintf(buf
, "step%sc", gmx_step_str(step
, buf2
));
316 write_constr_pdb(buf
, "coordinates after constraining", mtop
, start
, homenr
, cr
, xprime
, box
);
319 fprintf(log
, "Wrote pdb files with previous and current coordinates\n");
321 fprintf(stderr
, "Wrote pdb files with previous and current coordinates\n");
324 bool Constraints::apply(bool bLog
,
329 ArrayRefWithPadding
<RVec
> x
,
330 ArrayRefWithPadding
<RVec
> xprime
,
331 ArrayRef
<RVec
> min_proj
,
335 ArrayRefWithPadding
<RVec
> v
,
337 ConstraintVariable econq
)
339 return impl_
->apply(bLog
, bEner
, step
, delta_step
, step_scaling
, std::move(x
), std::move(xprime
),
340 min_proj
, box
, lambda
, dvdlambda
, std::move(v
), vir
, econq
);
343 bool Constraints::Impl::apply(bool bLog
,
348 ArrayRefWithPadding
<RVec
> x
,
349 ArrayRefWithPadding
<RVec
> xprime
,
350 ArrayRef
<RVec
> min_proj
,
354 ArrayRefWithPadding
<RVec
> v
,
356 ConstraintVariable econq
)
362 real invdt
, vir_fac
= 0, t
;
364 t_pbc pbc
, *pbc_null
;
368 wallcycle_start(wcycle
, ewcCONSTR
);
370 if (econq
== ConstraintVariable::ForceDispl
&& !EI_ENERGY_MINIMIZATION(ir
.eI
))
373 "constrain called for forces displacements while not doing energy minimization, "
374 "can not do this while the LINCS and SETTLE constraint connection matrices are "
384 scaled_delta_t
= step_scaling
* ir
.delta_t
;
386 /* Prepare time step for use in constraint implementations, and
387 avoid generating inf when ir.delta_t = 0. */
394 invdt
= 1.0 / scaled_delta_t
;
397 if (ir
.efep
!= efepNO
&& EI_DYNAMICS(ir
.eI
))
399 /* Set the constraint lengths for the step at which this configuration
400 * is meant to be. The invmasses should not be changed.
402 lambda
+= delta_step
* ir
.fepvals
->delta_lambda
;
407 clear_mat(vir_r_m_dr
);
409 const InteractionList
& settle
= idef
->il
[F_SETTLE
];
410 nsettle
= settle
.size() / (1 + NRAL(F_SETTLE
));
414 nth
= gmx_omp_nthreads_get(emntSETTLE
);
421 /* We do not need full pbc when constraints do not cross update groups
422 * i.e. when dd->constraint_comm==NULL.
423 * Note that PBC for constraints is different from PBC for bondeds.
424 * For constraints there is both forward and backward communication.
426 if (ir
.pbcType
!= PbcType::No
&& (cr
->dd
|| pbcHandlingRequired_
)
427 && !(cr
->dd
&& cr
->dd
->constraint_comm
== nullptr))
429 /* With pbc=screw the screw has been changed to a shift
430 * by the constraint coordinate communication routine,
431 * so that here we can use normal pbc.
433 pbc_null
= set_pbc_dd(&pbc
, ir
.pbcType
, DOMAINDECOMP(cr
) ? cr
->dd
->numCells
: nullptr, FALSE
, box
);
440 /* Communicate the coordinates required for the non-local constraints
441 * for LINCS and/or SETTLE.
445 dd_move_x_constraints(cr
->dd
, box
, x
.unpaddedArrayRef(), xprime
.unpaddedArrayRef(),
446 econq
== ConstraintVariable::Positions
);
450 /* We need to initialize the non-local components of v.
451 * We never actually use these values, but we do increment them,
452 * so we should avoid uninitialized variables and overflows.
454 clear_constraint_quantity_nonlocal(cr
->dd
, v
.unpaddedArrayRef());
458 if (lincsd
!= nullptr)
460 bOK
= constrain_lincs(bLog
|| bEner
, ir
, step
, lincsd
, md
.invmass
, cr
, ms
, x
, xprime
,
461 min_proj
, box
, pbc_null
, md
.nMassPerturbed
!= 0, lambda
, dvdlambda
,
462 invdt
, v
.unpaddedArrayRef(), vir
!= nullptr, vir_r_m_dr
, econq
, nrnb
,
463 maxwarn
, &warncount_lincs
);
464 if (!bOK
&& maxwarn
< INT_MAX
)
468 fprintf(log
, "Constraint error in algorithm %s at step %s\n",
469 econstr_names
[econtLINCS
], gmx_step_str(step
, buf
));
475 if (shaked
!= nullptr)
477 bOK
= constrain_shake(log
, shaked
.get(), md
.invmass
, *idef
, ir
, x
.unpaddedArrayRef(),
478 xprime
.unpaddedArrayRef(), min_proj
, pbc_null
, nrnb
, lambda
,
479 dvdlambda
, invdt
, v
.unpaddedArrayRef(), vir
!= nullptr, vir_r_m_dr
,
480 maxwarn
< INT_MAX
, econq
);
482 if (!bOK
&& maxwarn
< INT_MAX
)
486 fprintf(log
, "Constraint error in algorithm %s at step %s\n",
487 econstr_names
[econtSHAKE
], gmx_step_str(step
, buf
));
495 bool bSettleErrorHasOccurred0
= false;
499 case ConstraintVariable::Positions
:
500 #pragma omp parallel for num_threads(nth) schedule(static)
501 for (int th
= 0; th
< nth
; th
++)
507 clear_mat(vir_r_m_dr_th
[th
]);
510 csettle(settled
, nth
, th
, pbc_null
, x
, xprime
, invdt
, v
, vir
!= nullptr,
511 th
== 0 ? vir_r_m_dr
: vir_r_m_dr_th
[th
],
512 th
== 0 ? &bSettleErrorHasOccurred0
: &bSettleErrorHasOccurred
[th
]);
514 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
516 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
519 inc_nrnb(nrnb
, eNR_CONSTR_V
, nsettle
* 3);
523 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, nsettle
* 3);
526 case ConstraintVariable::Velocities
:
527 case ConstraintVariable::Derivative
:
528 case ConstraintVariable::Force
:
529 case ConstraintVariable::ForceDispl
:
530 #pragma omp parallel for num_threads(nth) schedule(static)
531 for (int th
= 0; th
< nth
; th
++)
535 int calcvir_atom_end
;
539 calcvir_atom_end
= 0;
543 calcvir_atom_end
= md
.homenr
;
548 clear_mat(vir_r_m_dr_th
[th
]);
551 int start_th
= (nsettle
* th
) / nth
;
552 int end_th
= (nsettle
* (th
+ 1)) / nth
;
554 if (start_th
>= 0 && end_th
- start_th
> 0)
556 settle_proj(settled
, econq
, end_th
- start_th
,
557 settle
.iatoms
.data() + start_th
* (1 + NRAL(F_SETTLE
)), pbc_null
,
558 x
.unpaddedArrayRef(), xprime
.unpaddedArrayRef(), min_proj
,
559 calcvir_atom_end
, th
== 0 ? vir_r_m_dr
: vir_r_m_dr_th
[th
]);
562 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
564 /* This is an overestimate */
565 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
567 case ConstraintVariable::Deriv_FlexCon
:
568 /* Nothing to do, since the are no flexible constraints in settles */
570 default: gmx_incons("Unknown constraint quantity for settle");
575 /* Reduce the virial contributions over the threads */
576 for (int th
= 1; th
< nth
; th
++)
578 m_add(vir_r_m_dr
, vir_r_m_dr_th
[th
], vir_r_m_dr
);
582 if (econq
== ConstraintVariable::Positions
)
584 for (int th
= 1; th
< nth
; th
++)
586 bSettleErrorHasOccurred0
= bSettleErrorHasOccurred0
|| bSettleErrorHasOccurred
[th
];
589 if (bSettleErrorHasOccurred0
)
595 ": One or more water molecules can not be settled.\n"
596 "Check for bad contacts and/or reduce the timestep if appropriate.\n",
600 fprintf(log
, "%s", buf
);
602 fprintf(stderr
, "%s", buf
);
604 if (warncount_settle
> maxwarn
)
606 too_many_constraint_warnings(-1, warncount_settle
);
617 /* The normal uses of constrain() pass step_scaling = 1.0.
618 * The call to constrain() for SD1 that passes step_scaling =
619 * 0.5 also passes vir = NULL, so cannot reach this
620 * assertion. This assertion should remain until someone knows
621 * that this path works for their intended purpose, and then
622 * they can use scaled_delta_t instead of ir.delta_t
624 assert(gmx_within_tol(step_scaling
, 1.0, GMX_REAL_EPS
));
627 case ConstraintVariable::Positions
: vir_fac
= 0.5 / (ir
.delta_t
* ir
.delta_t
); break;
628 case ConstraintVariable::Velocities
: vir_fac
= 0.5 / ir
.delta_t
; break;
629 case ConstraintVariable::Force
:
630 case ConstraintVariable::ForceDispl
: vir_fac
= 0.5; break;
631 default: gmx_incons("Unsupported constraint quantity for virial");
636 vir_fac
*= 2; /* only constraining over half the distance here */
638 for (int i
= 0; i
< DIM
; i
++)
640 for (int j
= 0; j
< DIM
; j
++)
642 (*vir
)[i
][j
] = vir_fac
* vir_r_m_dr
[i
][j
];
649 dump_confs(log
, step
, mtop
, start
, homenr
, cr
, x
.unpaddedArrayRef(),
650 xprime
.unpaddedArrayRef(), box
);
653 if (econq
== ConstraintVariable::Positions
)
655 if (ir
.bPull
&& pull_have_constraint(pull_work
))
657 if (EI_DYNAMICS(ir
.eI
))
659 t
= ir
.init_t
+ (step
+ delta_step
) * ir
.delta_t
;
665 set_pbc(&pbc
, ir
.pbcType
, box
);
666 pull_constraint(pull_work
, &md
, &pbc
, cr
, ir
.delta_t
, t
,
667 as_rvec_array(x
.unpaddedArrayRef().data()),
668 as_rvec_array(xprime
.unpaddedArrayRef().data()),
669 as_rvec_array(v
.unpaddedArrayRef().data()), *vir
);
671 if (ed
&& delta_step
> 0)
673 /* apply the essential dynamics constraints here */
674 do_edsam(&ir
, step
, cr
, as_rvec_array(xprime
.unpaddedArrayRef().data()),
675 as_rvec_array(v
.unpaddedArrayRef().data()), box
, ed
);
678 wallcycle_stop(wcycle
, ewcCONSTR
);
680 if (!v
.empty() && md
.cFREEZE
)
682 /* Set the velocities of frozen dimensions to zero */
683 ArrayRef
<RVec
> vRef
= v
.unpaddedArrayRef();
685 int gmx_unused numThreads
= gmx_omp_nthreads_get(emntUpdate
);
687 #pragma omp parallel for num_threads(numThreads) schedule(static)
688 for (int i
= 0; i
< md
.homenr
; i
++)
690 int freezeGroup
= md
.cFREEZE
[i
];
692 for (int d
= 0; d
< DIM
; d
++)
694 if (ir
.opts
.nFreeze
[freezeGroup
][d
])
705 ArrayRef
<real
> Constraints::rmsdData() const
709 return lincs_rmsdData(impl_
->lincsd
);
717 real
Constraints::rmsd() const
721 return lincs_rmsd(impl_
->lincsd
);
729 int Constraints::numConstraintsTotal()
731 return impl_
->ncon_tot
;
734 FlexibleConstraintTreatment
flexibleConstraintTreatment(bool haveDynamicsIntegrator
)
736 if (haveDynamicsIntegrator
)
738 return FlexibleConstraintTreatment::Include
;
742 return FlexibleConstraintTreatment::Exclude
;
746 /*! \brief Returns a block struct to go from atoms to constraints
748 * The block struct will contain constraint indices with lower indices
749 * directly matching the order in F_CONSTR and higher indices matching
750 * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
752 * \param[in] numAtoms The number of atoms to construct the list for
753 * \param[in] ilists The interaction lists, size F_NRE
754 * \param[in] iparams Interaction parameters, can be null when
755 * \p flexibleConstraintTreatment==Include
756 * \param[in] flexibleConstraintTreatment The flexible constraint treatment,
759 * \returns a block struct with all constraints for each atom
761 static ListOfLists
<int> makeAtomsToConstraintsList(int numAtoms
,
762 ArrayRef
<const InteractionList
> ilists
,
763 ArrayRef
<const t_iparams
> iparams
,
764 FlexibleConstraintTreatment flexibleConstraintTreatment
)
766 GMX_ASSERT(flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
|| !iparams
.empty(),
767 "With flexible constraint detection we need valid iparams");
769 std::vector
<int> count(numAtoms
);
771 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
773 const InteractionList
& ilist
= ilists
[ftype
];
774 const int stride
= 1 + NRAL(ftype
);
775 for (int i
= 0; i
< ilist
.size(); i
+= stride
)
777 if (flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
778 || !isConstraintFlexible(iparams
, ilist
.iatoms
[i
]))
780 for (int j
= 1; j
< 3; j
++)
782 int a
= ilist
.iatoms
[i
+ j
];
789 std::vector
<int> listRanges(numAtoms
+ 1);
790 for (int a
= 0; a
< numAtoms
; a
++)
792 listRanges
[a
+ 1] = listRanges
[a
] + count
[a
];
795 std::vector
<int> elements(listRanges
[numAtoms
]);
797 /* The F_CONSTRNC constraints have constraint numbers
798 * that continue after the last F_CONSTR constraint.
800 int numConstraints
= 0;
801 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
803 const InteractionList
& ilist
= ilists
[ftype
];
804 const int stride
= 1 + NRAL(ftype
);
805 for (int i
= 0; i
< ilist
.size(); i
+= stride
)
807 if (flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
808 || !isConstraintFlexible(iparams
, ilist
.iatoms
[i
]))
810 for (int j
= 1; j
< 3; j
++)
812 const int a
= ilist
.iatoms
[i
+ j
];
813 elements
[listRanges
[a
] + count
[a
]++] = numConstraints
;
820 return ListOfLists
<int>(std::move(listRanges
), std::move(elements
));
823 ListOfLists
<int> make_at2con(int numAtoms
,
824 ArrayRef
<const InteractionList
> ilist
,
825 ArrayRef
<const t_iparams
> iparams
,
826 FlexibleConstraintTreatment flexibleConstraintTreatment
)
828 return makeAtomsToConstraintsList(numAtoms
, ilist
, iparams
, flexibleConstraintTreatment
);
831 ListOfLists
<int> make_at2con(const gmx_moltype_t
& moltype
,
832 gmx::ArrayRef
<const t_iparams
> iparams
,
833 FlexibleConstraintTreatment flexibleConstraintTreatment
)
835 return makeAtomsToConstraintsList(moltype
.atoms
.nr
, makeConstArrayRef(moltype
.ilist
), iparams
,
836 flexibleConstraintTreatment
);
839 //! Return the number of flexible constraints in the \c ilist and \c iparams.
840 int countFlexibleConstraints(ArrayRef
<const InteractionList
> ilist
, ArrayRef
<const t_iparams
> iparams
)
843 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
845 const int numIatomsPerConstraint
= 3;
846 for (int i
= 0; i
< ilist
[ftype
].size(); i
+= numIatomsPerConstraint
)
848 const int type
= ilist
[ftype
].iatoms
[i
];
849 if (iparams
[type
].constr
.dA
== 0 && iparams
[type
].constr
.dB
== 0)
859 //! Returns the index of the settle to which each atom belongs.
860 static std::vector
<int> make_at2settle(int natoms
, const InteractionList
& ilist
)
862 /* Set all to no settle */
863 std::vector
<int> at2s(natoms
, -1);
865 const int stride
= 1 + NRAL(F_SETTLE
);
867 for (int s
= 0; s
< ilist
.size(); s
+= stride
)
869 at2s
[ilist
.iatoms
[s
+ 1]] = s
/ stride
;
870 at2s
[ilist
.iatoms
[s
+ 2]] = s
/ stride
;
871 at2s
[ilist
.iatoms
[s
+ 3]] = s
/ stride
;
877 void Constraints::Impl::setConstraints(gmx_localtop_t
* top
, const t_mdatoms
& md
)
883 /* With DD we might also need to call LINCS on a domain no constraints for
884 * communicating coordinates to other nodes that do have constraints.
886 if (ir
.eConstrAlg
== econtLINCS
)
888 set_lincs(*idef
, md
.nr
, md
.invmass
, md
.lambda
, EI_DYNAMICS(ir
.eI
), cr
, lincsd
);
890 if (ir
.eConstrAlg
== econtSHAKE
)
894 // We are using the local topology, so there are only
895 // F_CONSTR constraints.
896 make_shake_sblock_dd(shaked
.get(), idef
->il
[F_CONSTR
]);
900 make_shake_sblock_serial(shaked
.get(), &top
->idef
, md
.nr
);
907 settle_set_constraints(settled
, idef
->il
[F_SETTLE
], md
);
910 /* Make a selection of the local atoms for essential dynamics */
913 dd_make_local_ed_indices(cr
->dd
, ed
);
917 void Constraints::setConstraints(gmx_localtop_t
* top
, const t_mdatoms
& md
)
919 impl_
->setConstraints(top
, md
);
922 /*! \brief Makes a per-moleculetype container of mappings from atom
923 * indices to constraint indices.
925 * Note that flexible constraints are only enabled with a dynamical integrator. */
926 static std::vector
<ListOfLists
<int>> makeAtomToConstraintMappings(const gmx_mtop_t
& mtop
,
927 FlexibleConstraintTreatment flexibleConstraintTreatment
)
929 std::vector
<ListOfLists
<int>> mapping
;
930 mapping
.reserve(mtop
.moltype
.size());
931 for (const gmx_moltype_t
& moltype
: mtop
.moltype
)
933 mapping
.push_back(make_at2con(moltype
, mtop
.ffparams
.iparams
, flexibleConstraintTreatment
));
938 Constraints::Constraints(const gmx_mtop_t
& mtop
,
939 const t_inputrec
& ir
,
944 const gmx_multisim_t
* ms
,
946 gmx_wallcycle
* wcycle
,
947 bool pbcHandlingRequired
,
950 impl_(new Impl(mtop
, ir
, pull_work
, log
, md
, cr
, ms
, nrnb
, wcycle
, pbcHandlingRequired
, numConstraints
, numSettles
))
954 Constraints::Impl::Impl(const gmx_mtop_t
& mtop_p
,
955 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
),
973 pull_work(pull_work
),
978 if (numConstraints
+ numSettles
> 0 && ir
.epc
== epcMTTK
)
980 gmx_fatal(FARGS
, "Constraints are not implemented with MTTK pressure control.");
984 if (numConstraints
> 0)
986 at2con_mt
= makeAtomToConstraintMappings(mtop
, flexibleConstraintTreatment(EI_DYNAMICS(ir
.eI
)));
988 for (const gmx_molblock_t
& molblock
: mtop
.molblock
)
990 int count
= countFlexibleConstraints(mtop
.moltype
[molblock
.type
].ilist
, mtop
.ffparams
.iparams
);
991 nflexcon
+= molblock
.nmol
* count
;
998 fprintf(log
, "There are %d flexible constraints\n", nflexcon
);
999 if (ir
.fc_stepsize
== 0)
1003 "WARNING: step size for flexible constraining = 0\n"
1004 " All flexible constraints will be rigid.\n"
1005 " Will try to keep all flexible constraints at their original "
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
, nflexcon
, at2con_mt
,
1020 DOMAINDECOMP(cr
) && ddHaveSplitConstraints(*cr
->dd
), ir
.nLincsIter
,
1024 if (ir
.eConstrAlg
== econtSHAKE
)
1026 if (DOMAINDECOMP(cr
) && ddHaveSplitConstraints(*cr
->dd
))
1029 "SHAKE is not supported with domain decomposition and constraint that "
1030 "cross domain boundaries, use LINCS");
1035 "For this system also velocities and/or forces need to be constrained, "
1036 "this can not be done with SHAKE, you should select LINCS");
1038 please_cite(log
, "Ryckaert77a");
1041 please_cite(log
, "Barth95a");
1044 shaked
= std::make_unique
<shakedata
>();
1050 please_cite(log
, "Miyamoto92a");
1052 settled
= settle_init(mtop
);
1054 /* Make an atom to settle index for use in domain decomposition */
1055 for (size_t mt
= 0; mt
< mtop
.moltype
.size(); mt
++)
1057 at2settle_mt
.emplace_back(
1058 make_at2settle(mtop
.moltype
[mt
].atoms
.nr
, mtop
.moltype
[mt
].ilist
[F_SETTLE
]));
1061 /* Allocate thread-local work arrays */
1062 int nthreads
= gmx_omp_nthreads_get(emntSETTLE
);
1063 if (nthreads
> 1 && vir_r_m_dr_th
== nullptr)
1065 snew(vir_r_m_dr_th
, nthreads
);
1066 snew(bSettleErrorHasOccurred
, nthreads
);
1071 char* env
= getenv("GMX_MAXCONSTRWARN");
1075 sscanf(env
, "%8d", &maxwarn
);
1082 fprintf(log
, "Setting the maximum number of constraint warnings to %d\n", maxwarn
);
1086 fprintf(stderr
, "Setting the maximum number of constraint warnings to %d\n", maxwarn
);
1089 warncount_lincs
= 0;
1090 warncount_settle
= 0;
1093 Constraints::Impl::~Impl()
1095 if (bSettleErrorHasOccurred
!= nullptr)
1097 sfree(bSettleErrorHasOccurred
);
1099 if (vir_r_m_dr_th
!= nullptr)
1101 sfree(vir_r_m_dr_th
);
1103 if (settled
!= nullptr)
1105 settle_free(settled
);
1110 void Constraints::saveEdsamPointer(gmx_edsam
* ed
)
1115 ArrayRef
<const ListOfLists
<int>> Constraints::atom2constraints_moltype() const
1117 return impl_
->at2con_mt
;
1120 ArrayRef
<const std::vector
<int>> Constraints::atom2settle_moltype() const
1122 return impl_
->at2settle_mt
;
1125 void do_constrain_first(FILE* fplog
,
1126 gmx::Constraints
* constr
,
1127 const t_inputrec
* ir
,
1128 const t_mdatoms
* md
,
1130 ArrayRefWithPadding
<RVec
> x
,
1131 ArrayRefWithPadding
<RVec
> v
,
1135 int i
, m
, start
, end
;
1137 real dt
= ir
->delta_t
;
1140 PaddedVector
<RVec
> savex(natoms
);
1147 fprintf(debug
, "vcm: start=%d, homenr=%d, end=%d\n", start
, md
->homenr
, end
);
1149 /* Do a first constrain to reset particles... */
1150 step
= ir
->init_step
;
1153 char buf
[STEPSTRSIZE
];
1154 fprintf(fplog
, "\nConstraining the starting coordinates (step %s)\n", gmx_step_str(step
, buf
));
1158 /* constrain the current position */
1159 constr
->apply(TRUE
, FALSE
, step
, 0, 1.0, x
, x
, {}, box
, lambda
, &dvdl_dum
, {}, nullptr,
1160 gmx::ConstraintVariable::Positions
);
1163 /* constrain the inital velocity, and save it */
1164 /* also may be useful if we need the ekin from the halfstep for velocity verlet */
1165 constr
->apply(TRUE
, FALSE
, step
, 0, 1.0, x
, v
, v
.unpaddedArrayRef(), box
, lambda
, &dvdl_dum
,
1166 {}, nullptr, gmx::ConstraintVariable::Velocities
);
1168 /* constrain the inital velocities at t-dt/2 */
1169 if (EI_STATE_VELOCITY(ir
->eI
) && ir
->eI
!= eiVV
)
1171 auto subX
= x
.paddedArrayRef().subArray(start
, end
);
1172 auto subV
= v
.paddedArrayRef().subArray(start
, end
);
1173 for (i
= start
; (i
< end
); i
++)
1175 for (m
= 0; (m
< DIM
); m
++)
1177 /* Reverse the velocity */
1178 subV
[i
][m
] = -subV
[i
][m
];
1179 /* Store the position at t-dt in buf */
1180 savex
[i
][m
] = subX
[i
][m
] + dt
* subV
[i
][m
];
1183 /* Shake the positions at t=-dt with the positions at t=0
1184 * as reference coordinates.
1188 char buf
[STEPSTRSIZE
];
1189 fprintf(fplog
, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step
, buf
));
1192 constr
->apply(TRUE
, FALSE
, step
, -1, 1.0, x
, savex
.arrayRefWithPadding(), {}, box
, lambda
,
1193 &dvdl_dum
, v
, nullptr, gmx::ConstraintVariable::Positions
);
1195 for (i
= start
; i
< end
; i
++)
1197 for (m
= 0; m
< DIM
; m
++)
1199 /* Re-reverse the velocities */
1200 subV
[i
][m
] = -subV
[i
][m
];