Exact exchange for PARSEC
Home page
This web page is the pseudocode of the algorithm done by Emmanuel Lorin de la Grandmaison, that can be used to solve the right-hand-side term for the EXX. It happens to be the same term as the coupling matrix appearing in TDDFT (1st term, second term is the second derivative of Exc. Don't touch that for now.)
Expand   Collapse
    int main(int argc, char* argv[])
  • Initialization
    • MPI Initialization
      • MPI_Init
      • MPI_Comm_size
      • MPI_Comm_rank
    • time (&Total_time_start, &ReadInput_time_start)
    • open file "input.dat", variable: infile
    • read using fortran routines in file "read_write.f", unformatted, la messe, quoi!
    • MPI_Bcast (&nstate, &ndim, &ngrid, &grid_step) [fortran parameters]
    • read_data (C++ method that calls the FORTRAN subroutine read_potold inside "read_write.f". Encore la messe, quoi! Defined in "../INC/f_header.h"
    • MPI_Bcast (eig, occp, kx, ky, kz, wf)
    • Defines variable "elements" from the eigen class (../INC/eigen.h): eigen* elements = new eigen(nstate).
    • elements->values = eig; elements->occupied = occp
    • Defines variable "psi" as class wave_function defined in ../INC/wave_function.h. psi->f[][][] = wf[j+tmp], where "wf" is in the fortran format and "psi->f" is in the C++ format.
    • "test.out" is the output file, instantiated by "outfile".
    • delete eig, occp, kx, ky, kz, wf; kx, ky, and kz don't seem to be used later???
    • define "potential" from class wave_function, don't understand why he does that ???
    • outputs min and max occupied levels: minfil, maxfil
    • End of reading input files, apparently
    • Prepares the output file "spectr.infl" and prints in it state occupations and eigenvalues. Dipole transitions are not yet printed, apparently.
    • Define variables "physics" as class macro_physics
    • The variable "physics->density" seems to be the density matrix
    • Apply a cutoff on the density: that should be according to equation just before section 5.2 of the article [CPC, 167 (2005), 7-22]; "density_cutoff". Then, only grid points > density_cutoff are kept; variable "eff_grid_pts"
    • Attention: "wf" is redefined here as a new double, then wf[] = psi->f[][][] only on the effective grid. That means that at this point the wave functions is cut off.
  • Calculation of dVxc/drho
    • skip that part for now, I don't need it for the EXX
  • Coupling matrix calculations (look for "cout << "construction of the coupling matrix" << endl);
    • Define psi_ij as wave_function class pointer; it looks like this variable is also the short range part
    • Define "long_range" as multipole class pointer
    • Define rowArray, cnts, value_global
    • Define PD as par_decomp for parallel decomposition set in ../LIB/par_decomp.h
    • PD->distribute_rows() on the children processors
    • Define psi_ij, phi_ij, psi_ij_bar, psi_ij_bar_k
    • restrict_pot() is part of ../LIB/product.h it puts pot = potential->f[][][]
    • Define psi_ij_kl: should be containing the product of the 4 wavefunctions, maybe, unclear
    • Initialize FFTw variables fftw_plan, etc. Variables are: *freal, *fmpole, *aux, *fgreen, *fcmplx (which is fftw_complex), also the plan where calculation is performed: fwd_plan (forward FFT), bck_plan (backward FFT).
    • setup_green_fourier(): defined in product.cpp. It looks like Figure 3 on the CPC article, the supercell method
    • Define the plan: fwd_plan = fftw_plan_dft_r2c_3d(real to complex) and bck_plan = fftw_plan_dft_c2r_3d(complex to real)
    • [page 10: the core]: Coupling matrix calculation is between time(&CoupMat_time_start); and time(&CoupMat_time_finish);
    • time(&CoupMat_time_start);
    • for (int I=0;I < maxfil;I++)
      • for (int J=max(minfil,I);J < nb_wave_functions;J++)
        • calculates occupation ("occuped") and eigenvalues differences: elements->values[J] - elements->values[I]
        • Verify that the differences in eigenvalues is relevant: "if(trlda0 ≤ trmax)"
          • k0 += 1; condition on k0; k1 += 1;
          • Some checks on the child processor, look for "k0" and "k1", etc.
          • time(&Solve_poisson_start);
          • method == 0: Fourier expansion without multipole FALR
            • dvem (C++ method for scalar product). FORTRAN: dot_product(vector_a, vector_b)
            • prolong_wf(); the offset = ngrid/2*(1- alpha), alpha is an extra parameter (1.7 check figure 3), ngrid is in the wavefunctions files
            • fftw_execute(fwd_plan): this is the forward FFT on the product of wavefunctions, see eq. (13)
            • apply_green_fourier(fcmplx, fgreen, nfft); this is the 1/k2 term, the Green's function of a pure Laplacian
            • fftw_execute(bck_plan): this is the backward FFT on the product of wavefunctions with the 1/k2 term, see eq. (13)
            • REMARK: could possibly use IBM FFT's such as ESSL's
            • restrict_wf(); use same offset as above in the prolong_wf();
            • prod_sum() just includes the calculated dVxc/drho (not useful for EXX)
          • method != 0: Fourier expansion with multipole FALR (this was all commented in the version received)
            • product_wf_shiv(fmpole, psi[I], psi[J]); defined in product.cpp, performs a scalar product with identical indices for the wave functions
            • long_range->compute_alm(), defined in multipole.cpp
            • long_range->compute_rho_short(), defined in multipole.cpp
            • fftw_execute(fwd_plan): this is the forward FFT on the product of wavefunctions, see eq. (13)
            • division_by_green_fourier(fcmplx, fgreen, nfft); this is the 1/k2 term, the Green's function of a pure Laplacian
            • fftw_execute(bck_plan): this is the backward FFT on the product of wavefunctions with the 1/k2 term, see eq. (13)
          • print out "sol.data": variable is outfile80
          • time(&Solve_poisson_finish);
          • for (int K=0; K
          • dvem (psi_ij_bar, psi_ij_bar_k). Could also use "product_wf" in product.cpp. This looks like the second integral just with a direct summation
          • for (int L=max(minfil,K); L
          • k1+=1;

This web page is regularly updated Summer 2009. This web page is for personal use; I don't see how this web page could be useful to anyone besides myself.