libvdwxc-b1e0dd854310410012d05daf4c6661b49f31b2ab/src/vdw_core.c

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <complex.h>

#include "vdw_kernel.h"
#include "vdwxc.h"
#include "vdw_q0.h"

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "vdw_core.h"

#if defined VDW_TIMING
#include "vdw_timing.h"
#endif

#define VDWXC_DEBUG 0

// Get name of functional from numeric code.
const char* vdwxc_funcname(int functional)
{
    static char* names[] = { "vdW-DF", "vdW-DF2", "vdW-DF-cx", "custom", "unknown" };
    switch(functional) {
    case FUNC_VDWDF:
        return names[0];
    case FUNC_VDWDF2:
        return names[1];
    case FUNC_VDWDFCX:
        return names[2];
    case FUNC_CUSTOM:
        return names[3];
    default:
        return names[4];
    }
}

// Get name of parallelization mode from numeric code.
const char* vdwxc_fftname(int fftcode)
{
    static char* names[] = {"fftw-serial", "fftw-mpi", "pfft", "unknown"};
    switch(fftcode) {
    case VDWXC_FFTW_SERIAL:
        return names[0];
    case VDWXC_FFTW_MPI:
        return names[1];
    case VDWXC_PFFT:
        return names[2];
    default:
        return names[3];
    }
}

const char* vdwxc_errstring(int errcode)
{
    static char* errstrs[] = {
        "No errorcode set.",
        "struct vdwxc_data uninitialized. Call vdwxc_initialize(...);",
        "Unkown functional code.",
        "Unit cell not set. Call vdwxc_set_unit_cell(...);",
        "Unit cell already set in call to vdwxc_set_unit_cell(...);",
        "Radial grid not set. Call vdwxc_set_radial_grid(...);",
        "Unknown error code or uninitialized vdwxc_data object."
    };
    if(errcode >= 0 && errcode <= 5)
        return errstrs[errcode];
    else
        return errstrs[6];
}

void vdwxc_print_error(vdwxc_data data)
{
    // Internal error: Print the error and abort program.
    // We should implement something nicer for user errors though.
    printf("%s (errcode=%d)\n",
           vdwxc_errstring(data->errorcode), data->errorcode);
    if(data->errorcode != VDWXC_ERROR_NOERROR) {
        assert(0);
    }
}

// Get the size of a vdwxc_data.
int vdwxc_get_struct_size(void)
{
    return sizeof(struct vdwxc_data_obj);
}

// Set all bytes of vdwxc_data to zero.
void vdwxc_nullify(vdwxc_data data)
{
    memset(data, 0, vdwxc_get_struct_size());
}

// Set functional type and some default values.
void vdwxc_set_defaults(vdwxc_data vdw, int functional, int nspin)
{
    switch (functional) {
    case FUNC_VDWDF:
        vdw->Z_ab = -0.8491;
        break;
    case FUNC_VDWDF2:
        vdw->Z_ab = -1.887;
        break;
    case FUNC_VDWDFCX:
        vdw->Z_ab = -0.8491;
        break;
    default:
        printf("unknown vdw-df functional %d\n", functional);
        assert(0);
        break;
    }
    vdw->functional = functional;
    assert((nspin == 1) || (nspin == 2));
    vdw->nspin = nspin;
    vdw->q_cut = 5.0; // XXX to be declared by kernel
    vdw->fft_type = VDWXC_FFTW_SERIAL;
    vdw->mpi_size = 1;
    vdw->initialized = VDWXC_INITIALIZED;
}

// Allocate a new vdwxc_data and return a pointer to it.

vdwxc_data vdwxc_new_anyspin(int functional, int nspin)
{
    vdwxc_data vdw = malloc(vdwxc_get_struct_size());
    vdwxc_nullify(vdw);
    vdwxc_set_defaults(vdw, functional, nspin);
    vdw->kernel = vdwxc_default_kernel();
    return vdw;
}

vdwxc_data vdwxc_new(int functional)
{
    return vdwxc_new_anyspin(functional, 1);
}

