diff --git a/CMakeLists.txt b/CMakeLists.txt index c96eb74b..126ba68b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -145,8 +145,6 @@ endif() # Python if(PythonSupport) add_subdirectory(python/${PROJECT_NAME}) - # elk binary i/o code for wrappers - add_subdirectory(python/${PROJECT_NAME}/converters/elktools/elkwrappers) endif() # Docs diff --git a/python/triqs_dft_tools/converters/elk.py b/python/triqs_dft_tools/converters/elk.py index e58fb219..166dd888 100644 --- a/python/triqs_dft_tools/converters/elk.py +++ b/python/triqs_dft_tools/converters/elk.py @@ -787,8 +787,8 @@ def convert_transport_input(self): # velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics # and each k-point. - #load fortran wrapper module - import triqs_dft_tools.converters.elktools.elkwrappers.getpmatelk as et + #load pure Python module to read PMAT.OUT + from triqs_dft_tools.converters.elktools import getpmatelk #elk velocities for all bands pmat=numpy.zeros([self.nstsv,self.nstsv,3],dtype=complex) @@ -802,10 +802,8 @@ def convert_transport_input(self): mpi.report("Reading PMAT.OUT") #read velocities for each k-point for ik in range(self.n_k): - #need to use a fortran array for wrapper - f_vkl = numpy.asfortranarray(self.vkl[ik,:]) - #read the ik velocity using the wrapper - pmat[:,:,:]=et.getpmatelk(ik+1,self.nstsv,f_vkl) + #read the ik velocity using pure Python + pmat[:,:,:]=getpmatelk(ik+1, self.nstsv, self.vkl[ik,:]) #loop over spin for isp in range(n_spin_blocks): #no. correlated bands at ik diff --git a/python/triqs_dft_tools/converters/elktools/__init__.py b/python/triqs_dft_tools/converters/elktools/__init__.py index 583b4fc6..53d98a01 100644 --- a/python/triqs_dft_tools/converters/elktools/__init__.py +++ b/python/triqs_dft_tools/converters/elktools/__init__.py @@ -25,6 +25,7 @@ from .readElkfiles import readElkfiles from .elk_converter_tools import ElkConverterTools +from .getpmatelk import getpmatelk -__all__ =['readElkfiles','ElkConverterTools'] +__all__ =['readElkfiles','ElkConverterTools','getpmatelk'] diff --git a/python/triqs_dft_tools/converters/elktools/elkwrappers/CMakeLists.txt b/python/triqs_dft_tools/converters/elktools/elkwrappers/CMakeLists.txt deleted file mode 100644 index baea83c0..00000000 --- a/python/triqs_dft_tools/converters/elktools/elkwrappers/CMakeLists.txt +++ /dev/null @@ -1,36 +0,0 @@ -# F2PY headers -execute_process( - COMMAND "${TRIQS_PYTHON_EXECUTABLE}" -c - "import numpy.f2py; print(numpy.f2py.get_include())" - OUTPUT_VARIABLE F2PY_INCLUDE_DIR - OUTPUT_STRIP_TRAILING_WHITESPACE) - -add_library(fortranobject OBJECT "${F2PY_INCLUDE_DIR}/fortranobject.c") -target_link_libraries(fortranobject PUBLIC Python::NumPy) -target_include_directories(fortranobject PUBLIC "${F2PY_INCLUDE_DIR}") -set_property(TARGET fortranobject PROPERTY POSITION_INDEPENDENT_CODE ON) - -add_custom_command( - OUTPUT getpmatelkmodule.c getpmatelk-f2pywrappers.f getpmatelk-f2pywrappers2.f90 - DEPENDS getpmatelk.f90 - VERBATIM - COMMAND "${CMAKE_COMMAND}" -E touch "getpmatelkmodule.c" "getpmatelk-f2pywrappers.f" "getpmatelk-f2pywrappers2.f90" - COMMAND "${TRIQS_PYTHON_EXECUTABLE}" -m numpy.f2py - "${CMAKE_CURRENT_SOURCE_DIR}/getpmatelk.f90" -m getpmatelk --lower) - -Python_add_library(getpmatelk MODULE - "${CMAKE_CURRENT_BINARY_DIR}/getpmatelkmodule.c" - "${CMAKE_CURRENT_BINARY_DIR}/getpmatelk-f2pywrappers.f" - "${CMAKE_CURRENT_BINARY_DIR}/getpmatelk-f2pywrappers2.f90" - "${CMAKE_CURRENT_SOURCE_DIR}/getpmatelk.f90") -set_property(TARGET getpmatelk PROPERTY SUFFIX "${TRIQS_PYTHON_MODULE_EXT}") -target_link_libraries(getpmatelk PRIVATE fortranobject) - -install(TARGETS getpmatelk DESTINATION ${TRIQS_PYTHON_LIB_DEST_ROOT}/${PROJECT_NAME}/converters/elktools/elkwrappers) - -# user warning -message(STATUS "-----------------------------------------------------------------------------") -message(STATUS " ******** USER NOTE ******** ") -message(STATUS " This version of DFTTools contains interface routines to read Elk's binary ") -message(STATUS " files. ") -message(STATUS "-----------------------------------------------------------------------------") diff --git a/python/triqs_dft_tools/converters/elktools/elkwrappers/getpmatelk.f90 b/python/triqs_dft_tools/converters/elktools/elkwrappers/getpmatelk.f90 deleted file mode 100644 index 64156c79..00000000 --- a/python/triqs_dft_tools/converters/elktools/elkwrappers/getpmatelk.f90 +++ /dev/null @@ -1,65 +0,0 @@ - -! Copyright (C) 2010 S. Sharma, J. K. Dewhurst and E. K. U. Gross. -! This file is distributed under the terms of the GNU General Public License. - -! version 6.2.8 file modified by A. D. N. James for interface with TRIQS - -subroutine getpmatelk(ik,nstsv,vkl,pmat) -!use modmain -implicit none -! arguments -integer, intent(in) :: ik !Need to read this in for the interface -integer, intent(in) :: nstsv !Need to read this in for the interface -real(8), intent(in) :: vkl(3) !TRIQS uses reduced kpts set -complex(8), intent(out) :: pmat(nstsv,nstsv,3) -! local variables -integer recl,nstsv_,i -real(8) vkl_(3),t1 - -!adnj - set up tolerance for lattice vectors, although this is an input in Elk, -! it is not advised to change this in Elk. Therefore, it should be fine to set -!it as a constant here. -real(8) epslat -epslat=1.d-6 - -! find the record length -inquire(iolength=recl) vkl_,nstsv_,pmat -!$OMP CRITICAL(u150) -do i=1,2 - open(150,file='PMAT.OUT',form='UNFORMATTED',access='DIRECT',recl=recl,err=10) - read(150,rec=ik,err=10) vkl_,nstsv_,pmat - exit -10 continue - if (i.eq.2) then - write(*,*) - write(*,'("Error(getpmat): unable to read from PMAT.OUT")') - write(*,*) - stop - end if - close(150) -end do -!$OMP END CRITICAL(u150) -!adnj edit - updated for vkl array from TRIQS -t1=abs(vkl(1)-vkl_(1))+abs(vkl(2)-vkl_(2))+abs(vkl(3)-vkl_(3)) -if (t1.gt.epslat) then - write(*,*) - write(*,'("Error(getpmat): differing vectors for k-point ",I8)') ik - !write(*,'(" current : ",3G18.10)') vkl(:,ik) - write(*,'(" current : ",3G18.10)') vkl(:) - write(*,'(" PMAT.OUT : ",3G18.10)') vkl_ - write(*,*) - stop -end if -if (nstsv.ne.nstsv_) then - write(*,*) - write(*,'("Error(getpmat): differing nstsv for k-point ",I8)') ik - write(*,'(" current : ",I8)') nstsv - write(*,'(" PMAT.OUT : ",I8)') nstsv_ - write(*,*) - stop -end if - -return - -end subroutine - diff --git a/python/triqs_dft_tools/converters/elktools/getpmatelk.py b/python/triqs_dft_tools/converters/elktools/getpmatelk.py new file mode 100644 index 00000000..a17ef6c5 --- /dev/null +++ b/python/triqs_dft_tools/converters/elktools/getpmatelk.py @@ -0,0 +1,116 @@ + +########################################################################## +# +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# +# Copyright (C) 2025 by N. Wentzell +# +# TRIQS is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# TRIQS. If not, see . +# +########################################################################## + +import numpy as np +import os + +def getpmatelk(ik, nstsv, vkl, filename='PMAT.OUT'): + """ + Pure Python implementation to read momentum matrix elements from Elk's PMAT.OUT file. + + This replaces the previous f2py compiled Fortran module. + + Parameters + ---------- + ik : int + K-point index (1-based, as in Fortran) + nstsv : int + Number of states (second variation) + vkl : array_like + K-point lattice coordinates (3 floats) + filename : str, optional + Path to PMAT.OUT file (default: 'PMAT.OUT') + + Returns + ------- + pmat : ndarray + Momentum matrix elements with shape (nstsv, nstsv, 3), dtype=complex128 + """ + + if not os.path.exists(filename): + raise IOError(f"Error(getpmat): unable to read from {filename}") + + # Tolerance for k-point comparison (from Elk's epslat) + epslat = 1.0e-6 + + # Data types matching Fortran: + # real(8) -> float64 (8 bytes) + # integer -> int32 (4 bytes, standard Fortran default integer) + # Note: If Elk is compiled with -fdefault-integer-8, use int64 + # complex(8) -> complex128 (16 bytes: 8-byte real + 8-byte imag) + + # Calculate record length in bytes + # Record contains: vkl(3) + nstsv + pmat(nstsv, nstsv, 3) + vkl_bytes = 3 * 8 # 3 real(8) + nstsv_bytes = 4 # 1 integer + pmat_bytes = nstsv * nstsv * 3 * 16 # nstsv*nstsv*3 complex(8) + recl = vkl_bytes + nstsv_bytes + pmat_bytes + + # Read the record at position ik (1-based indexing) + # In direct access, record ik starts at byte position (ik-1) * recl + try: + with open(filename, 'rb') as f: + # Seek to the correct record + f.seek((ik - 1) * recl) + + # Read vkl (3 float64) + vkl_read = np.fromfile(f, dtype=np.float64, count=3) + if vkl_read.size != 3: + raise IOError(f"Error(getpmat): unable to read from {filename}") + + # Read nstsv (1 int32) + nstsv_array = np.fromfile(f, dtype=np.int32, count=1) + if nstsv_array.size != 1: + raise IOError(f"Error(getpmat): unable to read from {filename}") + nstsv_read = nstsv_array[0] + + # Read pmat (nstsv*nstsv*3 complex128) + # Fortran stores arrays in column-major order + # pmat(nstsv, nstsv, 3) in Fortran means first index varies fastest + pmat_flat = np.fromfile(f, dtype=np.complex128, count=nstsv*nstsv*3) + if pmat_flat.size != nstsv*nstsv*3: + raise IOError(f"Error(getpmat): unable to read from {filename}") + + # Reshape directly to (nstsv, nstsv, 3) with Fortran order + pmat = pmat_flat.reshape((nstsv, nstsv, 3), order='F') + except (OSError, IOError) as e: + raise IOError(f"Error(getpmat): unable to read from {filename}") from e + + # Validate k-point + vkl = np.asarray(vkl) + t1 = np.sum(np.abs(vkl - vkl_read)) + if t1 > epslat: + raise ValueError( + f"Error(getpmat): differing vectors for k-point {ik}\n" + f" current : {vkl}\n" + f" PMAT.OUT : {vkl_read}" + ) + + # Validate nstsv + if nstsv != nstsv_read: + raise ValueError( + f"Error(getpmat): differing nstsv for k-point {ik}\n" + f" current : {nstsv}\n" + f" PMAT.OUT : {nstsv_read}" + ) + + return pmat diff --git a/test/python/elk/elk_transport_convert/elk_transport_convert.out.h5 b/test/python/elk/elk_transport_convert/elk_transport_convert.out.h5 deleted file mode 100644 index a32aaef9..00000000 Binary files a/test/python/elk/elk_transport_convert/elk_transport_convert.out.h5 and /dev/null differ