diff --git a/.gitignore b/.gitignore index 11ec0f59..c4e78f0b 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ controllers/options/options_driver.c controllers/defaults/DEFAULTS_driver.F services/controllers/SERVICES_setup.F services/units/UNITS_defaults.F +tests/ diff --git a/compile/compile_me.mk b/compile/compile_me.mk new file mode 100644 index 00000000..d209e66a --- /dev/null +++ b/compile/compile_me.mk @@ -0,0 +1,14 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2025 The Yambo Team +# +# Authors (see AUTHORS file for details): MN +# +# Co-Authors (see AUTHORS file for details): RR AM +# LetzElPhC ep plugin - yambo build system stub +# Integration handled via SERVICES_LIBS in libraries.mk +# +ifneq (,$(findstring yambo_ep,$(MAKECMDGOALS))) + Y_PRECMP= -D_Y6_LETZ +endif diff --git a/compile/goals.mk b/compile/goals.mk new file mode 100644 index 00000000..9fa43d29 --- /dev/null +++ b/compile/goals.mk @@ -0,0 +1,9 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2025 The Yambo Team +# +# Authors (see AUTHORS file for details): AM MN RR +# + +.PHONY: lelphc clean_lelphc diff --git a/compile/headers_and_libs.mk b/compile/headers_and_libs.mk new file mode 100644 index 00000000..b5d57158 --- /dev/null +++ b/compile/headers_and_libs.mk @@ -0,0 +1,15 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2025 The Yambo Team +# +# Authors (see AUTHORS file for details): RR AM +# +# LetzElPhC ep plugin - extra include paths for the vendored C source. +# C files under plugins/ep/services//*.c use relative includes +# like "common/dtypes.h", so the plugin services root must be on the +# include path. +# +iheaders += $(IFLAG)$(srcdir)/plugins/letz/services +include_headers += $(IFLAG)$(srcdir)/plugins/letz/services +llibs += -lcblas diff --git a/compile/help.mk b/compile/help.mk new file mode 100644 index 00000000..3a7ee157 --- /dev/null +++ b/compile/help.mk @@ -0,0 +1,16 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2025 The Yambo Team +# +# Authors (see AUTHORS file for details): RR AM +# +# Help blurb for the ep (LetzElPhC) plugin, called from +# compile/global/functions/help.mk when plugins/ep.pulled is present. +# +define ep_help + $(ECHO) " ep plugin = LetzElPhC electron-phonon coupling plugin";\ + $(ECHO) " Trigger: yambo -collisions ep";\ + $(ECHO) " Source: plugins/ep/services/, integration: packages/el-ph/";\ + $(ECHO) +endef diff --git a/compile/libraries.mk b/compile/libraries.mk index bc64ef5c..3e5115dd 100644 --- a/compile/libraries.mk +++ b/compile/libraries.mk @@ -9,8 +9,31 @@ ifeq ($(wildcard compile/global/defs.mk),compile/global/defs.mk) include compile/defs.mk endif # -EP_S1= [Sep,nonloc] [Sep,io/ezxml] [Sep,parser/inih] [Sep,io] -EP_S2= [Sep,io/qe] [Sep,common] [Sep,dvloc] [Sep,elph] [Sep,wfc] [Sep,symmetries] [Sep,fft] -EP_S3= [Sep,common/cwalk] [Sep,preprocessor] [Sep,interpolation] [Sep,common/ELPH_hash_map] [Sep,phonon] [Sep,common/kdtree] [Sep,parser] +# Library link order is INTRINSIC to LetzElPhC's dependency. +# Static linker is left-to-right: symbol-users before symbol-providers. # -SERVICES_LIBS += ${EP_S1} ${EP_S2} ${EP_S3} +# Dependency tiers (each tier depends only on equal/lower tiers): +# T6 drivers : interpolation, elph +# T5 QE I/O : io/qe +# T4 main I/O : io +# T3 kernels+pp : dvloc, symmetries, nonloc, preprocessor +# T2 primitives : fft, parser, phonon, wfc +# T1 base+vendors : common, common/kdtree, io/ezxml, common/ELPH_hash_map, +# parser/inih, common/cwalk +# +# elph (T5) hosts the Fortran-to-C ABI boundary (ep_f2c_bridge.c): +# elph_driver_f2c / elph_driver_cb_f2c / elph_driver_cb2_f2c +# These are the ONLY entry points yambo_ep calls into LetzElPhC. +# +# Order generated by check_dep.py topological sort (most-dependent first). +# +EP_S1= [Sletz,interpolation] +EP_S2= [Sletz,io/qe] [Sletz,common/ELPH_hash_map] [Sletz,common] [Sletz,symmetries] [Sletz,nonloc] +EP_S3= [Sletz,dvloc] [Sletz,io] [Sletz,io/ezxml] +EP_S4= [Sletz,wfc] [Sletz,fft] [Sletz,dvloc] [Sletz,common/cwalk] [Sletz,preprocessor] +EP_S5= [Sletz,phonon] [Sletz,parser/inih] [Sletz,parser] +EP_S6= [Sletz,common/kdtree] [Sletz,elph] +# +ifneq (,$(filter yambo_ep ypp_ep,$(MAKECMDGOALS))) + SERVICES_LIBS += ${EP_S1} ${EP_S2} ${EP_S3} ${EP_S4} ${EP_S5} ${EP_S6} +endif diff --git a/compile/targets.mk b/compile/targets.mk new file mode 100644 index 00000000..5868b9c9 --- /dev/null +++ b/compile/targets.mk @@ -0,0 +1,12 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2025 The Yambo Team +# +# Authors (see AUTHORS file for details): AM MN RR +# +lelphc: + $(MAKE) -C $(srcdir)/plugins/letz/services + +clean_lelphc: + -$(MAKE) -C $(srcdir)/plugins/letz/services clean diff --git a/info/compilation.info b/info/compilation.info new file mode 100644 index 00000000..f82a95c8 --- /dev/null +++ b/info/compilation.info @@ -0,0 +1,7 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2026 The Yambo Team +# +# Authors (see AUTHORS file for details): AM RR +# diff --git a/info/intro.info b/info/intro.info new file mode 100644 index 00000000..211a647b --- /dev/null +++ b/info/intro.info @@ -0,0 +1,8 @@ +# +# License-Identifier: GPL +# +# Copyright (C) 2026 The Yambo Team +# +# Authors (see AUTHORS file for details): AM RR +# + [letz-plugin] This plugin connects Yambo with the LetzElPhC code diff --git a/services/Makefile b/services/Makefile index 53edd6b4..2382396e 100644 --- a/services/Makefile +++ b/services/Makefile @@ -4,12 +4,15 @@ TARGET := lelphc #### make start here -SUBDIR = nonloc io/ezxml parser/inih io \ - io/qe common dvloc elph wfc \ - symmetries fft common/cwalk \ - preprocessor interpolation \ - common/ELPH_hash_map phonon \ - common/kdtree parser +# SUBDIR order mirrors the link-tier (T0 base → T5 drivers). +# For library use, compile/libraries.mk enforces the same tier ordering. +SUBDIR = common common/cwalk common/ELPH_hash_map common/kdtree \ + io/ezxml parser/inih \ + wfc fft phonon parser \ + dvloc nonloc symmetries preprocessor \ + io \ + io/qe \ + elph interpolation .SUFFIXES : diff --git a/services/common/dtypes.h b/services/common/dtypes.h index 4010569e..dfa7ed0b 100644 --- a/services/common/dtypes.h +++ b/services/common/dtypes.h @@ -56,6 +56,8 @@ struct calc_details // name of the input file. // can be elph input or DFT-Phonon input // base on calc_type + char ph_save_dir[512]; + // phonon save directory path for QE preprocessor }; struct symmetry @@ -141,7 +143,9 @@ struct Lattice // array of symmetries (operated on cartisian coordinates) bool is_soc_present; // system has soc if true else false - // ------------------------------------------------------- + // Yambo parallel workload + int *K_par; + int NK_par; }; struct Phonon @@ -178,6 +182,9 @@ struct Phonon // (natom,3,3,3) // Quadrapole tensor. for now not used is kept to NULL. // ------------------------------------------------------- + // Yambo parallel workload + int *Q_par; + int NQ_par; }; struct Vloc_table diff --git a/services/common/free_dtypes.c b/services/common/free_dtypes.c index 3b6fc66e..c12c0f56 100644 --- a/services/common/free_dtypes.c +++ b/services/common/free_dtypes.c @@ -5,7 +5,6 @@ #include "dtypes.h" #include "elphC.h" -#include "nonloc/fcoeff.h" void free_wfc_type(struct WFC* Wfc) { @@ -109,33 +108,4 @@ void free_local_pseudo_type(struct local_pseudo* lpseudo) lpseudo->rab_grid = NULL; return; -} - -void free_Pseudo_type(struct Pseudo* pseudo) -{ - if (!pseudo) - { - return; - } - - free_f_Coeff(pseudo); - free(pseudo->fCoeff); - free_Vloc_table_type(pseudo->vloc_table); - free(pseudo->Fsign); - free(pseudo->PP_table); - - if (pseudo->loc_pseudo) - { - for (ND_int itype = 0; itype < pseudo->ntype; ++itype) - { - free_local_pseudo_type(pseudo->loc_pseudo + itype); - } - } - free(pseudo->loc_pseudo); - - pseudo->loc_pseudo = NULL; - pseudo->PP_table = NULL; - pseudo->Fsign = NULL; - pseudo->fCoeff = NULL; - return; -} +} \ No newline at end of file diff --git a/services/common/free_dtypes.h b/services/common/free_dtypes.h index 3549574f..a1eab0a1 100644 --- a/services/common/free_dtypes.h +++ b/services/common/free_dtypes.h @@ -6,5 +6,4 @@ void free_lattice_type(struct Lattice* lattice); void free_phonon_type(struct Phonon* phonon); void free_Vloc_table_type(struct Vloc_table* vloc_tab); void free_local_pseudo_type(struct local_pseudo* lpseudo); -void free_Pseudo_type(struct Pseudo* pseudo); -void free_wfc_type(struct WFC* wfc); +void free_wfc_type(struct WFC* wfc); \ No newline at end of file diff --git a/services/common/parallel.c b/services/common/parallel.c index d4a62402..9306a76b 100644 --- a/services/common/parallel.c +++ b/services/common/parallel.c @@ -9,6 +9,7 @@ This file contains function which distributes cpus #include "dtypes.h" #include "elphC.h" #include "error.h" +#include /*get block size and starting idx of dimension that is distrbuted amoung cpus*/ ND_int get_mpi_local_size_idx(const ND_int n, ND_int* start_idx, MPI_Comm Comm) diff --git a/services/common/print_info.c b/services/common/print_info.c index bfdd0549..8e825a32 100644 --- a/services/common/print_info.c +++ b/services/common/print_info.c @@ -12,18 +12,25 @@ #include "dtypes.h" +static FILE* elph_log_fp = NULL; + +void elph_set_log_file(FILE* fp) { elph_log_fp = fp; } + +FILE* elph_get_log_file(void) { return elph_log_fp ? elph_log_fp : stdout; } + void print_info_msg(int mpi_rank, const char* fmt, ...) { if (mpi_rank) { return; } + FILE* out = elph_get_log_file(); va_list args; va_start(args, fmt); - vfprintf(stdout, fmt, args); + vfprintf(out, fmt, args); va_end(args); - fprintf(stdout, "\n"); - fflush(stdout); + fprintf(out, "\n"); + fflush(out); } void print_input_info(const char* save_dir, const char* ph_save_dir, diff --git a/services/common/print_info.h b/services/common/print_info.h index 85a25b24..804aaae6 100644 --- a/services/common/print_info.h +++ b/services/common/print_info.h @@ -8,6 +8,9 @@ void print_ELPH_logo(int mpi_rank, FILE* output); void print_info_msg(int mpi_rank, const char* fmt, ...); +void elph_set_log_file(FILE* fp); +FILE* elph_get_log_file(void); + void print_input_info(const char* save_dir, const char* ph_save_dir, const char* kernel, const bool kminusq, const enum ELPH_dft_code dft_code, diff --git a/services/elph.h b/services/elph.h new file mode 100644 index 00000000..30350216 --- /dev/null +++ b/services/elph.h @@ -0,0 +1,68 @@ +#pragma once +#include +#include + +#include "common/dtypes.h" +#include "elphC.h" + +/* + * Per-(iq_BZ, ik_BZ) fill callback for yambo COLL integration. + * Called once per BZ (q,k) pair on commK_rank==0 processes. + * iq_BZ, ik_BZ : 0-based BZ indices + * data : ELPH_cmplx buffer, C row-major (nmodes, nspin, nbnds, nbnds) + * nq..nb_start : full-BZ dimensions (constant across calls); nb_start is 1-based + */ +typedef void (*elph_fill_fn)(int iq_BZ, int ik_BZ, + const void* data, + int nq, int nk, int nmodes, int nspin, + int nbnds, int nb_start); + +/* + * Callback for dV_q^nu(G) in reciprocal space, called once per iBZ q-point + * from commK rank 0. + * iq_iBZ : 0-based global iBZ q-point index + * dVG : C row-major (nmodes, nmag, nfft_x, nfft_y, nfft_z) after forward FFT + * ph_freqs : phonon frequencies omega_ph[nmodes] in Hartree (ELPH_float*) + * nq_iBZ : total number of iBZ q-points + */ +typedef void (*elph_dvG_fill_fn)(int iq_iBZ, + const void* dVG, + const void* ph_freqs, + int nq_iBZ, int nmodes, int nmag, + int nfft_x, int nfft_y, int nfft_z); + +void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, + MPI_Comm comm_world); + +/* Callback-enabled variant: skips ndb.elph; calls fill_fn per (q,k) instead. */ +void elph_driver_cb(const char* ELPH_input_file, enum ELPH_dft_code dft_code, + MPI_Comm comm_world, elph_fill_fn fill_fn); + +/* + * Extended callback variant: like elph_driver_cb plus calls dvG_fill_fn once per + * iBZ q-point with the full dV_q^nu(G) potential in reciprocal space. + * Either callback may be NULL to skip that output. + */ +void elph_driver_cb2(const char* ELPH_input_file, enum ELPH_dft_code dft_code, + MPI_Comm comm_world, elph_fill_fn fill_fn, + elph_dvG_fill_fn dvG_fill_fn); + +void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, + struct Pseudo* pseudo, struct Phonon* phonon, + const ND_int iqpt, ELPH_cmplx* eigVec, + ELPH_cmplx* dVscfq, const int ncid_elph, + const int varid_elph, const int ncid_dmat, + const int varid_dmat, const bool non_loc, + const bool kminusq, + const struct ELPH_MPI_Comms* Comm, + elph_fill_fn fill_fn); + +void compute_and_write_dmats(const char* file_name, const struct WFC* wfcs, + const struct Lattice* lattice, + const ND_int nph_sym, + const struct symmetry* sym_data, + const struct ELPH_MPI_Comms* Comm); + +void init_kernel(struct kernel_info* kernel); + +void set_kernel(const char* kernel_str, struct kernel_info* kernel); diff --git a/services/elph/.objects b/services/elph/.objects index dae7f16b..6eadfd85 100644 --- a/services/elph/.objects +++ b/services/elph/.objects @@ -1 +1,2 @@ -objs = set_kernel.o elph_driver.o compute_elph_matq.o compute_dmats.o +objs = elph_driver_cb.o set_kernel.o elph_driver.o ep_f2c_bridge.o dvG_utils.o compute_elph_matq.o compute_dmats.o \ + elph_driver_cb2.o diff --git a/services/elph/compute_dmats.c b/services/elph/compute_dmats.c index 2880e8ee..608c27df 100644 --- a/services/elph/compute_dmats.c +++ b/services/elph/compute_dmats.c @@ -94,10 +94,16 @@ void compute_and_write_dmats(const char* file_name, const struct WFC* wfcs, struct progress_bar pbar[1]; start_progressbar(pbar, Comm->commW_rank, ndmats); - for (ND_int idmat = 0; idmat < ndmats; ++idmat) + for (ND_int i_par = 0; i_par < lattice->NK_par; ++i_par) { - ND_int isym = (idmat + dmat_shift) / nk_totalBZ; - ND_int ikBZ = (idmat + dmat_shift) % nk_totalBZ; + + ND_int ikBZ=lattice->K_par[i_par]; + + for (ND_int isym = 0; isym < nph_sym; ++isym) + { + + //ND_int isym = (idmat + dmat_shift) / nk_totalBZ; + //ND_int ikBZ = (idmat + dmat_shift) % nk_totalBZ; startp[0] = isym; startp[1] = ikBZ; @@ -116,8 +122,9 @@ void compute_and_write_dmats(const char* file_name, const struct WFC* wfcs, } } - // update the progress bar - print_progressbar(pbar); + } + // update the progress bar + print_progressbar(pbar); } if (Comm->commK_rank == 0) { diff --git a/services/elph/compute_elph_matq.c b/services/elph/compute_elph_matq.c index 82b72650..75749002 100644 --- a/services/elph/compute_elph_matq.c +++ b/services/elph/compute_elph_matq.c @@ -24,7 +24,9 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, const int varid_elph, const int ncid_dmat, const int varid_dmat, const bool non_loc, const bool kminusq, - const struct ELPH_MPI_Comms* Comm) + const struct ELPH_MPI_Comms* Comm, + elph_fill_fn fill_fn, + const ND_int iqpt_iBZ) { /* dVscf -> (nmodes,nmag,Nx,Ny,Nz) @@ -53,6 +55,8 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, distribute_to_grps(nk_totalBZ, Comm->nkpools, Comm->commQ_rank / Comm->commK_size, &kshift); + nk_this_pool=lattice->NK_par; + if (nk_this_pool < 1) { error_msg( @@ -113,10 +117,10 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, struct progress_bar pbar[1]; start_progressbar(pbar, Comm->commW_rank, nk_this_pool); // compute electron-phonon matrix elements for each kpoint - for (ND_int ii = 0; ii < nk_this_pool; ++ii) + for (ND_int ii = 0; ii < lattice->NK_par; ++ii) { /* compute the global k index */ - ND_int i = kshift + ii; + ND_int i = lattice->K_par[ii]; int idx_k = i; int idx_kq = KplusQidxs[i]; @@ -146,8 +150,18 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, int nc_err; if (Comm->commK_rank == 0) { - if ((nc_err = nc_put_vara(ncid_elph, varid_elph, startp, countp, - elph_kq_mn))) + //fprintf(stderr,"\n CPU %i k %i %i", Comm->commW_rank,i,startp[1]); + if (fill_fn != NULL) + { + fill_fn((int)startp[0], (int)startp[1], elph_kq_mn, + (int)phonon->nq_BZ, (int)nk_totalBZ, (int)nmodes, + (int)lattice->nspin, (int)nbnds, + (int)lattice->start_band, + (int)iqpt_iBZ, + (const void*)(phonon->qpts_BZ + 3 * startp[0])); + } + else if ((nc_err = nc_put_vara(ncid_elph, varid_elph, startp, + countp, elph_kq_mn))) { ERR(nc_err); } @@ -209,8 +223,17 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, startp[0] = qpos_star; startp[1] = idx_Sk; // Write it for Sq and Sk point - if ((nc_err = nc_put_vara(ncid_elph, varid_elph, startp, countp, - gSq_buff))) + if (fill_fn != NULL) + { + fill_fn((int)startp[0], (int)startp[1], gSq_buff, + (int)phonon->nq_BZ, (int)nk_totalBZ, (int)nmodes, + (int)lattice->nspin, (int)nbnds, + (int)lattice->start_band, + (int)iqpt_iBZ, + (const void*)(phonon->qpts_BZ + 3 * startp[0])); + } + else if ((nc_err = nc_put_vara(ncid_elph, varid_elph, startp, + countp, gSq_buff))) { ERR(nc_err); } diff --git a/services/elph/dvG_utils.c b/services/elph/dvG_utils.c new file mode 100644 index 00000000..1db548ce --- /dev/null +++ b/services/elph/dvG_utils.c @@ -0,0 +1,141 @@ +/* + * License-Identifier: GPL + * + * Copyright (C) 2025 The Yambo Team + * + * Authors (see AUTHORS file for details): RR AM + * + * Gather distributed dVscf z-slabs to commK rank 0 and apply a forward 3D + * FFT, producing dV_q^nu(G) ready for the elph_dvG_fill_fn callback. + */ + +#include "dvG_utils.h" +#include "fft/fft.h" +#include "common/error.h" +#include +#include + +ELPH_cmplx* gather_dVscf_and_fft(const ELPH_cmplx* dVscf, + const struct Lattice* lattice, + MPI_Comm commK) +{ + ND_int nmodes = lattice->nmodes; + ND_int nmag = lattice->nmag; + ND_int Nx = lattice->fft_dims[0]; + ND_int Ny = lattice->fft_dims[1]; + ND_int Nz = lattice->fft_dims[2]; + ND_int Nz_loc = lattice->nfftz_loc; + ND_int shift = lattice->nfftz_loc_shift; + ND_int Nxy = Nx * Ny; + ND_int NmodMag = nmodes * nmag; + + int commK_rank, commK_size; + MPI_Comm_rank(commK, &commK_rank); + MPI_Comm_size(commK, &commK_size); + + /* Collect per-rank Nz_loc and shifts. */ + int* all_Nz_loc = malloc(commK_size * sizeof(int)); + int* all_shifts = malloc(commK_size * sizeof(int)); + CHECK_ALLOC(all_Nz_loc); + CHECK_ALLOC(all_shifts); + int my_Nz_loc = (int)Nz_loc; + int my_shift = (int)shift; + MPI_Allgather(&my_Nz_loc, 1, MPI_INT, all_Nz_loc, 1, MPI_INT, commK); + MPI_Allgather(&my_shift, 1, MPI_INT, all_shifts, 1, MPI_INT, commK); + + /* + * Pack local slab: rearrange from + * dVscf[nmodes][nmag][Nx][Ny][Nz_loc] (C row-major) + * to + * packed[Nz_loc][nmodes][nmag][Nx][Ny] + * so that the z-extent is the leading dimension for contiguous Gatherv. + */ + ND_int local_count = Nz_loc * NmodMag * Nxy; + ELPH_cmplx* packed = malloc(local_count * sizeof(ELPH_cmplx)); + CHECK_ALLOC(packed); + + for (ND_int iv = 0; iv < nmodes; ++iv) + for (ND_int is = 0; is < nmag; ++is) + for (ND_int ix = 0; ix < Nx; ++ix) + for (ND_int iy = 0; iy < Ny; ++iy) + for (ND_int iz = 0; iz < Nz_loc; ++iz) + { + packed[iz * NmodMag*Nxy + iv * nmag*Nxy + is * Nxy + ix * Ny + iy] = + dVscf[iv * nmag*Nxy*Nz_loc + is * Nxy*Nz_loc + + ix * Ny*Nz_loc + iy * Nz_loc + iz]; + } + + /* Build Gatherv send/recv counts and displacements (only on rank 0). */ + int* sendcounts = NULL; + int* displs = NULL; + ELPH_cmplx* gathered = NULL; + ND_int NmodMagNxy = NmodMag * Nxy; + + if (commK_rank == 0) + { + sendcounts = malloc(commK_size * sizeof(int)); + displs = malloc(commK_size * sizeof(int)); + CHECK_ALLOC(sendcounts); + CHECK_ALLOC(displs); + for (int r = 0; r < commK_size; ++r) + { + sendcounts[r] = all_Nz_loc[r] * (int)NmodMagNxy; + displs[r] = all_shifts[r] * (int)NmodMagNxy; + } + gathered = malloc((size_t)Nz * NmodMag * Nxy * sizeof(ELPH_cmplx)); + CHECK_ALLOC(gathered); + } + + MPI_Gatherv(packed, (int)local_count, ELPH_MPI_cmplx, + gathered, sendcounts, displs, ELPH_MPI_cmplx, + 0, commK); + + free(packed); + free(all_Nz_loc); + free(all_shifts); + if (commK_rank == 0) + { + free(sendcounts); + free(displs); + } + + if (commK_rank != 0) + return NULL; + + /* + * Transpose gathered[Nz][nmodes][nmag][Nx][Ny] back to + * dVG[nmodes][nmag][Nx][Ny][Nz]. + */ + ELPH_cmplx* dVG = malloc((size_t)NmodMag * Nxy * Nz * sizeof(ELPH_cmplx)); + CHECK_ALLOC(dVG); + + for (ND_int iz = 0; iz < Nz; ++iz) + for (ND_int iv = 0; iv < nmodes; ++iv) + for (ND_int is = 0; is < nmag; ++is) + for (ND_int ix = 0; ix < Nx; ++ix) + for (ND_int iy = 0; iy < Ny; ++iy) + { + dVG[iv * nmag*Nxy*Nz + is * Nxy*Nz + ix * Ny*Nz + iy * Nz + iz] = + gathered[iz * NmodMagNxy + iv * nmag*Nxy + is * Nxy + ix * Ny + iy]; + } + free(gathered); + + /* + * Forward 3D FFT for each (mode, spin) slice. + * fftw plan_dft_3d uses row-major [n0][n1][n2] = [Nx][Ny][Nz]. + */ + ND_int slice_sz = Nxy * Nz; + for (ND_int iv = 0; iv < nmodes; ++iv) + for (ND_int is = 0; is < nmag; ++is) + { + ELPH_cmplx* sl = dVG + (iv * nmag + is) * slice_sz; + fftw_generic_plan plan = fftw_fun(plan_dft_3d)( + (int)Nx, (int)Ny, (int)Nz, + sl, sl, + FFTW_FORWARD, FFTW_ESTIMATE); + fftw_fun(execute)(plan); + fftw_fun(destroy_plan)(plan); + } + + return dVG; +} diff --git a/services/elph/dvG_utils.h b/services/elph/dvG_utils.h new file mode 100644 index 00000000..fc4029e4 --- /dev/null +++ b/services/elph/dvG_utils.h @@ -0,0 +1,31 @@ +#pragma once +/* + * License-Identifier: GPL + * + * Copyright (C) 2025 The Yambo Team + * + * Authors (see AUTHORS file for details): RR AM + * + * Helpers for gathering the distributed dVscf real-space potential and + * transforming it to reciprocal space so that the elph_dvG_fill_fn callback + * can be invoked from commK rank 0. + */ + +#include "common/dtypes.h" +#include "elph.h" +#include + +/* + * Gather z-slabs of dVscf from all commK ranks onto rank 0, perform a + * forward 3D FFT for each (mode, spin) slice, and return the result. + * + * dVscf layout (C row-major, distributed): + * (nmodes, nmag, fft_dims[0], fft_dims[1], nfftz_loc) + * + * Return value (only meaningful on commK rank 0; NULL on other ranks): + * malloc'd buffer of shape (nmodes, nmag, Nx, Ny, Nz) in G-space (complex). + * Caller must free(). + */ +ELPH_cmplx* gather_dVscf_and_fft(const ELPH_cmplx* dVscf, + const struct Lattice* lattice, + MPI_Comm commK); diff --git a/services/elph/elph.h b/services/elph/elph.h index d0513c5d..622d5137 100644 --- a/services/elph/elph.h +++ b/services/elph/elph.h @@ -2,12 +2,58 @@ #include #include +#include "parser/parser.h" #include "common/dtypes.h" #include "elphC.h" +#include "yambo.h" + +/* + * Per-(iq_BZ, ik_BZ) fill callback for yambo COLL integration. + * Called once per BZ (q,k) pair on commK_rank==0 processes. + * iq_BZ, ik_BZ : 0-based BZ indices + * data : ELPH_cmplx buffer, C row-major (nmodes, nspin, nbnds, nbnds) + * nq..nb_start : full-BZ dimensions (constant across calls); nb_start is 1-based + * iq_iBZ : 0-based iBZ q-point index for this iq_BZ star + * qpt_BZ_crys : crystal-coordinate q-point for iq_BZ (ELPH_float[3]) + */ +typedef void (*elph_fill_fn)(int iq_BZ, int ik_BZ, + const void* data, + int nq, int nk, int nmodes, int nspin, + int nbnds, int nb_start, + int iq_iBZ, const void* qpt_BZ_crys); + +/* + * Callback for dV_q^nu(G) in reciprocal space, called once per iBZ q-point + * from commK rank 0. + * iq_iBZ : 0-based global iBZ q-point index + * dVG : C row-major (nmodes, nmag, nfft_x, nfft_y, nfft_z) after forward FFT + * ph_freqs : phonon frequencies omega_ph[nmodes] in Hartree (ELPH_float*) + * nq_iBZ : total number of iBZ q-points + */ +typedef void (*elph_dvG_fill_fn)(int iq_iBZ, + const void* dVG, + const void* ph_freqs, + int nq_iBZ, int nmodes, int nmag, + int nfft_x, int nfft_y, int nfft_z); void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, MPI_Comm comm_world); +/* Callback-enabled variant: skips ndb.elph; calls fill_fn per (q,k) instead. */ +void elph_driver_cb(const char* ELPH_input_file, enum ELPH_dft_code dft_code, + MPI_Comm comm_world, elph_fill_fn fill_fn); + +/* + * Extended callback variant: like elph_driver_cb plus calls dvG_fill_fn once per + * iBZ q-point with the full dV_q^nu(G) potential in reciprocal space. + * Either callback may be NULL to skip that output. + * comm_q, comm_k: Y6 PAR communicators for q,k distribution. + */ +void elph_driver_cb2(struct elph_usr_input* input_data, struct Y6_info* y6_data, struct Y6_parallel_work* y6_work, enum ELPH_dft_code dft_code, + elph_fill_fn fill_fn, + elph_dvG_fill_fn dvG_fill_fn,int i_control, + MPI_Comm comm_world ); + void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, struct Pseudo* pseudo, struct Phonon* phonon, const ND_int iqpt, ELPH_cmplx* eigVec, @@ -15,7 +61,9 @@ void compute_and_write_elphq(struct WFC* wfcs, struct Lattice* lattice, const int varid_elph, const int ncid_dmat, const int varid_dmat, const bool non_loc, const bool kminusq, - const struct ELPH_MPI_Comms* Comm); + const struct ELPH_MPI_Comms* Comm, + elph_fill_fn fill_fn, + const ND_int iqpt_iBZ); void compute_and_write_dmats(const char* file_name, const struct WFC* wfcs, const struct Lattice* lattice, diff --git a/services/elph/elph_driver.c b/services/elph/elph_driver.c index e09baeff..fb82cac8 100644 --- a/services/elph/elph_driver.c +++ b/services/elph/elph_driver.c @@ -17,6 +17,7 @@ THe starting point for the entire code #include "common/parallel.h" #include "common/print_info.h" #include "dvloc/dvloc.h" +#include "dvG_utils.h" #include "elph.h" #include "elphC.h" #include "fft/fft.h" @@ -24,6 +25,7 @@ THe starting point for the entire code #include "io/qe/qe_io.h" #include "parser/parser.h" #include "symmetries/symmetries.h" +#include "common/string_func.h" void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, MPI_Comm comm_world) @@ -51,7 +53,7 @@ void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, mpi_comms); // print logo and stated message - print_ELPH_logo(mpi_comms->commW_rank, stdout); + print_ELPH_logo(mpi_comms->commW_rank, elph_get_log_file()); print_info_msg(mpi_comms->commW_rank, "********** Program started **********"); print_input_info(input_data->save_dir, input_data->ph_save_dir, @@ -73,10 +75,18 @@ void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, struct WFC* wfcs; // read the SAVE data and phonon related data. + if (dft_code == DFT_CODE_QE) + { + get_data_from_qe(lattice, phonon, pseudo, input_data->ph_save_dir, NULL, + mpi_comms); + } + else + { + error_msg("Only QE supported"); + } read_and_alloc_save_data(input_data->save_dir, mpi_comms, input_data->start_bnd, input_data->end_bnd, &wfcs, - input_data->ph_save_dir, lattice, pseudo, phonon, - dft_code); + input_data->ph_save_dir, lattice, pseudo, phonon); // print info about lattice and phonons print_lattice_info(mpi_comms, lattice); @@ -303,7 +313,7 @@ void elph_driver(const char* ELPH_input_file, enum ELPH_dft_code dft_code, compute_and_write_elphq(wfcs, lattice, pseudo, phonon, iqpt_iBZg, eigVec, dVscf, ncid_elph, varid_elph, ncid_dmat, varid_dmat, kernel->non_loc, - input_data->kminusq, mpi_comms); + input_data->kminusq, mpi_comms, NULL, iqpt_iBZg); } free(eig_Sq); diff --git a/services/elph/elph_driver_cb.c b/services/elph/elph_driver_cb.c new file mode 100644 index 00000000..a28f71c0 --- /dev/null +++ b/services/elph/elph_driver_cb.c @@ -0,0 +1,226 @@ +/* + Yambo calling point +*/ +#include +#include +#include +#include +#include +#include +#include +#include + +#include "common/ELPH_timers.h" +#include "common/dtypes.h" +#include "common/error.h" +#include "common/init_dtypes.h" +#include "common/parallel.h" +#include "common/print_info.h" +#include "dvloc/dvloc.h" +#include "dvG_utils.h" +#include "elph.h" +#include "elphC.h" +#include "fft/fft.h" +#include "io/io.h" +#include "io/qe/qe_io.h" +#include "parser/parser.h" +#include "symmetries/symmetries.h" +#include "common/string_func.h" + +/* + * Callback-enabled driver: identical flow to elph_driver but writes gkkp + * data via fill_fn(iq_BZ, ik_BZ, data, ...) instead of ndb.elph. + * ndb.Dmats is still created and used internally for BZ symmetry expansion. + * Eigenvectors / frequencies are NOT written (yambo reads them elsewhere). + */ +void elph_driver_cb(const char* ELPH_input_file, enum ELPH_dft_code dft_code, + MPI_Comm comm_world, elph_fill_fn fill_fn) +{ + struct elph_usr_input* input_data; + init_ELPH_clocks(); + + read_elph_input_file(ELPH_input_file, &input_data, comm_world); + + struct kernel_info* kernel = malloc(sizeof(struct kernel_info)); + init_kernel(kernel); + set_kernel(input_data->kernel_str, kernel); + + struct ELPH_MPI_Comms* mpi_comms = malloc(sizeof(struct ELPH_MPI_Comms)); + CHECK_ALLOC(mpi_comms); + + create_parallel_comms(input_data->nqpool, input_data->nkpool, comm_world, + mpi_comms); + + print_ELPH_logo(mpi_comms->commW_rank, elph_get_log_file()); + print_info_msg(mpi_comms->commW_rank, + "********** Program started **********"); + print_input_info(input_data->save_dir, input_data->ph_save_dir, + input_data->kernel_str, input_data->kminusq, dft_code, + mpi_comms); + + struct Lattice* lattice = malloc(sizeof(struct Lattice)); + CHECK_ALLOC(lattice); + init_lattice_type(lattice); + + struct Pseudo* pseudo = malloc(sizeof(struct Pseudo)); + CHECK_ALLOC(pseudo); + init_Pseudo_type(pseudo); + + struct Phonon* phonon = malloc(sizeof(struct Phonon)); + CHECK_ALLOC(phonon); + init_phonon_type(phonon); + + struct WFC* wfcs; + + if (dft_code == DFT_CODE_QE) + { + get_data_from_qe(lattice, phonon, pseudo, input_data->ph_save_dir, NULL, + mpi_comms); + } + else + { + error_msg("Only QE supported"); + } + read_and_alloc_save_data(input_data->save_dir, mpi_comms, + input_data->start_bnd, input_data->end_bnd, &wfcs, + input_data->ph_save_dir, lattice, pseudo, phonon); + char DM_name[100]; + strlcpy_custom(DM_name, input_data->ph_save_dir,100); + strlcat(DM_name,"/ndb.Dmats", 100); + + print_lattice_info(mpi_comms, lattice); + print_phonon_info(mpi_comms, phonon); + + print_info_msg(mpi_comms->commW_rank, ""); + print_info_msg(mpi_comms->commW_rank, + "=== Computing Dmats for phonon symmetries ==="); + compute_and_write_dmats(DM_name, wfcs, lattice, phonon->nph_sym, + phonon->ph_syms, mpi_comms); + + ND_int nmodes = lattice->nmodes; + ND_int nfft_loc = + lattice->fft_dims[0] * lattice->fft_dims[1] * lattice->nfftz_loc; + + ELPH_cmplx* eigVec = malloc(sizeof(ELPH_cmplx) * nmodes * nmodes); + CHECK_ALLOC(eigVec); + + ELPH_cmplx* dVscf = + malloc(sizeof(ELPH_cmplx) * nmodes * lattice->nmag * nfft_loc); + CHECK_ALLOC(dVscf); + + ELPH_float* omega_ph = malloc(sizeof(ELPH_float) * nmodes); + CHECK_ALLOC(omega_ph); + + int ncid_dmat, nc_err; + int varid_dmat; + + if (mpi_comms->commK_rank == 0) + { + if ((nc_err = nc_open_par(DM_name, NC_NOWRITE, mpi_comms->commR, + MPI_INFO_NULL, &ncid_dmat))) + { + ERR(nc_err); + } + + if ((nc_err = nc_inq_varid(ncid_dmat, "Dmats", &varid_dmat))) + { + ERR(nc_err); + } + + size_t Dmat_counts[6] = {0, 0, 0, 0, 0, 0}; + if ((nc_err = nc_get_vara(ncid_dmat, varid_dmat, Dmat_counts, + Dmat_counts, NULL))) + { + ERR(nc_err); + } + } + + print_info_msg(mpi_comms->commW_rank, + "=== Computing Electron-phonon matrix elements ==="); + print_info_msg(mpi_comms->commW_rank, ""); + + for (ND_int iqpt = 0; iqpt < phonon->nq_iBZ_loc; ++iqpt) + { + print_info_msg(mpi_comms->commW_rank, "### q-point : %d/%d", + (int)(iqpt + 1), (int)phonon->nq_iBZ_loc); + + ND_int iqpt_iBZg = iqpt + phonon->nq_shift; + + if (dft_code == DFT_CODE_QE) + { + ELPH_cmplx* dvscf_read = NULL; + if (kernel->screening == ELPH_DFPT_SCREENING) + { + dvscf_read = dVscf; + } + else + { + ND_int dvscf_num = nmodes * lattice->nmag * nfft_loc; + for (ND_int ix = 0; ix < dvscf_num; ++ix) + { + dVscf[ix] = 0.0; + } + } + get_dvscf_dyn_qe(input_data->ph_save_dir, lattice, iqpt_iBZg, + eigVec, dvscf_read, omega_ph, mpi_comms); + } + else + { + error_msg("Currently only quantum espresso supported"); + } + + ELPH_cmplx* Vlocr = malloc(sizeof(ELPH_cmplx) * nmodes * nfft_loc); + CHECK_ALLOC(Vlocr); + dVlocq(phonon->qpts_iBZ + iqpt_iBZg * 3, lattice, pseudo, eigVec, Vlocr, + mpi_comms->commK); + + if (dft_code == DFT_CODE_QE) + { + add_dvscf_qe(dVscf, Vlocr, lattice); + } + else + { + error_msg("Currently only quantum espresso supported"); + } + free(Vlocr); + + compute_and_write_elphq(wfcs, lattice, pseudo, phonon, iqpt_iBZg, + eigVec, dVscf, 0, 0, ncid_dmat, + varid_dmat, kernel->non_loc, + input_data->kminusq, mpi_comms, fill_fn, + iqpt_iBZg); + } + + if (mpi_comms->commK_rank == 0) + { + if ((nc_err = nc_close(ncid_dmat))) + { + ERR(nc_err); + } + } + + free(omega_ph); + free(eigVec); + free(dVscf); + + int World_rank_tmp = mpi_comms->commW_rank; + + free(kernel); + free_elph_usr_input(input_data); + free_save_data(wfcs, lattice, pseudo, phonon); + free(lattice); + free(pseudo); + free(phonon); + free_parallel_comms(mpi_comms); + free(mpi_comms); + fftw_fun(cleanup)(); + + if (0 == World_rank_tmp) + { + print_ELPH_clock_summary(); + } + cleanup_ELPH_clocks(); + print_info_msg(World_rank_tmp, "********** Program ended **********"); +} + + diff --git a/services/elph/elph_driver_cb2.c b/services/elph/elph_driver_cb2.c new file mode 100644 index 00000000..8503a26b --- /dev/null +++ b/services/elph/elph_driver_cb2.c @@ -0,0 +1,346 @@ +/* + Yambo calling point +*/ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "common/ELPH_timers.h" +#include "common/constants.h" +#include "common/cwalk/cwalk.h" +#include "common/dtypes.h" +#include "common/error.h" +#include "common/init_dtypes.h" +#include "common/numerical_func.h" +#include "common/parallel.h" +#include "common/print_info.h" +#include "dvloc/dvloc.h" +#include "dvG_utils.h" +#include "yambo.h" +#include "elph.h" +#include "elphC.h" +#include "fft/fft.h" +#include "io/io.h" +#include "io/nc_utils.h" +#include "io/qe/qe_io.h" +#include "parser/parser.h" +#include "symmetries/symmetries.h" +#include "common/string_func.h" + +/* + * Extended callback variant: like elph_driver_cb plus dvG_fill_fn called + * once per iBZ q-point from commK rank 0 with dV_q^nu(G) in G-space. + * Either callback may be NULL to skip that output. + * comm_q, comm_k: Y6 PAR communicators for q,k distribution (nqpool/nkpool derived from these). + */ +void elph_driver_cb2(struct elph_usr_input* input_data,struct Y6_info* y6_data, + struct Y6_parallel_work* y6_work, enum ELPH_dft_code dft_code, + elph_fill_fn fill_fn, + elph_dvG_fill_fn dvG_fill_fn,int i_control, + MPI_Comm comm_world) +{ + /*struct elph_usr_input* input_data;*/ + init_ELPH_clocks(); + + struct kernel_info* kernel = malloc(sizeof(struct kernel_info)); + init_kernel(kernel); + set_kernel(input_data->kernel_str, kernel); + + struct ELPH_MPI_Comms* mpi_comms = malloc(sizeof(struct ELPH_MPI_Comms)); + CHECK_ALLOC(mpi_comms); + +#ifdef _Y6_LETZ + /* Y6 mode: comm_world is actually comm_q from yambo PAR scheme. + Create full MPI structure from scratch (will be reconciled with Y6 later). + For now, use Y6-provided communicators for proper coordination. */ + create_parallel_comms(input_data->nqpool, input_data->nkpool, comm_world, mpi_comms); +#else + /* Standalone mode: create 2D MPI structure from comm_world */ + create_parallel_comms(input_data->nqpool, input_data->nkpool, comm_world, mpi_comms); +#endif + + + if (i_control > 0 ) + { + print_ELPH_logo(mpi_comms->commW_rank, elph_get_log_file()); + print_info_msg(mpi_comms->commW_rank, + "********** Program started **********"); + print_input_info(input_data->save_dir, input_data->ph_save_dir, + input_data->kernel_str, input_data->kminusq, dft_code, + mpi_comms); + } + + struct Lattice* lattice = malloc(sizeof(struct Lattice)); + CHECK_ALLOC(lattice); + init_lattice_type(lattice); + + lattice->NK_par=y6_work->NK; + lattice->K_par = malloc(lattice->NK_par * sizeof(int)); + memcpy(lattice->K_par, y6_work->K, lattice->NK_par * sizeof(int)); + + struct Pseudo* pseudo = malloc(sizeof(struct Pseudo)); + CHECK_ALLOC(pseudo); + init_Pseudo_type(pseudo); + + struct Phonon* phonon = malloc(sizeof(struct Phonon)); + CHECK_ALLOC(phonon); + init_phonon_type(phonon); + + phonon->NQ_par=y6_work->NQ; + phonon->Q_par = malloc(phonon->NQ_par * sizeof(int)); + memcpy(phonon->Q_par, y6_work->Q, phonon->NQ_par * sizeof(int)); + + /* + fprintf(stderr,"\n"); + for (int i = 0; i < y6_work->NK; i++) { + fprintf(stderr," ID %i K %i\n ",mpi_comms->commW_rank, y6_work->K[i]); + } + fprintf(stderr,"\n"); + for (int i = 0; i < y6_work->NQ; i++) { + fprintf(stderr," R %i Q %i\n ",mpi_comms->commW_rank, y6_work->Q[i]); + } + fprintf(stderr,"\n"); + */ + + if (dft_code == DFT_CODE_QE) + { + get_data_from_qe(lattice, phonon, pseudo, input_data->ph_save_dir, NULL, + mpi_comms); + if (mpi_comms->commW_rank == 0) { + fprintf(stderr, "[ELPH DEBUG] LetzElPhC read nq_iBZ=%d, nq_BZ=%d from dyn files\n", + (int)phonon->nq_iBZ, (int)phonon->nq_BZ); + } + } + else + { + error_msg("Only QE supported"); + } + + /* Populate y6_data (elphC_info in Fortran) with lattice/phonon metadata. + This happens early so query-only mode (i_control < 0) can return with + nmodes, natom, nsym, etc. already set for Fortran-side allocation. */ + + y6_data->natom=lattice->natom; + y6_data->nsym=lattice->nsym; + y6_data->timerev=lattice->timerev; + y6_data->nspin=lattice->nspin; + y6_data->nspinor=lattice->nspinor; + y6_data->total_bands=lattice->total_bands; + y6_data->start_band=lattice->start_band; + y6_data->end_band=lattice->end_band; + y6_data->nbnds=lattice->nbnds; + y6_data->lattice_dim=lattice->dimension; + y6_data->nmag=lattice->nmag; + y6_data->nmodes=lattice->nmodes; + y6_data->nq_ibz=phonon->nq_iBZ; + y6_data->nq_bz=phonon->nq_BZ; + + /* Load k-points and other data needed for both query and computation modes */ + struct WFC* wfcs = NULL; + read_and_alloc_save_data(input_data->save_dir, mpi_comms, + input_data->start_bnd, input_data->end_bnd, &wfcs, + input_data->ph_save_dir, lattice, pseudo, phonon); + + /* Populate k/q metadata for Yambo */ + y6_data->nkpts_ibz = lattice->nkpts_iBZ; + y6_data->nkpts_bz = lattice->nkpts_BZ; + y6_data->nph_sym = phonon->nph_sym; + y6_data->kminusq = (int)input_data->kminusq; + + /* Point to internal Letz k/q arrays */ + y6_data->kpt_fullBZ_crys = lattice->kpt_fullBZ_crys; + y6_data->kmap = lattice->kmap; + y6_data->qpts_iBZ = phonon->qpts_iBZ; + y6_data->qpts_BZ = phonon->qpts_BZ; + y6_data->qmap = phonon->qmap; + y6_data->nqstar = (int*)phonon->nqstar; + + /* Compute k+q indices for each q-point and expose to Yambo */ + int* kplusq_all = malloc(phonon->nq_BZ * lattice->nkpts_BZ * sizeof(int)); + CHECK_ALLOC(kplusq_all); + + for (ND_int iq_bz = 0; iq_bz < phonon->nq_BZ; ++iq_bz) { + ND_int iq_ibz = phonon->qmap[2 * iq_bz]; + ELPH_float* qpt = phonon->qpts_iBZ + iq_ibz * 3; + + ELPH_float qpt_tmp[3]; + memcpy(qpt_tmp, qpt, 3 * sizeof(ELPH_float)); + if (input_data->kminusq) { + for (int xi = 0; xi < 3; ++xi) qpt_tmp[xi] = -qpt_tmp[xi]; + } + + int* kplusq_q = kplusq_all + iq_bz * lattice->nkpts_BZ; + get_KplusQ_idxs(lattice->nkpts_BZ, lattice->kpt_fullBZ_crys, qpt_tmp, kplusq_q); + } + y6_data->kplusq_idxs = kplusq_all; + + /* Query-only mode: return early without computing gkkp */ + if (i_control < 0) { + return; + } + + char DM_name[100]; + strlcpy_custom(DM_name, input_data->ph_save_dir,100); + strlcat(DM_name,"/ndb.Dmats", 100); + + print_lattice_info(mpi_comms, lattice); + print_phonon_info(mpi_comms, phonon); + + print_info_msg(mpi_comms->commW_rank, ""); + print_info_msg(mpi_comms->commW_rank, + "=== Computing Dmats for phonon symmetries ==="); + compute_and_write_dmats(DM_name, wfcs, lattice, phonon->nph_sym, + phonon->ph_syms, mpi_comms); + + ND_int nmodes = lattice->nmodes; + ND_int nfft_loc = + lattice->fft_dims[0] * lattice->fft_dims[1] * lattice->nfftz_loc; + + ELPH_cmplx* eigVec = malloc(sizeof(ELPH_cmplx) * nmodes * nmodes); + CHECK_ALLOC(eigVec); + + ELPH_cmplx* dVscf = + malloc(sizeof(ELPH_cmplx) * nmodes * lattice->nmag * nfft_loc); + CHECK_ALLOC(dVscf); + + ELPH_float* omega_ph = malloc(sizeof(ELPH_float) * nmodes); + CHECK_ALLOC(omega_ph); + + int ncid_dmat, nc_err; + int varid_dmat; + + if (mpi_comms->commK_rank == 0) + { + if ((nc_err = nc_open_par(DM_name, NC_NOWRITE, mpi_comms->commR, + MPI_INFO_NULL, &ncid_dmat))) + { + ERR(nc_err); + } + + if ((nc_err = nc_inq_varid(ncid_dmat, "Dmats", &varid_dmat))) + { + ERR(nc_err); + } + + size_t Dmat_counts[6] = {0, 0, 0, 0, 0, 0}; + if ((nc_err = nc_get_vara(ncid_dmat, varid_dmat, Dmat_counts, + Dmat_counts, NULL))) + { + ERR(nc_err); + } + } + + print_info_msg(mpi_comms->commW_rank, + "=== Computing Electron-phonon matrix elements ==="); + print_info_msg(mpi_comms->commW_rank, ""); + + for (ND_int iqpt = 0; iqpt < phonon->NQ_par; ++iqpt) + { + print_info_msg(mpi_comms->commW_rank, "### q-point : %d/%d", + (int)(iqpt + 1), (int)phonon->nq_iBZ_loc); + + ND_int iqpt_iBZg = phonon->Q_par[iqpt]; + + if (dft_code == DFT_CODE_QE) + { + ELPH_cmplx* dvscf_read = NULL; + if (kernel->screening == ELPH_DFPT_SCREENING) + { + dvscf_read = dVscf; + } + else + { + ND_int dvscf_num = nmodes * lattice->nmag * nfft_loc; + for (ND_int ix = 0; ix < dvscf_num; ++ix) + { + dVscf[ix] = 0.0; + } + } + get_dvscf_dyn_qe(input_data->ph_save_dir, lattice, iqpt_iBZg, + eigVec, dvscf_read, omega_ph, mpi_comms); + } + else + { + error_msg("Currently only quantum espresso supported"); + } + + ELPH_cmplx* Vlocr = malloc(sizeof(ELPH_cmplx) * nmodes * nfft_loc); + CHECK_ALLOC(Vlocr); + dVlocq(phonon->qpts_iBZ + iqpt_iBZg * 3, lattice, pseudo, eigVec, Vlocr, + mpi_comms->commK); + + if (dft_code == DFT_CODE_QE) + { + add_dvscf_qe(dVscf, Vlocr, lattice); + } + else + { + error_msg("Currently only quantum espresso supported"); + } + free(Vlocr); + + ND_int calc_dvG=-1; + if (dvG_fill_fn != NULL && i_control>=2 ) { calc_dvG=1; } + + /* dvG callback: gather dVscf z-slabs to rank 0, FFT, call handler. */ + if (calc_dvG ==1 ) + { + ELPH_cmplx* dVG = gather_dVscf_and_fft(dVscf, lattice, + mpi_comms->commK); + if (mpi_comms->commK_rank == 0) + { + dvG_fill_fn((int)iqpt_iBZg, (const void*)dVG, + (const void*)omega_ph, + (int)phonon->nq_iBZ, + (int)nmodes, (int)lattice->nmag, + (int)lattice->fft_dims[0], + (int)lattice->fft_dims[1], + (int)lattice->fft_dims[2]); + free(dVG); + } + } + + compute_and_write_elphq(wfcs, lattice, pseudo, phonon, iqpt_iBZg, + eigVec, dVscf, 0, 0, ncid_dmat, + varid_dmat, kernel->non_loc, + input_data->kminusq, mpi_comms, fill_fn, + iqpt_iBZg); + } + + if (mpi_comms->commK_rank == 0) + { + if ((nc_err = nc_close(ncid_dmat))) + { + ERR(nc_err); + } + } + free(omega_ph); + free(eigVec); + free(dVscf); + + int World_rank_tmp = mpi_comms->commW_rank; + + free(kernel); + /*free_elph_usr_input(input_data);*/ + free_save_data(wfcs, lattice, pseudo, phonon); + free(lattice); + free(pseudo); + free(phonon); + //free_parallel_comms(mpi_comms); + //free(mpi_comms); + fftw_fun(cleanup)(); + + if (0 == World_rank_tmp) + { + print_ELPH_clock_summary(); + } + cleanup_ELPH_clocks(); + print_info_msg(World_rank_tmp, "********** Program ended **********"); +} + diff --git a/services/elph/ep_f2c_bridge.c b/services/elph/ep_f2c_bridge.c new file mode 100644 index 00000000..80cafebf --- /dev/null +++ b/services/elph/ep_f2c_bridge.c @@ -0,0 +1,126 @@ +/* + * License-Identifier: GPL + * + * Copyright (C) 2025 The Yambo Team + * + * Authors (see AUTHORS file for details): RR AM + * + * Thin bridge: called from Fortran via ISO_C_BINDING. + * Converts a Fortran MPI communicator handle (MPI_Fint) to a C MPI_Comm + * and forwards the call to LetzElPhC's elph_driver(). + * + */ + +#include "elph.h" /* already pulls in common/dtypes.h and elphC.h */ +#include "common/print_info.h" +#include +#include +#include +#include +#include +#include +#include "parser/parser.h" + +/* Fortran bind(C) callbacks — linked from the ep packages library. */ +extern void elph_coll_fill_gkkp(int, int, const void*, int, int, int, int, int, int, int, const void*); +extern void elph_coll_fill_dvg(int, const void*, const void*, int, int, int, int, int, int); + +/* + * Redirect stdout to the log file for the duration of the LetzElPhC call. + */ +static int saved_stdout_fd = -1; + +static void open_letz_log(MPI_Comm comm, const char* log_path) +{ + int rank; + MPI_Comm_rank(comm, &rank); + if (rank == 0 && log_path != NULL && log_path[0] != '\0') + { + FILE* fp = fopen(log_path, "w"); + if (fp) + { + elph_set_log_file(fp); + fflush(stdout); + saved_stdout_fd = dup(STDOUT_FILENO); + dup2(fileno(fp), STDOUT_FILENO); + } + } +} + +static void close_letz_log(void) +{ + FILE* fp = elph_get_log_file(); + if (fp != NULL) + { + fflush(fp); + if (saved_stdout_fd >= 0) + { + dup2(saved_stdout_fd, STDOUT_FILENO); + close(saved_stdout_fd); + saved_stdout_fd = -1; + } + fclose(fp); + elph_set_log_file(NULL); + } +} + +void elph_driver_f2c(const char* input_file, int dft_code, MPI_Fint f_comm, + const char* log_path) +{ + MPI_Comm c_comm = MPI_Comm_f2c(f_comm); + open_letz_log(c_comm, log_path); + elph_driver(input_file, (enum ELPH_dft_code)dft_code, c_comm); + close_letz_log(); +} + +/* + * Callback-enabled bridge: fill_fn_ptr argument is kept for API compatibility + * but ignored — elph_coll_fill_gkkp is resolved directly by the linker. + */ +void elph_driver_cb_f2c(const char* input_file, int dft_code, MPI_Fint f_comm, + void* fill_fn_ptr, const char* log_path) +{ + MPI_Comm c_comm = MPI_Comm_f2c(f_comm); + open_letz_log(c_comm, log_path); + elph_driver_cb(input_file, (enum ELPH_dft_code)dft_code, c_comm, + (elph_fill_fn)elph_coll_fill_gkkp); + close_letz_log(); +} + +/* + * Extended callback bridge: callbacks passed from Fortran to handle Linux linking. + * Communicators (f_comm_q, f_comm_k) passed from Y6 PAR schemes for MPI distribution. + */ +void elph_driver_cb2_f2c(struct elph_usr_input* input_data, struct Y6_info* y6_data, + int dft_code, void* fill_fn_ptr, void* dvG_fill_fn_ptr, + const char* log_path, int i_control, + int NQ_todo, int* Q_todo, int NK_todo , int* K_todo, + MPI_Fint f_comm_world, MPI_Fint f_comm_q, MPI_Fint f_comm_k) +{ + MPI_Comm c_comm_world = MPI_Comm_f2c(f_comm_world); + MPI_Comm c_comm_q = MPI_Comm_f2c(f_comm_q); + MPI_Comm c_comm_k = MPI_Comm_f2c(f_comm_k); + + open_letz_log(c_comm_world, log_path); + + struct Y6_parallel_work* y6_work = malloc(sizeof(struct Y6_parallel_work)); + + y6_work->Q = malloc(NQ_todo * sizeof(int)); + y6_work->NQ = NQ_todo; + for (int i = 0; i < y6_work->NQ; i++) { + y6_work->Q[i] = Q_todo[i]-1; + } + + y6_work->K = malloc(NK_todo * sizeof(int)); + y6_work->NK = NK_todo; + for (int i = 0; i < y6_work->NK; i++) { + y6_work->K[i] = K_todo[i]-1; + } + + /* Y6 mode: use provided communicators and callbacks from Fortran */ + elph_driver_cb2(input_data,y6_data,y6_work,(enum ELPH_dft_code)dft_code, + (elph_fill_fn)fill_fn_ptr, + (elph_dvG_fill_fn)dvG_fill_fn_ptr,i_control,c_comm_world); + + close_letz_log(); +} diff --git a/services/elph/yambo.h b/services/elph/yambo.h new file mode 100644 index 00000000..0bc566e2 --- /dev/null +++ b/services/elph/yambo.h @@ -0,0 +1,45 @@ +/* +This file contains Yambo specific data-structures used in the Code. +*/ +#pragma once +#include + +struct Y6_parallel_work +{ + int *Q; + int NQ; + int *K; + int NK; +}; + +// Pieces of Lattice and Phonon to be passed to Yambo// +struct Y6_info +{ + int natom; + int nsym; + int timerev; + int nspin; + int nspinor; + int total_bands; + int start_band; + int end_band; + int nbnds; + int lattice_dim; + int nmag; + int nmodes; + int nq_ibz; + int nq_bz; + int nkpts_ibz; + int nkpts_bz; + int nph_sym; + int kminusq; + // Pointers to k/q lists (C-side allocated, valid for query lifetime) + float* kpt_fullBZ_crys; // (nkpts_bz, 3) full BZ k-points in crystal coords + int* kmap; // (nkpts_bz, 2) maps full BZ k to iBZ+sym + float* qpts_iBZ; // (nq_ibz, 3) iBZ q-points + float* qpts_BZ; // (nq_bz, 3) full BZ q-points + int* qmap; // (nq_bz, 2) maps full BZ q to iBZ+sym + int* nqstar; // (nq_ibz,) star multiplicity + int* kplusq_idxs; // (nq_ibz, nkpts_bz) k+q indices for each q,k pair +}; + diff --git a/services/interpolation/interpolation_driver.c b/services/interpolation/interpolation_driver.c index e5b3eae3..35347746 100644 --- a/services/interpolation/interpolation_driver.c +++ b/services/interpolation/interpolation_driver.c @@ -19,6 +19,7 @@ #include "common/print_info.h" #include "common/progess_bar.h" #include "dvloc/dvloc.h" +#include "nonloc/fcoeff.h" #include "elphC.h" #include "interpolation.h" #include "interpolation_utilities.h" diff --git a/services/io/.objects b/services/io/.objects index 7659cbda..c959e87b 100644 --- a/services/io/.objects +++ b/services/io/.objects @@ -1 +1 @@ -objs = mpi_bcast.o write_basic_data.o get_save_data.o nc4_def_var.o +objs = mpi_bcast.o write_basic_data.o get_save_data.o nc4_def_var.o nc_utils.o diff --git a/services/io/get_save_data.c b/services/io/get_save_data.c index 699b8554..6df00f71 100644 --- a/services/io/get_save_data.c +++ b/services/io/get_save_data.c @@ -22,24 +22,18 @@ #include "dvloc/dvloc.h" #include "elphC.h" #include "io.h" -#include "io/qe/qe_io.h" #include "mpi_bcast.h" +#include "nc_utils.h" #include "nonloc/fcoeff.h" #include "symmetries/symmetries.h" /*static functions */ -static void quick_read(const int ncid, char* var_name, void* data_out); -static void quick_read_float(const int ncid, char* var_name, - ELPH_float* data_out); static void alloc_and_set_Gvec( ELPH_float* gvec, const ND_int ik, const ELPH_float* totalGvecs, const ND_int ng_total, const ELPH_float* Gvecidxs, const ND_int ng_shell, const ELPH_float* lat_param, const ND_int nG, const ND_int nG_shift); -static void quick_read_sub(const int ncid, char* var_name, const size_t* startp, - const size_t* countp, void* data_out); - static void get_wfc_from_save(ND_int spin_stride_len, ND_int ik, ND_int nkiBZ, ND_int nspin, ND_int nspinor, ND_int start_band, ND_int nbnds, ND_int nG, ND_int G_shift, @@ -49,13 +43,121 @@ static void get_wfc_from_save(ND_int spin_stride_len, ND_int ik, ND_int nkiBZ, static inline ND_int get_wf_io_pool(ND_int ik, ND_int q, ND_int r); +/* Read lattice dimensions (nspin, nsym, etc.) from SAVE netCDF files.*/ +void read_lattice_dimensions(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, + struct Lattice* lattice, struct Phonon* phonon) +{ + int mpi_error; + int nsELid, nsWFid, nsLATid, tempid, retval; + + size_t temp_str_len = strlen(SAVEdir) + 100; + char* temp_str = malloc(temp_str_len); + CHECK_ALLOC(temp_str); + + int nkBZ = 0; + if (Comm->commW_rank == 0) + { + cwk_path_join(SAVEdir, "ndb.KPT_indexes", temp_str, temp_str_len); + if ((retval = nc_open(temp_str, NC_NOWRITE, &tempid))) + { + ERR(retval); + } + int header_id; + nc_type compile_prec; + if ((retval = nc_inq_varid(tempid, "HEAD_VERSION", &header_id))) + { + ERR(retval); + } + if ((retval = nc_inq_vartype(tempid, header_id, &compile_prec))) + { + ERR(retval); + } + if (compile_prec != ELPH_NC4_IO_FLOAT) + { + error_msg("Yambo and LetzElPhC compiled with different precision."); + } +#if defined(YAMBO_LT_5_1) + ELPH_float kindx_pars[7]; + quick_read(tempid, "PARS", kindx_pars); + nkBZ = (int)rint(kindx_pars[0]); +#else + int nkBZ_read; + quick_read(tempid, "nXkbz", &nkBZ_read); + nkBZ = nkBZ_read; +#endif + if ((retval = nc_close(tempid))) + { + ERR(retval); + } + } + mpi_error = MPI_Bcast(&nkBZ, 1, MPI_INT, 0, Comm->commW); + MPI_error_msg(mpi_error); + + lattice->nkpts_BZ = nkBZ; + + ELPH_float dimensions[18]; + memset(dimensions, 0, sizeof(dimensions)); + if (Comm->commW_rank == 0) + { + cwk_path_join(SAVEdir, "ns.electrons", temp_str, temp_str_len); + if ((retval = nc_open(temp_str, NC_NOWRITE, &nsELid))) + { + ERR(retval); + } + quick_read_float(nsELid, "number_of_bands", dimensions + 5); + quick_read_float(nsELid, "number_of_k-points", dimensions + 6); + quick_read_float(nsELid, "spinor_components", dimensions + 11); + quick_read_float(nsELid, "spin_polarizations", dimensions + 12); + if ((retval = nc_close(nsELid))) + { + ERR(retval); + } + + cwk_path_join(SAVEdir, "ns.lattices", temp_str, temp_str_len); + if ((retval = nc_open(temp_str, NC_NOWRITE, &nsLATid))) + { + ERR(retval); + } + quick_read_float(nsLATid, "Time_reversal", dimensions + 9); + quick_read_float(nsLATid, "N_of_symmetries", dimensions + 10); + quick_read_float(nsLATid, "N_of_RL_vectors", dimensions + 7); + if ((retval = nc_close(nsLATid))) + { + ERR(retval); + } + + cwk_path_join(SAVEdir, "ns.wf", temp_str, temp_str_len); + if ((retval = nc_open(temp_str, NC_NOWRITE, &nsWFid))) + { + ERR(retval); + } + quick_read_float(nsWFid, "WF_COMPONENTS", dimensions + 8); + if ((retval = nc_close(nsWFid))) + { + ERR(retval); + } + } + mpi_error = MPI_Bcast(dimensions, 18, ELPH_MPI_float, 0, Comm->commW); + MPI_error_msg(mpi_error); + + lattice->nspinor = rint(dimensions[11]); + lattice->nspin = rint(dimensions[12]); + lattice->timerev = rint(dimensions[9]); + lattice->total_bands = rint(dimensions[5]); + lattice->nsym = rint(dimensions[10]); + lattice->nkpts_iBZ = rint(dimensions[6]); + + int nibz = lattice->nkpts_iBZ; + + free(temp_str); +} + /* Function body */ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, ND_int start_band, ND_int end_band, struct WFC** wfcs, char* ph_save_dir, struct Lattice* lattice, struct Pseudo* pseudo, - struct Phonon* phonon, - enum ELPH_dft_code dft_code) + struct Phonon* phonon) { /* This function allocates and reads data from SAVE dir. The following data is read : wfcs(in iBZ), lattice and pseudo @@ -75,18 +177,7 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, int mpi_error; - // first get the basic dft/dfpt data from dft code (code specific) before - // anything char* pp_head = "ns.kb_pp_pwscf"; // Change this accordingly - if (dft_code == DFT_CODE_QE) - { - // char* pp_head = "ns.kb_pp_pwscf"; - get_data_from_qe(lattice, phonon, pseudo, ph_save_dir, NULL, Comm); - } - else - { - error_msg("Only QE supported"); - } lattice->nfftz_loc = get_mpi_local_size_idx( lattice->fft_dims[2], &(lattice->nfftz_loc_shift), Comm->commK); @@ -104,19 +195,15 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, char* temp_str = malloc(temp_str_len); CHECK_ALLOC(temp_str); - int nkBZ = 0; // total kpoints in BZ - + int nkBZ = 0; /*****/ if (Comm->commW_rank == 0) { cwk_path_join(SAVEdir, "ndb.KPT_indexes", temp_str, temp_str_len); - if ((retval = nc_open(temp_str, NC_NOWRITE, &tempid))) { ERR(retval); } - // Before doing anything crazy, check if SAVE and - // letzElph are compiled with same precission int header_id; nc_type compile_prec; if ((retval = nc_inq_varid(tempid, "HEAD_VERSION", &header_id))) @@ -127,29 +214,24 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, { ERR(retval); } - if (compile_prec != ELPH_NC4_IO_FLOAT) { - error_msg( - "Yambo and LetzElPhC are compiled with different precision."); + error_msg("Yambo and LetzElPhC are compiled with different precision."); } - #if defined(YAMBO_LT_5_1) ELPH_float kindx_pars[7]; quick_read(tempid, "PARS", kindx_pars); - nkBZ = (int)rint(kindx_pars[0]); // FIX ME !! or kindx_pars[5] ? + nkBZ = (int)rint(kindx_pars[0]); #else int nkBZ_read; quick_read(tempid, "nXkbz", &nkBZ_read); nkBZ = nkBZ_read; #endif - if ((retval = nc_close(tempid))) { ERR(retval); } } - /* broad cast (int nkBZ) */ mpi_error = MPI_Bcast(&nkBZ, 1, MPI_INT, 0, Comm->commW); MPI_error_msg(mpi_error); /*******/ @@ -159,9 +241,8 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, "There are no kpoints in some cpus, Make sure nkpool < # of " "kpoints in full BZ."); } - // set nBZ lattice->nkpts_BZ = nkBZ; - // printf("Debug-%d \n",1); + ELPH_float dimensions[18]; memset(dimensions, 0, sizeof(dimensions)); if (Comm->commW_rank == 0) @@ -192,24 +273,21 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, } quick_read_float(nsWFid, "WF_COMPONENTS", dimensions + 8); } - /* bcast ELPH_float dimensions[18] */ mpi_error = MPI_Bcast(dimensions, 18, ELPH_MPI_float, 0, Comm->commW); MPI_error_msg(mpi_error); - // over write nspinor. This will allow for do so non-SOC ph calcs with soc - // wfcs. lattice->nspinor = rint(dimensions[11]); lattice->nspin = rint(dimensions[12]); lattice->timerev = rint(dimensions[9]); lattice->total_bands = rint(dimensions[5]); lattice->nsym = rint(dimensions[10]); lattice->nkpts_iBZ = rint(dimensions[6]); - ; int nibz = lattice->nkpts_iBZ; if (start_band < 1 || end_band < 1) { +#if !defined _Y6_LETZ if (Comm->commW_rank == 0) { fprintf(stderr, @@ -217,12 +295,14 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, "matrix elements for all bands, " "Bands index belong to [1,nbnds] \n"); } +#endif start_band = 1; end_band = lattice->total_bands; } if (start_band > lattice->total_bands || end_band > lattice->total_bands || start_band >= end_band) { +#if !defined _Y6_LETZ if (Comm->commW_rank == 0) { fprintf( @@ -230,6 +310,7 @@ void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, "Warning : invalid bands used in calculation. computing matrix " "elements for all bands \n"); } +#endif start_band = 1; end_band = lattice->total_bands; } @@ -759,72 +840,6 @@ static void get_wfc_from_save(ND_int spin_stride_len, ND_int ik, ND_int nkiBZ, } } -static void quick_read(const int ncid, char* var_name, void* data_out) -{ /* Serial read - load the entire varible data to data_out pointer - */ - int varid, retval; - - if ((retval = nc_inq_varid(ncid, var_name, &varid))) - { - ERR(retval); // get the varible id of the file - } - - if ((retval = nc_get_var(ncid, varid, data_out))) - { - ERR(retval); // get data in floats - } -} - -static void quick_read_float(const int ncid, char* var_name, - ELPH_float* data_out) -{ /* Serial read - load the entire varible data to data_out pointer - */ - int varid, retval; - - if ((retval = nc_inq_varid(ncid, var_name, &varid))) - { - ERR(retval); // get the varible id of the file - } - -#if defined(COMPILE_ELPH_DOUBLE) - if ((retval = nc_get_var_double(ncid, varid, data_out))) - { - ERR(retval); - } -#else - if ((retval = nc_get_var_float(ncid, varid, data_out))) - { - ERR(retval); - } -#endif -} - -static void quick_read_sub(const int ncid, char* var_name, const size_t* startp, - const size_t* countp, void* data_out) -{ /* Serial read - load the slice of varible data to data_out pointer - */ - int varid, retval; - - if ((retval = nc_inq_varid(ncid, var_name, &varid))) - { - ERR(retval); // get the varible id of the file - } - - // collective IO - if ((retval = nc_var_par_access(ncid, varid, NC_COLLECTIVE))) - { - ERR(retval); // NC_COLLECTIVE or NC_INDEPENDENT - } - - if ((retval = nc_get_vara(ncid, varid, startp, countp, data_out))) - { - ERR(retval); // get data in floats - } -} - static inline ND_int get_wf_io_pool(ND_int ik, ND_int q, ND_int r) { // a tiny helper function to find which pool should read the wfc diff --git a/services/io/io.h b/services/io/io.h index 0214d72b..f2809fd8 100644 --- a/services/io/io.h +++ b/services/io/io.h @@ -11,12 +11,14 @@ #define NC4_DEFAULT_CHUCK_KB 2048 // default chunking for large nc varaibles (in Kilobytes) +void read_lattice_dimensions(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, + struct Lattice* lattice, struct Phonon* phonon); + void read_and_alloc_save_data(char* SAVEdir, const struct ELPH_MPI_Comms* Comm, ND_int start_band, ND_int end_band, struct WFC** wfcs, char* ph_save_dir, struct Lattice* lattice, struct Pseudo* pseudo, - struct Phonon* phonon, - enum ELPH_dft_code dft_code); + struct Phonon* phonon); void free_save_data(struct WFC* wfcs, struct Lattice* lattice, struct Pseudo* pseudo, struct Phonon* phonon); diff --git a/services/io/nc_utils.c b/services/io/nc_utils.c new file mode 100644 index 00000000..25716153 --- /dev/null +++ b/services/io/nc_utils.c @@ -0,0 +1,44 @@ +#include +#include +#include +#include "elphC.h" +#include "common/error.h" +#include "nc_utils.h" + +void quick_read(const int ncid, char* var_name, void* data_out) +{ + int varid, retval; + if ((retval = nc_inq_varid(ncid, var_name, &varid))) + ERR(retval); + if ((retval = nc_get_var(ncid, varid, data_out))) + ERR(retval); +} + +void quick_read_float(const int ncid, char* var_name, ELPH_float* data_out) +{ + int varid, retval; + if ((retval = nc_inq_varid(ncid, var_name, &varid))) + ERR(retval); + +#if defined(COMPILE_ELPH_DOUBLE) + if ((retval = nc_get_var_double(ncid, varid, data_out))) + ERR(retval); +#else + if ((retval = nc_get_var_float(ncid, varid, data_out))) + ERR(retval); +#endif +} + +void quick_read_sub(const int ncid, char* var_name, const size_t* startp, + const size_t* countp, void* data_out) +{ + int varid, retval; + if ((retval = nc_inq_varid(ncid, var_name, &varid))) + ERR(retval); + + if ((retval = nc_var_par_access(ncid, varid, NC_COLLECTIVE))) + ERR(retval); + + if ((retval = nc_get_vara(ncid, varid, startp, countp, data_out))) + ERR(retval); +} diff --git a/services/io/nc_utils.h b/services/io/nc_utils.h new file mode 100644 index 00000000..c5de19a9 --- /dev/null +++ b/services/io/nc_utils.h @@ -0,0 +1,17 @@ +/* + * NetCDF utility functions for reading variables + */ +#pragma once + +#include +#include "elphC.h" + +/* Read entire variable from NetCDF file */ +void quick_read(const int ncid, char* var_name, void* data_out); + +/* Read entire ELPH_float variable from NetCDF file */ +void quick_read_float(const int ncid, char* var_name, ELPH_float* data_out); + +/* Read slice of variable from NetCDF file (collective I/O) */ +void quick_read_sub(const int ncid, char* var_name, const size_t* startp, + const size_t* countp, void* data_out); diff --git a/services/io/qe/get_data_from_qe.c b/services/io/qe/get_data_from_qe.c index 80cb7810..0f45290d 100644 --- a/services/io/qe/get_data_from_qe.c +++ b/services/io/qe/get_data_from_qe.c @@ -52,8 +52,9 @@ void get_data_from_qe(struct Lattice* lattice, struct Phonon* phonon, phonon->nq_iBZ_loc = distribute_to_grps(phonon->nq_iBZ, Comm->nqpools, Comm->commW_rank / Comm->commQ_size, &phonon->nq_shift); + phonon->nq_iBZ_loc = phonon->NQ_par; - if (phonon->nq_iBZ_loc < 1) + if (phonon->NQ_par < 1) { error_msg( "There are no qpoints in some qpools, Make sure nqpool < # of " diff --git a/services/io/qe/get_dvscf_dyn_qe.c b/services/io/qe/get_dvscf_dyn_qe.c index ec3776f9..1d6c8bfc 100644 --- a/services/io/qe/get_dvscf_dyn_qe.c +++ b/services/io/qe/get_dvscf_dyn_qe.c @@ -34,7 +34,7 @@ void get_dvscf_dyn_qe(const char* ph_save_dir, struct Lattice* lattice, pat_vecs = malloc(sizeof(ELPH_cmplx) * nmodes * nmodes); CHECK_ALLOC(pat_vecs); } - + // Note : It is very important to note that we must read eigenvectors by a // single cpu and then broadcast to all others. This is because // small numerical accuracies can lead different diff --git a/services/main_elphC.c b/services/main_elphC.c index 27f4fbf2..325c6126 100644 --- a/services/main_elphC.c +++ b/services/main_elphC.c @@ -76,7 +76,7 @@ int main(int argc, char* argv[]) { if (calc_info->code == DFT_CODE_QE) { - create_ph_save_dir_pp_qe(calc_info->input_file); + create_ph_save_dir_pp_qe(calc_info->input_file,calc_info->ph_save_dir,"./"); } } } diff --git a/services/make.inc b/services/make.inc index 368d2924..3358ebea 100644 --- a/services/make.inc +++ b/services/make.inc @@ -1,38 +1,39 @@ -CC := mpicc +# Library paths and compiler are taken directly from yambo's configure output +# (config/setup, two levels above this plugin directory). +YAMBO_ROOT := $(abspath $(dir $(lastword $(MAKEFILE_LIST)))/../../..) +include $(YAMBO_ROOT)/config/setup -### if you are using yambo <= 5.1.2, you need to add "-DYAMBO_LT_5_1" to cflags -### add -DCOMPILE_ELPH_DOUBLE if you want to compile the code in double precession to cflags -#CFLAGS := -O3 -Wall -Wextra -march=native -fopenmp -#LD_FLAGS := -fopenmp +CC := $(cc) -##Debug flags, incase of testing always turn on the address sanitizer ! (clang is recommanded for this) -CFLAGS := -g -Wall -Wextra # -fsanitize=address,undefined -fno-omit-frame-pointer # for mkl add -LD_FLAGS := -g # -fsanitize=address,undefined -Wl,-no_compact_unwind +### add -DCOMPILE_ELPH_DOUBLE to compile in double precision +CFLAGS := -O2 -Wall -Wextra +LD_FLAGS := +OPENMP_FLAGS := +#-DELPH_OMP_PARALLEL_BUILD ## uncomment for openmp build -OPENMP_FLAGS := -#-DELPH_OMP_PARALLEL_BUILD ## uncomment for openmp build +# FFTW: yambo's lfft/ifft cover fftw3 (double); LetzElPhC defaults to single precision (fftw3f). +# Set -DCOMPILE_ELPH_DOUBLE in CFLAGS above to switch LetzElPhC to double precision. +FFTW_INC := $(ifft) +ifeq (,$(findstring -DCOMPILE_ELPH_DOUBLE,$(CFLAGS))) +# lfft from yambo already contains -lfftw3f; strip before appending to avoid duplicate +FFTW3_LIB := $(filter-out -lfftw3f,$(lfft)) -lfftw3f +else +FFTW3_LIB := $(lfft) -lfftw3 +endif -FFTW_INC := -I${YAMBO_EXT_LIBS}/gfortran/mpifort/include/ -FFTW3_LIB := -L${YAMBO_EXT_LIBS}/gfortran/mpifort/lib -lfftw3f -lfftw3 -#FFTW_INC := -I/opt/homebrew/include -#FFTW3_LIB := -L/opt/homebrew/lib -lfftw3_threads -lfftw3f -lfftw3f_omp -lfftw3_omp -lfftw3 +BLAS_LIB := $(lblas) -#BLAS_LIB := -L/opt/homebrew/opt/openblas/lib -lopenblas -BLAS_LIB := -L${YAMBO_EXT_LIBS}/gfortran/mpifort/lib -lopenblas +# NetCDF: yambo's lnetcdf/inetcdf point to the parallel-enabled build (nc_open_par) +NETCDF_INC := $(inetcdf) +NETCDF_LIB := $(lnetcdf) -NETCDF_INC := -I${YAMBO_EXT_LIBS}/gfortran/mpifort/v4/parallel/include -NETCDF_LIB := -L${YAMBO_EXT_LIBS}/gfortran/mpifort/v4/parallel/lib -lnetcdf -HDF5_LIB := -L${YAMBO_EXT_LIBS}/gfortran/mpifort/v4/parallel/lib -lhdf5_hl -lhdf5 -lz -#NETCDF_INC := -I/Users/murali/softwares/core/include -#NETCDF_LIB := -L/Users/murali/softwares/core/lib -lnetcdf -#HDF5_LIB := -L/opt/homebrew/lib -lhdf5 +# HDF5: yambo links -lhdf5_fortran; LetzElPhC needs the C high-level API (-lhdf5_hl). +# Extract the -L flags from yambo's hdf5 config and substitute C-appropriate libraries. +HDF5_LIB := $(filter -L%,$(lhdf5)) -lhdf5_hl -lhdf5 -lz +INC_DIRS := +LIBS := -INC_DIRS := -LIBS := - -#### Notes Extra CFLAGS -### add -DCOMPILE_ELPH_DOUBLE if you want to compile the code in double precession -### if you are using yambo <= 5.1.2, you need to add "-DYAMBO_LT_5_1" to cflags -### for openmp use -DELPH_OMP_PARALLEL_BUILD in CFLAGS and set -fopenmp in LD_FLAGS and CFLAGS +#### Notes / Extra CFLAGS +### for openmp: add -DELPH_OMP_PARALLEL_BUILD to CFLAGS and -fopenmp to LD_FLAGS and CFLAGS diff --git a/services/nonloc/fcoeff.c b/services/nonloc/fcoeff.c index 18ef289f..72582c40 100644 --- a/services/nonloc/fcoeff.c +++ b/services/nonloc/fcoeff.c @@ -14,6 +14,7 @@ These routines contain functions related the pseudo potentials. #include "common/constants.h" #include "common/dtypes.h" #include "common/error.h" +#include "common/free_dtypes.h" #include "elphC.h" /** Static declarations*/ @@ -281,3 +282,31 @@ static ELPH_float CGCoeff(bool jp1, bool spin_down, int l, int twomj) } } } + +void free_Pseudo_type(struct Pseudo* pseudo) +{ + if (!pseudo) + { + return; + } + + free_f_Coeff(pseudo); + free(pseudo->fCoeff); + free_Vloc_table_type(pseudo->vloc_table); + free(pseudo->Fsign); + free(pseudo->PP_table); + + if (pseudo->loc_pseudo) + { + for (ND_int itype = 0; itype < pseudo->ntype; ++itype) + { + free_local_pseudo_type(pseudo->loc_pseudo + itype); + } + } + free(pseudo->loc_pseudo); + + pseudo->loc_pseudo = NULL; + pseudo->PP_table = NULL; + pseudo->Fsign = NULL; + pseudo->fCoeff = NULL; +} diff --git a/services/nonloc/fcoeff.h b/services/nonloc/fcoeff.h index 5745ff1e..707cfb29 100644 --- a/services/nonloc/fcoeff.h +++ b/services/nonloc/fcoeff.h @@ -3,3 +3,4 @@ void alloc_and_Compute_f_Coeff(struct Lattice* lattice, struct Pseudo* pseudo); void free_f_Coeff(struct Pseudo* pseudo); +void free_Pseudo_type(struct Pseudo* pseudo); diff --git a/services/parser/parser.h b/services/parser/parser.h index ab3f85ff..4f4fad65 100644 --- a/services/parser/parser.h +++ b/services/parser/parser.h @@ -14,7 +14,9 @@ struct elph_usr_input int start_bnd; // starting band int end_bnd; // last band char* save_dir; // save dir - char* ph_save_dir; // ph_save directory + char* scf_path; // SCF run directory + char* ph_path; // ph run directory + char* ph_save_dir; // ph_path/ph_save directory char* kernel_str; // level of screening to include bool kminusq; // true if convention is "yambo" else false }; diff --git a/services/preprocessor/cli.c b/services/preprocessor/cli.c index 091eec45..25a5ba70 100644 --- a/services/preprocessor/cli.c +++ b/services/preprocessor/cli.c @@ -17,6 +17,7 @@ void ELPH_cli_parser(int argc, char* argv[], struct calc_details* calc_info) calc_info->calc = CALC_ELPH; calc_info->code = DFT_CODE_QE; strcpy(calc_info->input_file, ""); + strcpy(calc_info->ph_save_dir, "ph_save"); /* * Here are the options * --help (help) diff --git a/services/preprocessor/pp_ph_save.c b/services/preprocessor/pp_ph_save.c index 6aa2a35e..3cc8006d 100644 --- a/services/preprocessor/pp_ph_save.c +++ b/services/preprocessor/pp_ph_save.c @@ -18,18 +18,22 @@ This file contains functions which are os dependent. #include "io/ezxml/ezxml.h" #include "preprocessor.h" -static ND_int find_nqpools(const char* out_dir, char* buffer_tmp, - ND_int buffer_size); +static ND_int find_nqpools(const char* out_dir, char* buffer_tmp,ND_int buffer_size); -void create_ph_save_dir_pp_qe(const char* inp_file) +#if defined _Y6_LETZ +void create_ph_save_dir_pp_qe(const char* inp_file,char* ph_path,char* scf_path ) +#else +void create_ph_save_dir_pp_qe(const char* inp_file,const char* ph_path, const char* scf_path ) +#endif { // parse the input file // open the qe ph.x input file char PH_SAVE_DIR_NAME[ELPH_MAX_ENV_SIZE]; - // check if env exits for ph_save char* env_var_tmp = getenv("ELPH_PH_SAVE_DIR"); +/* + // check if env exits for ph_save if (env_var_tmp && strlen(env_var_tmp) > 0) { if (strlen(env_var_tmp) > (ELPH_MAX_ENV_SIZE - 1)) @@ -46,6 +50,15 @@ void create_ph_save_dir_pp_qe(const char* inp_file) strlcpy_custom(PH_SAVE_DIR_NAME, PH_SAVE_DIR_NAME_DEFAULT, ELPH_MAX_ENV_SIZE); } +*/ + + strlcpy_custom(PH_SAVE_DIR_NAME, ph_path,ELPH_MAX_ENV_SIZE); + strlcat(PH_SAVE_DIR_NAME, "/", ELPH_MAX_ENV_SIZE); + strlcat(PH_SAVE_DIR_NAME, PH_SAVE_DIR_NAME_DEFAULT, ELPH_MAX_ENV_SIZE); + + char ph_workdir_tmp[ELPH_MAX_ENV_SIZE]; + strlcpy_custom(ph_workdir_tmp, ph_path,ELPH_MAX_ENV_SIZE); + // FILE* fp = fopen(inp_file, "r"); if (fp == NULL) @@ -85,6 +98,8 @@ void create_ph_save_dir_pp_qe(const char* inp_file) { strcpy(out_dir, "./"); } + strlcat(ph_workdir_tmp,"/",ELPH_MAX_ENV_SIZE); + strlcat(ph_workdir_tmp,out_dir,ELPH_MAX_ENV_SIZE); strcpy(dyn_prefix, "matdyn"); strcpy(dvscf_prefix, ""); @@ -281,6 +296,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) // from now we use read buffer as file_name_buf // now create the ph_save directory + if (ELPH_mkdir(PH_SAVE_DIR_NAME)) { if (errno != EEXIST) @@ -294,7 +310,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) snprintf(prefix_dir, 100, "%s.save", scf_prefix); cwk_path_join_multiple( - (const char*[]){out_dir, prefix_dir, "data-file-schema.xml", NULL}, + (const char*[]){scf_path,out_dir, prefix_dir, "data-file-schema.xml", NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); cwk_path_join_multiple( @@ -340,7 +356,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) char* tmp_str = pseudo_file_xml->txt; cwk_path_join_multiple( - (const char*[]){out_dir, prefix_dir, tmp_str, NULL}, src_file_tmp, + (const char*[]){scf_path,out_dir, prefix_dir, tmp_str, NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); cwk_path_join_multiple((const char*[]){PH_SAVE_DIR_NAME, tmp_str, NULL}, @@ -356,7 +372,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) // 3) copy dyn files - snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s0", dyn_prefix); + snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s/%s0",ph_path,dyn_prefix); // first copy the dyn0 file cwk_path_join_multiple((const char*[]){PH_SAVE_DIR_NAME, "dyn0", NULL}, dest_file_tmp, PH_X_INP_READ_BUF_SIZE); @@ -396,13 +412,11 @@ void create_ph_save_dir_pp_qe(const char* inp_file) { if (xml_dyn_prefix) { - snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s%d.xml", - dyn_prefix, idyn + 1); + snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s/%s%d.xml",ph_path,dyn_prefix, idyn + 1); } else { - snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s%d", dyn_prefix, - idyn + 1); + snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s/%s%d", ph_path, dyn_prefix,idyn + 1); } char dyn_tmp_str[32]; @@ -419,8 +433,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) } // Find the number of qpools used in ph.x calculation - ND_int ph_x_qpools = - find_nqpools(out_dir, tmp_file_buf, PH_X_INP_READ_BUF_SIZE); + ND_int ph_x_qpools = find_nqpools(ph_workdir_tmp, tmp_file_buf, PH_X_INP_READ_BUF_SIZE); // 3) copy dvscf files for (int idyn = 0; idyn < ndyn; ++idyn) @@ -450,7 +463,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) scf_prefix, dvscf_prefix, dvscf_tmp_str); cwk_path_join_multiple( - (const char*[]){out_dir, ph0_tmp, dest_file_tmp, + (const char*[]){ph_path,out_dir, ph0_tmp, dest_file_tmp, tmp_file_buf, NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); } @@ -460,7 +473,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) scf_prefix, dvscf_prefix, dvscf_tmp_str); cwk_path_join_multiple( - (const char*[]){out_dir, ph0_tmp, tmp_file_buf, NULL}, + (const char*[]){ph_path,out_dir, ph0_tmp, tmp_file_buf, NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); } @@ -472,9 +485,9 @@ void create_ph_save_dir_pp_qe(const char* inp_file) (const char*[]){PH_SAVE_DIR_NAME, tmp_file_buf, NULL}, dest_file_tmp, PH_X_INP_READ_BUF_SIZE); - if (0 == copy_files(src_file_tmp, dest_file_tmp)) + if (0 != copy_files(src_file_tmp, dest_file_tmp)) { - break; + printf("Warning : Error copying file %s.\n", src_file_tmp); } } } @@ -487,14 +500,17 @@ void create_ph_save_dir_pp_qe(const char* inp_file) { char dyn_tmp_str[32]; - snprintf(dest_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s.phsave", - scf_prefix); + snprintf(src_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s/%s/_ph0/%s.phsave",ph_path,out_dir,scf_prefix); snprintf(dyn_tmp_str, 32, "patterns.%d.xml", idyn + 1); + strlcat(src_file_tmp,"/",PH_X_INP_READ_BUF_SIZE); + strlcat(src_file_tmp,dyn_tmp_str,PH_X_INP_READ_BUF_SIZE); +/* cwk_path_join_multiple( (const char*[]){out_dir, "_ph0", dest_file_tmp, dyn_tmp_str, NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); +*/ cwk_path_join_multiple( (const char*[]){PH_SAVE_DIR_NAME, dyn_tmp_str, NULL}, dest_file_tmp, PH_X_INP_READ_BUF_SIZE); @@ -509,7 +525,7 @@ void create_ph_save_dir_pp_qe(const char* inp_file) snprintf(dest_file_tmp, PH_X_INP_READ_BUF_SIZE, "%s.phsave", scf_prefix); cwk_path_join_multiple( - (const char*[]){out_dir, "_ph0", dest_file_tmp, "tensors.xml", NULL}, + (const char*[]){ph_path,out_dir, "_ph0", dest_file_tmp, "tensors.xml", NULL}, src_file_tmp, PH_X_INP_READ_BUF_SIZE); cwk_path_join_multiple( diff --git a/services/preprocessor/preprocessor.h b/services/preprocessor/preprocessor.h index ea2adc0f..0c2b077e 100644 --- a/services/preprocessor/preprocessor.h +++ b/services/preprocessor/preprocessor.h @@ -9,7 +9,11 @@ void ELPH_cli_parser(int argc, char* argv[], struct calc_details* calc_info); -void create_ph_save_dir_pp_qe(const char* inp_file); +#if defined _Y6_LETZ +void create_ph_save_dir_pp_qe(const char* inp_file,char* ph_path,char* scf_path); +#else +void create_ph_save_dir_pp_qe(const char* inp_file,const char* ph_path, const char* scf_path); +#endif void ELPH_print_version(void); void ELPH_print_help(void);