vdwxc_data vdwxc_new_spin(int functional)
{
    return vdwxc_new_anyspin(functional, 2);
}


// Write unit cell to string
int vdwxc_sprint_cell(char* str, struct vdwxc_unit_cell* cell)
{
    char* ch = str;
    ch += sprintf(ch, "    Nglobal = %d x %d x %d\n", cell->Nglobal[0], cell->Nglobal[1], cell->Nglobal[2]);
    ch += sprintf(ch, "    Nlocal  = %d x %d x %d\n", cell->Nlocal[0], cell->Nlocal[1], cell->Nlocal[2]);
    ch += sprintf(ch, "    offset  = %d, %d, %d\n", cell->offset[0], cell->offset[1], cell->offset[2]);

    double lengths[3];
    int i, j;
    for(i=0; i < 3; i++) {
        lengths[i] = 0.0;
        for(j=0; j < 3; j++) {
            lengths[i] += cell->vec[3 * i + j] * cell->vec[3 * i + j];
        }
        lengths[i] = sqrt(lengths[i]);
    }

    ch += sprintf(ch, "    Cell lengths %f, %f, %f\n", lengths[0], lengths[1], lengths[2]);
    ch += sprintf(ch, "      %f, %f, %f\n", cell->vec[0], cell->vec[1], cell->vec[2]);
    ch += sprintf(ch, "      %f, %f, %f\n", cell->vec[3], cell->vec[4], cell->vec[5]);
    ch += sprintf(ch, "      %f, %f, %f\n", cell->vec[6], cell->vec[7], cell->vec[8]);
    ch += sprintf(ch, "    dV %e\n", cell->dV);
    int nbytes = ch - str;
    return nbytes;
}

#ifdef HAVE_MPI
int vdwxc_has_mpi(void){ return 1; }
#else
int vdwxc_has_mpi(void){ return 0; }
#endif

#ifdef HAVE_PFFT
int vdwxc_has_pfft(void){ return 1; }
#else
int vdwxc_has_pfft(void){  return 0; }
#endif

// Print representation of vdwxc_data to stdout.
void vdwxc_print(vdwxc_data data)
{
    // Buffer can hold 200 lines of full-width standard-size terminal.
    // We will probably not print that much, ever.
    int maxsize = 80 * 200;
    char str[maxsize];
    vdwxc_tostring(data, maxsize, str);
    // printf(str) yields a warning
    printf("%s", str);
}

// Write representation of vdwxc_data to string.  Returns number of bytes written.
// Truncates to maxsize if necessary.
int vdwxc_tostring(vdwxc_data data, int maxsize, char *str)
{
    // We build within our own buffer, then copy into user-supplied str.
    int maxlen = 200 * 80; // max 200 lines of length 80
    if(maxsize > maxlen)
        maxlen = maxsize;
    char buf[maxlen];
    char* ch = buf;
    ch += sprintf(ch, "  === vdW-DF data at %p ===\n", data);
    ch += sprintf(ch, "  functional %s [%d]\n", vdwxc_funcname(data->functional), data->functional);
    if (data->nspin == 1) {
        ch += sprintf(ch, "  spins: %d (spin-paired)\n", data->nspin);
    } else {
        ch += sprintf(ch, "  spins: %d (spin-polarized)\n", data->nspin);
    }
    ch += sprintf(ch, "  Z=%f :: qcut=%f\n", data->Z_ab, data->q_cut);
    ch += sprintf(ch, "  cell\n");
    ch += vdwxc_sprint_cell(ch, &data->cell);
    ch += sprintf(ch, "  icell\n");
    ch += vdwxc_sprint_cell(ch, &data->icell);
    ch += sprintf(ch, "  FFT: %s [%d]\n", vdwxc_fftname(data->fft_type), data->fft_type);
    if(data->fft_type == VDWXC_PFFT) {
        ch += sprintf(ch, "  pfft grid %d x %d\n", data->pfft_grid[0], data->pfft_grid[1]);
    }
#ifdef HAVE_MPI
    ch += sprintf(ch, "  MPI :: rank=%d :: size=%d :: ptr=<%p>\n",
                  data->mpi_rank, data->mpi_size, data->mpi_comm);
#endif
    ch += sprintf(ch, "  alloc Ng=%d Nglobal=%d gLDA=%d\n",
                  data->Ng, data->Nglobal, data->gLDA);
    ch += sprintf(ch, "  kLDA=%d\n", data->kLDA);
    ch += sprintf(ch, "  ptrs1 q0=<%p> work_ka=<%p>\n",
                  data->q0_g, data->work_ka);
    ch += sprintf(ch, "  ptrs2 rho_dq0drho_sg=<%p> rho_dq0dsigma_sg=<%p>\n",
                  data->rho_dq0drho_sg, data->rho_dq0dsigma_sg);
    ch += sprintf(ch, "  fftw ptrs  r2c=<%p> c2r=<%p>\n", data->plan_r2c, data->plan_c2r);
#ifdef HAVE_PFFT
    ch += sprintf(ch, "  pfft ptrs r2c=<%p> c2r=<%p>\n", data->pfft_plan_r2c, data->pfft_plan_c2r);
#endif
    ch += sprintf(ch, "  Error state %d: %s\n", data->errorcode, vdwxc_errstring(data->errorcode));
    ch += sprintf(ch, "  =============================\n");
    int len = strlen(buf);
    assert(ch - buf == len);
    assert(len < maxlen);
    assert(len < maxsize);
    strcpy(str, buf);
    return len;
}

