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/state.h"
72 #include "gromacs/pbcutil/pbc.h"
73 #include "gromacs/pulling/pull.h"
74 #include "gromacs/timing/wallcycle.h"
75 #include "gromacs/topology/ifunc.h"
76 #include "gromacs/topology/mtop_lookup.h"
77 #include "gromacs/topology/mtop_util.h"
78 #include "gromacs/utility/exceptions.h"
79 #include "gromacs/utility/fatalerror.h"
80 #include "gromacs/utility/gmxassert.h"
81 #include "gromacs/utility/listoflists.h"
82 #include "gromacs/utility/pleasecite.h"
83 #include "gromacs/utility/txtdump.h"
88 /* \brief Impl class for Constraints
90 * \todo Members like md, idef are valid only for the lifetime of a
91 * domain, which would be good to make clearer in the structure of the
92 * code. It should not be possible to call apply() if setConstraints()
93 * has not been called. For example, this could be achieved if
94 * setConstraints returned a valid object with such a method. That
95 * still requires that we manage the lifetime of that object
96 * correctly, however. */
97 class Constraints::Impl
100 Impl(const gmx_mtop_t
& mtop_p
,
101 const t_inputrec
& ir_p
,
104 const t_commrec
* cr_p
,
105 const gmx_multisim_t
* ms
,
107 gmx_wallcycle
* wcycle_p
,
108 bool pbcHandlingRequired
,
112 void setConstraints(gmx_localtop_t
* top
,
117 bool hasMassPerturbedAtoms
,
119 unsigned short* cFREEZE
);
120 bool apply(bool bLog
,
125 ArrayRefWithPadding
<RVec
> x
,
126 ArrayRefWithPadding
<RVec
> xprime
,
127 ArrayRef
<RVec
> min_proj
,
131 ArrayRefWithPadding
<RVec
> v
,
133 tensor constraintsVirial
,
134 ConstraintVariable econq
);
135 //! The total number of constraints.
137 //! The number of flexible constraints.
139 //! A list of atoms to constraints for each moleculetype.
140 std::vector
<ListOfLists
<int>> at2con_mt
;
141 //! A list of atoms to settles for each moleculetype
142 std::vector
<std::vector
<int>> at2settle_mt
;
144 Lincs
* lincsd
= nullptr; // TODO this should become a unique_ptr
146 std::unique_ptr
<shakedata
> shaked
;
148 std::unique_ptr
<SettleData
> settled
;
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
* threadConstraintsVirial
= { 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 InteractionDefinitions
* idef
= nullptr;
167 //! Total number of atoms.
169 //! Number of local atoms.
170 int numHomeAtoms_
= 0;
174 real
* inverseMasses_
;
175 //! If there are atoms with perturbed mass (for FEP).
176 bool hasMassPerturbedAtoms_
;
177 //! FEP lambda value.
179 //! Freeze groups data
180 unsigned short* cFREEZE_
;
181 //! Whether we need to do pbc for handling bonds.
182 bool pbcHandlingRequired_
= false;
186 //! Communication support.
187 const t_commrec
* cr
= nullptr;
188 //! Multi-sim support.
189 const gmx_multisim_t
* ms
= nullptr;
190 //! Pulling code object, if any.
191 pull_t
* pull_work
= nullptr;
192 /*!\brief Input options.
194 * \todo Replace with IMdpOptions */
195 const t_inputrec
& ir
;
196 //! Flop counting support.
197 t_nrnb
* nrnb
= nullptr;
198 //! Tracks wallcycle usage.
199 gmx_wallcycle
* wcycle
;
202 Constraints::~Constraints() = default;
204 int Constraints::numFlexibleConstraints() const
206 return impl_
->nflexcon
;
209 bool Constraints::havePerturbedConstraints() const
211 const gmx_ffparams_t
& ffparams
= impl_
->mtop
.ffparams
;
213 for (size_t i
= 0; i
< ffparams
.functype
.size(); i
++)
215 if ((ffparams
.functype
[i
] == F_CONSTR
|| ffparams
.functype
[i
] == F_CONSTRNC
)
216 && ffparams
.iparams
[i
].constr
.dA
!= ffparams
.iparams
[i
].constr
.dB
)
225 //! Clears constraint quantities for atoms in nonlocal region.
226 static void clear_constraint_quantity_nonlocal(gmx_domdec_t
* dd
, ArrayRef
<RVec
> q
)
228 int nonlocal_at_start
, nonlocal_at_end
, at
;
230 dd_get_constraint_range(dd
, &nonlocal_at_start
, &nonlocal_at_end
);
232 for (at
= nonlocal_at_start
; at
< nonlocal_at_end
; at
++)
238 void too_many_constraint_warnings(int eConstrAlg
, int warncount
)
242 "Too many %s warnings (%d)\n"
243 "If you know what you are doing you can %s"
244 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
245 "but normally it is better to fix the problem",
246 (eConstrAlg
== econtLINCS
) ? "LINCS" : "SETTLE", warncount
,
247 (eConstrAlg
== econtLINCS
) ? "adjust the lincs warning threshold in your mdp file\nor " : "\n");
250 //! Writes out coordinates.
251 static void write_constr_pdb(const char* fn
,
253 const gmx_mtop_t
& mtop
,
257 ArrayRef
<const RVec
> x
,
262 int dd_ac0
= 0, dd_ac1
= 0, i
, ii
, resnr
;
264 const char * anm
, *resnm
;
267 if (DOMAINDECOMP(cr
))
270 dd_get_constraint_range(dd
, &dd_ac0
, &dd_ac1
);
277 sprintf(fname
, "%s_n%d.pdb", fn
, cr
->sim_nodeid
);
281 sprintf(fname
, "%s.pdb", fn
);
284 out
= gmx_fio_fopen(fname
, "w");
286 fprintf(out
, "TITLE %s\n", title
);
287 gmx_write_pdb_box(out
, PbcType::Unset
, box
);
289 for (i
= start
; i
< start
+ homenr
; i
++)
293 if (i
>= dd_numHomeAtoms(*dd
) && i
< dd_ac0
)
297 ii
= dd
->globalAtomIndices
[i
];
303 mtopGetAtomAndResidueName(mtop
, ii
, &molb
, &anm
, &resnr
, &resnm
, nullptr);
304 gmx_fprintf_pdb_atomline(out
, epdbATOM
, ii
+ 1, anm
, ' ', resnm
, ' ', resnr
, ' ',
305 10 * x
[i
][XX
], 10 * x
[i
][YY
], 10 * x
[i
][ZZ
], 1.0, 0.0, "");
307 fprintf(out
, "TER\n");
312 //! Writes out domain contents to help diagnose crashes.
313 static void dump_confs(FILE* log
,
315 const gmx_mtop_t
& mtop
,
319 ArrayRef
<const RVec
> x
,
320 ArrayRef
<const RVec
> xprime
,
323 char buf
[STRLEN
], buf2
[22];
325 char* env
= getenv("GMX_SUPPRESS_DUMP");
331 sprintf(buf
, "step%sb", gmx_step_str(step
, buf2
));
332 write_constr_pdb(buf
, "initial coordinates", mtop
, start
, homenr
, cr
, x
, box
);
333 sprintf(buf
, "step%sc", gmx_step_str(step
, buf2
));
334 write_constr_pdb(buf
, "coordinates after constraining", mtop
, start
, homenr
, cr
, xprime
, box
);
337 fprintf(log
, "Wrote pdb files with previous and current coordinates\n");
339 fprintf(stderr
, "Wrote pdb files with previous and current coordinates\n");
342 bool Constraints::apply(bool bLog
,
347 ArrayRefWithPadding
<RVec
> x
,
348 ArrayRefWithPadding
<RVec
> xprime
,
349 ArrayRef
<RVec
> min_proj
,
353 ArrayRefWithPadding
<RVec
> v
,
355 tensor constraintsVirial
,
356 ConstraintVariable econq
)
358 return impl_
->apply(bLog
, bEner
, step
, delta_step
, step_scaling
, std::move(x
),
359 std::move(xprime
), min_proj
, box
, lambda
, dvdlambda
, std::move(v
),
360 computeVirial
, constraintsVirial
, econq
);
363 bool Constraints::Impl::apply(bool bLog
,
368 ArrayRefWithPadding
<RVec
> x
,
369 ArrayRefWithPadding
<RVec
> xprime
,
370 ArrayRef
<RVec
> min_proj
,
374 ArrayRefWithPadding
<RVec
> v
,
376 tensor constraintsVirial
,
377 ConstraintVariable econq
)
382 real invdt
, vir_fac
= 0, t
;
384 t_pbc pbc
, *pbc_null
;
388 wallcycle_start(wcycle
, ewcCONSTR
);
390 if (econq
== ConstraintVariable::ForceDispl
&& !EI_ENERGY_MINIMIZATION(ir
.eI
))
393 "constrain called for forces displacements while not doing energy minimization, "
394 "can not do this while the LINCS and SETTLE constraint connection matrices are "
403 scaled_delta_t
= step_scaling
* ir
.delta_t
;
405 /* Prepare time step for use in constraint implementations, and
406 avoid generating inf when ir.delta_t = 0. */
413 invdt
= 1.0 / scaled_delta_t
;
416 if (ir
.efep
!= efepNO
&& EI_DYNAMICS(ir
.eI
))
418 /* Set the constraint lengths for the step at which this configuration
419 * is meant to be. The invmasses should not be changed.
421 lambda
+= delta_step
* ir
.fepvals
->delta_lambda
;
426 clear_mat(constraintsVirial
);
428 const InteractionList
& settle
= idef
->il
[F_SETTLE
];
429 nsettle
= settle
.size() / (1 + NRAL(F_SETTLE
));
433 nth
= gmx_omp_nthreads_get(emntSETTLE
);
440 /* We do not need full pbc when constraints do not cross update groups
441 * i.e. when dd->constraint_comm==NULL.
442 * Note that PBC for constraints is different from PBC for bondeds.
443 * For constraints there is both forward and backward communication.
445 if (ir
.pbcType
!= PbcType::No
&& (cr
->dd
|| pbcHandlingRequired_
)
446 && !(cr
->dd
&& cr
->dd
->constraint_comm
== nullptr))
448 /* With pbc=screw the screw has been changed to a shift
449 * by the constraint coordinate communication routine,
450 * so that here we can use normal pbc.
452 pbc_null
= set_pbc_dd(&pbc
, ir
.pbcType
, DOMAINDECOMP(cr
) ? cr
->dd
->numCells
: nullptr, FALSE
, box
);
459 /* Communicate the coordinates required for the non-local constraints
460 * for LINCS and/or SETTLE.
464 dd_move_x_constraints(cr
->dd
, box
, x
.unpaddedArrayRef(), xprime
.unpaddedArrayRef(),
465 econq
== ConstraintVariable::Positions
);
469 /* We need to initialize the non-local components of v.
470 * We never actually use these values, but we do increment them,
471 * so we should avoid uninitialized variables and overflows.
473 clear_constraint_quantity_nonlocal(cr
->dd
, v
.unpaddedArrayRef());
477 if (lincsd
!= nullptr)
479 bOK
= constrain_lincs(bLog
|| bEner
, ir
, step
, lincsd
, inverseMasses_
, cr
, ms
, x
, xprime
,
480 min_proj
, box
, pbc_null
, hasMassPerturbedAtoms_
, lambda
, dvdlambda
,
481 invdt
, v
.unpaddedArrayRef(), computeVirial
, constraintsVirial
, econq
,
482 nrnb
, maxwarn
, &warncount_lincs
);
483 if (!bOK
&& maxwarn
< INT_MAX
)
487 fprintf(log
, "Constraint error in algorithm %s at step %s\n",
488 econstr_names
[econtLINCS
], gmx_step_str(step
, buf
));
494 if (shaked
!= nullptr)
496 bOK
= constrain_shake(log
, shaked
.get(), inverseMasses_
, *idef
, ir
, x
.unpaddedArrayRef(),
497 xprime
.unpaddedArrayRef(), min_proj
, pbc_null
, nrnb
, lambda
,
498 dvdlambda
, invdt
, v
.unpaddedArrayRef(), computeVirial
,
499 constraintsVirial
, maxwarn
< INT_MAX
, econq
);
501 if (!bOK
&& maxwarn
< INT_MAX
)
505 fprintf(log
, "Constraint error in algorithm %s at step %s\n",
506 econstr_names
[econtSHAKE
], gmx_step_str(step
, buf
));
514 bool bSettleErrorHasOccurred0
= false;
518 case ConstraintVariable::Positions
:
519 #pragma omp parallel for num_threads(nth) schedule(static)
520 for (int th
= 0; th
< nth
; th
++)
526 clear_mat(threadConstraintsVirial
[th
]);
529 csettle(*settled
, nth
, th
, pbc_null
, x
, xprime
, invdt
, v
, computeVirial
,
530 th
== 0 ? constraintsVirial
: threadConstraintsVirial
[th
],
531 th
== 0 ? &bSettleErrorHasOccurred0
: &bSettleErrorHasOccurred
[th
]);
533 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
535 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
538 inc_nrnb(nrnb
, eNR_CONSTR_V
, nsettle
* 3);
542 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, nsettle
* 3);
545 case ConstraintVariable::Velocities
:
546 case ConstraintVariable::Derivative
:
547 case ConstraintVariable::Force
:
548 case ConstraintVariable::ForceDispl
:
549 #pragma omp parallel for num_threads(nth) schedule(static)
550 for (int th
= 0; th
< nth
; th
++)
554 int calcvir_atom_end
;
558 calcvir_atom_end
= 0;
562 calcvir_atom_end
= numHomeAtoms_
;
567 clear_mat(threadConstraintsVirial
[th
]);
570 int start_th
= (nsettle
* th
) / nth
;
571 int end_th
= (nsettle
* (th
+ 1)) / nth
;
573 if (start_th
>= 0 && end_th
- start_th
> 0)
575 settle_proj(*settled
, econq
, end_th
- start_th
,
576 settle
.iatoms
.data() + start_th
* (1 + NRAL(F_SETTLE
)),
577 pbc_null
, x
.unpaddedArrayRef(), xprime
.unpaddedArrayRef(),
578 min_proj
, calcvir_atom_end
,
579 th
== 0 ? constraintsVirial
: threadConstraintsVirial
[th
]);
582 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
584 /* This is an overestimate */
585 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
587 case ConstraintVariable::Deriv_FlexCon
:
588 /* Nothing to do, since the are no flexible constraints in settles */
590 default: gmx_incons("Unknown constraint quantity for settle");
595 /* Reduce the virial contributions over the threads */
596 for (int th
= 1; th
< nth
; th
++)
598 m_add(constraintsVirial
, threadConstraintsVirial
[th
], constraintsVirial
);
602 if (econq
== ConstraintVariable::Positions
)
604 for (int th
= 1; th
< nth
; th
++)
606 bSettleErrorHasOccurred0
= bSettleErrorHasOccurred0
|| bSettleErrorHasOccurred
[th
];
609 if (bSettleErrorHasOccurred0
)
615 ": One or more water molecules can not be settled.\n"
616 "Check for bad contacts and/or reduce the timestep if appropriate.\n",
620 fprintf(log
, "%s", buf
);
622 fprintf(stderr
, "%s", buf
);
624 if (warncount_settle
> maxwarn
)
626 too_many_constraint_warnings(-1, warncount_settle
);
637 /* The normal uses of constrain() pass step_scaling = 1.0.
638 * The call to constrain() for SD1 that passes step_scaling =
639 * 0.5 also passes vir = NULL, so cannot reach this
640 * assertion. This assertion should remain until someone knows
641 * that this path works for their intended purpose, and then
642 * they can use scaled_delta_t instead of ir.delta_t
644 assert(gmx_within_tol(step_scaling
, 1.0, GMX_REAL_EPS
));
647 case ConstraintVariable::Positions
: vir_fac
= 0.5 / (ir
.delta_t
* ir
.delta_t
); break;
648 case ConstraintVariable::Velocities
: vir_fac
= 0.5 / ir
.delta_t
; break;
649 case ConstraintVariable::Force
:
650 case ConstraintVariable::ForceDispl
: vir_fac
= 0.5; break;
651 default: gmx_incons("Unsupported constraint quantity for virial");
656 vir_fac
*= 2; /* only constraining over half the distance here */
658 for (int i
= 0; i
< DIM
; i
++)
660 for (int j
= 0; j
< DIM
; j
++)
662 constraintsVirial
[i
][j
] *= vir_fac
;
669 dump_confs(log
, step
, mtop
, start
, numHomeAtoms_
, cr
, x
.unpaddedArrayRef(),
670 xprime
.unpaddedArrayRef(), box
);
673 if (econq
== ConstraintVariable::Positions
)
675 if (ir
.bPull
&& pull_have_constraint(pull_work
))
677 if (EI_DYNAMICS(ir
.eI
))
679 t
= ir
.init_t
+ (step
+ delta_step
) * ir
.delta_t
;
685 set_pbc(&pbc
, ir
.pbcType
, box
);
686 pull_constraint(pull_work
, masses_
, &pbc
, cr
, ir
.delta_t
, t
,
687 as_rvec_array(x
.unpaddedArrayRef().data()),
688 as_rvec_array(xprime
.unpaddedArrayRef().data()),
689 as_rvec_array(v
.unpaddedArrayRef().data()), constraintsVirial
);
691 if (ed
&& delta_step
> 0)
693 /* apply the essential dynamics constraints here */
694 do_edsam(&ir
, step
, cr
, as_rvec_array(xprime
.unpaddedArrayRef().data()),
695 as_rvec_array(v
.unpaddedArrayRef().data()), box
, ed
);
698 wallcycle_stop(wcycle
, ewcCONSTR
);
700 if (!v
.empty() && cFREEZE_
)
702 /* Set the velocities of frozen dimensions to zero */
703 ArrayRef
<RVec
> vRef
= v
.unpaddedArrayRef();
705 int gmx_unused numThreads
= gmx_omp_nthreads_get(emntUpdate
);
707 #pragma omp parallel for num_threads(numThreads) schedule(static)
708 for (int i
= 0; i
< numHomeAtoms_
; i
++)
710 int freezeGroup
= cFREEZE_
[i
];
712 for (int d
= 0; d
< DIM
; d
++)
714 if (ir
.opts
.nFreeze
[freezeGroup
][d
])
725 ArrayRef
<real
> Constraints::rmsdData() const
729 return lincs_rmsdData(impl_
->lincsd
);
737 real
Constraints::rmsd() const
741 return lincs_rmsd(impl_
->lincsd
);
749 int Constraints::numConstraintsTotal()
751 return impl_
->ncon_tot
;
754 FlexibleConstraintTreatment
flexibleConstraintTreatment(bool haveDynamicsIntegrator
)
756 if (haveDynamicsIntegrator
)
758 return FlexibleConstraintTreatment::Include
;
762 return FlexibleConstraintTreatment::Exclude
;
766 /*! \brief Returns a block struct to go from atoms to constraints
768 * The block struct will contain constraint indices with lower indices
769 * directly matching the order in F_CONSTR and higher indices matching
770 * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
772 * \param[in] numAtoms The number of atoms to construct the list for
773 * \param[in] ilists The interaction lists, size F_NRE
774 * \param[in] iparams Interaction parameters, can be null when
775 * \p flexibleConstraintTreatment==Include
776 * \param[in] flexibleConstraintTreatment The flexible constraint treatment,
779 * \returns a block struct with all constraints for each atom
781 static ListOfLists
<int> makeAtomsToConstraintsList(int numAtoms
,
782 ArrayRef
<const InteractionList
> ilists
,
783 ArrayRef
<const t_iparams
> iparams
,
784 FlexibleConstraintTreatment flexibleConstraintTreatment
)
786 GMX_ASSERT(flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
|| !iparams
.empty(),
787 "With flexible constraint detection we need valid iparams");
789 std::vector
<int> count(numAtoms
);
791 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
793 const InteractionList
& ilist
= ilists
[ftype
];
794 const int stride
= 1 + NRAL(ftype
);
795 for (int i
= 0; i
< ilist
.size(); i
+= stride
)
797 if (flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
798 || !isConstraintFlexible(iparams
, ilist
.iatoms
[i
]))
800 for (int j
= 1; j
< 3; j
++)
802 int a
= ilist
.iatoms
[i
+ j
];
809 std::vector
<int> listRanges(numAtoms
+ 1);
810 for (int a
= 0; a
< numAtoms
; a
++)
812 listRanges
[a
+ 1] = listRanges
[a
] + count
[a
];
815 std::vector
<int> elements(listRanges
[numAtoms
]);
817 /* The F_CONSTRNC constraints have constraint numbers
818 * that continue after the last F_CONSTR constraint.
820 int numConstraints
= 0;
821 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
823 const InteractionList
& ilist
= ilists
[ftype
];
824 const int stride
= 1 + NRAL(ftype
);
825 for (int i
= 0; i
< ilist
.size(); i
+= stride
)
827 if (flexibleConstraintTreatment
== FlexibleConstraintTreatment::Include
828 || !isConstraintFlexible(iparams
, ilist
.iatoms
[i
]))
830 for (int j
= 1; j
< 3; j
++)
832 const int a
= ilist
.iatoms
[i
+ j
];
833 elements
[listRanges
[a
] + count
[a
]++] = numConstraints
;
840 return ListOfLists
<int>(std::move(listRanges
), std::move(elements
));
843 ListOfLists
<int> make_at2con(int numAtoms
,
844 ArrayRef
<const InteractionList
> ilist
,
845 ArrayRef
<const t_iparams
> iparams
,
846 FlexibleConstraintTreatment flexibleConstraintTreatment
)
848 return makeAtomsToConstraintsList(numAtoms
, ilist
, iparams
, flexibleConstraintTreatment
);
851 ListOfLists
<int> make_at2con(const gmx_moltype_t
& moltype
,
852 gmx::ArrayRef
<const t_iparams
> iparams
,
853 FlexibleConstraintTreatment flexibleConstraintTreatment
)
855 return makeAtomsToConstraintsList(moltype
.atoms
.nr
, makeConstArrayRef(moltype
.ilist
), iparams
,
856 flexibleConstraintTreatment
);
859 //! Return the number of flexible constraints in the \c ilist and \c iparams.
860 int countFlexibleConstraints(ArrayRef
<const InteractionList
> ilist
, ArrayRef
<const t_iparams
> iparams
)
863 for (int ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
865 const int numIatomsPerConstraint
= 3;
866 for (int i
= 0; i
< ilist
[ftype
].size(); i
+= numIatomsPerConstraint
)
868 const int type
= ilist
[ftype
].iatoms
[i
];
869 if (iparams
[type
].constr
.dA
== 0 && iparams
[type
].constr
.dB
== 0)
879 //! Returns the index of the settle to which each atom belongs.
880 static std::vector
<int> make_at2settle(int natoms
, const InteractionList
& ilist
)
882 /* Set all to no settle */
883 std::vector
<int> at2s(natoms
, -1);
885 const int stride
= 1 + NRAL(F_SETTLE
);
887 for (int s
= 0; s
< ilist
.size(); s
+= stride
)
889 at2s
[ilist
.iatoms
[s
+ 1]] = s
/ stride
;
890 at2s
[ilist
.iatoms
[s
+ 2]] = s
/ stride
;
891 at2s
[ilist
.iatoms
[s
+ 3]] = s
/ stride
;
897 void Constraints::Impl::setConstraints(gmx_localtop_t
* top
,
902 const bool hasMassPerturbedAtoms
,
904 unsigned short* cFREEZE
)
906 numAtoms_
= numAtoms
;
907 numHomeAtoms_
= numHomeAtoms
;
909 inverseMasses_
= inverseMasses
;
910 hasMassPerturbedAtoms_
= hasMassPerturbedAtoms
;
918 /* With DD we might also need to call LINCS on a domain no constraints for
919 * communicating coordinates to other nodes that do have constraints.
921 if (ir
.eConstrAlg
== econtLINCS
)
923 set_lincs(*idef
, numAtoms_
, inverseMasses_
, lambda_
, EI_DYNAMICS(ir
.eI
), cr
, lincsd
);
925 if (ir
.eConstrAlg
== econtSHAKE
)
929 // We are using the local topology, so there are only
930 // F_CONSTR constraints.
931 GMX_RELEASE_ASSERT(idef
->il
[F_CONSTRNC
].empty(),
932 "Here we should not have no-connect constraints");
933 make_shake_sblock_dd(shaked
.get(), idef
->il
[F_CONSTR
]);
937 make_shake_sblock_serial(shaked
.get(), &top
->idef
, numAtoms_
);
944 settled
->setConstraints(idef
->il
[F_SETTLE
], numHomeAtoms_
, masses_
, inverseMasses_
);
947 /* Make a selection of the local atoms for essential dynamics */
950 dd_make_local_ed_indices(cr
->dd
, ed
);
954 void Constraints::setConstraints(gmx_localtop_t
* top
,
956 const int numHomeAtoms
,
959 const bool hasMassPerturbedAtoms
,
961 unsigned short* cFREEZE
)
963 impl_
->setConstraints(top
, numAtoms
, numHomeAtoms
, masses
, inverseMasses
, hasMassPerturbedAtoms
,
967 /*! \brief Makes a per-moleculetype container of mappings from atom
968 * indices to constraint indices.
970 * Note that flexible constraints are only enabled with a dynamical integrator. */
971 static std::vector
<ListOfLists
<int>> makeAtomToConstraintMappings(const gmx_mtop_t
& mtop
,
972 FlexibleConstraintTreatment flexibleConstraintTreatment
)
974 std::vector
<ListOfLists
<int>> mapping
;
975 mapping
.reserve(mtop
.moltype
.size());
976 for (const gmx_moltype_t
& moltype
: mtop
.moltype
)
978 mapping
.push_back(make_at2con(moltype
, mtop
.ffparams
.iparams
, flexibleConstraintTreatment
));
983 Constraints::Constraints(const gmx_mtop_t
& mtop
,
984 const t_inputrec
& ir
,
988 const gmx_multisim_t
* ms
,
990 gmx_wallcycle
* wcycle
,
991 bool pbcHandlingRequired
,
994 impl_(new Impl(mtop
, ir
, pull_work
, log
, cr
, ms
, nrnb
, wcycle
, pbcHandlingRequired
, numConstraints
, numSettles
))
998 Constraints::Impl::Impl(const gmx_mtop_t
& mtop_p
,
999 const t_inputrec
& ir_p
,
1002 const t_commrec
* cr_p
,
1003 const gmx_multisim_t
* ms_p
,
1005 gmx_wallcycle
* wcycle_p
,
1006 bool pbcHandlingRequired
,
1009 ncon_tot(numConstraints
),
1011 pbcHandlingRequired_(pbcHandlingRequired
),
1015 pull_work(pull_work
),
1020 if (numConstraints
+ numSettles
> 0 && ir
.epc
== epcMTTK
)
1022 gmx_fatal(FARGS
, "Constraints are not implemented with MTTK pressure control.");
1026 if (numConstraints
> 0)
1028 at2con_mt
= makeAtomToConstraintMappings(mtop
, flexibleConstraintTreatment(EI_DYNAMICS(ir
.eI
)));
1030 for (const gmx_molblock_t
& molblock
: mtop
.molblock
)
1032 int count
= countFlexibleConstraints(mtop
.moltype
[molblock
.type
].ilist
, mtop
.ffparams
.iparams
);
1033 nflexcon
+= molblock
.nmol
* count
;
1040 fprintf(log
, "There are %d flexible constraints\n", nflexcon
);
1041 if (ir
.fc_stepsize
== 0)
1045 "WARNING: step size for flexible constraining = 0\n"
1046 " All flexible constraints will be rigid.\n"
1047 " Will try to keep all flexible constraints at their original "
1049 " but the lengths may exhibit some drift.\n\n");
1055 please_cite(log
, "Hess2002");
1059 if (ir
.eConstrAlg
== econtLINCS
)
1061 lincsd
= init_lincs(log
, mtop
, nflexcon
, at2con_mt
,
1062 DOMAINDECOMP(cr
) && ddHaveSplitConstraints(*cr
->dd
), ir
.nLincsIter
,
1066 if (ir
.eConstrAlg
== econtSHAKE
)
1068 if (DOMAINDECOMP(cr
) && ddHaveSplitConstraints(*cr
->dd
))
1071 "SHAKE is not supported with domain decomposition and constraint that "
1072 "cross domain boundaries, use LINCS");
1077 "For this system also velocities and/or forces need to be constrained, "
1078 "this can not be done with SHAKE, you should select LINCS");
1080 please_cite(log
, "Ryckaert77a");
1083 please_cite(log
, "Barth95a");
1086 shaked
= std::make_unique
<shakedata
>();
1092 please_cite(log
, "Miyamoto92a");
1094 settled
= std::make_unique
<SettleData
>(mtop
);
1096 /* Make an atom to settle index for use in domain decomposition */
1097 for (size_t mt
= 0; mt
< mtop
.moltype
.size(); mt
++)
1099 at2settle_mt
.emplace_back(
1100 make_at2settle(mtop
.moltype
[mt
].atoms
.nr
, mtop
.moltype
[mt
].ilist
[F_SETTLE
]));
1103 /* Allocate thread-local work arrays */
1104 int nthreads
= gmx_omp_nthreads_get(emntSETTLE
);
1105 if (nthreads
> 1 && threadConstraintsVirial
== nullptr)
1107 snew(threadConstraintsVirial
, nthreads
);
1108 snew(bSettleErrorHasOccurred
, nthreads
);
1113 char* env
= getenv("GMX_MAXCONSTRWARN");
1117 sscanf(env
, "%8d", &maxwarn
);
1124 fprintf(log
, "Setting the maximum number of constraint warnings to %d\n", maxwarn
);
1128 fprintf(stderr
, "Setting the maximum number of constraint warnings to %d\n", maxwarn
);
1131 warncount_lincs
= 0;
1132 warncount_settle
= 0;
1135 Constraints::Impl::~Impl()
1137 if (bSettleErrorHasOccurred
!= nullptr)
1139 sfree(bSettleErrorHasOccurred
);
1141 if (threadConstraintsVirial
!= nullptr)
1143 sfree(threadConstraintsVirial
);
1148 void Constraints::saveEdsamPointer(gmx_edsam
* ed
)
1153 ArrayRef
<const ListOfLists
<int>> Constraints::atom2constraints_moltype() const
1155 return impl_
->at2con_mt
;
1158 ArrayRef
<const std::vector
<int>> Constraints::atom2settle_moltype() const
1160 return impl_
->at2settle_mt
;
1163 void do_constrain_first(FILE* fplog
,
1164 gmx::Constraints
* constr
,
1165 const t_inputrec
* ir
,
1168 ArrayRefWithPadding
<RVec
> x
,
1169 ArrayRefWithPadding
<RVec
> v
,
1173 int i
, m
, start
, end
;
1175 real dt
= ir
->delta_t
;
1178 PaddedVector
<RVec
> savex(numAtoms
);
1185 fprintf(debug
, "vcm: start=%d, homenr=%d, end=%d\n", start
, numHomeAtoms
, end
);
1187 /* Do a first constrain to reset particles... */
1188 step
= ir
->init_step
;
1191 char buf
[STEPSTRSIZE
];
1192 fprintf(fplog
, "\nConstraining the starting coordinates (step %s)\n", gmx_step_str(step
, buf
));
1196 bool needsLogging
= true;
1197 bool computeEnergy
= false;
1198 bool computeVirial
= false;
1199 /* constrain the current position */
1200 constr
->apply(needsLogging
, computeEnergy
, step
, 0, 1.0, x
, x
, {}, box
, lambda
, &dvdl_dum
, {},
1201 computeVirial
, nullptr, gmx::ConstraintVariable::Positions
);
1204 /* constrain the inital velocity, and save it */
1205 /* also may be useful if we need the ekin from the halfstep for velocity verlet */
1206 constr
->apply(needsLogging
, computeEnergy
, step
, 0, 1.0, x
, v
, v
.unpaddedArrayRef(), box
, lambda
,
1207 &dvdl_dum
, {}, computeVirial
, nullptr, gmx::ConstraintVariable::Velocities
);
1209 /* constrain the inital velocities at t-dt/2 */
1210 if (EI_STATE_VELOCITY(ir
->eI
) && ir
->eI
!= eiVV
)
1212 auto subX
= x
.paddedArrayRef().subArray(start
, end
);
1213 auto subV
= v
.paddedArrayRef().subArray(start
, end
);
1214 for (i
= start
; (i
< end
); i
++)
1216 for (m
= 0; (m
< DIM
); m
++)
1218 /* Reverse the velocity */
1219 subV
[i
][m
] = -subV
[i
][m
];
1220 /* Store the position at t-dt in buf */
1221 savex
[i
][m
] = subX
[i
][m
] + dt
* subV
[i
][m
];
1224 /* Shake the positions at t=-dt with the positions at t=0
1225 * as reference coordinates.
1229 char buf
[STEPSTRSIZE
];
1230 fprintf(fplog
, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step
, buf
));
1233 constr
->apply(needsLogging
, computeEnergy
, step
, -1, 1.0, x
, savex
.arrayRefWithPadding(),
1234 {}, box
, lambda
, &dvdl_dum
, v
, computeVirial
, nullptr,
1235 gmx::ConstraintVariable::Positions
);
1237 for (i
= start
; i
< end
; i
++)
1239 for (m
= 0; m
< DIM
; m
++)
1241 /* Re-reverse the velocities */
1242 subV
[i
][m
] = -subV
[i
][m
];
1248 void constrain_velocities(gmx::Constraints
* constr
,
1254 gmx_bool computeVirial
,
1255 tensor constraintsVirial
)
1257 if (constr
!= nullptr)
1259 constr
->apply(do_log
, do_ene
, step
, 1, 1.0, state
->x
.arrayRefWithPadding(),
1260 state
->v
.arrayRefWithPadding(), state
->v
.arrayRefWithPadding().unpaddedArrayRef(),
1261 state
->box
, state
->lambda
[efptBONDED
], dvdlambda
, ArrayRefWithPadding
<RVec
>(),
1262 computeVirial
, constraintsVirial
, ConstraintVariable::Velocities
);
1266 void constrain_coordinates(gmx::Constraints
* constr
,
1271 ArrayRefWithPadding
<RVec
> xp
,
1273 gmx_bool computeVirial
,
1274 tensor constraintsVirial
)
1276 if (constr
!= nullptr)
1278 constr
->apply(do_log
, do_ene
, step
, 1, 1.0, state
->x
.arrayRefWithPadding(), std::move(xp
),
1279 ArrayRef
<RVec
>(), state
->box
, state
->lambda
[efptBONDED
], dvdlambda
,
1280 state
->v
.arrayRefWithPadding(), computeVirial
, constraintsVirial
,
1281 ConstraintVariable::Positions
);