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.
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/essentialdynamics/edsam.h"
51 #include "gromacs/fileio/confio.h"
52 #include "gromacs/fileio/gmxfio.h"
53 #include "gromacs/fileio/pdbio.h"
54 #include "gromacs/gmxlib/nrnb.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/gmx_omp_nthreads.h"
58 #include "gromacs/mdlib/lincs.h"
59 #include "gromacs/mdlib/mdrun.h"
60 #include "gromacs/mdlib/settle.h"
61 #include "gromacs/mdlib/shake.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/mdtypes/mdatom.h"
66 #include "gromacs/mdtypes/state.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/pulling/pull.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/pleasecite.h"
76 #include "gromacs/utility/smalloc.h"
77 #include "gromacs/utility/txtdump.h"
79 typedef struct gmx_constr
{
80 int ncon_tot
; /* The total number of constraints */
81 int nflexcon
; /* The number of flexible constraints */
82 int n_at2con_mt
; /* The size of at2con = #moltypes */
83 t_blocka
*at2con_mt
; /* A list of atoms to constraints */
84 int n_at2settle_mt
; /* The size of at2settle = #moltypes */
85 int **at2settle_mt
; /* A list of atoms to settles */
86 gmx_bool bInterCGsettles
;
87 gmx_lincsdata_t lincsd
; /* LINCS data */
88 gmx_shakedata_t shaked
; /* SHAKE data */
89 gmx_settledata_t settled
; /* SETTLE data */
90 int maxwarn
; /* The maximum number of warnings */
93 gmx_edsam_t ed
; /* The essential dynamics data */
95 /* Thread local working data */
96 tensor
*vir_r_m_dr_th
; /* Thread virial contribution */
97 bool *bSettleErrorHasOccurred
; /* Did a settle error occur? */
99 /* Only used for printing warnings */
100 const gmx_mtop_t
*warn_mtop
; /* Pointer to the global topology */
103 int n_flexible_constraints(struct gmx_constr
*constr
)
109 nflexcon
= constr
->nflexcon
;
119 static void clear_constraint_quantity_nonlocal(gmx_domdec_t
*dd
, rvec
*q
)
121 int nonlocal_at_start
, nonlocal_at_end
, at
;
123 dd_get_constraint_range(dd
, &nonlocal_at_start
, &nonlocal_at_end
);
125 for (at
= nonlocal_at_start
; at
< nonlocal_at_end
; at
++)
131 void too_many_constraint_warnings(int eConstrAlg
, int warncount
)
134 "Too many %s warnings (%d)\n"
135 "If you know what you are doing you can %s"
136 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
137 "but normally it is better to fix the problem",
138 (eConstrAlg
== econtLINCS
) ? "LINCS" : "SETTLE", warncount
,
139 (eConstrAlg
== econtLINCS
) ?
140 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
143 static void write_constr_pdb(const char *fn
, const char *title
,
144 const gmx_mtop_t
*mtop
,
145 int start
, int homenr
, const t_commrec
*cr
,
146 rvec x
[], matrix box
)
150 int dd_ac0
= 0, dd_ac1
= 0, i
, ii
, resnr
;
152 const char *anm
, *resnm
;
155 if (DOMAINDECOMP(cr
))
158 dd_get_constraint_range(dd
, &dd_ac0
, &dd_ac1
);
165 sprintf(fname
, "%s_n%d.pdb", fn
, cr
->sim_nodeid
);
169 sprintf(fname
, "%s.pdb", fn
);
172 out
= gmx_fio_fopen(fname
, "w");
174 fprintf(out
, "TITLE %s\n", title
);
175 gmx_write_pdb_box(out
, -1, box
);
177 for (i
= start
; i
< start
+homenr
; i
++)
181 if (i
>= dd
->nat_home
&& i
< dd_ac0
)
185 ii
= dd
->gatindex
[i
];
191 mtopGetAtomAndResidueName(mtop
, ii
, &molb
, &anm
, &resnr
, &resnm
, nullptr);
192 gmx_fprintf_pdb_atomline(out
, epdbATOM
, ii
+1, anm
, ' ', resnm
, ' ', resnr
, ' ',
193 10*x
[i
][XX
], 10*x
[i
][YY
], 10*x
[i
][ZZ
], 1.0, 0.0, "");
195 fprintf(out
, "TER\n");
200 static void dump_confs(FILE *fplog
, gmx_int64_t step
, const gmx_mtop_t
*mtop
,
201 int start
, int homenr
, const t_commrec
*cr
,
202 rvec x
[], rvec xprime
[], matrix box
)
204 char buf
[STRLEN
], buf2
[22];
206 char *env
= getenv("GMX_SUPPRESS_DUMP");
212 sprintf(buf
, "step%sb", gmx_step_str(step
, buf2
));
213 write_constr_pdb(buf
, "initial coordinates",
214 mtop
, start
, homenr
, cr
, x
, box
);
215 sprintf(buf
, "step%sc", gmx_step_str(step
, buf2
));
216 write_constr_pdb(buf
, "coordinates after constraining",
217 mtop
, start
, homenr
, cr
, xprime
, box
);
220 fprintf(fplog
, "Wrote pdb files with previous and current coordinates\n");
222 fprintf(stderr
, "Wrote pdb files with previous and current coordinates\n");
225 gmx_bool
constrain(FILE *fplog
, gmx_bool bLog
, gmx_bool bEner
,
226 struct gmx_constr
*constr
,
227 t_idef
*idef
, const t_inputrec
*ir
,
229 const gmx_multisim_t
*ms
,
230 gmx_int64_t step
, int delta_step
,
233 rvec
*x
, rvec
*xprime
, rvec
*min_proj
,
234 gmx_bool bMolPBC
, matrix box
,
235 real lambda
, real
*dvdlambda
,
236 rvec
*v
, tensor
*vir
,
237 t_nrnb
*nrnb
, int econq
)
243 real invdt
, vir_fac
= 0, t
;
246 t_pbc pbc
, *pbc_null
;
250 if (econq
== econqForceDispl
&& !EI_ENERGY_MINIMIZATION(ir
->eI
))
252 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");
261 scaled_delta_t
= step_scaling
* ir
->delta_t
;
263 /* Prepare time step for use in constraint implementations, and
264 avoid generating inf when ir->delta_t = 0. */
265 if (ir
->delta_t
== 0)
271 invdt
= 1.0/scaled_delta_t
;
274 if (ir
->efep
!= efepNO
&& EI_DYNAMICS(ir
->eI
))
276 /* Set the constraint lengths for the step at which this configuration
277 * is meant to be. The invmasses should not be changed.
279 lambda
+= delta_step
*ir
->fepvals
->delta_lambda
;
284 clear_mat(vir_r_m_dr
);
289 settle
= &idef
->il
[F_SETTLE
];
290 nsettle
= settle
->nr
/(1+NRAL(F_SETTLE
));
294 nth
= gmx_omp_nthreads_get(emntSETTLE
);
301 /* We do not need full pbc when constraints do not cross charge groups,
302 * i.e. when dd->constraint_comm==NULL.
303 * Note that PBC for constraints is different from PBC for bondeds.
304 * For constraints there is both forward and backward communication.
306 if (ir
->ePBC
!= epbcNONE
&&
307 (cr
->dd
|| bMolPBC
) && !(cr
->dd
&& cr
->dd
->constraint_comm
== nullptr))
309 /* With pbc=screw the screw has been changed to a shift
310 * by the constraint coordinate communication routine,
311 * so that here we can use normal pbc.
313 pbc_null
= set_pbc_dd(&pbc
, ir
->ePBC
,
314 DOMAINDECOMP(cr
) ? cr
->dd
->nc
: nullptr,
322 /* Communicate the coordinates required for the non-local constraints
323 * for LINCS and/or SETTLE.
327 dd_move_x_constraints(cr
->dd
, box
, x
, xprime
, econq
== econqCoord
);
331 /* We need to initialize the non-local components of v.
332 * We never actually use these values, but we do increment them,
333 * so we should avoid uninitialized variables and overflows.
335 clear_constraint_quantity_nonlocal(cr
->dd
, v
);
339 if (constr
->lincsd
!= nullptr)
341 bOK
= constrain_lincs(fplog
, bLog
, bEner
, ir
, step
, constr
->lincsd
, md
, cr
, ms
,
343 box
, pbc_null
, lambda
, dvdlambda
,
344 invdt
, v
, vir
!= nullptr, vir_r_m_dr
,
346 constr
->maxwarn
, &constr
->warncount_lincs
);
347 if (!bOK
&& constr
->maxwarn
< INT_MAX
)
349 if (fplog
!= nullptr)
351 fprintf(fplog
, "Constraint error in algorithm %s at step %s\n",
352 econstr_names
[econtLINCS
], gmx_step_str(step
, buf
));
358 if (constr
->shaked
!= nullptr)
360 bOK
= constrain_shake(fplog
, constr
->shaked
,
362 idef
, ir
, x
, xprime
, min_proj
, nrnb
,
364 invdt
, v
, vir
!= nullptr, vir_r_m_dr
,
365 constr
->maxwarn
< INT_MAX
, econq
);
367 if (!bOK
&& constr
->maxwarn
< INT_MAX
)
369 if (fplog
!= nullptr)
371 fprintf(fplog
, "Constraint error in algorithm %s at step %s\n",
372 econstr_names
[econtSHAKE
], gmx_step_str(step
, buf
));
380 bool bSettleErrorHasOccurred
= false;
385 #pragma omp parallel for num_threads(nth) schedule(static)
386 for (th
= 0; th
< nth
; th
++)
392 clear_mat(constr
->vir_r_m_dr_th
[th
]);
395 csettle(constr
->settled
,
399 invdt
, v
? v
[0] : nullptr,
401 th
== 0 ? vir_r_m_dr
: constr
->vir_r_m_dr_th
[th
],
402 th
== 0 ? &bSettleErrorHasOccurred
: &constr
->bSettleErrorHasOccurred
[th
]);
404 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
406 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
409 inc_nrnb(nrnb
, eNR_CONSTR_V
, nsettle
*3);
413 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, nsettle
*3);
419 case econqForceDispl
:
420 #pragma omp parallel for num_threads(nth) schedule(static)
421 for (th
= 0; th
< nth
; th
++)
425 int calcvir_atom_end
;
429 calcvir_atom_end
= 0;
433 calcvir_atom_end
= md
->homenr
;
438 clear_mat(constr
->vir_r_m_dr_th
[th
]);
441 int start_th
= (nsettle
* th
)/nth
;
442 int end_th
= (nsettle
*(th
+1))/nth
;
444 if (start_th
>= 0 && end_th
- start_th
> 0)
446 settle_proj(constr
->settled
, econq
,
448 settle
->iatoms
+start_th
*(1+NRAL(F_SETTLE
)),
451 xprime
, min_proj
, calcvir_atom_end
,
452 th
== 0 ? vir_r_m_dr
: constr
->vir_r_m_dr_th
[th
]);
455 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
457 /* This is an overestimate */
458 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
460 case econqDeriv_FlexCon
:
461 /* Nothing to do, since the are no flexible constraints in settles */
464 gmx_incons("Unknown constraint quantity for settle");
469 /* Reduce the virial contributions over the threads */
470 for (int th
= 1; th
< nth
; th
++)
472 m_add(vir_r_m_dr
, constr
->vir_r_m_dr_th
[th
], vir_r_m_dr
);
476 if (econq
== econqCoord
)
478 for (int th
= 1; th
< nth
; th
++)
480 bSettleErrorHasOccurred
= bSettleErrorHasOccurred
|| constr
->bSettleErrorHasOccurred
[th
];
483 if (bSettleErrorHasOccurred
)
487 "\nstep " "%" GMX_PRId64
": One or more water molecules can not be settled.\n"
488 "Check for bad contacts and/or reduce the timestep if appropriate.\n",
492 fprintf(fplog
, "%s", buf
);
494 fprintf(stderr
, "%s", buf
);
495 constr
->warncount_settle
++;
496 if (constr
->warncount_settle
> constr
->maxwarn
)
498 too_many_constraint_warnings(-1, constr
->warncount_settle
);
509 /* The normal uses of constrain() pass step_scaling = 1.0.
510 * The call to constrain() for SD1 that passes step_scaling =
511 * 0.5 also passes vir = NULL, so cannot reach this
512 * assertion. This assertion should remain until someone knows
513 * that this path works for their intended purpose, and then
514 * they can use scaled_delta_t instead of ir->delta_t
516 assert(gmx_within_tol(step_scaling
, 1.0, GMX_REAL_EPS
));
520 vir_fac
= 0.5/(ir
->delta_t
*ir
->delta_t
);
523 vir_fac
= 0.5/ir
->delta_t
;
526 case econqForceDispl
:
530 gmx_incons("Unsupported constraint quantity for virial");
535 vir_fac
*= 2; /* only constraining over half the distance here */
537 for (int i
= 0; i
< DIM
; i
++)
539 for (int j
= 0; j
< DIM
; j
++)
541 (*vir
)[i
][j
] = vir_fac
*vir_r_m_dr
[i
][j
];
548 dump_confs(fplog
, step
, constr
->warn_mtop
, start
, homenr
, cr
, x
, xprime
, box
);
551 if (econq
== econqCoord
)
553 if (ir
->bPull
&& pull_have_constraint(ir
->pull_work
))
555 if (EI_DYNAMICS(ir
->eI
))
557 t
= ir
->init_t
+ (step
+ delta_step
)*ir
->delta_t
;
563 set_pbc(&pbc
, ir
->ePBC
, box
);
564 pull_constraint(ir
->pull_work
, md
, &pbc
, cr
, ir
->delta_t
, t
, x
, xprime
, v
, *vir
);
566 if (constr
->ed
&& delta_step
> 0)
568 /* apply the essential dynamics constraints here */
569 do_edsam(ir
, step
, cr
, xprime
, v
, box
, constr
->ed
);
576 real
*constr_rmsd_data(struct gmx_constr
*constr
)
580 return lincs_rmsd_data(constr
->lincsd
);
588 real
constr_rmsd(struct gmx_constr
*constr
)
592 return lincs_rmsd(constr
->lincsd
);
600 t_blocka
make_at2con(int start
, int natoms
,
601 const t_ilist
*ilist
, const t_iparams
*iparams
,
602 gmx_bool bDynamics
, int *nflexiblecons
)
604 int *count
, ncon
, con
, con_tot
, nflexcon
, ftype
, i
, a
;
611 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
613 ncon
= ilist
[ftype
].nr
/3;
614 ia
= ilist
[ftype
].iatoms
;
615 for (con
= 0; con
< ncon
; con
++)
617 bFlexCon
= (iparams
[ia
[0]].constr
.dA
== 0 &&
618 iparams
[ia
[0]].constr
.dB
== 0);
623 if (bDynamics
|| !bFlexCon
)
625 for (i
= 1; i
< 3; i
++)
634 *nflexiblecons
= nflexcon
;
637 at2con
.nalloc_index
= at2con
.nr
+1;
638 snew(at2con
.index
, at2con
.nalloc_index
);
640 for (a
= 0; a
< natoms
; a
++)
642 at2con
.index
[a
+1] = at2con
.index
[a
] + count
[a
];
645 at2con
.nra
= at2con
.index
[natoms
];
646 at2con
.nalloc_a
= at2con
.nra
;
647 snew(at2con
.a
, at2con
.nalloc_a
);
649 /* The F_CONSTRNC constraints have constraint numbers
650 * that continue after the last F_CONSTR constraint.
653 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
655 ncon
= ilist
[ftype
].nr
/3;
656 ia
= ilist
[ftype
].iatoms
;
657 for (con
= 0; con
< ncon
; con
++)
659 bFlexCon
= (iparams
[ia
[0]].constr
.dA
== 0 &&
660 iparams
[ia
[0]].constr
.dB
== 0);
661 if (bDynamics
|| !bFlexCon
)
663 for (i
= 1; i
< 3; i
++)
666 at2con
.a
[at2con
.index
[a
]+count
[a
]++] = con_tot
;
679 static int *make_at2settle(int natoms
, const t_ilist
*ilist
)
685 /* Set all to no settle */
686 for (a
= 0; a
< natoms
; a
++)
691 stride
= 1 + NRAL(F_SETTLE
);
693 for (s
= 0; s
< ilist
->nr
; s
+= stride
)
695 at2s
[ilist
->iatoms
[s
+1]] = s
/stride
;
696 at2s
[ilist
->iatoms
[s
+2]] = s
/stride
;
697 at2s
[ilist
->iatoms
[s
+3]] = s
/stride
;
703 void set_constraints(struct gmx_constr
*constr
,
704 gmx_localtop_t
*top
, const t_inputrec
*ir
,
705 const t_mdatoms
*md
, const t_commrec
*cr
)
707 t_idef
*idef
= &top
->idef
;
709 if (constr
->ncon_tot
> 0)
711 /* With DD we might also need to call LINCS on a domain no constraints for
712 * communicating coordinates to other nodes that do have constraints.
714 if (ir
->eConstrAlg
== econtLINCS
)
716 set_lincs(idef
, md
, EI_DYNAMICS(ir
->eI
), cr
, constr
->lincsd
);
718 if (ir
->eConstrAlg
== econtSHAKE
)
722 // We are using the local topology, so there are only
723 // F_CONSTR constraints.
724 make_shake_sblock_dd(constr
->shaked
, &idef
->il
[F_CONSTR
], &top
->cgs
, cr
->dd
);
728 make_shake_sblock_serial(constr
->shaked
, idef
, md
);
735 settle_set_constraints(constr
->settled
,
736 &idef
->il
[F_SETTLE
], md
);
739 /* Make a selection of the local atoms for essential dynamics */
740 if (constr
->ed
&& cr
->dd
)
742 dd_make_local_ed_indices(cr
->dd
, constr
->ed
);
746 gmx_constr_t
init_constraints(FILE *fplog
,
747 const gmx_mtop_t
*mtop
, const t_inputrec
*ir
,
748 bool doEssentialDynamics
,
752 gmx_mtop_ftype_count(mtop
, F_CONSTR
) +
753 gmx_mtop_ftype_count(mtop
, F_CONSTRNC
);
755 gmx_mtop_ftype_count(mtop
, F_SETTLE
);
757 GMX_RELEASE_ASSERT(!ir
->bPull
|| ir
->pull_work
!= nullptr, "init_constraints called with COM pulling before/without initializing the pull code");
759 if (nconstraints
+ nsettles
== 0 &&
760 !(ir
->bPull
&& pull_have_constraint(ir
->pull_work
)) &&
761 !doEssentialDynamics
)
766 struct gmx_constr
*constr
;
770 constr
->ncon_tot
= nconstraints
;
771 constr
->nflexcon
= 0;
772 if (nconstraints
> 0)
774 constr
->n_at2con_mt
= mtop
->moltype
.size();
775 snew(constr
->at2con_mt
, constr
->n_at2con_mt
);
776 for (int mt
= 0; mt
< static_cast<int>(mtop
->moltype
.size()); mt
++)
779 constr
->at2con_mt
[mt
] = make_at2con(0, mtop
->moltype
[mt
].atoms
.nr
,
780 mtop
->moltype
[mt
].ilist
,
781 mtop
->ffparams
.iparams
,
782 EI_DYNAMICS(ir
->eI
), &nflexcon
);
783 for (const gmx_molblock_t
&molblock
: mtop
->molblock
)
785 if (molblock
.type
== mt
)
787 constr
->nflexcon
+= molblock
.nmol
*nflexcon
;
792 if (constr
->nflexcon
> 0)
796 fprintf(fplog
, "There are %d flexible constraints\n",
798 if (ir
->fc_stepsize
== 0)
801 "WARNING: step size for flexible constraining = 0\n"
802 " All flexible constraints will be rigid.\n"
803 " Will try to keep all flexible constraints at their original length,\n"
804 " but the lengths may exhibit some drift.\n\n");
805 constr
->nflexcon
= 0;
808 if (constr
->nflexcon
> 0)
810 please_cite(fplog
, "Hess2002");
814 if (ir
->eConstrAlg
== econtLINCS
)
816 constr
->lincsd
= init_lincs(fplog
, mtop
,
817 constr
->nflexcon
, constr
->at2con_mt
,
818 DOMAINDECOMP(cr
) && cr
->dd
->bInterCGcons
,
819 ir
->nLincsIter
, ir
->nProjOrder
);
822 if (ir
->eConstrAlg
== econtSHAKE
)
824 if (DOMAINDECOMP(cr
) && cr
->dd
->bInterCGcons
)
826 gmx_fatal(FARGS
, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
828 if (constr
->nflexcon
)
830 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");
832 please_cite(fplog
, "Ryckaert77a");
835 please_cite(fplog
, "Barth95a");
838 constr
->shaked
= shake_init();
844 please_cite(fplog
, "Miyamoto92a");
846 constr
->bInterCGsettles
= inter_charge_group_settles(mtop
);
848 constr
->settled
= settle_init(mtop
);
850 /* Make an atom to settle index for use in domain decomposition */
851 constr
->n_at2settle_mt
= mtop
->moltype
.size();
852 snew(constr
->at2settle_mt
, constr
->n_at2settle_mt
);
853 for (size_t mt
= 0; mt
< mtop
->moltype
.size(); mt
++)
855 constr
->at2settle_mt
[mt
] =
856 make_at2settle(mtop
->moltype
[mt
].atoms
.nr
,
857 &mtop
->moltype
[mt
].ilist
[F_SETTLE
]);
860 /* Allocate thread-local work arrays */
861 int nthreads
= gmx_omp_nthreads_get(emntSETTLE
);
862 if (nthreads
> 1 && constr
->vir_r_m_dr_th
== nullptr)
864 snew(constr
->vir_r_m_dr_th
, nthreads
);
865 snew(constr
->bSettleErrorHasOccurred
, nthreads
);
869 if (nconstraints
+ nsettles
> 0 && ir
->epc
== epcMTTK
)
871 gmx_fatal(FARGS
, "Constraints are not implemented with MTTK pressure control.");
874 constr
->maxwarn
= 999;
875 char *env
= getenv("GMX_MAXCONSTRWARN");
879 sscanf(env
, "%8d", &constr
->maxwarn
);
880 if (constr
->maxwarn
< 0)
882 constr
->maxwarn
= INT_MAX
;
887 "Setting the maximum number of constraint warnings to %d\n",
893 "Setting the maximum number of constraint warnings to %d\n",
897 constr
->warncount_lincs
= 0;
898 constr
->warncount_settle
= 0;
900 constr
->warn_mtop
= mtop
;
905 /* Put a pointer to the essential dynamics constraints into the constr struct */
906 void saveEdsamPointer(gmx_constr_t constr
, gmx_edsam_t ed
)
911 const t_blocka
*atom2constraints_moltype(gmx_constr_t constr
)
913 return constr
->at2con_mt
;
916 const int **atom2settle_moltype(gmx_constr_t constr
)
918 return (const int **)constr
->at2settle_mt
;
922 gmx_bool
inter_charge_group_constraints(const gmx_mtop_t
*mtop
)
924 const gmx_moltype_t
*molt
;
927 int *at2cg
, cg
, a
, ftype
, i
;
931 for (size_t mb
= 0; mb
< mtop
->molblock
.size() && !bInterCG
; mb
++)
933 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
935 if (molt
->ilist
[F_CONSTR
].nr
> 0 ||
936 molt
->ilist
[F_CONSTRNC
].nr
> 0 ||
937 molt
->ilist
[F_SETTLE
].nr
> 0)
940 snew(at2cg
, molt
->atoms
.nr
);
941 for (cg
= 0; cg
< cgs
->nr
; cg
++)
943 for (a
= cgs
->index
[cg
]; a
< cgs
->index
[cg
+1]; a
++)
949 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
951 il
= &molt
->ilist
[ftype
];
952 for (i
= 0; i
< il
->nr
&& !bInterCG
; i
+= 1+NRAL(ftype
))
954 if (at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+2]])
968 gmx_bool
inter_charge_group_settles(const gmx_mtop_t
*mtop
)
970 const gmx_moltype_t
*molt
;
973 int *at2cg
, cg
, a
, ftype
, i
;
977 for (size_t mb
= 0; mb
< mtop
->molblock
.size() && !bInterCG
; mb
++)
979 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
981 if (molt
->ilist
[F_SETTLE
].nr
> 0)
984 snew(at2cg
, molt
->atoms
.nr
);
985 for (cg
= 0; cg
< cgs
->nr
; cg
++)
987 for (a
= cgs
->index
[cg
]; a
< cgs
->index
[cg
+1]; a
++)
993 for (ftype
= F_SETTLE
; ftype
<= F_SETTLE
; ftype
++)
995 il
= &molt
->ilist
[ftype
];
996 for (i
= 0; i
< il
->nr
&& !bInterCG
; i
+= 1+NRAL(F_SETTLE
))
998 if (at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+2]] ||
999 at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+3]])