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 by the GROMACS development team.
7 * Copyright (c) 2018,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.
40 #include "enerdata_utils.h"
42 #include "gromacs/mdtypes/enerdata.h"
43 #include "gromacs/mdtypes/inputrec.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/utility/smalloc.h"
47 gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups
, int numFepLambdas
) :
48 grpp(numEnergyGroups
),
49 enerpart_lambda(numFepLambdas
== 0 ? 0 : numFepLambdas
+ 1),
50 foreign_grpp(numEnergyGroups
)
54 static real
sum_v(int n
, gmx::ArrayRef
<const real
> v
)
60 for (i
= 0; (i
< n
); i
++)
68 void sum_epot(gmx_grppairener_t
* grpp
, real
* epot
)
72 /* Accumulate energies */
73 epot
[F_COUL_SR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULSR
]);
74 epot
[F_LJ
] = sum_v(grpp
->nener
, grpp
->ener
[egLJSR
]);
75 epot
[F_LJ14
] = sum_v(grpp
->nener
, grpp
->ener
[egLJ14
]);
76 epot
[F_COUL14
] = sum_v(grpp
->nener
, grpp
->ener
[egCOUL14
]);
78 /* lattice part of LR doesnt belong to any group
79 * and has been added earlier
81 epot
[F_BHAM
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMSR
]);
84 for (i
= 0; (i
< F_EPOT
); i
++)
86 if (i
!= F_DISRESVIOL
&& i
!= F_ORIRESDEV
)
88 epot
[F_EPOT
] += epot
[i
];
93 void sum_dhdl(gmx_enerdata_t
* enerd
, gmx::ArrayRef
<const real
> lambda
, const t_lambda
& fepvals
)
97 enerd
->dvdl_lin
[efptVDW
] += enerd
->term
[F_DVDL_VDW
]; /* include dispersion correction */
98 enerd
->term
[F_DVDL
] = 0.0;
99 for (int i
= 0; i
< efptNR
; i
++)
101 if (fepvals
.separate_dvdl
[i
])
103 /* could this be done more readably/compactly? */
106 case (efptMASS
): index
= F_DKDL
; break;
107 case (efptCOUL
): index
= F_DVDL_COUL
; break;
108 case (efptVDW
): index
= F_DVDL_VDW
; break;
109 case (efptBONDED
): index
= F_DVDL_BONDED
; break;
110 case (efptRESTRAINT
): index
= F_DVDL_RESTRAINT
; break;
111 default: index
= F_DVDL
; break;
113 enerd
->term
[index
] = enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
116 fprintf(debug
, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n", efpt_names
[i
], i
,
117 enerd
->term
[index
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
122 enerd
->term
[F_DVDL
] += enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
125 fprintf(debug
, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n", efpt_names
[0], i
,
126 enerd
->term
[F_DVDL
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
131 if (fepvals
.separate_dvdl
[efptBONDED
])
133 enerd
->term
[F_DVDL_BONDED
] += enerd
->term
[F_DVDL_CONSTR
];
137 enerd
->term
[F_DVDL
] += enerd
->term
[F_DVDL_CONSTR
];
140 for (int i
= 0; i
< fepvals
.n_lambda
; i
++)
142 /* note we are iterating over fepvals here!
143 For the current lam, dlam = 0 automatically,
144 so we don't need to add anything to the
145 enerd->enerpart_lambda[0] */
147 /* we don't need to worry about dvdl_lin contributions to dE at
148 current lambda, because the contributions to the current
149 lambda are automatically zeroed */
151 double& enerpart_lambda
= enerd
->enerpart_lambda
[i
+ 1];
153 for (gmx::index j
= 0; j
< lambda
.ssize(); j
++)
155 /* Note that this loop is over all dhdl components, not just the separated ones */
156 const double dlam
= fepvals
.all_lambda
[j
][i
] - lambda
[j
];
158 enerpart_lambda
+= dlam
* enerd
->dvdl_lin
[j
];
160 /* Constraints can not be evaluated at foreign lambdas, so we add
161 * a linear extrapolation. This is an approximation, but usually
162 * quite accurate since constraints change little between lambdas.
164 if ((j
== efptBONDED
&& fepvals
.separate_dvdl
[efptBONDED
])
165 || (j
== efptFEP
&& !fepvals
.separate_dvdl
[efptBONDED
]))
167 enerpart_lambda
+= dlam
* enerd
->term
[F_DVDL_CONSTR
];
170 if (j
== efptMASS
&& !fepvals
.separate_dvdl
[j
])
172 enerpart_lambda
+= dlam
* enerd
->term
[F_DKDL
];
177 fprintf(debug
, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
178 fepvals
.all_lambda
[j
][i
], efpt_names
[j
],
179 enerpart_lambda
- enerd
->enerpart_lambda
[0], dlam
, enerd
->dvdl_lin
[j
]);
184 /* The constrain contribution is now included in other terms, so clear it */
185 enerd
->term
[F_DVDL_CONSTR
] = 0;
189 void reset_foreign_enerdata(gmx_enerdata_t
* enerd
)
193 /* First reset all foreign energy components. Foreign energies always called on
194 neighbor search steps */
195 for (i
= 0; (i
< egNR
); i
++)
197 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
199 enerd
->foreign_grpp
.ener
[i
][j
] = 0.0;
203 /* potential energy components */
204 for (i
= 0; (i
<= F_EPOT
); i
++)
206 enerd
->foreign_term
[i
] = 0.0;
210 void reset_enerdata(gmx_enerdata_t
* enerd
)
214 /* First reset all energy components. */
215 for (i
= 0; (i
< egNR
); i
++)
217 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
219 enerd
->grpp
.ener
[i
][j
] = 0.0;
222 for (i
= 0; i
< efptNR
; i
++)
224 enerd
->dvdl_lin
[i
] = 0.0;
225 enerd
->dvdl_nonlin
[i
] = 0.0;
228 /* Normal potential energy components */
229 for (i
= 0; (i
<= F_EPOT
); i
++)
231 enerd
->term
[i
] = 0.0;
233 enerd
->term
[F_DVDL
] = 0.0;
234 enerd
->term
[F_DVDL_COUL
] = 0.0;
235 enerd
->term
[F_DVDL_VDW
] = 0.0;
236 enerd
->term
[F_DVDL_BONDED
] = 0.0;
237 enerd
->term
[F_DVDL_RESTRAINT
] = 0.0;
238 enerd
->term
[F_DKDL
] = 0.0;
239 std::fill(enerd
->enerpart_lambda
.begin(), enerd
->enerpart_lambda
.end(), 0);
240 /* reset foreign energy data - separate function since we also call it elsewhere */
241 reset_foreign_enerdata(enerd
);