The reciprocal energy is calculated in the scalar_sum() function in charge_grid.c (called by do_pmesh_kspace() in pmesh_kspace.c). The Ewald reciprocal energy is calculated according to the equation 4.7 :
The scalar_sum() function goes through all grid points using the counter variable named “ind”. For each grid point, it calculates its contribution to the reciprocal energy by applying equation 4.7 . The following terms in equation 4.7 are calculated or used: m2, the exponential term, the B(m1, m2, m3), and the IFFT transformed value of the charge array element. Also, in this step, the Q charge array is overwritten by the product of itself with the arrays B and C which are defined in equation 3.9 and 4.8 respectively. This new Q array which is equivalent to B*C*F-1(Q) is used in the reciprocal force calculation step. In the software implementation, the viral is also calculated; however, it is beyond the discussion of this appendix.
One thing is worthwhile to notice is that the indexes (m1, m2, and m3) for the C array is shifted according to the statement in :: “with m defined by m1’a1* + m2’a2* + m3’a3* , where mi’= mi for 0<mi<K/2 and mi’ = mi - Ki otherwise”.
In scalar_sum() function in charge_grid.c… int scalar_sum(double *q, double *ewaldcof, double *volume, double *recip,
double *bsp_mod1, double *bsp_mod2, double *bsp_mod3, int *nfft1, int *nfft2,
int *nfft3, int *nfftdim1, int *nfftdim2, int *nfftdim3, double *eer, double *vir)
3.2.10 Computation of F(Q) using 3D-FFT and Q Array Update
Similar to the IFFT operation described in the section 3.2.8, the 3D FFT is done in fft_forward() function in the fftcall.c. It invokes the dynamic library libpubfft.a to perform the 3D-FFT operation. The transformed elements are stored in the original charge array; hence, this is termed in-place FFT operation. The 3D FFT is performed on the updated Q charge array (in step 3.2.9) to obtain the convolution (θrec * Q) which is necessary to calculate the reciprocal forces (as shown in equation 4.9).
The following calculation shows the derivation of the convolution result. As described in the section 3.2.9, the updated Q array contains the value of B∙C∙F-1(Q).
(θrec * Q)
In the fft_forward() function in fftcalls.c… int fft_forward(double *array, double *fftable, double *ffwork, int *nfft1, int *nfft2, int *nfft3, int *nfftdim1, int *nfftdim2, int *nfftdim3, int *nfftable, int *nffwork)
The reciprocal forces for all particles are calculated in the grad_sum() function in charge_grid.c (called by do_pmesh_kspace() in pmesh_kspace.c).
This function goes through all particles. For each particle, it identifies all the grid points that this particle has been interpolated to. Then, it calculates the contribution of the reciprocal force exerted on this charged particle by each of these grid points. For a 3-D simulation space, there are three components of the forces, namely, the x, y and z directional forces. The convolution (θrec * Q) is already resided in the Q charge array which is calculated in the operation described in the section 3.2.10. The derivates of the B-Spline coefficients (dtheta array) which are needed in the calculation of the ∂Q/∂rαi are already computed in the program step described in the section 3.2.6.
For x direction: ∂Q/∂rαi = q*dthetax*thetay*thetaz.
For y direction, ∂Q/∂rαi = q*thetax*dthetay*thetaz.
For z direction, ∂Q/∂rαi = q*thetax*thetay*dthetaz.
In the grad_sum( ) function in charge_grid.c… int grad_sum( int *numatoms, PmeParticlePtr ParticlePtr,