void vdwxc_set_unit_cell(vdwxc_data data,
                         int Nx, int Ny, int Nz,
                         double C00, double C10, double C20,
                         double C01, double C11, double C21,
                         double C02, double C12, double C22)
{
    double det = (C00 * (C11 * C22 - C21 * C12) \
                  - C01 * (C10 * C22 - C12 * C20) \
                  + C02 * (C10 * C21 - C11 * C20));
    data->cell.dV = det / (Nx * Ny * Nz);
    data->cell.vec[0] = C00;
    data->cell.vec[1] = C01;
    data->cell.vec[2] = C02;
    data->cell.vec[3] = C10;
    data->cell.vec[4] = C11;
    data->cell.vec[5] = C12;
    data->cell.vec[6] = C20;
    data->cell.vec[7] = C21;
    data->cell.vec[8] = C22;

    data->cell.Nglobal[0] = Nx;
    data->cell.Nglobal[1] = Ny;
    data->cell.Nglobal[2] = Nz;

    // Calculate reciprocal cell
    double idet = 1.0 / (C00 * (C11 * C22 - C21 * C12) \
                         - C01 * (C10 * C22 - C12 * C20) \
                         + C02 * (C10 * C21 - C11 * C20)) * 2 * VDWXC_PI;
    data->icell.Nglobal[0] = Nx;
    data->icell.Nglobal[1] = Ny;
    data->icell.Nglobal[2] = Nz / 2 + 1;
    data->icell.dV = 1e100;
    data->icell.vec[0] =  (C11 * C22 - C21 * C12) * idet;
    data->icell.vec[1] = -(C01 * C22 - C02 * C21) * idet;
    data->icell.vec[2] =  (C01 * C12 - C02 * C11) * idet;
    data->icell.vec[3] = -(C10 * C22 - C12 * C20) * idet;
    data->icell.vec[4] =  (C00 * C22 - C02 * C20) * idet;
    data->icell.vec[5] = -(C00 * C12 - C10 * C02) * idet;
    data->icell.vec[6] =  (C10 * C21 - C20 * C11) * idet;
    data->icell.vec[7] = -(C00 * C21 - C20 * C01) * idet;
    data->icell.vec[8] =  (C00 * C11 - C10 * C01) * idet;

    data->Nglobal = data->cell.Nglobal[0]
                  * data->cell.Nglobal[1]
                  * data->cell.Nglobal[2];
}

