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-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020,2021, 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.
44 #include "gromacs/mdtypes/atominfo.h"
45 #include "gromacs/topology/ifunc.h"
46 #include "gromacs/topology/topology.h"
47 #include "gromacs/utility/arrayref.h"
48 #include "gromacs/utility/fatalerror.h"
49 #include "gromacs/utility/smalloc.h"
51 static gmx_bool ip_pert(int ftype, const t_iparams* ip)
53 if (NRFPB(ftype) == 0)
67 bPert = (ip->harmonic.rA != ip->harmonic.rB || ip->harmonic.krA != ip->harmonic.krB);
70 bPert = (ip->morse.b0A != ip->morse.b0B || ip->morse.cbA != ip->morse.cbB
71 || ip->morse.betaA != ip->morse.betaB);
74 bPert = (ip->restraint.lowA != ip->restraint.lowB || ip->restraint.up1A != ip->restraint.up1B
75 || ip->restraint.up2A != ip->restraint.up2B
76 || ip->restraint.kA != ip->restraint.kB);
79 bPert = (ip->u_b.thetaA != ip->u_b.thetaB || ip->u_b.kthetaA != ip->u_b.kthetaB
80 || ip->u_b.r13A != ip->u_b.r13B || ip->u_b.kUBA != ip->u_b.kUBB);
86 bPert = (ip->pdihs.phiA != ip->pdihs.phiB || ip->pdihs.cpA != ip->pdihs.cpB);
90 for (int i = 0; i < NR_RBDIHS; i++)
92 if (ip->rbdihs.rbcA[i] != ip->rbdihs.rbcB[i])
101 case F_TABDIHS: bPert = (ip->tab.kA != ip->tab.kB); break;
104 for (int i = 0; i < DIM; i++)
106 if (ip->posres.pos0A[i] != ip->posres.pos0B[i] || ip->posres.fcA[i] != ip->posres.fcB[i])
113 bPert = ((ip->dihres.phiA != ip->dihres.phiB) || (ip->dihres.dphiA != ip->dihres.dphiB)
114 || (ip->dihres.kfacA != ip->dihres.kfacB));
117 bPert = (ip->lj14.c6A != ip->lj14.c6B || ip->lj14.c12A != ip->lj14.c12B);
119 case F_CMAP: bPert = FALSE; break;
124 "Function type %s does not support currentely free energy calculations",
125 interaction_function[ftype].longname);
128 "Function type %s not implemented in ip_pert",
129 interaction_function[ftype].longname);
136 //! Return whether the two atom indices in a 1-4 interaction have perturbed charges per \c atomInfo
137 static bool hasPerturbedChargeIn14Interaction(int atom0, int atom1, gmx::ArrayRef<const int64_t> atomInfo)
139 return (bool(atomInfo[atom0] & gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction)
140 || bool(atomInfo[atom1] & gmx::sc_atomInfo_HasPerturbedChargeIn14Interaction));
143 gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t* mtop)
145 const gmx_ffparams_t* ffparams = &mtop->ffparams;
147 /* Loop over all the function types and compare the A/B parameters */
148 gmx_bool bPert = FALSE;
149 for (int i = 0; i < ffparams->numTypes(); i++)
151 int ftype = ffparams->functype[i];
152 if (interaction_function[ftype].flags & IF_BOND)
154 if (ip_pert(ftype, &ffparams->iparams[i]))
161 /* Check perturbed charges for 1-4 interactions */
162 for (const gmx_molblock_t& molb : mtop->molblock)
164 const t_atom* atom = mtop->moltype[molb.type].atoms.atom;
165 const InteractionList& il = mtop->moltype[molb.type].ilist[F_LJ14];
166 gmx::ArrayRef<const int> ia = il.iatoms;
167 for (int i = 0; i < il.size(); i += 3)
169 if (atom[ia[i + 1]].q != atom[ia[i + 1]].qB || atom[ia[i + 2]].q != atom[ia[i + 2]].qB)
179 void gmx_sort_ilist_fe(InteractionDefinitions* idef, gmx::ArrayRef<const int64_t> atomInfo)
181 bool havePerturbedInteractions = false;
183 int iabuf_nalloc = 0;
184 t_iatom* iabuf = nullptr;
186 for (int ftype = 0; ftype < F_NRE; ftype++)
188 if (interaction_function[ftype].flags & IF_BOND)
190 InteractionList* ilist = &idef->il[ftype];
191 int* iatoms = ilist->iatoms.data();
192 const int nral = NRAL(ftype);
196 while (i < ilist->size())
198 /* Check if this interaction is perturbed */
199 if (ip_pert(ftype, idef->iparams.data() + iatoms[i])
201 && hasPerturbedChargeIn14Interaction(iatoms[i + 1], iatoms[i + 2], atomInfo)))
203 /* Copy to the perturbed buffer */
204 if (ib + 1 + nral > iabuf_nalloc)
206 iabuf_nalloc = over_alloc_large(ib + 1 + nral);
207 srenew(iabuf, iabuf_nalloc);
209 for (int a = 0; a < 1 + nral; a++)
211 iabuf[ib++] = iatoms[i++];
214 havePerturbedInteractions = true;
219 for (int a = 0; a < 1 + nral; a++)
221 iatoms[ic++] = iatoms[i++];
225 /* Now we know the number of non-perturbed interactions */
226 idef->numNonperturbedInteractions[ftype] = ic;
228 /* Copy the buffer with perturbed interactions to the ilist */
229 for (int a = 0; a < ib; a++)
231 iatoms[ic++] = iabuf[a];
236 const int numNonperturbed = idef->numNonperturbedInteractions[ftype];
238 "%s non-pert %d pert %d\n",
239 interaction_function[ftype].longname,
241 ilist->size() - numNonperturbed);
248 idef->ilsort = (havePerturbedInteractions ? ilsortFE_SORTED : ilsortNO_FE);