From 9e7ad090444afa460dfe7d1b09da00d89d6421ba Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 16 Jun 2026 17:05:30 +0200 Subject: [PATCH] Preserve double precision in native dmftproj and test symmetry unitarity --- fortran/dmftproj/SRC_templates/case.cf_f_mm2 | 24 +++++------ .../dmftproj/SRC_templates/case.cf_p_cubic | 8 ++-- fortran/dmftproj/dmftproj.f | 2 +- fortran/dmftproj/orthogonal.f | 2 +- fortran/dmftproj/outputqmc.f | 32 +++++++------- fortran/dmftproj/rot_dens.f | 8 ++-- fortran/dmftproj/set_ang_trans.f | 10 ++--- fortran/dmftproj/set_rotloc.f | 2 +- fortran/dmftproj/setsym.f | 12 +++--- fortran/dmftproj/symmetrize_mat.f | 8 ++-- test/python/wien2k/CMakeLists.txt | 4 ++ test/python/wien2k/wien2k_native_precision.py | 43 +++++++++++++++++++ 12 files changed, 101 insertions(+), 54 deletions(-) create mode 100644 test/python/wien2k/wien2k_native_precision.py diff --git a/fortran/dmftproj/SRC_templates/case.cf_f_mm2 b/fortran/dmftproj/SRC_templates/case.cf_f_mm2 index 3dea83d..188beaf 100644 --- a/fortran/dmftproj/SRC_templates/case.cf_f_mm2 +++ b/fortran/dmftproj/SRC_templates/case.cf_f_mm2 @@ -1,14 +1,14 @@ 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -*0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. A1 -*0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. -.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. A2 - 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -*0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. B1 - 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -*0. 0. 0. 0. 0.70710678 0. 0. 0. -.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. B2 +*0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. A1 +*0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. -.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. A2 + 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. +*0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. B1 + 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. +*0. 0. 0. 0. 0.70710678118654752 0. 0. 0. -.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. B2 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. -*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. A1 -*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. -.70710678 0. 0. A2 - 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 -*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0.70710678 0. 0. 0. 0. B1 - 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -.70710678 0. -*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678 0. 0. 0. -.70710678 0. 0. 0. 0. 0. B2 +*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. A1 +*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. -.70710678118654752 0. 0. A2 + 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 +*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. B1 + 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -.70710678118654752 0. +*0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.70710678118654752 0. 0. 0. -.70710678118654752 0. 0. 0. 0. 0. B2 diff --git a/fortran/dmftproj/SRC_templates/case.cf_p_cubic b/fortran/dmftproj/SRC_templates/case.cf_p_cubic index 69e90cd..c3eb51c 100644 --- a/fortran/dmftproj/SRC_templates/case.cf_p_cubic +++ b/fortran/dmftproj/SRC_templates/case.cf_p_cubic @@ -1,7 +1,7 @@ - 0.707106781 0. 0. 0. -0.707106781 0. 0. 0. 0. 0. 0. 0. - 0. 0.707106781 0. 0. 0. 0.707106781 0. 0. 0. 0. 0. 0. + 0.707106781186547521 0. 0. 0. -0.707106781186547521 0. 0. 0. 0. 0. 0. 0. + 0. 0.707106781186547521 0. 0. 0. 0.707106781186547521 0. 0. 0. 0. 0. 0. *0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. - 0. 0. 0. 0. 0. 0. 0.707106781 0. 0. 0. -0.707106781 0. - 0. 0. 0. 0. 0. 0. 0. 0.707106781 0. 0. 0. 0.707106781 + 0. 0. 0. 0. 0. 0. 0.707106781186547521 0. 0. 0. -0.707106781186547521 0. + 0. 0. 0. 0. 0. 0. 0. 0.707106781186547521 0. 0. 0. 0.707106781186547521 *0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. diff --git a/fortran/dmftproj/dmftproj.f b/fortran/dmftproj/dmftproj.f index 7106526..5704738 100644 --- a/fortran/dmftproj/dmftproj.f +++ b/fortran/dmftproj/dmftproj.f @@ -649,7 +649,7 @@ PROGRAM dmftproj ENDIF DO ib=kp(ik,is)%nbmin,kp(ik,is)%nbmax READ(iualmblm,*)rtetr,kp(ik,is)%eband(ib) - kp(ik,is)%tetrweight(ib)=CMPLX(rtetr,0d0) + kp(ik,is)%tetrweight(ib)=DCMPLX(rtetr,0d0) ENDDO C rtetr = tetrahedron weights of the band ib at this kpoint C the field kp(ik,is)%eband(ib) = eigenvalues of the ib band at this kpoint diff --git a/fortran/dmftproj/orthogonal.f b/fortran/dmftproj/orthogonal.f index 7d24858..479743a 100644 --- a/fortran/dmftproj/orthogonal.f +++ b/fortran/dmftproj/orthogonal.f @@ -156,7 +156,7 @@ SUBROUTINE sqrt_eigenvec(cmat,D1,m,inv) STOP ENDIF C W contains the eigenvalues of cmat. - W_comp=CMPLX(W,0d0) + W_comp=DCMPLX(W,0d0) C C Checking of the validity of the computation : C --------------------------------------------- diff --git a/fortran/dmftproj/outputqmc.f b/fortran/dmftproj/outputqmc.f index f0df75d..83f35d1 100644 --- a/fortran/dmftproj/outputqmc.f +++ b/fortran/dmftproj/outputqmc.f @@ -220,22 +220,22 @@ SUBROUTINE outqmc(elecn,qbbot) C The spinor-rotation matrix is directly calculated from the Euler angles a,b and c. IF (rotloc(iatom)%timeinv) THEN factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - spinrot(2,1)=EXP(CMPLX(0.d0,factor))* + spinrot(2,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) spinrot(1,2)=-CONJG(spinrot(2,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - spinrot(2,2)=-EXP(CMPLX(0.d0,factor))* + spinrot(2,2)=-EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) spinrot(1,1)=CONJG(spinrot(2,2)) ELSE factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - spinrot(1,1)=EXP(CMPLX(0.d0,factor))* + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) spinrot(2,2)=CONJG(spinrot(1,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - spinrot(1,2)=EXP(CMPLX(0.d0,factor))* + spinrot(1,2)=EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) spinrot(2,1)=-CONJG(spinrot(1,2)) ENDIF @@ -674,7 +674,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (g-a) if beta=Pi. C Up/up and Dn/dn terms - spinrot(1,1)=EXP(CMPLX(0.d0,factor)) + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(2,2)=CONJG(spinrot(1,1)) C spinrot(1,1) = -exp(+i(alpha-gamma)/2) ; spinrot(2,2) = -exp(-i(alpha-gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. @@ -685,7 +685,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (a+g) if beta=0. C Up/up and Dn/dn terms - spinrot(1,1)=EXP(CMPLX(0.d0,factor)) + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(2,2)=CONJG(spinrot(1,1)) C the field phase is 2pi-(alpha+gamma) in this case. C spinrot(1,1) = -exp(-i(alpha+gamma)/2) ; spinrot(2,2) = -exp(i(alpha-gamma)/2) @@ -755,7 +755,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (g-a) in this case. C Up/up block : - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) C As a result, ephase = -exp(i(alpha-gamma)/2) spinrot(1:2*l+1,1:2*l+1)= = ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l) @@ -771,7 +771,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (a+g) in this case. C Up/up block : - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) C As a result, ephase = -exp(-i(alpha+gamma)/2) spinrot(1:2*l+1,1:2*l+1)= = ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l) @@ -1085,22 +1085,22 @@ SUBROUTINE outqmc(elecn,qbbot) C The spinor-rotation matrix is directly calculated from the Euler angles a,b and c. IF (rotloc(iatom)%timeinv) THEN factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - spinrot(2,1)=EXP(CMPLX(0.d0,factor))* + spinrot(2,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) spinrot(1,2)=-CONJG(spinrot(2,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - spinrot(2,2)=-EXP(CMPLX(0.d0,factor))* + spinrot(2,2)=-EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) spinrot(1,1)=CONJG(spinrot(2,2)) ELSE factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - spinrot(1,1)=EXP(CMPLX(0.d0,factor))* + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) spinrot(2,2)=CONJG(spinrot(1,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - spinrot(1,2)=EXP(CMPLX(0.d0,factor))* + spinrot(1,2)=EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) spinrot(2,1)=-CONJG(spinrot(1,2)) ENDIF @@ -1242,7 +1242,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (g-a) if beta=Pi. C Up/up and Dn/dn terms - spinrot(1,1)=EXP(CMPLX(0.d0,factor)) + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(2,2)=CONJG(spinrot(1,1)) C spinrot(1,1) = -exp(i(alpha-gamma)/2) ; spinrot(2,2) = -exp(-i(alpha-gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. @@ -1253,7 +1253,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (a+g) if beta=0. C Up/up and Dn/dn terms - spinrot(1,1)=EXP(CMPLX(0.d0,factor)) + spinrot(1,1)=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(2,2)=CONJG(spinrot(1,1)) C spinrot(1,1) = -exp(-i(alpha+gamma)/2) ; spinrot(2,2) = -exp(i(alpha-gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. @@ -1302,7 +1302,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is (g-a) in this case. C Up/up block : - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) C AS a result, ephase = -exp(i(alpha-gamma)/2) spinrot(1:2*l+1,1:2*l+1)= = ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l) @@ -1318,7 +1318,7 @@ SUBROUTINE outqmc(elecn,qbbot) factor=srot(isym)%phase/2.d0 C We remind that the field phase is 2pi-(alpha+gamma) in this case. C Up/up block : - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) C As a result, ephase = -exp(-i(alpha+gamma)/2) spinrot(1:2*l+1,1:2*l+1)= = ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l) diff --git a/fortran/dmftproj/rot_dens.f b/fortran/dmftproj/rot_dens.f index 32e6261..34b708c 100644 --- a/fortran/dmftproj/rot_dens.f +++ b/fortran/dmftproj/rot_dens.f @@ -64,12 +64,12 @@ SUBROUTINE rotdens_mat(Dmat,orbit,norbit) ALLOCATE(rot_dmat(1:2,1:2)) IF (rotloc(iatom)%timeinv) THEN factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - tmp_mat(2,1)=EXP(CMPLX(0.d0,factor))* + tmp_mat(2,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) tmp_mat(1,2)=-CONJG(tmp_mat(2,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - tmp_mat(2,2)=-EXP(CMPLX(0.d0,factor))* + tmp_mat(2,2)=-EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) tmp_mat(1,1)=CONJG(tmp_mat(2,2)) C definition of the total density matrix @@ -85,12 +85,12 @@ SUBROUTINE rotdens_mat(Dmat,orbit,norbit) & rot_dmat(1:2,1:2)) ELSE factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0 - tmp_mat(1,1)=EXP(CMPLX(0.d0,factor))* + tmp_mat(1,1)=EXP(CMPLX(0.d0,factor,KIND=8))* & DCOS(rotloc(iatom)%b/2.d0) tmp_mat(2,2)=CONJG(tmp_mat(1,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0 - tmp_mat(1,2)=EXP(CMPLX(0.d0,factor))* + tmp_mat(1,2)=EXP(CMPLX(0.d0,factor,KIND=8))* & DSIN(rotloc(iatom)%b/2.d0) tmp_mat(2,1)=-CONJG(tmp_mat(1,2)) C definition of the total density matrix diff --git a/fortran/dmftproj/set_ang_trans.f b/fortran/dmftproj/set_ang_trans.f index 361893e..bf2bb4a 100644 --- a/fortran/dmftproj/set_ang_trans.f +++ b/fortran/dmftproj/set_ang_trans.f @@ -37,7 +37,7 @@ SUBROUTINE set_ang_trans USE prnt IMPLICIT NONE CHARACTER(len=150) :: fullpath - CHARACTER(len=250) :: buf1 + CHARACTER(len=2000) :: buf1 CHARACTER(len=25) :: basis_file CHARACTER(len=1) :: repsign INTEGER, DIMENSION(2*(2*lmax+1)) :: degrep @@ -140,10 +140,10 @@ SUBROUTINE set_ang_trans degrep(irep)=m-ind+1 ind=m+1 ENDIF - READ(buf1(2:250),*)(rtrans(m1),itrans(m1),m1=-l,l) + READ(buf1(2:),*)(rtrans(m1),itrans(m1),m1=-l,l) C The line of the file is stored in the column of reptrans, which is temporarly "P". reptrans(l,isrt)%transmat(-l:l,m)= - & CMPLX(rtrans(-l:l),itrans(-l:l)) + & CMPLX(rtrans(-l:l),itrans(-l:l),KIND=8) m = m + 1 ENDIF ENDDO @@ -201,10 +201,10 @@ SUBROUTINE set_ang_trans degrep(irep)=m-ind+1 ind=m+1 ENDIF - READ(buf1(2:250),*)(rtrans(m1),itrans(m1), + READ(buf1(2:),*)(rtrans(m1),itrans(m1), & m1=1,2*(2*l+1)) tempmat(1:2*(2*l+1),m)= - = CMPLX(rtrans(1:2*(2*l+1)),itrans(1:2*(2*l+1))) + = CMPLX(rtrans(1:2*(2*l+1)),itrans(1:2*(2*l+1)),KIND=8) C The lines of the read matrix are stored in the column of tempmat, which is then P. ENDDO C diff --git a/fortran/dmftproj/set_rotloc.f b/fortran/dmftproj/set_rotloc.f index d2cbec1..478382c 100644 --- a/fortran/dmftproj/set_rotloc.f +++ b/fortran/dmftproj/set_rotloc.f @@ -89,7 +89,7 @@ SUBROUTINE set_rotloc C Whatever the value of beta (0 or Pi), the spinor rotation matrix of isym is block-diagonal. C because the time-reversal operation have been applied if necessary. factor=srot(isym)%phase/2.d0 - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) C We remind that the field phase is (g-a) if beta=Pi. As a result, ephase = exp(+i(g-a)/2) = -exp(+i(alpha-gamma)/2) C We remind that the field phase is (a+g) if beta=0. As a result, ephase = exp(+i(a+g)/2)=-exp(-i(alpha+gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. diff --git a/fortran/dmftproj/setsym.f b/fortran/dmftproj/setsym.f index 877a200..0e6afa1 100644 --- a/fortran/dmftproj/setsym.f +++ b/fortran/dmftproj/setsym.f @@ -505,12 +505,12 @@ SUBROUTINE setsym C in agreement with Wien conventions used for the definition of spmt (in SRC_lapwdm/sym.f) C Up/up and Dn/dn terms factor=(rotloc(iatomref)%a+rotloc(iatomref)%g)/2.d0 - spmt(1,1)=EXP(CMPLX(0.d0,factor)) + spmt(1,1)=EXP(CMPLX(0.d0,factor,KIND=8)) & *DCOS(rotloc(iatomref)%b/2.d0) spmt(2,2)=CONJG(spmt(1,1)) C Up/dn and Dn/up terms factor=-(rotloc(iatomref)%a-rotloc(iatomref)%g)/2.d0 - spmt(1,2)=EXP(CMPLX(0.d0,factor)) + spmt(1,2)=EXP(CMPLX(0.d0,factor,KIND=8)) & *DSIN(rotloc(iatomref)%b/2.d0) spmt(2,1)=-CONJG(spmt(1,2)) C Up/up block : @@ -826,7 +826,7 @@ SUBROUTINE spinrotmat(spinrot,isym,l) C We remind that the field phase is (g-a) in this case. C as a result, ephase = exp(+i(g-a)/2) = -exp(+i(alpha-gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(1:2*l+1,2*l+2:2*(2*l+1))= = ephase*srot(isym)%rotl(-l:l,-l:l,l) C Dn/up block : @@ -842,7 +842,7 @@ SUBROUTINE spinrotmat(spinrot,isym,l) C We remind that the field phase is (a+g) in this case. C as a result, ephase = exp(+i(a+g)/2)=-exp(-i(alpha+gamma)/2) C in good agreement with Wien conventions for the definition of this phase factor. - ephase=EXP(CMPLX(0.d0,factor)) + ephase=EXP(CMPLX(0.d0,factor,KIND=8)) spinrot(1:2*l+1,1:2*l+1)= = ephase*srot(isym)%rotl(-l:l,-l:l,l) C Dn/dn block : @@ -861,11 +861,11 @@ SUBROUTINE spinrotmat(spinrot,isym,l) C in agreement with Wien conventions used for the definition of spmt (in SRC_lapwdm/sym.f) C Up/up and Dn/dn terms factor=(srot(isym)%a+srot(isym)%g)/2.d0 - spmt(1,1)=EXP(CMPLX(0.d0,factor))*DCOS(srot(isym)%b/2.d0) + spmt(1,1)=EXP(CMPLX(0.d0,factor,KIND=8))*DCOS(srot(isym)%b/2.d0) spmt(2,2)=CONJG(spmt(1,1)) C Up/dn and Dn/up terms factor=-(srot(isym)%a-srot(isym)%g)/2.d0 - spmt(1,2)=EXP(CMPLX(0.d0,factor))*DSIN(srot(isym)%b/2.d0) + spmt(1,2)=EXP(CMPLX(0.d0,factor,KIND=8))*DSIN(srot(isym)%b/2.d0) spmt(2,1)=-CONJG(spmt(1,2)) C Up/up block : spinrot(1:2*l+1,1:2*l+1)= diff --git a/fortran/dmftproj/symmetrize_mat.f b/fortran/dmftproj/symmetrize_mat.f index 694f75e..774c590 100644 --- a/fortran/dmftproj/symmetrize_mat.f +++ b/fortran/dmftproj/symmetrize_mat.f @@ -99,13 +99,13 @@ SUBROUTINE symmetrize_mat(Dmat,orbit,norbit) ephase=1.d0 C For the up/dn block, initialisation of the phase factor IF(is==3) THEN - ephase=EXP(CMPLX(0d0,srot(isym)%phase)) + ephase=EXP(CMPLX(0d0,srot(isym)%phase,KIND=8)) C if srot%timeinv = .TRUE. , phase= g-a = 2pi+(alpha-gamma) and ephase = exp(+i(g-a)) = exp(+i(alpha-gamma)) C if srot%timeinv = .FALSE., phase= a+g = 2pi-(alpha+gamma) and ephase = exp(+i(a+g)) = exp(-i(alpha+gamma)) ENDIF C For the dn/up block, initialisation of the phase factor IF(is==4) THEN - ephase=EXP(CMPLX(0d0,-srot(isym)%phase)) + ephase=EXP(CMPLX(0d0,-srot(isym)%phase,KIND=8)) C if srot%timeinv = .TRUE. , phase= g-a = 2pi+(alpha-gamma) and ephase = exp(-i(g-a)) = exp(-i(alpha-gamma)) C if srot%timeinv = .FALSE., phase= a+g = 2pi-(alpha+gamma) and ephase = exp(-i(a+g)) = exp(+i(alpha+gamma)) ENDIF @@ -231,13 +231,13 @@ SUBROUTINE symmetrize_mat(Dmat,orbit,norbit) ephase=1.d0 C For the up/dn block, initialisation of the phase factor IF(is==3) THEN - ephase=EXP(CMPLX(0d0,srot(isym)%phase)) + ephase=EXP(CMPLX(0d0,srot(isym)%phase,KIND=8)) C if srot%timeinv = .TRUE. , phase= g-a = 2pi+(alpha-gamma) and ephase = exp(+i(g-a)) = exp(+i(alpha-gamma)) C if srot%timeinv = .FALSE., phase= a+g = 2pi-(alpha+gamma) and ephase = exp(+i(a+g)) = exp(-i(alpha+gamma)) ENDIF C For the dn/up block, initialisation of the phase factor IF(is==4) THEN - ephase=EXP(CMPLX(0d0,-srot(isym)%phase)) + ephase=EXP(CMPLX(0d0,-srot(isym)%phase,KIND=8)) C if srot%timeinv = .TRUE. , phase= g-a = 2pi+(alpha-gamma) and ephase = exp(-i(g-a)) = exp(-i(alpha-gamma)) C if srot%timeinv = .FALSE., phase= a+g = 2pi-(alpha+gamma) and ephase = exp(-i(a+g)) = exp(+i(alpha+gamma)) ENDIF diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index c43cbed..f6ed02c 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/test/python/wien2k/CMakeLists.txt @@ -26,6 +26,10 @@ set_property(TEST Py_wien2k_soc_convert APPEND PROPERTY ENVIRONMENT PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD}) # Pure-Python case.symqmc generator vs the shipped dmftproj Fortran output. +# Native precision is checked against unitarity, without regenerating fixtures. +add_test(NAME Py_wien2k_native_precision + COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_native_precision.py $) + # The generator tests read source fixtures and write only to scratch directories. add_test(NAME Py_wien2k_symqmc_python COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_python.py diff --git a/test/python/wien2k/wien2k_native_precision.py b/test/python/wien2k/wien2k_native_precision.py new file mode 100644 index 0000000..2a1f521 --- /dev/null +++ b/test/python/wien2k/wien2k_native_precision.py @@ -0,0 +1,43 @@ +"""Check native precision against unitarity, independently of golden outputs.""" + +import gzip +import os +from pathlib import Path +import shutil +import subprocess +import sys +import tempfile + +import numpy as np + +here = Path(__file__).resolve().parent +with tempfile.TemporaryDirectory() as tmp: + case_dir = Path(tmp) / 'CaOs2' + case_dir.mkdir() + for suffix in ('struct', 'dmftsym', 'indmftpr'): + shutil.copyfile(here / f'CaOs2.{suffix}', case_dir / f'CaOs2.{suffix}') + for spin in ('up', 'dn'): + with gzip.open(here / f'CaOs2.almblm{spin}.gz', 'rb') as src: + (case_dir / f'CaOs2.almblm{spin}').write_bytes(src.read()) + templates = Path(tmp) / 'SRC_templates' + templates.mkdir() + # Wien2k supplies the d template externally. Use full precision here to + # distinguish constructor rounding from truncation in the input template. + basis = (here / 'case.cf_d_eg_t2g').read_text() + basis = basis.replace('.70710678', '.7071067811865475244') + (templates / 'case.cf_d_eg_t2g').write_text(basis) + run = subprocess.run([sys.argv[1], '-so', '-sp'], cwd=case_dir, + env=dict(os.environ, WIENROOT=tmp), + stdout=subprocess.PIPE, stderr=subprocess.STDOUT, + text=True) + assert run.returncode == 0, run.stdout + tokens = (case_dir / 'CaOs2.symqmc').read_text().split() + header_size = 2 + 16 * 3 + 16 + blocks = np.array(tokens[header_size:], dtype=float).reshape(16, 2, 2, 10, 10) + matrices = blocks[:, :, 0] + 1j * blocks[:, :, 1] + for operation in matrices: + for matrix in operation: + np.testing.assert_allclose(matrix @ matrix.conj().T, np.eye(10), + atol=1e-12, rtol=0) + +print('wien2k_native_precision: ok')