void vdwxc_allocate_buffers(vdwxc_data data)
{
    data->Ng = data->cell.Nlocal[0] * data->cell.Nlocal[1] * data->cell.Nlocal[2];

    size_t gridsize = data->Ng * sizeof(double);

    data->q0_g = (double *) malloc(gridsize);
    data->rho_dq0drho_sg = (double *) malloc(gridsize * data->nspin);
    data->rho_dq0dsigma_sg = (double *) malloc(gridsize * data->nspin);
}

void vdwxc_init_serial(vdwxc_data data)
{
    data->kLDA = data->icell.Nglobal[2];
    data->gLDA = 2 * data->kLDA;
    int nktotal = data->cell.Nglobal[0] * data->cell.Nglobal[1] * data->kLDA * data->kernel.nalpha;
    data->work_ka = (double complex*) fftw_malloc(nktotal * sizeof(double complex));
    int i;
    for(i = 0; i < 3; i++) {
        data->cell.Nlocal[i] = data->cell.Nglobal[i];
        data->icell.Nlocal[i]= data->cell.Nglobal[i];
    }
    data->icell.Nlocal[2] = data->icell.Nlocal[2] / 2 + 1;

    data->plan_r2c = fftw_plan_many_dft_r2c(3,
                                            data->cell.Nglobal,
                                            data->kernel.nalpha,
                                            (double*)data->work_ka,
                                            NULL,
                                            data->kernel.nalpha,
                                            1,
                                            data->work_ka,
                                            NULL,
                                            data->kernel.nalpha,
                                            1,
                                            FFTW_ESTIMATE);
    data->plan_c2r = fftw_plan_many_dft_c2r(3,
                                            data->cell.Nglobal,
                                            data->kernel.nalpha,
                                            data->work_ka,
                                            NULL,
                                            data->kernel.nalpha,
                                            1,
                                            (double*)data->work_ka,
                                            NULL,
                                            data->kernel.nalpha,
                                            1,
                                            FFTW_ESTIMATE);
    assert(data->plan_r2c != NULL);
    assert(data->plan_c2r != NULL);
    vdwxc_allocate_buffers(data);
}

void vdwxc_compute_thetas(struct vdwxc_kernel *kernel,
                        int nx, int ny, int nz, int leading_dim,
                        double* rho, double* q0, double* output)
{
    // We have to translate between the contiguous array rho
    // and the padded array output.
    int ix, iy, iz;
    int ixyz_fftw;
    for(ix = 0; ix < nx; ix++) {
        for(iy = 0; iy < ny; iy++) {
            for(iz = 0; iz < nz; iz++) {
                //ijk = iz + nz * (iy + (ny * ix));
                ixyz_fftw = (iz + leading_dim * (iy + ny * ix));
                vdwxc_evaluate_palpha_splines(kernel, *rho, *q0, output + ixyz_fftw * kernel->nalpha);
                rho++;
                q0++;
            }
        }
    }
}


void vdwxc_write_workbuffer(vdwxc_data data, char *fname)
{
    FILE *fd;
    int ix, iy, iz;
    int ixyz_fftw;
    int alpha;
    int nx = data->cell.Nlocal[0];
    int ny = data->cell.Nlocal[1];
    int nz = data->cell.Nlocal[2];
    int leading_dim = data->gLDA;
    double *work_ka = (double*)data->work_ka;
    fd = fopen(fname, "w");
    for(ix = 0; ix < nx; ix++) {
        for(iy = 0; iy < ny; iy++) {
            for(iz = 0; iz < nz; iz++) {
                //ijk = iz + nz * (iy + (ny * ix));
                ixyz_fftw = (iz + leading_dim * (iy + ny * ix));
                for(alpha=0; alpha < data->kernel.nalpha; alpha++) {
                    fprintf(fd, "%.14e\n", work_ka[alpha + ixyz_fftw * data->kernel.nalpha]);
                }
            }
        }
    }
    fclose(fd);
}

