From b70f70a3f29dcea4a8b21ff9e827aff680706db9 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 15 Jun 2026 23:37:31 +0200 Subject: [PATCH] [dmftproj] build fromfile basis representation in double precision The fromfile angular basis (case with a "fromfile" harmonics choice, e.g. a |j,m_j> basis) is read and represented in single precision, so case.symqmc and case.sympar carry a ~1e-7 error there. Two causes in set_ang_trans.f / setsym.f: - CMPLX(rtrans, itrans) and EXP(CMPLX(0.d0, factor)) construct complex values from REAL(KIND=8) arguments without a KIND, so the constructor evaluates in default (single) precision before widening to the COMPLEX(KIND=8) targets. Add KIND=8 to the basis-read and spinor-phase casts. - the basis lines are read through a 250-column buffer (buf1) with an explicit buf1(2:250) slice, which truncates a higher-precision basis file mid-line. Widen the buffer and read the full line. Everything around these casts (the transmat/rotrep/rotl storage, the MATMULs, the read buffers) is already double precision; only the constructor intermediates leaked. The change matches them to the surrounding code. Effect: with a full-precision |j,m_j> basis, the d-shell case.symqmc now agrees with an independent double-precision reference to 1.9e-14 (was 6.0e-8), and the spinor matrices are unitary to 4e-16 (was 9e-9). The single-precision path is the one discussed on #148. No test regenerates from dmftproj (the converter tests read committed fixtures), so existing output is unaffected. ## Verification cmake + make + ctest (TRIQS, Build_Tests=ON): 100% tests passed, 0 tests failed out of 18 dmftproj on a full-precision |j,m_j> d-basis (16 ops, 8 time-reversal), vs an independent double-precision generator: before: max|fortran - ref| = 6.0e-8, unitarity dev = 9e-9 (single precision) after: max|fortran - ref| = 1.9e-14, unitarity dev = 4e-16 (double precision) Relates to #148. --- fortran/dmftproj/set_ang_trans.f | 10 +++++----- fortran/dmftproj/setsym.f | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) 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/setsym.f b/fortran/dmftproj/setsym.f index 877a200..daeb14c 100644 --- a/fortran/dmftproj/setsym.f +++ b/fortran/dmftproj/setsym.f @@ -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 :