TrinkleGroup · ydanuiuc · Jan 10, 2025
diff --git a/README.md b/README.md
@@ -1,3 +1,40 @@
+# Charge Smoothing
+This is the charge smoothing branch of the Bader Integration repository. Although it is not directly related to Bader integration or the weight method, it was developed from some of the code pieces in this repository, especially the CHGCAR manipulation routines in directory ```chgcar-manipulation``` written by Prof. Trinkle.
+
+Here, charge smoothing refers to suppressing the high-frequency components of a charge density in the reciprocal space. Because these high-frequency components could lead to numerical noise in the real space, such operation is possible to make the charge density a smoother quantity. In this code, the smoothing procedure is decribed as
+
+$$\rho(\mathbf{r}) \xrightarrow{\text{FFT}} \rho(\mathbf{G}) \xrightarrow{\text{smoothing}} \rho'(\mathbf{G}) \xrightarrow{\text{IFFT}} \rho'(\mathbf{r}) $$
+
+where we take the original charge density $\rho(\mathbf{r})$, do a fast Fourier transform (FFT), apply the smoothing method or a low-pass filter, and get the smoothed charge density $\rho'(\mathbf{r})$ after an inverse fast fourier transform (IFFT).
+
+Two smoothing methods or low-pass filters are implemented in file ```smooth.H```: a power-law truncation filter and a Gaussian filter. The code extracts the parameters of the chosen filter by fitting a power-law/Gaussian curve to a reference (expected to be smooth and well-behaved) charge density, and then apply the filter to the target charge density to smooth it. See functions ```extract_filter()``` and ```apply_filter()```, or ```extract_filter_gaussian()``` and ```apply_filter_gaussian()``` in file ```smooth.H``` for more details.
+
+Two parameters are currently hard-coded in this code: an integer ```GAUSSIAN``` which controls whether the Gaussian filter or the power-law truncation filter is to be used, and a double-type ```FILTER_MULT``` which manually scales the parameter ```l0``` of the Gaussian filter. (Yang found in his own research that the extracted Gaussian filter was not the optimal one, so he made manual adjustments and chose to determine the parameter with other methods.) 
+
+Please read the source code in ```chgcar_smoothing.C``` and ```smoothing.H```, understand what it is doing before reusing it.
+
+The charge smoothing source code is put in directory ```chgcar-smoothing```. To compile the code, run command
+
+    make
+
+A cpp compiler such as g++ is required for compilation. This code uses the GNU Scientific Library (gsl library). In Ubuntu, this library can be acquired by running ```sudo apt-get install libgsl-dev```.
+
+To run the compiled binary ```chgcar_smoothing```, run command
+
+    ./chgcar_smoothing CHGCAR_ref CHGCAR_tgt > CHGCAR_smoothed
+
+where ```CHGCAR_ref``` is the reference charge density for filter extraction, ```CHGCAR_tgt``` is the target charge density to be smoothed, and ```CHGCAR_smoothed``` is the smoothed charge density. It is assumed that the charge density files are in the CHGCAR-format of VASP 4. For files generated by VASP 5 or newer versions, please remove the line in the CHGCAR file header that specifies the ion species; alternatively, consider updating the code (one can replace the ```chgcar.H``` file here with the updated ```chgcar.H``` file in the master branch, and add a line in file ```chgcar_smoothing.C``` to output the ion species).
+
+This code was developed by Yang Dan, and used in the publication below. Details about charge smoothing are provided in the supplementary material.
+
+* Y. Dan and D. R. Trinkle, "First-principles core energies of isolated basal and prism screw dislocations in magnesium." *Mater. Res. Lett.* **10**, 360–368 (2022). [doi](https://doi.org/10.1080/21663831.2022.2051763)
+
+***
+
+Below is the original README.
+
+***
+
 # Bader integration using weighted integration
 This code takes the grid-based output from [VASP](https://www.vasp.at/) as from CHGCAR or AECCAR files and performs *Bader integration* over the basins of attraction. The call
 

diff --git a/chgcar-smoothing/3dfft.H b/chgcar-smoothing/3dfft.H
@@ -0,0 +1,184 @@
+#ifndef __3DFFT_H
+#define __3DFFT_H
+#include <math.h>
+#include <gsl/gsl_fft_complex.h>
+
+#define REAL(z,i) ((z)[2*(i)])
+#define IMAG(z,i) ((z)[2*(i)+1])
+
+int forward_3dfft (double* Z, int N0, int N1, int N2, double* Xr, double* Xi);
+
+int inverse_3dfft (double* Xr, double* Xi, int N0, int N1, int N2, double* Z);
+
+// A grid N0 x N1 x N2.  Our indexing of Z is (fortran?) compatible with
+// CHGCAR: Z[i + N0*(j + N1*k)] == Z[i][j][k]
+// Assume EVERYTHING is already allocated.
+int forward_3dfft (double* Z, int N0, int N1, int N2, double* Xr, double* Xi) 
+{
+  // Wavetables and workspace:
+  gsl_fft_complex_wavetable* wave[3];
+  wave[0] = gsl_fft_complex_wavetable_alloc(N0);
+  wave[1] = gsl_fft_complex_wavetable_alloc(N1);
+  wave[2] = gsl_fft_complex_wavetable_alloc(N2);
+
+  gsl_fft_complex_workspace* work[3];
+  work[0] = gsl_fft_complex_workspace_alloc(N0);
+  work[1] = gsl_fft_complex_workspace_alloc(N1);
+  work[2] = gsl_fft_complex_workspace_alloc(N2);
+
+  // dataspace: ("packed" gsl complex arrays)
+  double* data[3], *d;
+  data[0] = new double[2*N0];
+  data[1] = new double[2*N1];
+  data[2] = new double[2*N2];
+
+  int i, j, k;
+  int nind, N01=N0*N1;
+
+  for (k=0; k<N2; ++k) {
+
+    d = data[0];
+    for (j=0; j<N1; ++j) {
+      // transform Z(i)[j,k] -> Z~(i)[j,k]
+      nind = N0*(j+N1*k);
+      for (i=0; i<N0; ++i) {
+	REAL(d, i) = Z[i+nind];
+	IMAG(d, i) = 0;
+      }
+      gsl_fft_complex_forward(d, 1, N0, wave[0], work[0]);
+      for (i=0; i<N0; ++i) {
+	Xr[i+nind] = REAL(d, i);
+	Xi[i+nind] = IMAG(d, i);
+      }
+    }
+
+    d = data[1];
+    for (i=0; i<N0; ++i) {
+      // transform X[i](j)[k] -> X~[i](j)[k]
+      nind = i+N0*N1*k;
+      for (j=0; j<N1; ++j) {
+	REAL(d, j) = Xr[j*N0+nind];
+	IMAG(d, j) = Xi[j*N0+nind];
+      }
+      gsl_fft_complex_forward(d, 1, N1, wave[1], work[1]);
+      for (j=0; j<N1; ++j) {
+	Xr[j*N0+nind] = REAL(d, j);
+	Xi[j*N0+nind] = IMAG(d, j);
+      }
+    }
+  }
+
+  d = data[2];
+  for (i=0; i<N0; ++i)
+    for (j=0; j<N1; ++j) {
+      nind = i + j*N0;
+      for (k=0; k<N2; ++k) {
+	REAL(d, k) = Xr[k*N01+nind];
+	IMAG(d, k) = Xi[k*N01+nind];
+      }
+      gsl_fft_complex_forward(d, 1, N2, wave[2], work[2]);
+      for (k=0; k<N2; ++k) {
+	Xr[k*N01+nind] = REAL(d, k);
+	Xi[k*N01+nind] = IMAG(d, k);
+      }
+    }
+
+  // Garbage collection:
+  for (i=0; i<3; ++i) delete[] data[i];
+  for (i=0; i<3; ++i) gsl_fft_complex_workspace_free(work[i]);
+  for (i=0; i<3; ++i) gsl_fft_complex_wavetable_free(wave[i]);
+
+  return 0;
+}
+
+// Note: this is *destructive* with Xr and Xi.
+int inverse_3dfft (double* Xr, double* Xi, int N0, int N1, int N2, double* Z) 
+{
+  // Wavetables and workspace:
+  gsl_fft_complex_wavetable* wave[3];
+  wave[0] = gsl_fft_complex_wavetable_alloc(N0);
+  wave[1] = gsl_fft_complex_wavetable_alloc(N1);
+  wave[2] = gsl_fft_complex_wavetable_alloc(N2);
+
+  gsl_fft_complex_workspace* work[3];
+  work[0] = gsl_fft_complex_workspace_alloc(N0);
+  work[1] = gsl_fft_complex_workspace_alloc(N1);
+  work[2] = gsl_fft_complex_workspace_alloc(N2);
+
+  // dataspace: ("packed" gsl complex arrays)
+  double* data[3], *d;
+  data[0] = new double[2*N0];
+  data[1] = new double[2*N1];
+  data[2] = new double[2*N2];
+
+  int i, j, k;
+  int nind, N01=N0*N1;
+
+  for (k=0; k<N2; ++k) {
+
+    d = data[0];
+    for (j=0; j<N1; ++j) {
+      // transform Z(i)[j,k] -> Z~(i)[j,k]
+      nind = N0*(j+N1*k);
+      for (i=0; i<N0; ++i) {
+	REAL(d, i) = Xr[i+nind];
+	IMAG(d, i) = Xi[i+nind];;
+      }
+      gsl_fft_complex_inverse(d, 1, N0, wave[0], work[0]);
+      for (i=0; i<N0; ++i) {
+	Xr[i+nind] = REAL(d, i);
+	Xi[i+nind] = IMAG(d, i);
+      }
+    }
+
+    d = data[1];
+    for (i=0; i<N0; ++i) {
+      // transform X[i](j)[k] -> X~[i](j)[k]
+      nind = i+N0*N1*k;
+      for (j=0; j<N1; ++j) {
+	REAL(d, j) = Xr[j*N0+nind];
+	IMAG(d, j) = Xi[j*N0+nind];
+      }
+      gsl_fft_complex_inverse(d, 1, N1, wave[1], work[1]);
+      for (j=0; j<N1; ++j) {
+	Xr[j*N0+nind] = REAL(d, j);
+	Xi[j*N0+nind] = IMAG(d, j);
+      }
+    }
+  }
+
+  // Last step--finally overwrite Z:
+  int is_real=1;
+  double sum_imag = 0;
+  d = data[2];
+  for (i=0; i<N0; ++i)
+    for (j=0; j<N1; ++j) {
+      nind = i + j*N0;
+      for (k=0; k<N2; ++k) {
+	REAL(d, k) = Xr[k*N01+nind];
+	IMAG(d, k) = Xi[k*N01+nind];
+      }
+      gsl_fft_complex_inverse(d, 1, N2, wave[2], work[2]);
+      for (k=0; k<N2; ++k) {
+	Z[k*N01+nind] = REAL(d, k);
+	sum_imag += IMAG(d,k) * IMAG(d,k);
+	is_real = is_real && (zero(IMAG(d, k)));
+      }
+    }
+
+  sum_imag = sqrt(sum_imag/(N01*N2));
+  if (! is_real)
+    fprintf(stderr, "## iFFT had rms imaginary component: %.8le\n", sum_imag);
+
+  // Garbage collection:
+  for (i=0; i<3; ++i) delete[] data[i];
+  for (i=0; i<3; ++i) gsl_fft_complex_workspace_free(work[i]);
+  for (i=0; i<3; ++i) gsl_fft_complex_wavetable_free(wave[i]);
+
+  return !is_real;
+}
+
+#undef REAL
+#undef IMAG
+
+#endif
diff --git a/chgcar-smoothing/Makefile b/chgcar-smoothing/Makefile
@@ -0,0 +1,20 @@
+CC:=g++
+INCLUDE:=-I.
+LIBS:=-L.
+LIBM:=-lm
+LIBGSL:=-lgsl -lgslcblas
+
+INCLUDES = *.H 
+
+TARGET = chgcar_smoothing
+
+all: $(TARGET)
+
+chgcar_smoothing.o: chgcar_smoothing.C $(INCLUDES)
+	$(CC) -Wall -c chgcar_smoothing.C
+
+chgcar_smoothing: chgcar_smoothing.o
+	$(CC) -o chgcar_smoothing chgcar_smoothing.o $(LIBGSL) $(LIBM)
+
+clean:
+	rm -f *.o chgcar_smoothing