// XXXXX this is a copy of calculate_thetas.  Fix this!
void vdwxc_compute_thetas_spin(struct vdwxc_kernel *kernel,
                             int nx, int ny, int nz, int leading_dim,
                             double *rho_up, double *rho_dn,
                             double *q0, double *output)
{
    // We have to translate between the contiguous array rho
    // and the padded array output.
    int ix, iy, iz;
    int ixyz_fftw;
    for(ix = 0; ix < nx; ix++) {
        for(iy = 0; iy < ny; iy++) {
            for(iz = 0; iz < nz; iz++) {
                //ijk = iz + nz * (iy + (ny * ix));
                double rho = (*rho_up) + (*rho_dn);
                ixyz_fftw = (iz + leading_dim * (iy + ny * ix));
                vdwxc_evaluate_palpha_splines(kernel, rho, *q0, output + ixyz_fftw * kernel->nalpha);
                rho_up++;
                rho_dn++;
                q0++;
            }
        }
    }
}

void vdwxc_calculate_thetas(vdwxc_data data)
{
    if(data->nspin == 1) {
        vdwxc_compute_thetas(&data->kernel,
                           data->cell.Nlocal[0],
                           data->cell.Nlocal[1],
                           data->cell.Nlocal[2],
                           data->gLDA,
                           data->rho, data->q0_g,
                           (double*)data->work_ka);
    } else if (data->nspin == 2) {
        vdwxc_compute_thetas_spin(&data->kernel,
                                data->cell.Nlocal[0],
                                data->cell.Nlocal[1],
                                data->cell.Nlocal[2],
                                data->gLDA,
                                data->rho_up, data->rho_dn, data->q0_g,
                                (double*)data->work_ka);
    } else {
        assert(0);
    }
}

// Calculate potential and derivatives from F_ka.
void potential(struct vdwxc_kernel *kernel, int nx, int ny, int nz, int gLDA, int Nglobal,
               double *rho, double *q0, double *rho_dq0drho_sg, double *rho_dq0dsigma_sg,
               double *F_ga, double *dedn, double *dedsigma)
{
    START_TIMER(potential, "potential");
    double prefactor = 1.0 / Nglobal;
    //double* F_ga = ; // XXXXX
    int ia_fftw;
    int ixyz, ixyz_fftw;
    int ix, iy, iz;
    int Ng = nx * ny * nz;
    //double integrals[2] = {0.0, 0.0};
    int alpha;
    // TODO: Move this to inner loop to save Memory

    double *p_a = (double *)malloc(sizeof(double) * kernel->nalpha);
    double *dpdq_a = (double *)malloc(sizeof(double) * kernel->nalpha);

    for(ix = 0; ix < nx; ix++) {
        for(iy = 0; iy < ny; iy++) {
            for(iz = 0; iz < nz; iz++) {
                ixyz = iz + nz * (iy + (ny * ix));
                ixyz_fftw = (iz + gLDA * (iy + ny * ix));
                double rho_dq0drho = rho_dq0drho_sg[ixyz];
                double rho_dq0dsigma = rho_dq0dsigma_sg[ixyz];
                assert(ixyz < Ng);
                dedn[ixyz] = 0.0; // zero or not?  From perspective of calling code...
                dedsigma[ixyz] = 0.0; // zero or not?
                vdwxc_evaluate_palpha_splines(kernel, 1.0, q0[ixyz], p_a);
                vdwxc_evaluate_palpha_splines_derivative(kernel, q0[ixyz], dpdq_a);
                for(alpha=0; alpha < kernel->nalpha; alpha++) {
                    ia_fftw = ixyz_fftw * kernel->nalpha + alpha;
                    dedn[ixyz] += prefactor * F_ga[ia_fftw] *
                        (p_a[alpha] + dpdq_a[alpha] * rho_dq0drho);
                    //printf("%.14e\n", prefactor * F_ga[ia_fftw] *
                    //       (p_a[alpha] + dpdq_a[alpha] * rho_dq0drho));
                    dedsigma[ixyz] += prefactor * F_ga[ia_fftw] *
                        dpdq_a[alpha] * rho_dq0dsigma;
                }
                // Here we can evaluate energy density as well:
                //
                //   epsilon(r) = (1/2) sum[alpha] n(r) p[alpha](r) F[alpha](r)
                //
                // n(r) times p[alpha](r) is actually theta[alpha](r),
                // which we do not have because we overwrote it.
                // But I don't know if anyone actually needs epsilon(r) for anything.

                //integrals[0] += dedn_g[ixyz] * rho_g[ixyz];
                //integrals[1] += dedsigma_g[ixyz] * rho_g[ixyz];
            }
        }
    }
    free(p_a);
    free(dpdq_a);
    //integrals[0] *= data->cell.dV;
    //integrals[1] *= data->cell.dV;
    END_TIMER(potential, "potential");
    // store integrals somehow for checking/debugging?
}

