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, 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/topology/ifunc.h"
45 #include "gromacs/topology/topology.h"
46 #include "gromacs/utility/arrayref.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/utility/smalloc.h"
50 static gmx_bool
ip_pert(int ftype
, const t_iparams
* ip
)
55 if (NRFPB(ftype
) == 0)
68 bPert
= (ip
->harmonic
.rA
!= ip
->harmonic
.rB
|| ip
->harmonic
.krA
!= ip
->harmonic
.krB
);
71 bPert
= (ip
->morse
.b0A
!= ip
->morse
.b0B
|| ip
->morse
.cbA
!= ip
->morse
.cbB
72 || ip
->morse
.betaA
!= ip
->morse
.betaB
);
75 bPert
= (ip
->restraint
.lowA
!= ip
->restraint
.lowB
|| ip
->restraint
.up1A
!= ip
->restraint
.up1B
76 || ip
->restraint
.up2A
!= ip
->restraint
.up2B
77 || ip
->restraint
.kA
!= ip
->restraint
.kB
);
80 bPert
= (ip
->u_b
.thetaA
!= ip
->u_b
.thetaB
|| ip
->u_b
.kthetaA
!= ip
->u_b
.kthetaB
81 || ip
->u_b
.r13A
!= ip
->u_b
.r13B
|| ip
->u_b
.kUBA
!= ip
->u_b
.kUBB
);
87 bPert
= (ip
->pdihs
.phiA
!= ip
->pdihs
.phiB
|| ip
->pdihs
.cpA
!= ip
->pdihs
.cpB
);
91 for (i
= 0; i
< NR_RBDIHS
; i
++)
93 if (ip
->rbdihs
.rbcA
[i
] != ip
->rbdihs
.rbcB
[i
])
102 case F_TABDIHS
: bPert
= (ip
->tab
.kA
!= ip
->tab
.kB
); break;
105 for (i
= 0; i
< DIM
; i
++)
107 if (ip
->posres
.pos0A
[i
] != ip
->posres
.pos0B
[i
] || ip
->posres
.fcA
[i
] != ip
->posres
.fcB
[i
])
114 bPert
= ((ip
->dihres
.phiA
!= ip
->dihres
.phiB
) || (ip
->dihres
.dphiA
!= ip
->dihres
.dphiB
)
115 || (ip
->dihres
.kfacA
!= ip
->dihres
.kfacB
));
118 bPert
= (ip
->lj14
.c6A
!= ip
->lj14
.c6B
|| ip
->lj14
.c12A
!= ip
->lj14
.c12B
);
120 case F_CMAP
: bPert
= FALSE
; break;
124 gmx_fatal(FARGS
, "Function type %s does not support currentely free energy calculations",
125 interaction_function
[ftype
].longname
);
127 gmx_fatal(FARGS
, "Function type %s not implemented in ip_pert",
128 interaction_function
[ftype
].longname
);
134 static gmx_bool
ip_q_pert(int ftype
, const t_iatom
* ia
, const t_iparams
* ip
, const real
* qA
, const real
* qB
)
136 /* 1-4 interactions do not have the charges stored in the iparams list,
137 * so we need a separate check for those.
139 return (ip_pert(ftype
, ip
+ ia
[0])
140 || (ftype
== F_LJ14
&& (qA
[ia
[1]] != qB
[ia
[1]] || qA
[ia
[2]] != qB
[ia
[2]])));
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(t_idef
* idef
, const real
* qA
, const real
* qB
)
181 int ftype
, nral
, i
, ic
, ib
, a
;
195 const t_iparams
* iparams
= idef
->iparams
;
197 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
199 if (interaction_function
[ftype
].flags
& IF_BOND
)
201 ilist
= &idef
->il
[ftype
];
202 iatoms
= ilist
->iatoms
;
207 while (i
< ilist
->nr
)
209 /* Check if this interaction is perturbed */
210 if (ip_q_pert(ftype
, iatoms
+ i
, iparams
, qA
, qB
))
212 /* Copy to the perturbed buffer */
213 if (ib
+ 1 + nral
> iabuf_nalloc
)
215 iabuf_nalloc
= over_alloc_large(ib
+ 1 + nral
);
216 srenew(iabuf
, iabuf_nalloc
);
218 for (a
= 0; a
< 1 + nral
; a
++)
220 iabuf
[ib
++] = iatoms
[i
++];
226 for (a
= 0; a
< 1 + nral
; a
++)
228 iatoms
[ic
++] = iatoms
[i
++];
232 /* Now we now the number of non-perturbed interactions */
233 idef
->numNonperturbedInteractions
[ftype
] = ic
;
235 /* Copy the buffer with perturbed interactions to the ilist */
236 for (a
= 0; a
< ib
; a
++)
238 iatoms
[ic
++] = iabuf
[a
];
243 const int numNonperturbed
= idef
->numNonperturbedInteractions
[ftype
];
244 fprintf(debug
, "%s non-pert %d pert %d\n", interaction_function
[ftype
].longname
,
245 numNonperturbed
, ilist
->nr
- numNonperturbed
);
252 idef
->ilsort
= ilsortFE_SORTED
;