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.
37 /* This file is completely threadsafe - keep it that way! */
40 #include "constraintrange.h"
46 #include "gromacs/mdlib/constr.h"
47 #include "gromacs/mdtypes/inputrec.h"
48 #include "gromacs/topology/block.h"
49 #include "gromacs/topology/mtop_util.h"
50 #include "gromacs/utility/basedefinitions.h"
51 #include "gromacs/utility/fatalerror.h"
52 #include "gromacs/utility/logger.h"
53 #include "gromacs/utility/real.h"
54 #include "gromacs/utility/smalloc.h"
59 //! Recursing function to help find all adjacent constraints.
60 static void constr_recur(const t_blocka *at2con,
61 const InteractionList *ilist, const t_iparams *iparams,
63 int at, int depth, int nc, int *path,
64 real r0, real r1, real *r2max,
73 gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
74 gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
76 /* Loop over all constraints connected to this atom */
77 for (c = at2con->index[at]; c < at2con->index[at+1]; c++)
80 /* Do not walk over already used constraints */
82 for (a1 = 0; a1 < depth; a1++)
91 const int *ia = constr_iatomptr(ia1, ia2, con);
92 /* Flexible constraints currently have length 0, which is incorrect */
95 len = iparams[ia[0]].constr.dA;
99 len = iparams[ia[0]].constr.dB;
101 /* In the worst case the bond directions alternate */
112 /* Assume angles of 120 degrees between all bonds */
113 if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max)
115 *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
118 fprintf(debug, "Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0, rn1, sqrt(*r2max));
119 for (a1 = 0; a1 < depth; a1++)
121 fprintf(debug, " %d %5.3f",
123 iparams[constr_iatomptr(ia1, ia2, con)[0]].constr.dA);
125 fprintf(debug, " %d %5.3f\n", con, len);
128 /* Limit the number of recursions to 1000*nc,
129 * so a call does not take more than a second,
130 * even for highly connected systems.
132 if (depth + 1 < nc && *count < 1000*nc)
144 constr_recur(at2con, ilist, iparams,
145 bTopB, a1, depth+1, nc, path, rn0, rn1, r2max, count);
152 //! Find the interaction radius needed for constraints for this molecule type.
153 static real constr_r_max_moltype(const gmx_moltype_t *molt,
154 const t_iparams *iparams,
155 const t_inputrec *ir)
157 int natoms, *path, at, count;
160 real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
162 if (molt->ilist[F_CONSTR].size() == 0 &&
163 molt->ilist[F_CONSTRNC].size() == 0)
168 natoms = molt->atoms.nr;
170 at2con = make_at2con(*molt, iparams,
171 flexibleConstraintTreatment(EI_DYNAMICS(ir->eI)));
172 snew(path, 1+ir->nProjOrder);
173 for (at = 0; at < 1+ir->nProjOrder; at++)
179 for (at = 0; at < natoms; at++)
185 constr_recur(&at2con, molt->ilist, iparams,
186 FALSE, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxA, &count);
188 if (ir->efep == efepNO)
195 for (at = 0; at < natoms; at++)
200 constr_recur(&at2con, molt->ilist, iparams,
201 TRUE, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxB, &count);
203 lam0 = ir->fepvals->init_lambda;
204 if (EI_DYNAMICS(ir->eI))
206 lam0 += ir->init_step*ir->fepvals->delta_lambda;
208 rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
209 if (EI_DYNAMICS(ir->eI))
211 lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
212 rmax = std::max(rmax, (1 - lam1)*std::sqrt(r2maxA) + lam1*std::sqrt(r2maxB));
216 done_blocka(&at2con);
222 real constr_r_max(const MDLogger &mdlog, const gmx_mtop_t *mtop, const t_inputrec *ir)
225 for (const gmx_moltype_t &molt : mtop->moltype)
227 rmax = std::max(rmax,
228 constr_r_max_moltype(&molt,
229 mtop->ffparams.iparams, ir));
232 GMX_LOG(mdlog.info).appendTextFormatted(
233 "Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm",
234 1+ir->nProjOrder, rmax);