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, 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.
39 * \brief This file contains function definitions necessary for
40 * computing energies and forces for the plain-Ewald long-ranged part,
41 * and the correction for overall system charge for all Ewald-family
44 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
45 * \author Mark Abraham <mark.j.abraham@gmail.com>
46 * \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/md_enums.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/smalloc.h"
73 struct gmx_ewald_tab_t
77 t_complex
*tab_xy
, *tab_qxyz
;
80 void init_ewald_tab(struct gmx_ewald_tab_t
**et
, const t_inputrec
*ir
, FILE *fp
)
85 fprintf(fp
, "Will do ordinary reciprocal space Ewald sum.\n");
88 (*et
)->nx
= ir
->nkx
+1;
89 (*et
)->ny
= ir
->nky
+1;
90 (*et
)->nz
= ir
->nkz
+1;
91 (*et
)->kmax
= std::max((*et
)->nx
, std::max((*et
)->ny
, (*et
)->nz
));
93 (*et
)->tab_xy
= nullptr;
94 (*et
)->tab_qxyz
= nullptr;
97 //! Calculates wave vectors.
98 static void calc_lll(const rvec box
, rvec lll
)
100 lll
[XX
] = 2.0*M_PI
/box
[XX
];
101 lll
[YY
] = 2.0*M_PI
/box
[YY
];
102 lll
[ZZ
] = 2.0*M_PI
/box
[ZZ
];
105 //! Make tables for the structure factor parts
106 static void tabulateStructureFactors(int natom
, rvec x
[], int kmax
, cvec
**eir
, rvec lll
)
112 printf("Go away! kmax = %d\n", kmax
);
116 for (i
= 0; (i
< natom
); i
++)
118 for (m
= 0; (m
< 3); m
++)
124 for (m
= 0; (m
< 3); m
++)
126 eir
[1][i
][m
].re
= cos(x
[i
][m
]*lll
[m
]);
127 eir
[1][i
][m
].im
= sin(x
[i
][m
]*lll
[m
]);
129 for (j
= 2; (j
< kmax
); j
++)
131 for (m
= 0; (m
< 3); m
++)
133 eir
[j
][i
][m
] = cmul(eir
[j
-1][i
][m
], eir
[1][i
][m
]);
139 real
do_ewald(t_inputrec
*ir
,
141 real chargeA
[], real chargeB
[],
143 t_commrec
*cr
, int natoms
,
144 matrix lrvir
, real ewaldcoeff
,
145 real lambda
, real
*dvdlambda
,
146 struct gmx_ewald_tab_t
*et
)
148 real factor
= -1.0/(4*ewaldcoeff
*ewaldcoeff
);
149 real
*charge
, energy_AB
[2], energy
;
151 int lowiy
, lowiz
, ix
, iy
, iz
, n
, q
;
152 real tmp
, cs
, ss
, ak
, akv
, mx
, my
, mz
, m2
, scale
;
153 gmx_bool bFreeEnergy
;
159 gmx_fatal(FARGS
, "No parallel Ewald. Use PME instead.\n");
163 /* Scale box with Ewald wall factor */
165 EwaldBoxZScaler
boxScaler(*ir
);
166 boxScaler
.scaleBox(box
, scaledBox
);
169 for (int i
= 0; (i
< DIM
); i
++)
171 boxDiag
[i
] = scaledBox
[i
][i
];
175 real scaleRecip
= 4.0*M_PI
/(boxDiag
[XX
]*boxDiag
[YY
]*boxDiag
[ZZ
])*ONE_4PI_EPS0
/ir
->epsilon_r
;
177 if (!et
->eir
) /* allocate if we need to */
179 snew(et
->eir
, et
->kmax
);
180 for (n
= 0; n
< et
->kmax
; n
++)
182 snew(et
->eir
[n
], natoms
);
184 snew(et
->tab_xy
, natoms
);
185 snew(et
->tab_qxyz
, natoms
);
188 bFreeEnergy
= (ir
->efep
!= efepNO
);
192 calc_lll(boxDiag
, lll
);
193 tabulateStructureFactors(natoms
, x
, et
->kmax
, et
->eir
, lll
);
195 for (q
= 0; q
< (bFreeEnergy
? 2 : 1); q
++)
205 scale
= 1.0 - lambda
;
215 for (ix
= 0; ix
< et
->nx
; ix
++)
218 for (iy
= lowiy
; iy
< et
->ny
; iy
++)
223 for (n
= 0; n
< natoms
; n
++)
225 et
->tab_xy
[n
] = cmul(et
->eir
[ix
][n
][XX
], et
->eir
[iy
][n
][YY
]);
230 for (n
= 0; n
< natoms
; n
++)
232 et
->tab_xy
[n
] = cmul(et
->eir
[ix
][n
][XX
], conjugate(et
->eir
[-iy
][n
][YY
]));
235 for (iz
= lowiz
; iz
< et
->nz
; iz
++)
238 m2
= mx
*mx
+my
*my
+mz
*mz
;
239 ak
= exp(m2
*factor
)/m2
;
240 akv
= 2.0*ak
*(1.0/m2
-factor
);
243 for (n
= 0; n
< natoms
; n
++)
245 et
->tab_qxyz
[n
] = rcmul(charge
[n
], cmul(et
->tab_xy
[n
],
246 et
->eir
[iz
][n
][ZZ
]));
251 for (n
= 0; n
< natoms
; n
++)
253 et
->tab_qxyz
[n
] = rcmul(charge
[n
], cmul(et
->tab_xy
[n
],
254 conjugate(et
->eir
[-iz
][n
][ZZ
])));
259 for (n
= 0; n
< natoms
; n
++)
261 cs
+= et
->tab_qxyz
[n
].re
;
262 ss
+= et
->tab_qxyz
[n
].im
;
264 energy_AB
[q
] += ak
*(cs
*cs
+ss
*ss
);
265 tmp
= scale
*akv
*(cs
*cs
+ss
*ss
);
266 lrvir
[XX
][XX
] -= tmp
*mx
*mx
;
267 lrvir
[XX
][YY
] -= tmp
*mx
*my
;
268 lrvir
[XX
][ZZ
] -= tmp
*mx
*mz
;
269 lrvir
[YY
][YY
] -= tmp
*my
*my
;
270 lrvir
[YY
][ZZ
] -= tmp
*my
*mz
;
271 lrvir
[ZZ
][ZZ
] -= tmp
*mz
*mz
;
272 for (n
= 0; n
< natoms
; n
++)
274 /*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/
275 tmp
= scale
*ak
*(cs
*et
->tab_qxyz
[n
].im
-ss
*et
->tab_qxyz
[n
].re
);
276 f
[n
][XX
] += tmp
*mx
*2*scaleRecip
;
277 f
[n
][YY
] += tmp
*my
*2*scaleRecip
;
278 f
[n
][ZZ
] += tmp
*mz
*2*scaleRecip
;
294 energy
= energy_AB
[0];
298 energy
= (1.0 - lambda
)*energy_AB
[0] + lambda
*energy_AB
[1];
299 *dvdlambda
+= scaleRecip
*(energy_AB
[1] - energy_AB
[0]);
302 lrvir
[XX
][XX
] = -0.5*scaleRecip
*(lrvir
[XX
][XX
]+energy
);
303 lrvir
[XX
][YY
] = -0.5*scaleRecip
*(lrvir
[XX
][YY
]);
304 lrvir
[XX
][ZZ
] = -0.5*scaleRecip
*(lrvir
[XX
][ZZ
]);
305 lrvir
[YY
][YY
] = -0.5*scaleRecip
*(lrvir
[YY
][YY
]+energy
);
306 lrvir
[YY
][ZZ
] = -0.5*scaleRecip
*(lrvir
[YY
][ZZ
]);
307 lrvir
[ZZ
][ZZ
] = -0.5*scaleRecip
*(lrvir
[ZZ
][ZZ
]+energy
);
309 lrvir
[YY
][XX
] = lrvir
[XX
][YY
];
310 lrvir
[ZZ
][XX
] = lrvir
[XX
][ZZ
];
311 lrvir
[ZZ
][YY
] = lrvir
[YY
][ZZ
];
313 energy
*= scaleRecip
;
318 real
ewald_charge_correction(t_commrec
*cr
, t_forcerec
*fr
, real lambda
,
320 real
*dvdlambda
, tensor vir
)
323 real vol
, fac
, qs2A
, qs2B
, vc
, enercorr
;
328 /* Apply charge correction */
329 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
331 fac
= M_PI
*ONE_4PI_EPS0
/(fr
->ic
->epsilon_r
*2.0*vol
*vol
*gmx::square(fr
->ic
->ewaldcoeff_q
));
333 qs2A
= fr
->qsum
[0]*fr
->qsum
[0];
334 qs2B
= fr
->qsum
[1]*fr
->qsum
[1];
336 vc
= (qs2A
*(1 - lambda
) + qs2B
*lambda
)*fac
;
340 *dvdlambda
+= -vol
*(qs2B
- qs2A
)*fac
;
342 for (d
= 0; d
< DIM
; d
++)
349 fprintf(debug
, "Total charge correction: Vcharge=%g\n", enercorr
);