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, 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.
43 #include "gromacs/legacyheaders/types/ifunc.h"
44 #include "gromacs/topology/topology.h"
45 #include "gromacs/utility/fatalerror.h"
46 #include "gromacs/utility/smalloc.h"
48 static gmx_bool
ip_pert(int ftype
, const t_iparams
*ip
)
53 if (NRFPB(ftype
) == 0)
66 bPert
= (ip
->harmonic
.rA
!= ip
->harmonic
.rB
||
67 ip
->harmonic
.krA
!= ip
->harmonic
.krB
);
70 bPert
= (ip
->morse
.b0A
!= ip
->morse
.b0B
||
71 ip
->morse
.cbA
!= ip
->morse
.cbB
||
72 ip
->morse
.betaA
!= ip
->morse
.betaB
);
75 bPert
= (ip
->restraint
.lowA
!= ip
->restraint
.lowB
||
76 ip
->restraint
.up1A
!= ip
->restraint
.up1B
||
77 ip
->restraint
.up2A
!= ip
->restraint
.up2B
||
78 ip
->restraint
.kA
!= ip
->restraint
.kB
);
81 bPert
= (ip
->u_b
.thetaA
!= ip
->u_b
.thetaB
||
82 ip
->u_b
.kthetaA
!= ip
->u_b
.kthetaB
||
83 ip
->u_b
.r13A
!= ip
->u_b
.r13B
||
84 ip
->u_b
.kUBA
!= ip
->u_b
.kUBB
);
90 bPert
= (ip
->pdihs
.phiA
!= ip
->pdihs
.phiB
||
91 ip
->pdihs
.cpA
!= ip
->pdihs
.cpB
);
95 for (i
= 0; i
< NR_RBDIHS
; i
++)
97 if (ip
->rbdihs
.rbcA
[i
] != ip
->rbdihs
.rbcB
[i
])
107 bPert
= (ip
->tab
.kA
!= ip
->tab
.kB
);
111 for (i
= 0; i
< DIM
; i
++)
113 if (ip
->posres
.pos0A
[i
] != ip
->posres
.pos0B
[i
] ||
114 ip
->posres
.fcA
[i
] != ip
->posres
.fcB
[i
])
121 bPert
= ((ip
->dihres
.phiA
!= ip
->dihres
.phiB
) ||
122 (ip
->dihres
.dphiA
!= ip
->dihres
.dphiB
) ||
123 (ip
->dihres
.kfacA
!= ip
->dihres
.kfacB
));
126 bPert
= (ip
->lj14
.c6A
!= ip
->lj14
.c6B
||
127 ip
->lj14
.c12A
!= ip
->lj14
.c12B
);
136 gmx_fatal(FARGS
, "Function type %s does not support currentely free energy calculations",
137 interaction_function
[ftype
].longname
);
140 gmx_fatal(FARGS
, "Function type %s not implemented in ip_pert",
141 interaction_function
[ftype
].longname
);
147 static gmx_bool
ip_q_pert(int ftype
, const t_iatom
*ia
,
148 const t_iparams
*ip
, const real
*qA
, const real
*qB
)
150 /* 1-4 interactions do not have the charges stored in the iparams list,
151 * so we need a separate check for those.
153 return (ip_pert(ftype
, ip
+ia
[0]) ||
154 (ftype
== F_LJ14
&& (qA
[ia
[1]] != qB
[ia
[1]] ||
155 qA
[ia
[2]] != qB
[ia
[2]])));
158 gmx_bool
gmx_mtop_bondeds_free_energy(const gmx_mtop_t
*mtop
)
160 const gmx_ffparams_t
*ffparams
;
168 ffparams
= &mtop
->ffparams
;
170 /* Loop over all the function types and compare the A/B parameters */
172 for (i
= 0; i
< ffparams
->ntypes
; i
++)
174 ftype
= ffparams
->functype
[i
];
175 if (interaction_function
[ftype
].flags
& IF_BOND
)
177 if (ip_pert(ftype
, &ffparams
->iparams
[i
]))
184 /* Check perturbed charges for 1-4 interactions */
185 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
187 atom
= mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
.atom
;
188 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[F_LJ14
];
190 for (i
= 0; i
< il
->nr
; i
+= 3)
192 if (atom
[ia
[i
+1]].q
!= atom
[ia
[i
+1]].qB
||
193 atom
[ia
[i
+2]].q
!= atom
[ia
[i
+2]].qB
)
203 void gmx_sort_ilist_fe(t_idef
*idef
, const real
*qA
, const real
*qB
)
205 int ftype
, nral
, i
, ic
, ib
, a
;
221 iparams
= idef
->iparams
;
223 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
225 if (interaction_function
[ftype
].flags
& IF_BOND
)
227 ilist
= &idef
->il
[ftype
];
228 iatoms
= ilist
->iatoms
;
233 while (i
< ilist
->nr
)
235 /* Check if this interaction is perturbed */
236 if (ip_q_pert(ftype
, iatoms
+i
, iparams
, qA
, qB
))
238 /* Copy to the perturbed buffer */
239 if (ib
+ 1 + nral
> iabuf_nalloc
)
241 iabuf_nalloc
= over_alloc_large(ib
+1+nral
);
242 srenew(iabuf
, iabuf_nalloc
);
244 for (a
= 0; a
< 1+nral
; a
++)
246 iabuf
[ib
++] = iatoms
[i
++];
252 for (a
= 0; a
< 1+nral
; a
++)
254 iatoms
[ic
++] = iatoms
[i
++];
258 /* Now we now the number of non-perturbed interactions */
259 ilist
->nr_nonperturbed
= ic
;
261 /* Copy the buffer with perturbed interactions to the ilist */
262 for (a
= 0; a
< ib
; a
++)
264 iatoms
[ic
++] = iabuf
[a
];
269 fprintf(debug
, "%s non-pert %d pert %d\n",
270 interaction_function
[ftype
].longname
,
271 ilist
->nr_nonperturbed
,
272 ilist
->nr
-ilist
->nr_nonperturbed
);
279 idef
->ilsort
= ilsortFE_SORTED
;