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,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.
40 * \brief This file contains function definitions necessary for
41 * computing energies and forces for the plain-Ewald long-ranged part,
42 * and the correction for overall system charge for all Ewald-family
45 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
46 * \author Mark Abraham <mark.j.abraham@gmail.com>
47 * \ingroup module_ewald
59 #include "gromacs/ewald/ewald_utils.h"
60 #include "gromacs/math/functions.h"
61 #include "gromacs/math/gmxcomplex.h"
62 #include "gromacs/math/units.h"
63 #include "gromacs/math/utilities.h"
64 #include "gromacs/math/vec.h"
65 #include "gromacs/math/vectypes.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/forcerec.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/interaction_const.h"
70 #include "gromacs/mdtypes/md_enums.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/smalloc.h"
74 struct gmx_ewald_tab_t
78 t_complex
*tab_xy
, *tab_qxyz
;
81 void init_ewald_tab(struct gmx_ewald_tab_t
** et
, const t_inputrec
* ir
, FILE* fp
)
86 fprintf(fp
, "Will do ordinary reciprocal space Ewald sum.\n");
89 (*et
)->nx
= ir
->nkx
+ 1;
90 (*et
)->ny
= ir
->nky
+ 1;
91 (*et
)->nz
= ir
->nkz
+ 1;
92 (*et
)->kmax
= std::max((*et
)->nx
, std::max((*et
)->ny
, (*et
)->nz
));
94 (*et
)->tab_xy
= nullptr;
95 (*et
)->tab_qxyz
= nullptr;
98 //! Calculates wave vectors.
99 static void calc_lll(const rvec box
, rvec lll
)
101 lll
[XX
] = 2.0 * M_PI
/ box
[XX
];
102 lll
[YY
] = 2.0 * M_PI
/ box
[YY
];
103 lll
[ZZ
] = 2.0 * M_PI
/ box
[ZZ
];
106 //! Make tables for the structure factor parts
107 static void tabulateStructureFactors(int natom
, const rvec x
[], int kmax
, cvec
** eir
, const rvec lll
)
113 printf("Go away! kmax = %d\n", kmax
);
117 for (i
= 0; (i
< natom
); i
++)
119 for (m
= 0; (m
< 3); m
++)
125 for (m
= 0; (m
< 3); m
++)
127 eir
[1][i
][m
].re
= std::cos(x
[i
][m
] * lll
[m
]);
128 eir
[1][i
][m
].im
= std::sin(x
[i
][m
] * lll
[m
]);
130 for (j
= 2; (j
< kmax
); j
++)
132 for (m
= 0; (m
< 3); m
++)
134 eir
[j
][i
][m
] = cmul(eir
[j
- 1][i
][m
], eir
[1][i
][m
]);
140 real
do_ewald(const t_inputrec
* ir
,
143 const real chargeA
[],
144 const real chargeB
[],
154 real factor
= -1.0 / (4 * ewaldcoeff
* ewaldcoeff
);
156 real energy_AB
[2], energy
;
158 int lowiy
, lowiz
, ix
, iy
, iz
, n
, q
;
159 real tmp
, cs
, ss
, ak
, akv
, mx
, my
, mz
, m2
, scale
;
160 gmx_bool bFreeEnergy
;
166 gmx_fatal(FARGS
, "No parallel Ewald. Use PME instead.\n");
170 /* Scale box with Ewald wall factor */
172 EwaldBoxZScaler
boxScaler(*ir
);
173 boxScaler
.scaleBox(box
, scaledBox
);
176 for (int i
= 0; (i
< DIM
); i
++)
178 boxDiag
[i
] = scaledBox
[i
][i
];
182 real scaleRecip
= 4.0 * M_PI
/ (boxDiag
[XX
] * boxDiag
[YY
] * boxDiag
[ZZ
]) * ONE_4PI_EPS0
/ ir
->epsilon_r
;
184 if (!et
->eir
) /* allocate if we need to */
186 snew(et
->eir
, et
->kmax
);
187 for (n
= 0; n
< et
->kmax
; n
++)
189 snew(et
->eir
[n
], natoms
);
191 snew(et
->tab_xy
, natoms
);
192 snew(et
->tab_qxyz
, natoms
);
195 bFreeEnergy
= (ir
->efep
!= efepNO
);
199 calc_lll(boxDiag
, lll
);
200 tabulateStructureFactors(natoms
, x
, et
->kmax
, et
->eir
, lll
);
202 for (q
= 0; q
< (bFreeEnergy
? 2 : 1); q
++)
212 scale
= 1.0 - lambda
;
222 for (ix
= 0; ix
< et
->nx
; ix
++)
225 for (iy
= lowiy
; iy
< et
->ny
; iy
++)
230 for (n
= 0; n
< natoms
; n
++)
232 et
->tab_xy
[n
] = cmul(et
->eir
[ix
][n
][XX
], et
->eir
[iy
][n
][YY
]);
237 for (n
= 0; n
< natoms
; n
++)
239 et
->tab_xy
[n
] = cmul(et
->eir
[ix
][n
][XX
], conjugate(et
->eir
[-iy
][n
][YY
]));
242 for (iz
= lowiz
; iz
< et
->nz
; iz
++)
245 m2
= mx
* mx
+ my
* my
+ mz
* mz
;
246 ak
= std::exp(m2
* factor
) / m2
;
247 akv
= 2.0 * ak
* (1.0 / m2
- factor
);
250 for (n
= 0; n
< natoms
; n
++)
252 et
->tab_qxyz
[n
] = rcmul(charge
[n
], cmul(et
->tab_xy
[n
], et
->eir
[iz
][n
][ZZ
]));
257 for (n
= 0; n
< natoms
; n
++)
259 et
->tab_qxyz
[n
] = rcmul(
260 charge
[n
], cmul(et
->tab_xy
[n
], conjugate(et
->eir
[-iz
][n
][ZZ
])));
265 for (n
= 0; n
< natoms
; n
++)
267 cs
+= et
->tab_qxyz
[n
].re
;
268 ss
+= et
->tab_qxyz
[n
].im
;
270 energy_AB
[q
] += ak
* (cs
* cs
+ ss
* ss
);
271 tmp
= scale
* akv
* (cs
* cs
+ ss
* ss
);
272 lrvir
[XX
][XX
] -= tmp
* mx
* mx
;
273 lrvir
[XX
][YY
] -= tmp
* mx
* my
;
274 lrvir
[XX
][ZZ
] -= tmp
* mx
* mz
;
275 lrvir
[YY
][YY
] -= tmp
* my
* my
;
276 lrvir
[YY
][ZZ
] -= tmp
* my
* mz
;
277 lrvir
[ZZ
][ZZ
] -= tmp
* mz
* mz
;
278 for (n
= 0; n
< natoms
; n
++)
280 /*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/
281 tmp
= scale
* ak
* (cs
* et
->tab_qxyz
[n
].im
- ss
* et
->tab_qxyz
[n
].re
);
282 f
[n
][XX
] += tmp
* mx
* 2 * scaleRecip
;
283 f
[n
][YY
] += tmp
* my
* 2 * scaleRecip
;
284 f
[n
][ZZ
] += tmp
* mz
* 2 * scaleRecip
;
300 energy
= energy_AB
[0];
304 energy
= (1.0 - lambda
) * energy_AB
[0] + lambda
* energy_AB
[1];
305 *dvdlambda
+= scaleRecip
* (energy_AB
[1] - energy_AB
[0]);
308 lrvir
[XX
][XX
] = -0.5 * scaleRecip
* (lrvir
[XX
][XX
] + energy
);
309 lrvir
[XX
][YY
] = -0.5 * scaleRecip
* (lrvir
[XX
][YY
]);
310 lrvir
[XX
][ZZ
] = -0.5 * scaleRecip
* (lrvir
[XX
][ZZ
]);
311 lrvir
[YY
][YY
] = -0.5 * scaleRecip
* (lrvir
[YY
][YY
] + energy
);
312 lrvir
[YY
][ZZ
] = -0.5 * scaleRecip
* (lrvir
[YY
][ZZ
]);
313 lrvir
[ZZ
][ZZ
] = -0.5 * scaleRecip
* (lrvir
[ZZ
][ZZ
] + energy
);
315 lrvir
[YY
][XX
] = lrvir
[XX
][YY
];
316 lrvir
[ZZ
][XX
] = lrvir
[XX
][ZZ
];
317 lrvir
[ZZ
][YY
] = lrvir
[YY
][ZZ
];
319 energy
*= scaleRecip
;
324 real
ewald_charge_correction(const t_commrec
* cr
,
325 const t_forcerec
* fr
,
332 real vol
, fac
, qs2A
, qs2B
, vc
, enercorr
;
337 /* Apply charge correction */
338 vol
= box
[XX
][XX
] * box
[YY
][YY
] * box
[ZZ
][ZZ
];
340 fac
= M_PI
* ONE_4PI_EPS0
341 / (fr
->ic
->epsilon_r
* 2.0 * vol
* vol
* gmx::square(fr
->ic
->ewaldcoeff_q
));
343 qs2A
= fr
->qsum
[0] * fr
->qsum
[0];
344 qs2B
= fr
->qsum
[1] * fr
->qsum
[1];
346 vc
= (qs2A
* (1 - lambda
) + qs2B
* lambda
) * fac
;
348 enercorr
= -vol
* vc
;
350 *dvdlambda
+= -vol
* (qs2B
- qs2A
) * fac
;
352 for (d
= 0; d
< DIM
; d
++)
359 fprintf(debug
, "Total charge correction: Vcharge=%g\n", enercorr
);