// Calculate potential and derivatives from F_ka.
void vdwxc_potential(vdwxc_data data)
{
    int Ng = data->Ng;
    if(data->nspin == 1) {
        potential(&data->kernel,
                  data->cell.Nlocal[0], data->cell.Nlocal[1], data->cell.Nlocal[2],
                  data->gLDA, data->Nglobal,
                  data->rho, data->q0_g, data->rho_dq0drho_sg, data->rho_dq0dsigma_sg, (double*)data->work_ka,
                  data->dedn, data->dedsigma);
    } else {
        assert(data->nspin == 2);
        potential(&data->kernel,
                  data->cell.Nlocal[0], data->cell.Nlocal[1], data->cell.Nlocal[2],
                  data->gLDA, data->Nglobal,
                  data->rho_up, data->q0_g,
                  data->rho_dq0drho_sg, data->rho_dq0dsigma_sg,
                  (double*)data->work_ka,
                  data->dedn_up, data->dedsigma_up);
        potential(&data->kernel,
                  data->cell.Nlocal[0], data->cell.Nlocal[1], data->cell.Nlocal[2],
                  data->gLDA, data->Nglobal,
                  data->rho_dn, data->q0_g,
                  data->rho_dq0drho_sg + Ng, data->rho_dq0dsigma_sg + Ng,
                  (double*)data->work_ka,
                  data->dedn_dn, data->dedsigma_dn);
    }
}

void vdwxc_writefile(char *name, int N, double *array)
{
    FILE *fd;
    int i;
    fd = fopen(name, "w");
    for(i=0; i < N; i++) {
        fprintf(fd, "%.14e\n", array[i]);
    }
    fclose(fd);
}

void vdwxc_dump(vdwxc_data data)
{
    if(data->nspin == 1) {
        vdwxc_writefile("libvdwxc.rho.dat", data->Ng, data->rho);
        vdwxc_writefile("libvdwxc.sigma.dat", data->Ng, data->sigma);
        vdwxc_writefile("libvdwxc.dedn.dat", data->Ng, data->dedn);
        vdwxc_writefile("libvdwxc.dedsigma.dat", data->Ng, data->dedsigma);
    } else {
        vdwxc_writefile("libvdwxc.rho_up.dat", data->Ng, data->rho_up);
        vdwxc_writefile("libvdwxc.rho_dn.dat", data->Ng, data->rho_dn);
        vdwxc_writefile("libvdwxc.sigma_up.dat", data->Ng, data->sigma_up);
        vdwxc_writefile("libvdwxc.sigma_dn.dat", data->Ng, data->sigma_dn);
        vdwxc_writefile("libvdwxc.dedn_up.dat", data->Ng, data->dedn_up);
        vdwxc_writefile("libvdwxc.dedn_dn.dat", data->Ng, data->dedn_dn);
        vdwxc_writefile("libvdwxc.dedsigma_up.dat", data->Ng, data->dedsigma_up);
        vdwxc_writefile("libvdwxc.dedsigma_dn.dat", data->Ng, data->dedsigma_dn);
    }
    vdwxc_writefile("libvdwxc.q0.dat", data->Ng, data->q0_g);
    vdwxc_writefile("libvdwxc.rho_dq0drho_sg.dat", data->Ng * data->nspin, data->rho_dq0drho_sg);
    vdwxc_writefile("libvdwxc.rho_dq0dsigma_sg.dat", data->Ng * data->nspin, data->rho_dq0dsigma_sg);
}


