The Coulomb interaction energy is the energy associated with the interaction of electrically charged particles. It is given by the equation:
\[U = \frac{1}{4π\epsilon_0} \sum_{i > j} \frac{q_i q_j}{|r_i-r_j|}\] where \(q_{i}\) and \(q_{j}\) are the charge of the two interacting particles, and \(\epsilon_0\) is the vacuum permittivity. The integration is performed over all space.

The problem with evaluating this expression arises from the fact that the interaction energy between two charged particles is infinite when the particles are at the same point in space. Also the sum is conditionally convergent for alternating charges. This divergence must be removed in order to obtain a meaningful result.

One approach to removing this divergence is the Ewald summation method. This method involves splitting the Coulomb interaction energy into a real space contribution, which is calculated using a finite cutoff distance, and a reciprocal space contribution, which is calculated using fast Fourier transforms. The real space contribution is responsible for the short-range interactions, while the reciprocal space contribution is responsible for the long-range interactions.

The Ewald summation method can provide accurate results for large systems, such as in molecular dynamics simulations of dusty plasmas or in the study of colloidal suspensions. The basic idea is to split the long-range interaction into two components, one that can be treated efficiently using Fourier transforms and one that can be treated using a short-range approximation. The Ewald sum is then the sum of these two components.

The method was first proposed by the German physicist Paul Peter Ewald in 1921, and it has become a standard tool for calculating the electrostatic interactions in many-particle systems. The technique has been extended to include other types of interactions, such as dipole-dipole interactions or magnetic interactions.

In molecular dynamics simulations, the Ewald summation method allows for efficient calculation of the electrostatic interactions between all of the atoms in a system, taking into account the periodic boundary conditions of the simulation cell. This makes it possible to study systems with thousands or even millions of atoms, which is important in fields such as biophysics or materials science..

Mathematical Formulation

