-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathewald.h
164 lines (138 loc) · 5.68 KB
/
ewald.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
// This file is part of the ESPResSo distribution (http://www.espresso.mpg.de).
// It is therefore subject to the ESPResSo license agreement which you accepted upon receiving the distribution
// and by which you are legally bound while utilizing this file in any form or way.
// There is NO WARRANTY, not even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// You should have received a copy of that license along with this program;
// if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
// write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
// Copyright (c) 2002-2004; all rights reserved unless otherwise stated.
#ifndef EWALD_H
#define EWALD_H
/** \file ewald.h Ewald algorithm for long range coulomb interaction.
*
* Implementation of the standard Ewald-Summation
*
* Further reading:
* <ul>
* <li> P.P. Ewald,
* <i>Die Berechnung optischer und elektrostatischer Gitterpotentiale</i>,
* Ann. Phys. (64) 253-287, 1921
* <li> M. Deserno and C. Holm,
* <i>How to mesh up {E}wald sums. I. + II.</i>,
* J. Chem. Phys. (109) 7678, 1998; (109) 7694, 1998
* <li> M. Deserno, C. Holm and H. J. Limbach,
* <i>How to mesh up {E}wald sums. </i>,
* in Molecular Dynamics on Parallel Computers,
* Ed. R. Esser et al., World Scientific, Singapore, 2000
* <li>
* </ul>
*
* For more information about the ewald algorithm,
* see \ref ewald.c "ewald.c"
*/
#include "config.h"
#include "debug.h"
#include "interaction_data.h"
#ifdef ELECTROSTATICS
/** This value for ewald.epsilon indicates metallic boundary conditions. */
#define EWALD_EPSILON_METALLIC 0.0
/************************************************
* data types
************************************************/
/** Structure to hold Ewald parameters and some dependend variables. */
typedef struct {
/** Ewald splitting parameter (0<alpha<1), rescaled to alpha_L = alpha * box_l. */
double alpha_L;
/** Cutoff radius for real space electrostatics (>0), rescaled to r_cut_iL = r_cut * box_l_i. */
double r_cut_iL;
/** epsilon of the "surrounding dielectric". */
double epsilon;
/** unscaled \ref alpha_L for use with fast inline functions only */
double alpha;
/** unscaled \ref r_cut_iL for use with fast inline functions only */
double r_cut;
/** Maximum KVec, e.g. maximum frequency in k-space*/
int kmax;
/** squared \ref kmax */
int kmaxsq;
} ewald_struct;
/** \name Exported Variables */
/************************************************************/
/*@{*/
/** Ewald parameters. */
extern ewald_struct ewald;
/*@}*/
/** \name Exported Functions */
/************************************************************/
/*@{*/
/// print the ewald parameters to the interpreters result
int printEWALDToResult(Tcl_Interp *interp);
/// parse the ewald parameters
int inter_parse_ewald(Tcl_Interp * interp, int argc, char ** argv);
/// sanity checks
int EWALD_sanity_checks();
/** Initialize all structures, parameters and arrays needed for the
* Ewald algorithm.
*/
void EWALD_init();
/** Calculate number of charged particles, the sum of the squared
charges and the squared sum of the charges. */
void EWALD_count_charged_particles();
/** Reallocate memory for k-space caches */
void EWALD_on_resort_particles();
/** Updates \ref ewald_struct::alpha and \ref ewald_struct::r_cut if \ref box_l changed. */
void EWALD_scaleby_box_l();
/** Calculate the k-space contribution to the coulomb interaction forces. */
double EWALD_calc_kspace_forces(int force_flag, int energy_flag);
/** Calculate real space contribution of coulomb pair forces.
If NPT is compiled in, it returns the energy, which is needed for NPT. */
MDINLINE double add_ewald_coulomb_pair_force(Particle *p1, Particle *p2,
double *d,double dist2,double dist,double force[3])
{
int j;
double fac1,fac2, adist, erfc_part_ri;
if(dist < ewald.r_cut) {
adist = ewald.alpha * dist;
#if USE_ERFC_APPROXIMATION
erfc_part_ri = AS_erfc_part(adist) / dist;
fac1 = coulomb.prefactor * p1->p.q * p2->p.q * exp(-adist*adist);
fac2 = fac1 * (erfc_part_ri + 2.0*ewald.alpha*wupii) / dist2;
#else
erfc_part_ri = erfc(adist) / dist;
fac1 = coulomb.prefactor * p1->p.q * p2->p.q;
fac2 = fac1 * (erfc_part_ri + 2.0*ewald.alpha*wupii*exp(-adist*adist)) / dist2;
#endif
for(j=0;j<3;j++)
force[j] += fac2 * d[j];
ESR_TRACE(fprintf(stderr,"%d: RSE: Pair (%d-%d) dist=%.3f: force (%.3e,%.3e,%.3e)\n",this_node,
p1->p.identity,p2->p.identity,dist,fac*d[0],fac*d[1],fac*d[2]));
ONEPART_TRACE(if(p1->p.identity==check_id) fprintf(stderr,"%d: OPT: ESR f = (%.3e,%.3e,%.3e) with part id=%d at dist %f fac %.3e\n",this_node,p1->f.f[0],p1->f.f[1],p1->f.f[2],p2->p.identity,dist,fac2));
ONEPART_TRACE(if(p2->p.identity==check_id) fprintf(stderr,"%d: OPT: ESR f = (%.3e,%.3e,%.3e) with part id=%d at dist %f fac %.3e\n",this_node,p2->f.f[0],p2->f.f[1],p2->f.f[2],p1->p.identity,dist,fac2));
#ifdef NPT
return fac1 * erfc_part_ri;
#endif
}
return 0.0;
}
/** Calculate real space contribution of coulomb pair energy. */
MDINLINE double ewald_coulomb_pair_energy(Particle *p1, Particle *p2,
double *d,double dist2,double dist)
{
double adist, erfc_part_ri;
if(dist < ewald.r_cut) {
adist = ewald.alpha * dist;
#if USE_ERFC_APPROXIMATION
erfc_part_ri = AS_erfc_part(adist) / dist;
return coulomb.prefactor*p1->p.q*p2->p.q *erfc_part_ri*exp(-adist*adist);
#else
erfc_part_ri = erfc(adist) / dist;
return coulomb.prefactor*p1->p.q*p2->p.q *erfc_part_ri;
#endif
}
return 0.0;
}
/** Clean up Ewald memory allocations. */
void Ewald_exit();
/*@}*/
#endif
#endif