void vdwxc_calculate_q0_anyspin(vdwxc_data data)
{
    if(data->nspin == 1) {
        vdwxc_calculate_q0(data->Ng, data->Z_ab, data->q_cut,
                           data->rho, data->sigma, data->q0_g,
                           data->rho_dq0drho_sg, data->rho_dq0dsigma_sg);
    } else if(data->nspin == 2) {
        vdwxc_calculate_q0_spin(data->Ng, data->Z_ab, //data->q_cut,
                                data->rho_up, data->rho_dn,
                                data->sigma_up, data->sigma_dn,
                                data->q0_g,
                                data->rho_dq0drho_sg, data->rho_dq0dsigma_sg);
    } else {
        assert(0);
    }
}

void vdwxc_get_q0(vdwxc_data data, double *q0)
{
    int i;
    for(i=0; i < data->Ng; i++) {
        q0[i] = data->q0_g[i];
    }
}

// Calculate q0 and theta from density and gradient.
void vdwxc_q0_and_theta(vdwxc_data data)
{
    if (data->initialized != VDWXC_INITIALIZED)
    {
        data->errorcode = VDWXC_ERROR_UNINITIALIZED;
        vdwxc_print_error(data);
        return;
    }

    START_TIMER(q0, "q0");
    vdwxc_calculate_q0_anyspin(data);
    END_TIMER(q0, "q0");
    START_TIMER(thetas, "thetas");
    vdwxc_calculate_thetas(data);
    END_TIMER(thetas, "thetas");
}

#define MODE_SERIAL
#define METH(name) name##_serial
#include "vdw_include.c"
#undef NAME
#undef METH
#undef MODE_SERIAL

#define MODE_MPI
#define METH(name) name##_mpi
#include "vdw_include.c"
#undef NAME
#undef METH
#undef MODE_MPI

#define MODE_PFFT
#define METH(name) name##_pfft
#include "vdw_include.c"
#undef NAME
#undef METH
#undef MODE_PFFT


void vdwxc_check_convenience_pointers(vdwxc_data data)
{
    if(data->nspin == 1) {
        assert(data->rho != NULL);
        assert(data->sigma != NULL);
        assert(data->dedn != NULL);
        assert(data->dedsigma != NULL);
    } else if(data->nspin == 2) {
        assert(data->rho_up != NULL);
        assert(data->rho_dn != NULL);
        assert(data->sigma_up != NULL);
        assert(data->sigma_dn != NULL);
        assert(data->dedn_up != NULL);
        assert(data->dedn_dn != NULL);
        assert(data->dedsigma_up != NULL);
        assert(data->dedsigma_dn != NULL);
    } else {
        assert(0);
    }
}


void vdwxc_nullify_convenience_pointers(vdwxc_data data)
{
    // Set all the pointers to input/output parameters to zero.
    // We have not allocated any of these arrays so no freeing.
    // This is standard cleanup after a single calculation.
    data->rho = NULL;
    data->rho_up = NULL;
    data->rho_dn = NULL;
    data->sigma = NULL;
    data->sigma_up = NULL;
    data->sigma_dn = NULL;
    data->dedn = NULL;
    data->dedn_up = NULL;
    data->dedn_dn = NULL;
    data->dedsigma = NULL;
    data->dedsigma_up = NULL;
    data->dedsigma_dn = NULL;
}


double vdwxc_calculate_anyspin(vdwxc_data data)
{
    // Internal calculate function.  The convenience pointers must be allocated
    // appropriately
    if (data->initialized != VDWXC_INITIALIZED)
    {
        data->errorcode = VDWXC_ERROR_UNINITIALIZED;
        vdwxc_print_error(data);
        return 0.0;
    }
    if (data->cell.Nglobal[0] * data->cell.Nglobal[1] * data->cell.Nglobal[2] <= 0)
    {
        data->errorcode = VDWXC_ERROR_UNIT_CELL_NOT_SET;
        vdwxc_print_error(data);
        return 0.0;
    }

    vdwxc_check_convenience_pointers(data);

    double energy;
    switch(data->fft_type) {
    case VDWXC_FFTW_SERIAL:
        energy = vdwxc_calculate_serial(data);
        break;
    case VDWXC_FFTW_MPI:
        energy = vdwxc_calculate_mpi(data);
        break;
    case VDWXC_PFFT:
        energy = vdwxc_calculate_pfft(data);
        break;
    default:
        energy = 0.0;
        assert(0);
    }
#if VDWXC_DEBUG
    vdwxc_dump(data);
#endif
    return energy;
}