The real part of the Ewald sum is \[U^{(r)} = \sum_{i=1}^N \sum_{j=1}^N \sum_{|\textbf n|=0}^{\infty}q_iq_j\frac{\text{erfc}(\alpha|\vec{r}_{ij}+\textbf nL|)}{|\vec{r}_{ij}+\textbf nL|}.\] The expression is for the electrostatic interaction energy, \(U^{(r)}\), between two charges \(q_i\) and \(q_j\) located at positions \(\vec{r}_i\) and \(\vec{r}_j\) in a periodic system with a periodic boundary condition of length L. The term \(\textbf n\) is a 3D vector of integers and the notation \(|\textbf n|\) indicates the magnitude or length of the vector. The function \(\text{erfc}(x)\) is the complementary error function, defined as \(\text{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_x^\infty e^{-t^2}dt\).

The k-part of the Ewald sum is
\[U^{(k)} = \frac{1}{2L^3} \sum_{\vec k\ne0} \sum_{i=1}^N \sum_{j=1}^N \frac{4 \pi q_i q_j}{k^2} \exp^{[i\vec k \cdot (\vec r_j - \vec r_i)]} \exp^{(-k^2 / 2\alpha)},\]
where, \(\vec{k} = 2\pi \textbf n/L.\)

The first sum is over all wave vectors \(\vec k \neq 0\), and the second and third sums are over all particle pairs.
In this expression, the exponential term \(\exp^{[i\vec k \cdot (\vec r_j - \vec r_i)]}\) represents the Fourier transform of the pairwise interaction between particles, and the exponential term \(\exp^{(-k^2 / 2\alpha)}\) is a Debye screening factor that accounts for the effect of the medium on the interaction between particles. The factor of \(1/2L^3\) is a normalization constant, where L is the length of the simulation box.
,
The self energy term is
\[U^{(self)} = -\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^Nq_i^2\] This expression is for the self-energy term of the charged particles, which is often included in the calculation of electrostatic potential energy to account for the interaction of a particle with its own electric field. The self-energy term is represented by \(U^{(self)}\). In this expression, the sum is over all N charged particles, and the term \(q_i^2\) represents the square of the charge of particle i. The factor of \(-\alpha/\sqrt{\pi}\) is a constant. This term accounts for the interaction of each particle with its own electric field, which is screened by the medium.

Now, the total force can be written as,
\[\vec F_i = \vec F_i^{(k)} +\vec F_i^{(r)}+\vec F_i^{(self)}.\]
We get the force as \(\vec F_i = -\nabla(U^{(total)})\).
\[\vec F_i^{(r)} = \sum_{i=1}^N \sum_{j=1}^N \sum_{|\textbf n|=0}^{\infty}q_iq_j \Big ( \frac{2\alpha}{\sqrt{\pi}} \exp^{(-\alpha^2|\vec{r}_{ij}+\textbf mL|^2)} +\frac{\text{erfc}(\alpha|\vec{r}_{ij}+\textbf nL|)}{|\vec{r}_{ij}+\textbf nL|}\Big ) . \frac{\vec{r}_{ij}+\textbf nL}{|\vec{r}_{ij}+\textbf nL|^2}.\] \[\vec F_i^{(k)} = q_i\frac{4\pi}{L^3} \sum_{j=0} q_j \sum_{\vec k\ne0} \frac{\vec k}{k^2} \exp^{(-k^2/{4\alpha^2})}\sin[\vec k \cdot \vec r_{ij}].\] \[\vec F_i^{(self)} = 0.\] For the efficient calculation, one can rewrite the Fourier part as \[\vec F_i^{(k)} = \frac{4\pi q_i}{2\pi L^3}\sum_{\vec{k}\ne 0}\frac{\vec{k} \exp^{(\frac{-k^2}{4\alpha^2})}}{k^2}\bigg(\sin(\vec k_i. \vec r_i) \text {Re}(S(\vec k)) + \cos(\vec k_i. \vec r_i) \text {Im}(S(\vec k))\bigg),\] where, \(S(\vec k) = \sum_{i=1}^N q_i\text{e}^{i\vec{k}_i.\vec{r}_i}.\)

C code Implementation

Below is the C function to calculate the real part energy and the forces.

double PME_Real(double **r_pos, double **prop, double **f_rewald, double *L, double *Ext_L){
    int i=0, j=0, k=0, nx=0, ny=0, nz=0;
    double rx=0.0, ry=0.0, rz=0.0;
    double alphafac = (2.0*alpha)/sqrt(M_PI);
    double alpha2 = alpha*alpha;
    double qi=0, qj=0, rsq=0.0, sqrt_rsq=0.0, erfc_sqrt_rsq=0.0, dcdr=0.0, e_rewald=0.0;
    double dij[dimen], Dij[dimen];

    // Initiate the force arrays with zero values //
    for (i=0;i<ParticleN;i++){
        for (j=0;j<dimen;j++){
            f_rewald[i][j] = 0.0; 
        }
    }
    
    // Initiate the distance calculation arrays //
    for (k=0;k<dimen;k++){
        dij[k] = 0.0; 
        Dij[k] = 0.0;
    }
    
    for (i=0;i<ParticleN-1;i++){
        qi = prop[i][0]; // charge of the particle with index i //
        for (j=i+1;j<ParticleN;j++){
            qj = prop[j][0]; // charge of the particle with index j //
            
            for (k=0;k<dimen;k++){
                dij[k] = r_pos[i][k]-r_pos[j][k]; // calculate the distance between particle i and j //
                //dij[k] = dij[k] - L[k]*round(dij[k]/L[k]); // Minimum Image Convension [0,L]
                dij[k] = dij[k] - 2*Ext_L[k]*round(0.5*dij[k]/Ext_L[k]); // Minimum Image Convension [-L,L] (boundary condition)
            }
            
            rx = dij[0]; // x component of the vector |r_i - r_j| //
            ry = dij[1]; // y component of the vector |r_i - r_j| //
            rz = dij[2]; // z component of the vector |r_i - r_j| //
            
            // Loop over the central and the image cells //
            for(nx = -PME_nmax; nx <= PME_nmax; nx++){
                for(ny = -PME_nmax; ny <= PME_nmax; ny++){
                    for(nz = -PME_nmax; nz <= PME_nmax; nz++){
                        
                        rsq = rx*rx + ry*ry + rz*rz;
                        sqrt_rsq = sqrt(rsq); // absolute distance //
                        
                        if(nx == 0 && ny == 0 && nz == 0){
                            // for the central cell: ignore i == j //
                            if(i != j){
                                erfc_sqrt_rsq = erfc(alpha*sqrt_rsq)/sqrt_rsq;
                                dcdr = ONE_PI_EPS0*qi*qj*(alphafac*exp(-alpha2*rsq) + erfc_sqrt_rsq)/rsq;

                                for (k=0;k<dimen;k++){
                                    f_rewald[i][k] = f_rewald[i][k] + dcdr*dij[k];
                                    f_rewald[j][k] = f_rewald[j][k] - dcdr*dij[k];
                                }
                                e_rewald += ONE_PI_EPS0*qi*qj*erfc_sqrt_rsq;
                            }
                            else{
                                // ignore, because of same particle //
                            }
                        }
                        
                        else{
                            // image cells (around the main cell) //
                            Dij[0] = rx + 2*nx*Ext_L[0]; // x component of the distance |r_i - r_j + n.L| //
                            Dij[1] = ry + 2*ny*Ext_L[1]; // y component of the distance |r_i - r_j + n.L| //
                            Dij[2] = rz + 2*nz*Ext_L[2]; // z component of the distance |r_i - r_j + n.L| //
                            
                            rsq = Dij[0]*Dij[0] + Dij[1]*Dij[1] + Dij[2]*Dij[2];
                            sqrt_rsq = sqrt(rsq); // absolute distance //
                            
                            erfc_sqrt_rsq = erfc(alpha*sqrt_rsq)/sqrt_rsq;
                            dcdr = ONE_PI_EPS0*qi*qj*(alphafac*exp(-alpha2*rsq) + erfc_sqrt_rsq)/rsq;
                            
                            for (k=0;k<dimen;k++){
                                f_rewald[i][k] = f_rewald[i][k] + dcdr*Dij[k]; // add the forces //
                                f_rewald[j][k] = f_rewald[j][k] - dcdr*Dij[k]; // Newton's Third Law: force are equal and opposite //
                            }
                            e_rewald += ONE_PI_EPS0*qi*qj*erfc_sqrt_rsq; // Ewald real part energy //
                        }
                    }
                }
            }								
        }	
    }
    
    printf("PME REwald: %lf\n", e_rewald);
    return e_rewald;
    
}

double PME_Fourier(double **r_pos, double **prop, double **f_kewald, double *L, double *Ext_L){
	int i=0, j=0, kx=0, ky=0,kz=0;
	double qi=0.0, ksq=0.0, kdotr=0.0, tmp=0.0, tmp2=0.0, e_kewald=0.0;
	double mx=0.0, my=0.0 ,mz=0.0, a=0.0 ,b=0.0, akak=0.0;
	double GAMMA = -0.25/(alpha*alpha);
	double recip = 2*M_PI*ONE_PI_EPS0/(8*Ext_L[0]*Ext_L[1]*Ext_L[2]);
	double ak[2];
	
	ak[0] = 0.0;
	ak[1] = 0.0;
	
	for (i=0;i<ParticleN;i++){
		for (j=0;j<dimen;j++){
			f_kewald[i][j] = 0.0; 
		}
	}
	
	for(kx = -PME_kmax; kx <= PME_kmax; kx++){
		mx = M_PI*kx/Ext_L[0];
		for(ky = -PME_kmax; ky <= PME_kmax; ky++){
			my = M_PI*ky/Ext_L[1];
			for(kz = -PME_kmax; kz <= PME_kmax; kz++){
				mz = M_PI*kz/Ext_L[2];
				ksq = mx*mx + my*my + mz*mz;
				
				if(ksq != 0){
					ak[0] = 0;
					ak[1] = 0;
					for(i=0;i<ParticleN;i++){
						qi = prop[i][0];
						kdotr = mx*r_pos[i][0] + my*r_pos[i][1] + mz*r_pos[i][2];
						
						ak[0] += qi*cos(kdotr);
						ak[1] -= qi*sin(kdotr);
					}
					
					a = ak[0];
					b = ak[1];
					akak = (a*a + b*b);
					tmp = recip * exp(GAMMA*ksq)/ksq;
					
					for(i=0;i<ParticleN;i++){
						kdotr = mx*r_pos[i][0] + my*r_pos[i][1] + mz*r_pos[i][2];
						qi = prop[i][0];
						tmp2 = 2*tmp*qi*(sin(kdotr) * a + cos(kdotr) * b);

						f_kewald[i][0] += tmp2 * mx;
						f_kewald[i][1] += tmp2 * my;
						f_kewald[i][2] += tmp2 * mz;
					}
					e_kewald += tmp * akak;
				}
			}
		}
	}
	
	printf("PME KEwald: %lf\n", e_kewald);
	return e_kewald;
}

double PME_Self(double **r_pos, double **prop){
	int i=0;
	double e_selfewald = 0.0;
	for(i=0;i<ParticleN;i++){
		e_selfewald += prop[i][0]*prop[i][0];
	}
	e_selfewald *= -ONE_PI_EPS0*alpha/sqrt(M_PI);
	
	printf("PME Self_Ewald: %lf\n", e_selfewald);
	return e_selfewald;
} 

// Number of Particles in the System //
int ParticleN=8;
int dimen=3; // dimensionality of the system

double default_mass=1.0;
double default_gama=0.7;

// Box boundaries [0,L]//
double Lx=32.0;
double Ly=32.0;
double Lz=32.0;

// Box boundaries [-L,L] //
double Ext_Lx=1.0;
double Ext_Ly=1.0;
double Ext_Lz=1.0;

// Ewald parameters //
double ONE_PI_EPS0=1.0;
double alpha=4.0;
int PME_nmax=5; // Real space cells
int PME_kmax=13; // K space max vectors

Full Implementation

For the full Implementation and the test results, please follow my Github link: 3D Ewald Summation in C.