// Calculate van der Waals energy plus energy/density derivative and
// energy/gradient-density derivative from density and gradient.
double vdwxc_calculate(vdwxc_data data,
                       double *rho_g, double *sigma_g,
                       double *dedn_g, double *dedsigma_g)
{
    assert(data->nspin == 1);
    data->rho = rho_g;
    data->sigma = sigma_g;
    data->dedn = dedn_g;
    data->dedsigma = dedsigma_g;
    double energy = vdwxc_calculate_anyspin(data);
    vdwxc_nullify_convenience_pointers(data);
    return energy;
}

double vdwxc_calculate_spin(vdwxc_data data,
                            double *rho_up_g, double *rho_dn_g,
                            double *sigma_up_g, double *sigma_dn_g,
                            double *dedn_up_g, double *dedn_dn_g,
                            double *dedsigma_up_g, double *dedsigma_dn_g)
{
    assert(data->nspin == 2);
    data->rho_up = rho_up_g;
    data->rho_dn = rho_dn_g;
    data->sigma_up = sigma_up_g;
    data->sigma_dn = sigma_dn_g;
    data->dedn_up = dedn_up_g;
    data->dedn_dn = dedn_dn_g;
    data->dedsigma_up = dedsigma_up_g;
    data->dedsigma_dn = dedsigma_dn_g;
    double energy = vdwxc_calculate_anyspin(data);
    vdwxc_nullify_convenience_pointers(data);
    return energy;
}


void vdwxc_fatal(char* msg)
{
    printf("Fatal: %s", msg);
    assert(0); // TODO: Do things well
}

void vdwxc_finalize(vdwxc_data* vdw)
{
    // We have to free all the arrays that we have allocated and zero
    // the pointers.  We will zero the pointers in the end by setting
    // all the data to zero before freeing the vdwxc_data.  We
    // actually take a pointer *to* a vdwxc_data, so we can zero the
    // vdwxc_data pointer as well.  (Like MPI_Finalize.)
    vdwxc_data data = vdw[0];
#ifdef HAVE_PFFT
    if(data->fft_type == VDWXC_PFFT) {
        if (data->pfft_plan_r2c != NULL) {
            pfft_destroy_plan(data->pfft_plan_r2c);
        }
        if (data->pfft_plan_c2r != NULL) {
            pfft_destroy_plan(data->pfft_plan_c2r);
        }
        if (data->mpi_cart_comm != NULL) {
            MPI_Comm_free(&data->mpi_cart_comm);
        }
        if (data->work_ka != NULL) {
            pfft_free(data->work_ka);
            // Set to null so we won't call fftw_free on it in a moment
            data->work_ka = NULL;
        }
    }
#else
    assert(data->fft_type != VDWXC_PFFT);
#endif
    if (data->plan_r2c != NULL) {
        fftw_destroy_plan(data->plan_r2c);
    }
    if (data->plan_c2r != NULL) {
        fftw_destroy_plan(data->plan_c2r);
    }
    if (data->work_ka != NULL) {
        fftw_free(data->work_ka);
    }
    if (data->q0_g != NULL) {
        free(data->q0_g);
    }
    if (data->rho_dq0drho_sg != NULL) {
        free(data->rho_dq0drho_sg);
    }
    if (data->rho_dq0dsigma_sg != NULL) {
        free(data->rho_dq0dsigma_sg);
    }
    vdwxc_nullify(data);
    free(data);
    vdw[0] = NULL;
}