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/python/triqs_dftkit/wien2k/_dmftproj.py b/python/triqs_dftkit/wien2k/_dmftproj.py
new file mode 100644
index 0000000..e6bb4d7
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/_dmftproj.py
@@ -0,0 +1,580 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola, C. Martins
+#
+# 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 .
+#
+################################################################################
+
+"""Shared dmftproj machinery for the pure-Python case.* generators.
+
+The five generators (symqmc, oubwin, ctqmcout, sympar, parproj) reproduce
+different outputs of the dmftproj Fortran executable but draw on the same
+underlying pieces: the list-directed almblm reader, the Fortran real/complex
+token parsers, the Wigner-D rotation in the dmftproj convention, the angular
+basis transforms (transmat), the case.indmftpr / case.dmftsym parsers, the
+proj_mode==0 band-window selection and the gfortran-style real formatters.
+This module owns each of those exactly once.
+
+All generators use the exact double-precision cubic harmonics. dmftproj reads
+the same coefficients from full-precision SRC_templates (the precision-fix PR),
+so the port reproduces its output to machine precision. `reptrans` keeps a
+`cast` flag only to reproduce the legacy single-precision template on demand.
+"""
+
+import math
+import os
+import numpy as np
+
+
+# --- list-directed token stream over a Fortran free-form text file -----------
+
+class Reader:
+ """Record-oriented reader mirroring Fortran list-directed and READ('()')
+ semantics on a free-form almblm file.
+
+ `record()` returns the whitespace tokens of the next physical line (the
+ common case: every almblm logical record fits on one line). `skip()`
+ consumes one physical line unconditionally, matching READ('()')."""
+
+ def __init__(self, path):
+ with open(path) as fh:
+ self._lines = fh.read().splitlines()
+ self._i = 0
+
+ def skip(self):
+ self._i += 1
+
+ def record(self):
+ toks = self._lines[self._i].split()
+ self._i += 1
+ return toks
+
+
+def to_float(tok):
+ """Parse a Fortran real token, accepting D/d exponents (1.0D-3) and the
+ list-directed 4.57E-002 form. Python float already handles E/e and the
+ embedded sign in the exponent; only D/d needs translation."""
+ return float(tok.replace('D', 'E').replace('d', 'e'))
+
+
+def to_complex(s):
+ s = s.strip().lstrip('(').rstrip(')')
+ re_s, im_s = s.split(',')
+ return complex(to_float(re_s), to_float(im_s))
+
+
+def read_complex(r):
+ """Read one list-directed complex `(re,im)` value, possibly split over
+ whitespace tokens."""
+ buf = list(r.record())
+ while ')' not in ''.join(buf):
+ buf += r.record()
+ return to_complex(''.join(buf))
+
+
+def read_two_complex(r):
+ """Read a record holding two list-directed complex values (Alm, Blm)."""
+ buf = list(r.record())
+ while ''.join(buf).count(')') < 2:
+ buf += r.record()
+ s = ''.join(buf)
+ cut = s.index(')') + 1
+ return to_complex(s[:cut]), to_complex(s[cut:])
+
+
+# --- angular bases: transmat = , m = -l..l --------------------------
+# Standard cubic harmonics, Wien2k convention, exact double precision. dmftproj
+# reads the same coefficients from full-precision SRC_templates (the precision-
+# fix PR); `reptrans(..., cast=True)` reproduces the legacy float32 template.
+
+_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
+
+_CUBIC = {
+ 1: np.array([[2 ** -0.5, 0, -(2 ** -0.5)],
+ [-1j * 2 ** -0.5, 0, -1j * 2 ** -0.5],
+ [0, 1, 0]], dtype=complex),
+ 2: np.array([
+ [0, 0, 1, 0, 0],
+ [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
+ [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
+ [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
+ [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
+ ], dtype=complex),
+}
+
+
+def reptrans(basis, l, cast=False):
+ """transmat = , exact double-precision cubic harmonics. dmftproj
+ reads these from full-precision SRC_templates (the precision-fix PR), so the
+ port uses the exact analytic values. A complex basis is the exact identity.
+ `cast` is retained for callers that still want the legacy float32 truncation."""
+ if basis == 'cubic' and l in _CUBIC:
+ c = _CUBIC[l]
+ return c.astype(np.complex64).astype(np.complex128) if cast else c
+ if basis == 'cubic' and l != 0:
+ raise NotImplementedError('cubic f harmonics are not implemented yet')
+ return _COMPLEX[l]
+
+
+def read_fromfile(path, l):
+ """Parse a dmftproj fromfile basis (each line: the coefficients of a new
+ basis vector in {|m,up>, |m,dn>}, m = -l..l, real/imag interleaved; '*'
+ marks the end of an irep). Return (P, mixing) where P = is the
+ transform matrix (reptrans.transmat), full 2(2l+1) for a spin-mixing basis
+ or the (2l+1) up/up block otherwise."""
+ n = 2 * (2 * l + 1)
+ rows = []
+ for line in open(path):
+ line = line.rstrip('\n')
+ if not line.strip():
+ continue
+ body = line[1:]
+ vals = [to_float(x) for x in body.split()[:2 * n]]
+ if len(vals) != 2 * n:
+ raise ValueError(f'{path}: each basis row requires {2 * n} real coefficients')
+ rows.append([vals[2 * k] + 1j * vals[2 * k + 1] for k in range(n)])
+ if len(rows) == n:
+ break
+ R = np.array(rows) # R[i, :] = |new_i> in old basis
+ if R.shape != (n, n) or not np.all(np.isfinite(R)):
+ raise ValueError(f'{path}: expected a finite {n} by {n} basis matrix')
+ d = 2 * l + 1
+ up_up, up_dn = R[:d, :d], R[:d, d:]
+ dn_up, dn_dn = R[d:, :d], R[d:, d:]
+ mixing = not (np.array_equal(dn_dn, up_up) and
+ np.all(up_dn == 0) and np.all(dn_up == 0))
+ P = np.conj(R) if mixing else np.conj(up_up) # = conj()
+ return P, mixing
+
+
+# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
+
+def small_d(l, m, n, b):
+ f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
+ (math.factorial(l + n) * math.factorial(l - n))
+ s = 0.0
+ for t in range(0, 2 * l + 1):
+ if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
+ f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
+ (math.factorial(l - m - t) * math.factorial(m + n + t) *
+ math.factorial(l - n - t) * math.factorial(t))
+ f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 \
+ else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
+ f4 = 1.0 if (2 * t + n + m) == 0 \
+ else math.cos(b / 2) ** (2 * t + n + m)
+ s += (-1) ** (l - m - t) * f2 * f3 * f4
+ return math.sqrt(f1) * s
+
+
+def dmat(l, a, b, c, det):
+ # det is the improper-rotation parity: pass op['iprop'] (+-1), not
+ # np.linalg.det(krotm). krotm in case.dmftsym is the proper part and its
+ # determinant is +1 even for improper ops, so it drops the (-1)^l parity
+ # factor. Harmless for even l, wrong sign for odd l (p, f).
+ D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ for n in range(-l, l + 1):
+ v = np.exp(1j * n * a) * np.exp(1j * m * c) * small_d(l, m, n, b)
+ if det < -0.5:
+ v *= (-1) ** l
+ D[m + l, n + l] = v
+ return D
+
+
+def tmat(l):
+ """Complex-conjugation operator in the spherical-harmonic basis,
+ T[m, -m] = (-1)^m (timeinv.f)."""
+ T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ T[-m + l, m + l] = (-1) ** m
+ return T
+
+
+def timeinv_orbital(l, mat):
+ return tmat(l) @ np.conj(mat)
+
+
+def mixing_timeinv_op(l, P):
+ """The spinor time-reversal operator -i sigma_y (x) T in the mixing basis:
+ tinv_{new} = P tinv_{lm} P^T (timeinv.f, ifmixing branch). Returned as the
+ bare operator; callers apply it as tinv @ conj(mat)."""
+ d = 2 * l + 1
+ tm = tmat(l)
+ tinv = np.zeros((2 * d, 2 * d), dtype=complex)
+ tinv[:d, d:] = -tm
+ tinv[d:, :d] = tm
+ return P @ tinv @ P.T
+
+
+def rotloc_rotl_so(l, ref, ops, iatom, iref):
+ """The composed 2(2l+1) Rloc rotation rotloc(iatom)%rotl(l) under SP+SO
+ (setsym.f:496-528 + set_rotloc.f): the representative spinor rotloc
+ spmt (x) D(rotloc_ref) composed with the first symmetry op R[isym] mapping
+ iref onto iatom, with the orbital time-reversal applied for the magnetic op.
+ Returns (rotl, timeinv); callers apply their own basis transform (the full
+ transmat for a mixing basis, blkdiag(transmat, transmat) otherwise)."""
+ d = 2 * l + 1
+ Dref = dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
+ f = (ref['a'] + ref['g']) / 2.0
+ spmt = np.zeros((2, 2), dtype=complex)
+ spmt[0, 0] = np.exp(1j * f) * math.cos(ref['b'] / 2.0)
+ spmt[1, 1] = np.conj(spmt[0, 0])
+ f = -(ref['a'] - ref['g']) / 2.0
+ spmt[0, 1] = np.exp(1j * f) * math.sin(ref['b'] / 2.0)
+ spmt[1, 0] = -np.conj(spmt[0, 1])
+ rotl = np.zeros((2 * d, 2 * d), dtype=complex)
+ rotl[:d, :d] = spmt[0, 0] * Dref
+ rotl[d:, d:] = spmt[1, 1] * Dref
+ rotl[:d, d:] = spmt[0, 1] * Dref
+ rotl[d:, :d] = spmt[1, 0] * Dref
+
+ op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+ det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+ - op['krotm'][0, 1] * op['krotm'][1, 0])
+ timeinv = det2 < 0.0
+ srot_phase = (op['g'] - op['a']) if timeinv else (op['a'] + op['g'])
+ rotl_sym = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ if timeinv:
+ rotl_sym = tmat(l) @ np.conj(rotl_sym)
+ ephase = np.exp(1j * srot_phase / 2.0)
+ tmp = np.zeros((2 * d, 2 * d), dtype=complex)
+ tmp[:d, :d] = ephase * rotl_sym
+ tmp[d:, d:] = np.conj(ephase) * rotl_sym
+ rotl = (tmp @ np.conj(rotl)) if timeinv else (tmp @ rotl)
+ return rotl, timeinv
+
+
+def rotloc_orbital(l, ref, ops, iatom, iref):
+ """Orbital part of the full local rotation R[isym] * rotloc_ref."""
+ op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+ return dmat(l, op['a'], op['b'], op['g'], float(op['iprop'])) @ \
+ dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
+
+
+def mixing_rotrep(op, l, P, ti):
+ """The full 2(2l+1) spinor representation D(R)_{new_i} = P spinrot P^dag of
+ one symmetry operation in a spin-coupling (mixing) basis, with the spinor
+ time-reversal operator applied for the magnetic (timeinv) operations
+ (setsym.f spinrotmat + timeinv_op). This is srot%rotrep(l,isrt)%mat, shared
+ by the symqmc and sympar shell matrices."""
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ phase = (op['c'] - op['a']) if ti else (op['a'] + op['c'])
+ e = np.exp(1j * phase / 2)
+ d = 2 * l + 1
+ spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
+ if ti: # beta = pi, block-antidiagonal
+ spinrot[:d, d:] = e * rotl
+ spinrot[d:, :d] = -np.conj(e) * rotl
+ else: # beta = 0, block-diagonal
+ spinrot[:d, :d] = e * rotl
+ spinrot[d:, d:] = np.conj(e) * rotl
+ rotrep = P @ spinrot @ np.conj(P.T)
+ if ti:
+ rotrep = mixing_timeinv_op(l, P) @ np.conj(rotrep)
+ return rotrep
+
+
+# --- case.indmftpr -----------------------------------------------------------
+
+def read_indmftpr(indmftpr):
+ """Full structured parse of case.indmftpr.
+
+ Returns a dict with nsort, mult, lmax, the spin-orbit flag `so`, the energy
+ window (e_bot, e_top, proj_mode) and one `sorts` entry per atomic sort with
+ its basis name, optional fromfile sourcefile path, and the correlated /
+ included l lists (l_inc==2 / l_inc in {1,2}). Every generator derives its
+ shell or orbital list from this one parse."""
+ raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+ raw = [l for l in raw if l != '']
+ nsort = int(raw[0].split()[0])
+ mult = [int(x) for x in raw[1].split()][:nsort]
+ lmax = int(raw[2].split()[0])
+ i = 3
+ so = 0
+ sorts = []
+ for isort in range(nsort):
+ basis = raw[i].split()[0]
+ i += 1
+ sourcefile = None
+ if basis == 'fromfile':
+ sourcefile = os.path.join(os.path.dirname(indmftpr), raw[i])
+ i += 1
+ l_inc = [int(x) for x in raw[i].split()]
+ i += 1
+ ireps = [int(x) for x in raw[i].split()]
+ i += 1
+ correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
+ included_ls = [l for l in range(len(l_inc)) if l_inc[l] in (1, 2)]
+ for l in correlated_ls:
+ if ireps[l] > 1:
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ if ireps[l] > 0:
+ if raw[i].strip() != '1':
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ i += 1 # each correlated l has its own correps line
+ if correlated_ls:
+ sort_so = int(raw[i].split()[0])
+ if sort_so not in (0, 1):
+ raise ValueError('spin-orbit flags must be 0 or 1')
+ if any(s['correlated_ls'] for s in sorts) and sort_so != so:
+ raise NotImplementedError('different spin-orbit flags across correlated sorts are not implemented')
+ so = sort_so # SO flag follows a correlated sort
+ i += 1
+ sorts.append(dict(basis=basis, sourcefile=sourcefile,
+ correlated_ls=correlated_ls, included_ls=included_ls))
+ last = raw[-1].split()
+ e_bot, e_top = to_float(last[0]), to_float(last[1])
+ proj_mode = int(last[2]) if len(last) >= 3 else 0
+ if proj_mode not in (0, 1, 2) or not np.isfinite([e_bot, e_top]).all() or e_bot > e_top:
+ raise ValueError('projection window requires finite ordered bounds and mode 0, 1 or 2')
+ return dict(nsort=nsort, mult=mult, lmax=lmax, sorts=sorts, so=so,
+ e_bot=e_bot, e_top=e_top, proj_mode=proj_mode)
+
+
+# --- case.dmftsym ------------------------------------------------------------
+
+def read_dmftsym(path, rotloc=False):
+ """Symmetry operations from case.dmftsym. Each op carries (perm, a, b, g,
+ iprop, krotm); the third Euler angle is named both `c` and `g` for the
+ callers that use either name. With rotloc=True also parse the
+ 'Global->local' representative rotloc per sort and return it as a third
+ value (setsym.f:437-455)."""
+ lines = open(path).read().split('\n')
+ nsym = int(lines[0].split()[0])
+ perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
+ starts = [i for i, l in enumerate(lines) if 'Sym. op.' in l]
+ ops = []
+ for k, s in enumerate(starts[:nsym]):
+ toks = lines[s + 1].split()
+ a, b, g = (math.radians(float(x)) for x in toks[:3])
+ iprop = int(toks[3])
+ krotm = np.array([[to_float(x) for x in lines[s + 2 + r].split()]
+ for r in range(3)])
+ ops.append(dict(perm=perms[k], a=a, b=b, c=g, g=g, iprop=iprop,
+ krotm=krotm))
+ if not rotloc:
+ return nsym, ops
+
+ gl = next(i for i, l in enumerate(lines) if 'Global->local' in l)
+ rotloc_ref = []
+ j = gl + 1
+ while j < len(lines):
+ if lines[j].strip() == '' or not lines[j].split()[0].lstrip('-').isdigit():
+ j += 1
+ continue
+ # sort index line, then 3 krotm rows, then the Euler/iprop line.
+ if len(lines[j].split()) == 1:
+ krotm = np.array([[to_float(x) for x in lines[j + 1 + r].split()]
+ for r in range(3)])
+ ang = lines[j + 4].split()
+ a, b, g = (math.radians(float(x)) for x in ang[:3])
+ iprop = int(ang[3])
+ rotloc_ref.append(dict(krotm=krotm, a=a, b=b, g=g, iprop=iprop))
+ j += 5
+ else:
+ j += 1
+ return nsym, ops, rotloc_ref
+
+
+# --- almblm parsing ----------------------------------------------------------
+
+def read_almblm(path, info, projectors=False):
+ """Read one spin's almblm file in the dmftproj.f read order.
+
+ With projectors=False keep only the header (elecn, nk, nloat, eferm), the
+ per-(sort,l) overlap block (u_dot_norm, nLO, ovl_LO_u, ovl_LO_udot) and the
+ per-k band data (nbmin/nbmax, Fermi-shifted eband, tetrahedron weights); the
+ per-(l,m) coefficient records are still consumed to keep the reader in sync
+ but their values are discarded (oubwin's window-only path).
+
+ With projectors=True also accumulate the per-(l,m) Alm/Blm/Clm coefficients
+ in the canonical Wien2k packing lm = l*l + (m+l), for every atom of every
+ sort (ctqmcout/parproj)."""
+ nsort = info['nsort']
+ lmax = info['lmax']
+ mult = info['mult']
+ nlm = (lmax + 1) ** 2
+ natom = sum(mult)
+
+ r = Reader(path)
+ elecn = to_float(r.record()[0])
+ nk = int(r.record()[0])
+ r.record() # nloat
+ eferm = to_float(r.record()[0])
+
+ nLO = {}
+ u_dot_norm = {}
+ ovl_LO_u = {}
+ ovl_LO_udot = {}
+ kp = [None] * nk
+
+ for isrt in range(1, nsort + 1):
+ for l in range(lmax + 1):
+ u_dot_norm[(l, isrt)] = to_float(r.record()[0])
+ n = int(r.record()[0])
+ if n > 1:
+ raise NotImplementedError('more than one local orbital per angular momentum is not implemented')
+ nLO[(l, isrt)] = n
+ for ilo in range(1, n + 1):
+ toks = r.record()
+ ovl_LO_u[(ilo, l, isrt)] = to_float(toks[0])
+ ovl_LO_udot[(ilo, l, isrt)] = to_float(toks[1])
+ for ik in range(nk):
+ r.skip() # "IK = .." banner
+ r.skip() # 3-int line
+ head = r.record()
+ nbmin, nbmax = int(head[1]), int(head[2])
+ nb = nbmax - nbmin + 1
+ if kp[ik] is None:
+ kp[ik] = dict(nbmin=nbmin, nbmax=nbmax, eband=None,
+ weight=None, tetr=None)
+ if projectors:
+ kp[ik].update(
+ Alm=np.zeros((nlm, natom + 1, nb), dtype=complex),
+ Blm=np.zeros((nlm, natom + 1, nb), dtype=complex),
+ Clm=np.zeros((4, nlm, natom + 1, nb), dtype=complex))
+ eband = np.empty(nb)
+ tetr = np.empty(nb)
+ for off in range(nb):
+ toks = r.record()
+ tetr[off] = to_float(toks[0])
+ eband[off] = to_float(toks[1])
+ eband = eband - eferm
+ if kp[ik]['eband'] is None:
+ kp[ik]['eband'] = eband
+ kp[ik]['weight'] = tetr[0]
+ kp[ik]['tetr'] = tetr
+ for imu in range(1, mult[isrt - 1] + 1):
+ iatom = sum(mult[:isrt - 1]) + imu
+ r.skip() # banner
+ r.record() # idum
+ for off in range(nb):
+ lm = 0
+ for l in range(lmax + 1):
+ for m in range(-l, l + 1):
+ alm, blm = read_two_complex(r)
+ for ilo in range(nLO[(l, isrt)]):
+ clm = read_complex(r)
+ if projectors:
+ kp[ik]['Clm'][ilo, lm, iatom, off] = clm
+ if projectors:
+ kp[ik]['Alm'][lm, iatom, off] = alm
+ kp[ik]['Blm'][lm, iatom, off] = blm
+ lm += 1
+
+ return dict(elecn=elecn, eferm=eferm, nk=nk, nlm=nlm, natom=natom,
+ nLO=nLO, u_dot_norm=u_dot_norm, ovl_LO_u=ovl_LO_u,
+ ovl_LO_udot=ovl_LO_udot, kp=kp)
+
+
+# --- band-window selection (set_projections.f, proj_mode==0) -----------------
+
+def select_window(nbmin, nbmax, eband, e1, e2):
+ opened = False
+ included = False
+ nb_bot = nb_top = 0
+ for off, ib in enumerate(range(nbmin, nbmax + 1)):
+ e = eband[off]
+ if not opened and e > e1 and e <= e2:
+ opened = True
+ nb_bot = ib
+ elif opened and e > e2:
+ nb_top = ib - 1
+ included = True
+ break
+ elif ib == nbmax and e > e1 and e <= e2:
+ nb_top = ib
+ included = True
+ if not included:
+ nb_bot = nb_top = 0
+ return included, nb_bot, nb_top
+
+
+def select_band_window(nbmin, nbmax, b_bot, b_top):
+ """proj_mode 1/2 band-index selection for one k-point (set_projections.f
+ 70-88). e1/e2 are band indices, not energies: every k-point is included,
+ nb_bot = b_bot clamped up to nbmin (strict INT(e1) > nbmin), nb_top = b_top
+ clamped down to nbmax. Returns (included=True, nb_bot, nb_top)."""
+ nb_bot = b_bot if b_bot > nbmin else nbmin
+ nb_top = b_top if b_top < nbmax else nbmax
+ return True, nb_bot, nb_top
+
+
+def band_index_window(info, spins):
+ """Resolve the (b_bot, b_top) band-index window for proj_mode 1 and 2.
+
+ proj_mode 2 (dmftproj.f:233-237): b_bot=INT(e_bot), b_top=INT(e_top) taken
+ directly from the indmftpr window line.
+
+ proj_mode 1 (dmftproj.f:704-722): e_bot/e_top are Fermi-shifted energies;
+ scan every spin and k-point for bands with e_bot < E <= e_top and take the
+ global min/max band index, seeded with b_bot=1000, b_top=1 so an empty scan
+ keeps that seed. `spins` is the list of read_almblm dicts (one per spin)."""
+ if info['proj_mode'] == 2:
+ return int(info['e_bot']), int(info['e_top'])
+ if info['proj_mode'] != 1:
+ raise ValueError('band_index_window is only valid for proj_mode 1 or 2')
+ e_bot, e_top = info['e_bot'], info['e_top']
+ b_bot, b_top = 1000, 1
+ for sp in spins:
+ for kp in sp['kp']:
+ for off, ib in enumerate(range(kp['nbmin'], kp['nbmax'] + 1)):
+ e = kp['eband'][off]
+ if e > e_bot and e <= e_top:
+ if ib > b_top:
+ b_top = ib
+ if ib < b_bot:
+ b_bot = ib
+ return b_bot, b_top
+
+
+# --- gfortran-style real formatters ------------------------------------------
+
+def fmt(x):
+ """One list-directed real, gfortran-style (leading sign space, ~17 sig)."""
+ return ' %.16E' % float(x)
+
+
+def write_row(f, arr):
+ f.write(''.join(fmt(x) for x in arr) + '\n')
+
+
+def fmt_scalar(x):
+ x = float(x)
+ return '%.16f' % x if abs(x) < 1e5 else '%.16E' % x
+
+
+def validate_so_windows(spins, info, windows, below=None):
+ """Require valid selected windows and the same absolute bands in both spins."""
+ if len(spins) != 2:
+ raise ValueError('two-spin projectors require both up and down almblm files')
+ bw = band_index_window(info, spins) if info['proj_mode'] else None
+ for spin in spins[1:]:
+ if spin['nk'] != spins[0]['nk']:
+ raise ValueError('up/down k-point counts differ')
+ for ik, kp in enumerate(spin['kp']):
+ selected = (select_window(kp['nbmin'], kp['nbmax'], kp['eband'], info['e_bot'], info['e_top'])
+ if bw is None else select_band_window(kp['nbmin'], kp['nbmax'], *bw))
+ if selected != windows[ik]:
+ raise ValueError('projector run requires identical up/dn band windows')
+ if below is not None:
+ lower = (select_window(kp['nbmin'], kp['nbmax'], kp['eband'], -1e6, info['e_bot'])
+ if bw is None else select_band_window(kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
+ if lower != below[ik]:
+ raise ValueError('projector run requires identical up/dn below-window bands')
+ if not all(incl and lo <= hi for incl, lo, hi in windows):
+ raise ValueError('projectors require a nonempty included window at every k-point')
diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py
new file mode 100644
index 0000000..59304a9
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/ctqmcout.py
@@ -0,0 +1,423 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# 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 .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj correlated-shell projector file
+(case.ctqmcout), replacing the projector path of the dmftproj Fortran
+executable (outputqmc.f subroutine outqmc).
+
+Given the alm/blm coefficient files (case.almblmup / case.almblmdn, one per
+spin), the projector definition (case.indmftpr), the structure
+(case.struct) and the symmetry input (case.dmftsym), this rebuilds the
+correlated (Wannier) projector pr_crorb%mat_rep and writes case.ctqmcout in
+the format the converter reads.
+
+It follows the Fortran path field by field:
+
+- almblm read order (dmftproj.f:559-689): header, per-sort per-l overlap
+ block, then per-k banners + idum/nbmin/nbmax + rtetr/eband lines, then the
+ per-(l,m) (Alm,Blm) / Clm complex coefficients in the canonical Wien2k lm
+ packing lm = l*l + (m+l) + 1.
+- band-window selection (set_projections.f:89-129, proj_mode==0).
+- raw correlated projector P(m,ib) = Alm(lm) + sum_ilo Clm(ilo,lm)*ovl_LO_u
+ (set_projections.f:243-267); Blm and the radial s12 do not enter pr_crorb.
+- local rotation rot_projectmat (Wigner-D of rotloc, identity here) then the
+ angular-basis transform mat_rep = transmat . (Rloc . P) where transmat =
+ is the cubic-harmonics transform (set_projections.f:450-478).
+- Loewdin orthonormalization of the stacked correlated projectors,
+ P <- O^{-1/2} P with O = D D^H over all correlated rows
+ (orthogonal_wannier_SO, dmftproj.f:847). This step is what the ctqmcout
+ numbers carry; it is applied per k over the full ndim x nbands stack.
+- the Rloc spinor block rotloc%rotrep (outputqmc.f:303-313), built from the
+ symmetry operation mapping the representative atom onto each equivalent one
+ (set_rotloc.f / setsym.f).
+
+ctqmcout layout is outputqmc.f:66-613. The cubic transform uses the exact
+analytic harmonics; with the precision-fixed dmftproj (full double precision
+templates and KIND=8 casts) the port reproduces case.ctqmcout to machine
+precision on a full-rank window.
+"""
+
+import os
+import numpy as np
+from scipy.linalg.lapack import zheev
+from scipy.linalg.blas import zgemm
+
+from ._dmftproj import (validate_so_windows, band_index_window, dmat, fmt_scalar, read_almblm,
+ read_dmftsym, read_fromfile, read_indmftpr, reptrans,
+ rotloc_orbital, rotloc_rotl_so, select_band_window, select_window,
+ tmat, write_row)
+
+
+def _sqrt_inv(O):
+ """O^{-1/2} of a Hermitian matrix, reproducing orthogonal.f sqrtm
+ (inv=.TRUE.: Z diag(w^{-1/2}) Z^H). Use the same LAPACK routine the Fortran
+ calls, ZHEEV('V','U'), via scipy so the eigenvectors match rather than the
+ divide-and-conquer ZHEEVD of numpy.linalg.eigh. dmftproj takes 1/sqrt of the
+ eigenvalue as a *complex* sqrt (W_comp = CMPLX(W,0)); reproduce that with a
+ complex power. The result D1 @ conj(Z).T matches the Fortran's ZGEMM('N','T').
+
+ Full-rank overlaps give a well-defined inverse square root. Singular
+ overlaps cannot define orthonormal correlated projectors and are rejected.
+ """
+ w, Z, info = zheev(O, compute_v=1, lower=0) # 'V', 'U'
+ if info != 0:
+ raise np.linalg.LinAlgError(f'Loewdin eigensolver failed (info={info})')
+ if w[0] <= np.finfo(float).eps * len(w) * abs(w[-1]):
+ raise ValueError('Loewdin overlap is rank deficient; widen the projection window')
+ D1 = Z * (w.astype(complex) ** -0.5) # Z @ diag(w^{-1/2})
+ return zgemm(1.0, D1, np.conj(Z), trans_b=1) # ZGEMM('N','T'): D1 @ conj(Z)^T
+
+
+# --- correlated / included orbital descriptors -------------------------------
+
+def _build_crorbs(info):
+ crorbs = []
+ for isort in range(1, info['nsort'] + 1):
+ sortinfo = info['sorts'][isort - 1]
+ for l in sortinfo['correlated_ls']:
+ for imu in range(1, info['mult'][isort - 1] + 1):
+ atom = sum(info['mult'][:isort - 1]) + imu
+ crorbs.append(dict(atom=atom, sort=isort, l=l,
+ basis=sortinfo['basis'], first=(imu == 1)))
+ return crorbs
+
+
+def _build_orbs(info):
+ orbs = []
+ for isort in range(1, info['nsort'] + 1):
+ sortinfo = info['sorts'][isort - 1]
+ for l in sortinfo['included_ls']:
+ for imu in range(1, info['mult'][isort - 1] + 1):
+ atom = sum(info['mult'][:isort - 1]) + imu
+ orbs.append(dict(atom=atom, sort=isort, l=l))
+ return orbs
+
+
+# --- Rloc spinor representation (set_rotloc.f / setsym.f) ---------------------
+
+def _rloc_rotrep(l, ref, ops, iatom, iref, transmat, ifSO):
+ """rotloc(iatom)%rotrep(l)%mat for a non-mixing shell: the composed Rloc
+ spinor rotation rotloc_ref then the symmetry op mapping iref onto iatom
+ (set_rotloc.f / setsym.f), put into the new basis with S = blkdiag(transmat,
+ transmat). Shares rotloc_rotl_so with the mixing path so the non-trivial
+ struct reference rotation rotloc_ref (e.g. Te, beta=pi/2) is not dropped.
+ Without SO it reduces to the bare (2l+1) D(op)D(ref) transmat sandwich
+ (set_rotloc.f non-SO branch). Returns (rotrep, timeinv)."""
+ if not ifSO:
+ D = rotloc_orbital(l, ref, ops, iatom, iref)
+ return transmat @ D @ np.conj(transmat.T), False
+ rotl, timeinv = rotloc_rotl_so(l, ref, ops, iatom, iref)
+ d = 2 * l + 1
+ S = np.zeros((2 * d, 2 * d), dtype=complex)
+ S[:d, :d] = transmat
+ S[d:, d:] = transmat
+ rotrep = (S @ rotl @ S.T) if timeinv else (S @ rotl @ np.conj(S.T))
+ return rotrep, timeinv
+
+
+def _rloc_rotrep_mixing(l, ref, ops, iatom, iref, P):
+ """rotloc(iatom)%rotrep(l)%mat for a mixing SO shell: the composed Rloc
+ spinor rotation (rotloc_rotl_so) put into the new basis with the full
+ 2(2l+1) transmat P (set_rotloc.f mixing branch). Returns (rotrep, timeinv)."""
+ rotl, timeinv = rotloc_rotl_so(l, ref, ops, iatom, iref)
+ rotrep = (P @ rotl @ P.T) if timeinv else (P @ rotl @ np.conj(P.T))
+ return rotrep, timeinv
+
+
+# --- ctqmcout writer ---------------------------------------------------------
+
+def write_ctqmcout(case):
+ """Read .almblm{up,dn}, .indmftpr, .struct,
+ .dmftsym and write .ctqmcout in the dmftproj format."""
+ info = read_indmftpr(case + '.indmftpr')
+ nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ if os.path.exists(up) and os.path.exists(dn):
+ spin_files = [up, dn]
+ else:
+ spin_files = [case + '.almblm']
+ ifSP = len(spin_files) == 2
+ ifSO = bool(info['so'])
+ ns = 2 if ifSP else 1
+
+ spins = [read_almblm(p, info, projectors=True) for p in spin_files]
+ nk = spins[0]['nk']
+ elecn = spins[0]['elecn']
+
+ crorbs = _build_crorbs(info)
+ orbs = _build_orbs(info)
+ ncrorb = len(crorbs)
+ norb = len(orbs)
+
+ # band-index window for proj_mode 1/2 (set_projections is called with band
+ # indices, not energies). proj_mode 1 scans all spins for the global window.
+ bw = band_index_window(info, spins) if info['proj_mode'] != 0 else None
+
+ # window in [e_bot, e_top] (proj_mode 0) or [b_bot, b_top] (mode 1/2)
+ windows = []
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ if info['proj_mode'] == 0:
+ windows.append(select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+ else:
+ windows.append(select_band_window(
+ kp['nbmin'], kp['nbmax'], bw[0], bw[1]))
+
+ # window below e_bot (mode 0) or below b_bot (mode 1/2), for qbbot.
+ # Mode 1/2: set_projections(1, b_bot-1) -> select_band_window with top=b_bot-1.
+ win_below = []
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ if info['proj_mode'] == 0:
+ win_below.append(select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'], -1e6, info['e_bot']))
+ else:
+ win_below.append(select_band_window(
+ kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
+
+ if ifSO or len(spins) == 2:
+ validate_so_windows(spins, info, windows, win_below)
+ elif not all(window[0] and window[1] <= window[2] for window in windows):
+ raise ValueError('projectors require a nonempty included window at every k-point')
+
+ # qbbot: point integration over bands below e_bot, is=1 only under SO
+ qbbot = 0.0
+ for ispin in range(ns):
+ for ik in range(nk):
+ incl, nb_bot, nb_top = win_below[ik]
+ if incl:
+ qbbot += (nb_top - nb_bot + 1) * spins[ispin]['kp'][ik]['weight']
+ if ifSO:
+ break
+
+ transmats = {}
+ mixing = {}
+ for cr in crorbs:
+ key = (cr['l'], cr['sort'])
+ if cr['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][cr['sort'] - 1]['sourcefile'], cr['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(cr['basis'], cr['l']), False
+ if not ifSO and any(mixing.values()):
+ raise NotImplementedError('spin-mixing bases require spin-orbit coupling')
+
+ # rotloc Euler angles per crorb: the symmetry op mapping the representative
+ # atom of the sort onto this atom (set_rotloc.f). With identity struct
+ # local rotation, rotloc%(a,b,g,iprop) are that op's Euler angles, and
+ # rot_projectmat applies dmat(l, a, b, g, iprop) in the |lm> basis
+ # (rot_projectmat.f:59-65). For atom 2 this is a non-trivial in-plane
+ # rotation (a=270) that mixes the m-rows; omitting it flips the m=1,2 rows.
+ # ---- raw correlated projector mat_rep, then Loewdin orthonormalize ----
+ # Non-mixing: mat_rep[(icr, ik, is)] -> (2l+1, nbsel) per spin block.
+ # Mixing: mat_rep[(icr, ik)] -> (2*(2l+1), nbsel), the stacked spinor block.
+ mat_rep = {}
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ atom = cr['atom']
+ sort = cr['sort']
+ d = 2 * l + 1
+ transmat = transmats[(l, sort)]
+ ismix = mixing[(l, sort)]
+ ref = rotloc_ref[sort - 1]
+ rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ nbsel = nb_top - nb_bot + 1
+ spin_blocks = []
+ for ispin in range(ns):
+ sp = spins[ispin]
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ nlo = sp['nLO'][(l, sort)]
+ P = np.zeros((d, nbsel), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l) # 0-based packed index
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ val = kp['Alm'][lm, atom, off]
+ for ilo in range(nlo):
+ val += kp['Clm'][ilo, lm, atom, off] * \
+ sp['ovl_LO_u'][(ilo + 1, l, sort)]
+ P[mi, j] = val
+ spin_blocks.append(rot @ P) # rot_projectmat per spin
+ if ismix:
+ stack = np.vstack(spin_blocks) # (2*(2l+1), nbsel), up then dn
+ mat_rep[(icr, ik)] = transmat @ stack
+ else:
+ for ispin in range(ns):
+ mat_rep[(icr, ik, ispin)] = transmat @ spin_blocks[ispin]
+
+ # Loewdin orthonormalization. Under SO (orthogonal_wannier_SO) the up+dn
+ # blocks of each crorb are stacked together and a single overlap is
+ # orthonormalized per k. Without SO (orthogonal_wannier) spin is a good
+ # quantum number, so each spin is orthonormalized independently over the
+ # stack of (2l+1)-wide crorb blocks of that spin.
+ if ifSO:
+ ortho_groups = [[(ik, None)] for ik in range(nk)]
+ else:
+ ortho_groups = [[(ik, ispin)] for ik in range(nk) for ispin in range(ns)]
+ for group in ortho_groups:
+ (ik, gspin) = group[0]
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ blocks = []
+ layout = [] # (key, nrows)
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if ifSO and mixing[(l, cr['sort'])]:
+ blocks.append(mat_rep[(icr, ik)])
+ layout.append(((icr, ik), 2 * (2 * l + 1)))
+ elif ifSO:
+ for ispin in range(ns):
+ blocks.append(mat_rep[(icr, ik, ispin)])
+ layout.append(((icr, ik, ispin), 2 * l + 1))
+ else:
+ blocks.append(mat_rep[(icr, ik, gspin)])
+ layout.append(((icr, ik, gspin), 2 * l + 1))
+ D = np.vstack(blocks) # ndim x nbnd
+ # match the Fortran's exact BLAS calls (orthogonal_wannier[_SO]): the
+ # near-singular O^{-1/2} amplifies any last-bit difference, so use the
+ # identical ZGEMM trans flags rather than numpy's conj-transpose copies.
+ O = zgemm(1.0, D, D, trans_b=2) # ZGEMM('N','C'): D @ D^H
+ S = _sqrt_inv(O)
+ D_orth = zgemm(1.0, S, D) # ZGEMM('N','N'): O^{-1/2} @ D
+ row = 0
+ for key, nrows in layout:
+ mat_rep[key] = D_orth[row:row + nrows, :]
+ row += nrows
+
+ # ---- Rloc rotrep per crorb ----
+ rloc_blocks = []
+ for cr in crorbs:
+ l = cr['l']
+ transmat = transmats[(l, cr['sort'])]
+ iref = sum(info['mult'][:cr['sort'] - 1]) + 1
+ if mixing[(l, cr['sort'])]:
+ rloc_blocks.append(_rloc_rotrep_mixing(
+ l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat))
+ else:
+ rloc_blocks.append(_rloc_rotrep(
+ l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat, ifSO))
+
+ # ---- write ----
+ with open(case + '.ctqmcout', 'w') as f:
+ f.write('13.605698\n')
+ f.write('%6d\n' % nk)
+ f.write('%6d\n' % (1 if ifSP else 0))
+ f.write('%6d\n' % (1 if ifSO else 0))
+ f.write(' %s\n' % fmt_scalar(qbbot))
+ f.write(' %s\n' % fmt_scalar(elecn))
+
+ f.write('%6d\n' % norb)
+ for o in orbs:
+ dim = 2 * (2 * o['l'] + 1) if ifSO else 2 * o['l'] + 1
+ f.write('%6d %6d %6d %6d \n' % (o['atom'], o['sort'], o['l'], dim))
+
+ f.write('%6d\n' % ncrorb)
+ for cr in crorbs:
+ l = cr['l']
+ size = 2 * (2 * l + 1) if ifSO else 2 * l + 1
+ f.write('%6d %6d %6d %6d %6d %6d \n' %
+ (cr['atom'], cr['sort'], l, size, 1 if ifSO else 0, 1))
+
+ # Rloc block per crorb. SO writes the 2*(2l+1) spinor whole shell;
+ # non-SO the bare (2l+1) shell. The time-reversal flag follows under
+ # SO always, under non-SO only when ifSP (outputqmc.f:347-351).
+ for rotrep, timeinv in rloc_blocks:
+ for m in range(rotrep.shape[0]):
+ write_row(f, rotrep[m, :].real)
+ for m in range(rotrep.shape[0]):
+ write_row(f, rotrep[m, :].imag)
+ if ifSO or ifSP:
+ f.write('%6d\n' % (1 if timeinv else 0))
+
+ # complex-harmonics -> basis transform block (crorb%first only). Mixing
+ # writes the full 2(2l+1) transmat; non-mixing the spin block-diagonal.
+ for cr in crorbs:
+ if not cr['first']:
+ continue
+ l = cr['l']
+ transmat = transmats[(l, cr['sort'])]
+ d = 2 * l + 1
+ if mixing[(l, cr['sort'])]:
+ spinrot = transmat
+ elif ifSO:
+ spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
+ spinrot[:d, :d] = transmat
+ spinrot[d:, d:] = transmat
+ else:
+ spinrot = transmat # bare (2l+1) transform, non-SO
+ dim = spinrot.shape[0]
+ f.write('%6d %6d \n' % (1, dim))
+ for m in range(dim):
+ write_row(f, spinrot[m, :].real)
+ for m in range(dim):
+ write_row(f, spinrot[m, :].imag)
+
+ # number of bands per k (skip is=2 under SO)
+ for ispin in range(ns):
+ if ifSP and ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ f.write('%6d\n' % abs(nb_top - nb_bot + 1))
+
+ # projector block: DO ik, DO icrorb. Mixing writes the full 2(2l+1)
+ # block from is=1; non-mixing the two (2l+1) spin blocks.
+ for ik in range(nk):
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ P = mat_rep[(icr, ik)]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # k-weights
+ for ik in range(nk):
+ f.write(' %s\n' % fmt_scalar(spins[0]['kp'][ik]['weight']))
+
+ # H(k) eigenvalues (skip is=2 under SO)
+ for ispin in range(ns):
+ if ifSP and ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ kp = spins[ispin]['kp'][ik]
+ for ib in range(nb_bot, nb_top + 1):
+ f.write(' %s\n' % fmt_scalar(kp['eband'][ib - kp['nbmin']]))
diff --git a/python/triqs_dftkit/wien2k/oubwin.py b/python/triqs_dftkit/wien2k/oubwin.py
new file mode 100644
index 0000000..a46e13c
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/oubwin.py
@@ -0,0 +1,127 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# 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 .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj band-window file (case.oubwin),
+replacing the corresponding output of the dmftproj Fortran executable.
+
+Given the alm/blm coefficient file (case.almblm, or the spin pair
+case.almblmup / case.almblmdn) and the projector definition (case.indmftpr),
+this selects, per k-point and spin, the contiguous band range whose Fermi-
+shifted Kohn-Sham eigenvalues fall in the energy window, and writes the
+result in the case.oubwin format read by the charge self-consistency.
+
+It reproduces the Fortran path: read elecn/nk/nloat/eferm from the almblm
+header (outbwin.f / dmftproj.f), parse the energy window (e_bot, e_top,
+proj_mode) from the last line of case.indmftpr, shift every band eigenvalue
+by eferm, and apply the set_projections.f selection. For proj_mode==0 this is
+the energy-window rule (strict lower bound e_bot < E, inclusive upper bound
+E <= e_top, yielding a single contiguous [nb_bot, nb_top] per k). For
+proj_mode 1 and 2 the window is a pair of band indices (b_bot, b_top): every
+k-point is included with nb_bot/nb_top those indices clamped to the local
+nbmin/nbmax (set_projections.f:70-88). In mode 2 the indices come straight
+from the window line; in mode 1 they are the global min/max band index over
+all spins/k whose Fermi-shifted energy lies in (e_bot, e_top]
+(dmftproj.f:704-722). The weight written per included k-point is the
+tetrahedron weight of the lowest band nbmin.
+
+Spin-polarization is detected from the input files: case.almblmup /
+case.almblmdn present => two spin files (case.oubwinup / case.oubwindn),
+otherwise the single case.oubwin. The spin-orbit flag (ifSO) comes from
+case.indmftpr; with -so the up and dn windows must coincide and dmftproj
+aborts otherwise, a check this generator also performs.
+"""
+
+import os
+
+from ._dmftproj import (band_index_window, read_almblm, read_indmftpr,
+ select_band_window, select_window)
+
+
+def _windows_from_spin(sp, info, band_window):
+ """Per-k (included, nb_bot, nb_top, weight) for one already-read spin dict.
+ proj_mode==0 uses the energy window; modes 1/2 use the band-index window
+ (b_bot, b_top) precomputed in band_window."""
+ out = []
+ for kp in sp['kp']:
+ if info['proj_mode'] == 0:
+ incl, nb_bot, nb_top = select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top'])
+ else:
+ incl, nb_bot, nb_top = select_band_window(
+ kp['nbmin'], kp['nbmax'], band_window[0], band_window[1])
+ out.append((incl, nb_bot, nb_top, kp['weight']))
+ return out
+
+
+def _write_oubwin_file(path, ifso, windows):
+ """Write one case.oubwin spin file (outbwin.f). Integer fields are i6;
+ the per-k weight is a list-directed real. Not-included k-points emit only
+ the flag line."""
+ with open(path, 'w') as f:
+ f.write('%6d\n' % len(windows))
+ f.write('%6d\n' % (1 if ifso else 0))
+ for incl, nb_bot, nb_top, weight in windows:
+ f.write('%6d\n' % (1 if incl else 0))
+ if incl:
+ f.write('%6d%6d\n' % (nb_bot, nb_top))
+ f.write(' %.16f \n' % weight)
+
+
+def write_oubwin(case):
+ """Read .almblm{up,dn} (or .almblm) and write the matching
+ .oubwin{up,dn} (or .oubwin) in the dmftproj format.
+
+ Spin-polarization is inferred from the input files. The spin-orbit flag
+ written into every file comes from case.indmftpr; under -sp+-so the up/dn
+ windows must coincide (dmftproj aborts otherwise)."""
+ info = read_indmftpr(case + '.indmftpr')
+ ifso = bool(info['so'])
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ spin_polarized = os.path.exists(up) and os.path.exists(dn)
+
+ if ifso and not spin_polarized:
+ raise ValueError('spin-orbit windows require both up and down almblm files')
+ if spin_polarized:
+ sp_up = read_almblm(up, info)
+ sp_dn = read_almblm(dn, info)
+ # mode 1 scans both spins for the global band-index window.
+ bw = (band_index_window(info, [sp_up, sp_dn])
+ if info['proj_mode'] != 0 else None)
+ win_up = _windows_from_spin(sp_up, info, bw)
+ win_dn = _windows_from_spin(sp_dn, info, bw)
+ if ifso:
+ if len(win_up) != len(win_dn):
+ raise ValueError('up/down k-point counts differ')
+ for u, d in zip(win_up, win_dn):
+ if u[0] != d[0] or u[1] != d[1] or u[2] != d[2]:
+ raise ValueError(
+ 'spin-orbit run requires identical up/dn band '
+ 'windows at every k-point')
+ _write_oubwin_file(case + '.oubwinup', ifso, win_up)
+ _write_oubwin_file(case + '.oubwindn', ifso, win_dn)
+ else:
+ sp = read_almblm(case + '.almblm', info)
+ bw = (band_index_window(info, [sp])
+ if info['proj_mode'] != 0 else None)
+ win = _windows_from_spin(sp, info, bw)
+ _write_oubwin_file(case + '.oubwin', ifso, win)
diff --git a/python/triqs_dftkit/wien2k/outband.py b/python/triqs_dftkit/wien2k/outband.py
new file mode 100644
index 0000000..4a8acd1
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/outband.py
@@ -0,0 +1,331 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# 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 .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj band-structure projector file
+(case.outband), replacing the -band path of the dmftproj Fortran executable
+(outband.f).
+
+case.outband feeds the converter method convert_bands_input: it is the
+band-structure analog of case.ctqmcout. The correlated-shell projectors are
+built exactly as in ctqmcout.py (raw Alm/Clm projector -> rot_projectmat local
+rotation -> cubic/fromfile transmat -> Loewdin orthonormalization
+orthogonal_wannier_SO), but evaluated on a band k-path and written in the
+outband layout. case.outband also carries the raw radial-normalized Theta
+projector pr_orb%matn_rep for every included shell (the parproj.py path, no
+Loewdin) and the k-path labels.
+
+Two band-mode differences from the ctqmcout path (dmftproj.f:557-581):
+
+- nkband, the number of k-points along the plotted k-path, is read from
+ .klist_band (read_k_list.f). It is written in the header and drives the
+ k-label block; the per-k projector/band data loop over the almblm nk, which is
+ the single representative k-point lapw2 -almd writes in band mode.
+- the Fermi energy is read from the LAST line of .indmftpr, not from the
+ almblm header; the band almblm keeps the eferm record as a placeholder that
+ the reader skips (READ(iualmblm,*) with no target).
+
+outband layout is outband.f:56-278: nkband, the per-spin per-k band counts, the
+Loewdin correlated projectors (DO ik, DO icrorb), H(k) eigenvalues, the
+included-orbital radial sizes, the raw Theta projectors, and the k-labels.
+"""
+
+import os
+import tempfile
+import numpy as np
+from scipy.linalg.blas import zgemm
+
+from ._dmftproj import (validate_so_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
+ read_indmftpr, reptrans, rotloc_orbital, select_window, to_float,
+ write_row)
+from .ctqmcout import _build_crorbs, _sqrt_inv
+from .parproj import (_build_matn_rep, _build_matn_rep_mixing, _build_orbs)
+
+
+def _read_indmftpr_band(indmftpr):
+ """case.indmftpr in band mode: the standard parse, plus the Fermi energy
+ appended as the last physical line after the (e_bot e_top proj_mode) window
+ line. The trailing scalar would be mistaken for the window line, so feed
+ read_indmftpr the file without that record and return eferm separately."""
+ raw = [l.split('!')[0].rstrip('\n') for l in open(indmftpr)]
+ nonblank = [i for i, l in enumerate(raw) if l.strip() != '']
+ eferm = to_float(raw[nonblank[-1]].split()[0])
+ tmp = tempfile.NamedTemporaryFile('w', suffix='.indmftpr', delete=False,
+ dir=os.path.dirname(os.path.abspath(indmftpr)))
+ tmp.write('\n'.join(raw[:nonblank[-1]]) + '\n')
+ tmp.close()
+ try:
+ info = read_indmftpr(tmp.name)
+ finally:
+ os.unlink(tmp.name)
+ return info, eferm
+
+
+def _read_almblm_band(path, info, eferm):
+ """Band-mode almblm read: identical to read_almblm but the header eferm
+ record is a placeholder the band path ignores; the Fermi shift uses the
+ eferm passed in (from case.indmftpr, dmftproj.f:576-581). read_almblm reads
+ the 4th record as eferm, so overwrite that record with the indmftpr value."""
+ text = open(path).read().splitlines()
+ text[3] = ' %.16E' % eferm # Reader record 4 (physical line index 3)
+ tmp = tempfile.NamedTemporaryFile('w', suffix='.almblm', delete=False)
+ tmp.write('\n'.join(text) + '\n')
+ tmp.close()
+ try:
+ return read_almblm(tmp.name, info, projectors=True)
+ finally:
+ os.unlink(tmp.name)
+
+
+def read_k_list(path):
+ """nkband and the k-path labels from .klist_band (read_k_list.f).
+
+ Every physical line up to the END marker is a k-point; a line whose first
+ character is not a blank carries a label, whose position is the k-point
+ index and whose name is the leading non-blank token. Returns
+ (nkband, [(pos, name), ...])."""
+ lines = open(path).read().splitlines()
+ nkband = 0
+ labels = []
+ for line in lines:
+ if line[:3] == 'END':
+ break
+ nkband += 1
+ if line[:1] != ' ':
+ labels.append((nkband, line.split()[0]))
+ return nkband, labels
+
+
+def write_outband(case):
+ """Read .almblm{up,dn}, .indmftpr, .struct,
+ .dmftsym and .klist_band and write .outband in the
+ dmftproj band format."""
+ info, eferm = _read_indmftpr_band(case + '.indmftpr')
+ if not info['so']:
+ raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
+ if info['proj_mode'] != 0:
+ raise NotImplementedError('this output currently supports projection mode 0')
+ nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ if os.path.exists(up) and os.path.exists(dn):
+ spin_files = [up, dn]
+ else:
+ spin_files = [case + '.almblm']
+ ifSP = len(spin_files) == 2
+ ifSO = bool(info['so'])
+ ns = 2 if ifSP else 1
+
+ spins = [_read_almblm_band(p, info, eferm) for p in spin_files]
+ nk = spins[0]['nk']
+ info['nk'] = nk
+
+ nkband, labels = read_k_list(case + '.klist_band')
+ if nkband != nk:
+ raise ValueError('klist_band and almblm k-point counts differ')
+
+ crorbs = _build_crorbs(info)
+ orbs = _build_orbs(info)
+ ncrorb = len(crorbs)
+ norb = len(orbs)
+
+ # proj_mode 0 energy window [e_bot, e_top] per k (band path: proj_mode 0).
+ windows = []
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ windows.append(select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+
+ validate_so_windows(spins, info, windows)
+
+ transmats = {}
+ mixing = {}
+ for o in orbs:
+ key = (o['l'], o['sort'])
+ if key in transmats:
+ continue
+ if o['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][o['sort'] - 1]['sourcefile'], o['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(o['basis'], o['l']), False
+
+ # rot_projectmat local rotation: the op mapping the sort representative onto
+ # this atom (set_projections.f via rot_projectmat). Same as ctqmcout/parproj.
+ # ---- correlated projector mat_rep, then Loewdin orthonormalize ----
+ mat_rep = {}
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ atom = cr['atom']
+ sort = cr['sort']
+ d = 2 * l + 1
+ transmat = transmats[(l, sort)]
+ ismix = mixing[(l, sort)]
+ ref = rotloc_ref[sort - 1]
+ rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ nbsel = nb_top - nb_bot + 1
+ spin_blocks = []
+ for ispin in range(ns):
+ sp = spins[ispin]
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ nlo = sp['nLO'][(l, sort)]
+ P = np.zeros((d, nbsel), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l)
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ val = kp['Alm'][lm, atom, off]
+ for ilo in range(nlo):
+ val += kp['Clm'][ilo, lm, atom, off] * \
+ sp['ovl_LO_u'][(ilo + 1, l, sort)]
+ P[mi, j] = val
+ spin_blocks.append(rot @ P)
+ if ismix:
+ stack = np.vstack(spin_blocks)
+ mat_rep[(icr, ik)] = transmat @ stack
+ else:
+ for ispin in range(ns):
+ mat_rep[(icr, ik, ispin)] = transmat @ spin_blocks[ispin]
+
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ blocks = []
+ layout = []
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ blocks.append(mat_rep[(icr, ik)])
+ layout.append(((icr, ik), 2 * (2 * l + 1)))
+ else:
+ for ispin in range(ns):
+ blocks.append(mat_rep[(icr, ik, ispin)])
+ layout.append(((icr, ik, ispin), 2 * l + 1))
+ D = np.vstack(blocks)
+ # match the Fortran's exact BLAS calls (orthogonal_wannier_SO): the
+ # near-singular O^{-1/2} amplifies any last-bit difference, so use the
+ # identical ZGEMM trans flags rather than numpy's conj-transpose copies.
+ O = zgemm(1.0, D, D, trans_b=2) # ZGEMM('N','C'): D @ D^H
+ S = _sqrt_inv(O)
+ D_orth = zgemm(1.0, S, D) # ZGEMM('N','N'): O^{-1/2} @ D
+ row = 0
+ for key, nrows in layout:
+ mat_rep[key] = D_orth[row:row + nrows, :]
+ row += nrows
+
+ # ---- raw Theta projectors matn_rep per included orbital (parproj path) ----
+ matn_reps = {}
+ for o in orbs:
+ l = o['l']
+ transmat = transmats[(l, o['sort'])]
+ ref = rotloc_ref[o['sort'] - 1]
+ iref = sum(info['mult'][:o['sort'] - 1]) + 1
+ rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
+ if mixing[(l, o['sort'])]:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, windows, transmat, rot)
+ else:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, ns, windows, transmat, rot)
+
+ # ---- write ----
+ with open(case + '.outband', 'w') as f:
+ f.write('%6d\n' % nkband)
+
+ # number of bands per k (skip is=2 under SP+SO)
+ for ispin in range(ns):
+ if ifSP and ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ f.write('%6d\n' % abs(nb_top - nb_bot + 1))
+
+ # correlated projector block: DO ik, DO icrorb (Loewdin orthonormalized).
+ for ik in range(nk):
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ P = mat_rep[(icr, ik)]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # H(k) eigenvalues (skip is=2 under SO)
+ for ispin in range(ns):
+ if ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ kp = spins[ispin]['kp'][ik]
+ for ib in range(nb_bot, nb_top + 1):
+ f.write(' %.16E\n' % kp['eband'][ib - kp['nbmin']])
+
+ # included-orbital radial sizes norm_radf%n
+ for o in orbs:
+ n = spins[0]['nLO'][(o['l'], o['sort'])] + 2
+ f.write('%6d\n' % n)
+
+ # Theta projector block: DO iorb, DO ik, DO ir. Mixing writes the full
+ # 2(2l+1) block from is=1; non-mixing the two (2l+1) spin blocks.
+ for o in orbs:
+ l = o['l']
+ atom = o['atom']
+ n = spins[0]['nLO'][(l, o['sort'])] + 2
+ ismix = mixing[(l, o['sort'])]
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ for ir in range(n):
+ if ismix:
+ P = matn_reps[(atom, l)][ik][:, :, ir]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(ns):
+ P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # k-labels
+ for i, (pos, name) in enumerate(labels, start=1):
+ f.write('%6d%6d%s\n' % (i, pos, name))
diff --git a/python/triqs_dftkit/wien2k/parproj.py b/python/triqs_dftkit/wien2k/parproj.py
new file mode 100644
index 0000000..aa04389
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/parproj.py
@@ -0,0 +1,643 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# 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 .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj partial-projector file
+(case.parproj), replacing the parproj path of the dmftproj Fortran executable
+(outputqmc.f).
+
+This is the sibling of ctqmcout.py. Where ctqmcout covers the correlated
+shells (crorb, l_inc==2) and runs them through the Loewdin orthonormalization,
+parproj covers ALL included shells (orb, l_inc in {1,2}) and writes the RAW
+radial-normalized Theta projector. The three structural differences from
+ctqmcout are:
+
+1. No Loewdin orthonormalization. pr_orb%matn_rep is never touched by
+ orthogonal_wannier; the only normalization is the per-orbital radial
+ transform s12 = O_radial^{+1/2} (orthogonal_r with inv=.FALSE.) applied as
+ coeff.s12 per (m, ib) in set_projections.f:273-332.
+
+2. The Theta projector keeps the full radial vector coeff = [Alm, Blm, Clm...]
+ of length n = nLO+2; s12 is built from u_dot_norm, ovl_LO_u and the second
+ LO-overlap token ovl_LO_udot (set_projections.f:163-197). All n radial
+ channels are written separately (the ir loop).
+
+3. Three extra per-orbital blocks: the BZ-symmetrized, tetrahedron-weighted,
+ Rloc-rotated density matrix densmat (density.f + symmetrize_mat + rotdens_mat),
+ the per-orbital Rloc spinor rotrep (set_rotloc.f), and a time-reversal flag.
+
+The reference case CaOs2 is SP+SO (ifSP=ifSO=1, ns=2), so every per-orbital
+matrix is the full 2*(2l+1)=10 spin+orbital block. The exact cubic transmat and
+the rot_projectmat local rotation are applied to matn_rep before writing,
+exactly as in ctqmcout (dmftproj reads the same full-precision templates after
+the precision fix).
+
+outputqmc.f parproj writer: 897-1183. set_projections.f s12/projector:
+163-197/273-622. density.f Theta path: 588-915. symmetrize_mat.f /
+rot_dens.f / setsym.f / set_rotloc.f for the symmetrization and Rloc rotation.
+"""
+
+import math
+import os
+import numpy as np
+
+from ._dmftproj import (validate_so_windows, band_index_window, dmat, mixing_rotrep, read_almblm, read_dmftsym,
+ read_fromfile, read_indmftpr, reptrans,
+ rotloc_orbital, rotloc_rotl_so, select_band_window, select_window, tmat, write_row)
+
+
+def _op_timeinv(op, ifSO):
+ det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+ - op['krotm'][0, 1] * op['krotm'][1, 0])
+ return bool(ifSO and det2 < 0.0)
+
+
+def _sqrtm_real_sym(O):
+ """O^{+1/2} of a real symmetric matrix (orthogonal_r with inv=.FALSE.:
+ Z diag(sqrt(w)) Z^T, real part kept). dmftproj evaluates sqrt of the
+ eigenvalue as a complex sqrt (W_comp = CMPLX(W,0)); reproduce that with a
+ complex power so a negative eigenvalue gives i*sqrt(|w|)."""
+ w, Z = np.linalg.eigh(O)
+ D1 = Z * np.sqrt(w.astype(complex))
+ return (D1 @ Z.T).real
+
+
+# --- included orbital descriptors --------------------------------------------
+
+def _build_orbs(info):
+ orbs = []
+ for isort in range(1, info['nsort'] + 1):
+ sortinfo = info['sorts'][isort - 1]
+ for l in sortinfo['included_ls']:
+ for imu in range(1, info['mult'][isort - 1] + 1):
+ atom = sum(info['mult'][:isort - 1]) + imu
+ orbs.append(dict(atom=atom, sort=isort, l=l,
+ basis=sortinfo['basis']))
+ return orbs
+
+
+# --- s12 radial transform (set_projections.f:163-197) ------------------------
+
+def _build_s12(l, isrt, n, spin):
+ """O_radial^{+1/2} for one (l, sort, spin): the n x n radial overlap built
+ from u_dot_norm + ovl_LO_u + ovl_LO_udot, then its matrix square root."""
+ s = np.zeros((n, n))
+ s[0, 0] = 1.0
+ s[1, 1] = spin['u_dot_norm'][(l, isrt)]
+ for ilo in range(1, spin['nLO'][(l, isrt)] + 1):
+ s[1 + ilo, 1 + ilo] = 1.0
+ s[1 + ilo, 0] = spin['ovl_LO_u'][(ilo, l, isrt)]
+ s[0, 1 + ilo] = spin['ovl_LO_u'][(ilo, l, isrt)]
+ s[1 + ilo, 1] = spin['ovl_LO_udot'][(ilo, l, isrt)]
+ s[1, 1 + ilo] = spin['ovl_LO_udot'][(ilo, l, isrt)]
+ return _sqrtm_real_sym(s)
+
+
+# --- srot rotrep (setsym.f) --------------------------------------------------
+
+def _srot_rotrep_nonmixing(op, l, transmat, ifSP, ifSO):
+ """srot(isym)%rotrep(l,isrt)%mat, the (2l+1) up/up block of the symmetry op
+ in the new basis: transmat . D(R)_{lm} . transmat^H, with the orbital
+ time-reversal operator applied for the magnetic (timeinv) operations.
+
+ Returns (rotrep[2l+1,2l+1], timeinv, phase)."""
+ a, b, g, iprop = op['a'], op['b'], op['g'], op['iprop']
+ krotm = op['krotm']
+ det2 = krotm[0, 0] * krotm[1, 1] - krotm[0, 1] * krotm[1, 0]
+ if ifSP and ifSO:
+ timeinv = det2 < 0.0
+ phase = (g - a) if timeinv else (a + g)
+ else:
+ timeinv = False
+ phase = 0.0
+ rotl = dmat(l, a, b, g, float(iprop))
+ rotrep = transmat @ rotl @ np.conj(transmat.T)
+ if timeinv:
+ # timeinv_op in the new basis: reptrans T reptrans^T applied to conj.
+ tinv = transmat @ tmat(l) @ transmat.T
+ rotrep = tinv @ np.conj(rotrep)
+ return rotrep, timeinv, phase
+
+
+# --- rotloc rotrep (set_rotloc.f) under SP+SO, non-mixing --------------------
+
+def _rotloc_rotrep_nonso(orb, ops, info, ref, transmat):
+ """rotloc(iatom)%rotrep(l)%mat without SO: the bare (2l+1) Rloc rotation in
+ the new basis, transmat D(Rloc)_{lm} transmat^H with D(Rloc) the symmetry op
+ mapping the sort representative onto this atom (set_rotloc.f non-SO branch).
+ srot%timeinv is always false under non-SO, so timeinv=False."""
+ l = orb['l']
+ iref = sum(info['mult'][:orb['sort'] - 1]) + 1
+ D = rotloc_orbital(l, ref, ops, orb['atom'], iref)
+ return transmat @ D @ np.conj(transmat.T), False
+
+
+def _rotloc_rotrep_so(orb, ops, info, ref, transmat, mixing):
+ """rotloc(iatom)%rotrep(l)%mat, the full 2*(2l+1) Rloc spinor rotation in
+ the new basis under SP+SO (set_rotloc.f), plus timeinv flag.
+
+ ref is the representative-sort rotloc; rotloc_rotl_so composes the
+ representative spinor rotloc with the first symmetry op mapping the sort
+ representative onto this atom. The new-basis transform is the full transmat
+ for a mixing (spin-coupling) basis, and blkdiag(transmat, transmat) for a
+ spin-diagonal one."""
+ l = orb['l']
+ iref = sum(info['mult'][:orb['sort'] - 1]) + 1
+ d = 2 * l + 1
+ rotl, timeinv = rotloc_rotl_so(l, ref, ops, orb['atom'], iref)
+ if mixing:
+ S = transmat
+ else:
+ S = np.zeros((2 * d, 2 * d), dtype=complex)
+ S[:d, :d] = transmat
+ S[d:, d:] = transmat
+ rotrep = (S @ rotl @ S.T) if timeinv else (S @ rotl @ np.conj(S.T))
+ return rotrep, timeinv
+
+
+# --- projector matn_rep (set_projections.f Theta path) -----------------------
+
+def _build_matn_rep(orb, info, spins, ns, windows, transmat, rot):
+ """matn_rep[(ik, is)] -> (2l+1, nbsel, n): the raw radial-normalized Theta
+ projector with the local rotation and the basis transform applied per
+ radial channel (non-mixing SP+SO path, set_projections.f:591-621)."""
+ l = orb['l']
+ atom = orb['atom']
+ sort = orb['sort']
+ out = {}
+ for ik in range(info['nk']):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ nbsel = nb_top - nb_bot + 1
+ for ispin in range(ns):
+ sp = spins[ispin]
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ n = sp['nLO'][(l, sort)] + 2
+ s12 = _build_s12(l, sort, n, sp)
+ matn = np.zeros((2 * l + 1, nbsel, n), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l)
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ coeff = np.zeros(n, dtype=complex)
+ coeff[0] = kp['Alm'][lm, atom, off]
+ coeff[1] = kp['Blm'][lm, atom, off]
+ for ilo in range(n - 2):
+ coeff[2 + ilo] = kp['Clm'][ilo, lm, atom, off]
+ matn[mi, j, :] = coeff @ s12
+ # rot_projectmat then the basis transform, per radial channel.
+ for ir in range(n):
+ matn[:, :, ir] = transmat @ (rot @ matn[:, :, ir])
+ out[(ik, ispin)] = matn
+ return out
+
+
+def _raw_matn(orb, sp, ik, nb_bot, nb_top, rot):
+ """The rot_projectmat'd (2l+1, nbsel, n) Theta projector in the |lm> basis
+ for one spin, before the angular transform (set_projections.f matn_rep)."""
+ l, atom, sort = orb['l'], orb['atom'], orb['sort']
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ nbsel = off_top - off_bot + 1
+ n = sp['nLO'][(l, sort)] + 2
+ s12 = _build_s12(l, sort, n, sp)
+ matn = np.zeros((2 * l + 1, nbsel, n), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l)
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ coeff = np.zeros(n, dtype=complex)
+ coeff[0] = kp['Alm'][lm, atom, off]
+ coeff[1] = kp['Blm'][lm, atom, off]
+ for ilo in range(n - 2):
+ coeff[2 + ilo] = kp['Clm'][ilo, lm, atom, off]
+ matn[mi, j, :] = coeff @ s12
+ for ir in range(n):
+ matn[:, :, ir] = rot @ matn[:, :, ir]
+ return matn
+
+
+def _build_matn_rep_mixing(orb, info, spins, windows, transmat, rot):
+ """matn_rep[ik] -> (2*(2l+1), nbsel, n): the mixing Theta projector. Each
+ spin's rot_projectmat'd |lm> block is stacked (up then dn) and multiplied by
+ the full 2(2l+1) transmat per radial channel (set_projections.f:515-585)."""
+ l = orb['l']
+ d = 2 * l + 1
+ out = {}
+ for ik in range(info['nk']):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ up = _raw_matn(orb, spins[0], ik, nb_bot, nb_top, rot)
+ dn = _raw_matn(orb, spins[1], ik, nb_bot, nb_top, rot)
+ nbsel, n = up.shape[1], up.shape[2]
+ matn = np.zeros((2 * d, nbsel, n), dtype=complex)
+ matn[:d] = up
+ matn[d:] = dn
+ for ir in range(n):
+ matn[:, :, ir] = transmat @ matn[:, :, ir]
+ out[ik] = matn
+ return out
+
+
+def _orbital_densmat_mixing(orb, info, spins, windows, matn_rep):
+ """The single 2*(2l+1) raw density block for a mixing orbital, point-
+ integrated with the geometric k-weight (density.f:786-812): D = sum over k
+ and radial channels of matn_rep matn_rep^H, already in the new basis."""
+ l = orb['l']
+ d = 2 * l + 1
+ n = spins[0]['nLO'][(l, orb['sort'])] + 2
+ dens = np.zeros((2 * d, 2 * d), dtype=complex)
+ for ik in range(info['nk']):
+ incl, _, _ = windows[ik]
+ if not incl:
+ continue
+ weight = spins[0]['kp'][ik]['weight']
+ for i in range(n):
+ mat = matn_rep[(orb['atom'], orb['l'])][ik][:, :, i]
+ dens += (mat @ np.conj(mat.T)) * weight
+ return dens
+
+
+def _symmetrize_densmat_mixing(orbs, ops, nsym, rotreps, dens_raw):
+ """Scatter each mixing shell to its symmetry-permuted atom before averaging."""
+ out = {o['atom']: np.zeros_like(dens_raw[o['atom']]) for o in orbs}
+ for o in orbs:
+ atom = o['atom']
+ for isym, op in enumerate(ops):
+ rotrep, timeinv = rotreps[isym]
+ raw = dens_raw[atom]
+ tmp = np.conj(raw) if timeinv else raw
+ target = op['perm'][atom - 1]
+ out[target] += rotrep @ tmp @ rotrep.conj().T
+ return {atom: value / nsym for atom, value in out.items()}
+
+
+def _rotdens_densmat_mixing(blk, rotrep_loc, timeinv):
+ """rotdens_mat for a mixing orbital (rot_dens.f:119-141): inverse(Rloc) D
+ Rloc on the single 2*(2l+1) block."""
+ if timeinv:
+ return rotrep_loc.T @ np.conj(blk @ rotrep_loc)
+ return np.conj(rotrep_loc.T) @ (blk @ rotrep_loc)
+
+
+# --- density matrix (density.f Theta path, non-mixing SP+SO) -----------------
+
+def _orbital_densmat_blocks(orb, info, spins, ns, windows, matn_reps):
+ """The nsp=4 raw density-matrix blocks (up/up, dn/dn, up/dn, dn/up) for one
+ orbital, point-integrated with the geometric k-weight (density.f:889-907,
+ tetr=.FALSE. path). The window is the bands below e_bot (dmftproj.f:798)."""
+ l = orb['l']
+ d = 2 * l + 1
+ blocks = [np.zeros((d, d), dtype=complex) for _ in range(4)]
+ for ik in range(info['nk']):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ kp0 = spins[0]['kp'][ik]
+ off_bot = nb_bot - kp0['nbmin']
+ off_top = nb_top - kp0['nbmin']
+ n = spins[0]['nLO'][(l, orb['sort'])] + 2
+ for iss in range(1, 5):
+ if iss <= 2:
+ is_, is1 = iss, iss
+ else:
+ is_ = iss - 2
+ is1 = 3 - is_
+ isp, isp1 = is_ - 1, is1 - 1
+ weight = spins[isp]['kp'][ik]['weight']
+ acc = np.zeros((d, d), dtype=complex)
+ for i in range(n):
+ mat = matn_reps[(orb['atom'], orb['l'])][(ik, isp)][:, :, i]
+ cmat = matn_reps[(orb['atom'], orb['l'])][(ik, isp1)][:, :, i]
+ acc += mat @ np.conj(cmat.T)
+ blocks[iss - 1] += acc * weight
+ return blocks
+
+
+def _symmetrize_densmat(orbs, info, ops, nsym, srot_rotreps, dens_raw, ns,
+ ifSP, ifSO):
+ """symmetrize_mat for the non-mixing SP+SO path (symmetrize_mat.f:196-276).
+ dens_raw[(atom, l)] -> [4 blocks]. Returns symmetrized blocks per shell.
+
+ The blocks are summed over symmetry ops with srot%rotrep, the spin block
+ phase ephase (up/dn, dn/up), and the orbital-time-reversal conjugation for
+ magnetic ops; equivalent atoms of a sort scatter into one another via the
+ permutation. Result divided by nsym."""
+ out = {}
+ # group orbitals by sort, in orb order (they are already sort-major).
+ isort_groups = {}
+ for o in orbs:
+ isort_groups.setdefault((o['sort'], o['l']), []).append(o)
+ for (isrt, l), group in isort_groups.items():
+ d = 2 * l + 1
+ mult = len(group)
+ sym = [[np.zeros((d, d), dtype=complex) for _ in range(4)]
+ for _ in range(mult)]
+ for imult in range(mult):
+ iatom = group[imult]['atom']
+ for isym in range(nsym):
+ op = ops[isym]
+ rotrep, timeinv, phase = srot_rotreps[(l, isrt, isym)]
+ jorb = op['perm'][iatom - 1] - iatom + imult # 0-based target
+ for iss in range(1, 5):
+ is_ = iss if iss <= 2 else iss - 2
+ isp = is_ - 1
+ tmp = dens_raw[(iatom, l)][iss - 1].copy()
+ if ifSP and timeinv:
+ tmp = np.conj(tmp)
+ ephase = 1.0
+ if iss == 3:
+ ephase = np.exp(1j * phase)
+ elif iss == 4:
+ ephase = np.exp(-1j * phase)
+ val = rotrep @ (tmp @ np.conj(rotrep.T)) * ephase
+ sym[jorb][iss - 1] += val
+ for imult in range(mult):
+ iatom = group[imult]['atom']
+ out[(iatom, l)] = [sym[imult][k] / nsym for k in range(4)]
+ return out
+
+
+def _symmetrize_densmat_nonso(orbs, info, ops, nsym, srot_rotreps, dens_raw):
+ """symmetrize_mat for the non-SO non-mixing path (symmetrize_mat.f:196-262,
+ nsp=2). Only the two diagonal blocks (up/up, dn/dn) exist; srot%timeinv is
+ always false so the ephase phase factors are 1. Each block is summed over
+ symmetry ops as rotrep block rotrep^H, equivalent atoms scattering via the
+ permutation, divided by nsym. Returns per-atom [up/up, dn/dn]."""
+ out = {}
+ isort_groups = {}
+ for o in orbs:
+ isort_groups.setdefault((o['sort'], o['l']), []).append(o)
+ for (isrt, l), group in isort_groups.items():
+ d = 2 * l + 1
+ mult = len(group)
+ sym = [[np.zeros((d, d), dtype=complex) for _ in range(2)]
+ for _ in range(mult)]
+ for imult in range(mult):
+ iatom = group[imult]['atom']
+ for isym in range(nsym):
+ op = ops[isym]
+ rotrep, _, _ = srot_rotreps[(l, isrt, isym)]
+ jorb = op['perm'][iatom - 1] - iatom + imult
+ for iss in range(2):
+ blk = dens_raw[(iatom, l)][iss]
+ sym[jorb][iss] += rotrep @ (blk @ np.conj(rotrep.T))
+ for imult in range(mult):
+ iatom = group[imult]['atom']
+ out[(iatom, l)] = [sym[imult][k] / nsym for k in range(2)]
+ return out
+
+
+def _rotdens_densmat_nonso(blocks, rotrep_loc):
+ """rotdens_mat for the non-SO non-mixing path (rot_dens.f:197-228, nsp=2):
+ each (2l+1) spin block is rotated as conj(Rloc^T) D Rloc with the bare
+ (2l+1) rotloc rotrep (timeinv always false). Returns [up/up, dn/dn]."""
+ return [np.conj(rotrep_loc.T) @ (blk @ rotrep_loc) for blk in blocks]
+
+
+def _rotdens_densmat(orb, blocks, rotrep_loc, timeinv):
+ """rotdens_mat for non-mixing SP+SO (rot_dens.f:151-193): assemble the four
+ blocks into a 2*(2l+1) matrix, apply inverse(Rloc) D Rloc, return it as the
+ full 2*(2l+1) densprint matrix."""
+ l = orb['l']
+ d = 2 * l + 1
+ rd = np.zeros((2 * d, 2 * d), dtype=complex)
+ rd[:d, :d] = blocks[0]
+ rd[d:, d:] = blocks[1]
+ rd[:d, d:] = blocks[2]
+ rd[d:, :d] = blocks[3]
+ if timeinv:
+ rd = np.conj(rd @ rotrep_loc)
+ rd = rotrep_loc.T @ rd
+ else:
+ rd = rd @ rotrep_loc
+ rd = np.conj(rotrep_loc.T) @ rd
+ return rd
+
+
+# --- parproj writer ----------------------------------------------------------
+
+def write_parproj(case):
+ """Read .almblm{up,dn}, .indmftpr, .struct,
+ .dmftsym and write .parproj in the dmftproj format."""
+ info = read_indmftpr(case + '.indmftpr')
+ nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ if os.path.exists(up) and os.path.exists(dn):
+ spin_files = [up, dn]
+ else:
+ spin_files = [case + '.almblm']
+ ifSP = len(spin_files) == 2
+ ifSO = bool(info['so'])
+ ns = 2 if ifSP else 1
+
+ spins = [read_almblm(p, info, projectors=True) for p in spin_files]
+ nk = spins[0]['nk']
+ info['nk'] = nk
+
+ orbs = _build_orbs(info)
+ norb = len(orbs)
+
+ # Two band ranges (dmftproj.f:798,836): the energy window [e_bot, e_top]
+ # for the written Theta projector (block A), and the bands below e_bot
+ # (-Elarge, e_bot) for the density matrix. dmftproj computes the density
+ # matrix over the below-e_bot range with point integration LAST (the third
+ # density call is correlated-only), so densmat in outputqmc is that one.
+ windows = []
+ below_windows = []
+ bw = band_index_window(info, spins) if info['proj_mode'] else None
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ if bw is None:
+ windows.append(select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+ below_windows.append(select_window(kp['nbmin'], kp['nbmax'],
+ kp['eband'], -1e6, info['e_bot']))
+ else:
+ windows.append(select_band_window(kp['nbmin'], kp['nbmax'], *bw))
+ below_windows.append(select_band_window(kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
+
+ if ifSO or len(spins) == 2:
+ validate_so_windows(spins, info, windows, below_windows)
+ elif not all(window[0] and window[1] <= window[2] for window in windows):
+ raise ValueError('projectors require a nonempty included window at every k-point')
+
+ transmats = {}
+ mixing = {}
+ for o in orbs:
+ key = (o['l'], o['sort'])
+ if o['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][o['sort'] - 1]['sourcefile'], o['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(o['basis'], o['l']), False
+ if not ifSO and any(mixing.values()):
+ raise NotImplementedError('spin-mixing bases require spin-orbit coupling')
+
+ # rot_projectmat local rotation: the op mapping the representative atom of
+ # the sort onto this atom (set_projections.f via rot_projectmat).
+ # ---- projectors matn_rep per orbital (energy window for block A, full
+ # band range for the density matrix) ----
+ matn_reps = {}
+ matn_reps_full = {}
+ for o in orbs:
+ l = o['l']
+ transmat = transmats[(l, o['sort'])]
+ ref = rotloc_ref[o['sort'] - 1]
+ iref = sum(info['mult'][:o['sort'] - 1]) + 1
+ rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
+ if mixing[(l, o['sort'])]:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, windows, transmat, rot)
+ matn_reps_full[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, below_windows, transmat, rot)
+ else:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, ns, windows, transmat, rot)
+ matn_reps_full[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, ns, below_windows, transmat, rot)
+
+ # ---- rotloc rotrep per orbital ----
+ rotloc_rotrep = {}
+ for o in orbs:
+ transmat = transmats[(o['l'], o['sort'])]
+ ref = rotloc_ref[o['sort'] - 1]
+ if not ifSO:
+ rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_nonso(
+ o, ops, info, ref, transmat)
+ else:
+ rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_so(
+ o, ops, info, ref, transmat, mixing[(o['l'], o['sort'])])
+
+ # ---- density matrices: raw -> symmetrize -> rotdens ----
+ # Non-mixing sorts use the block path; mixing sorts the single-block path.
+ # Under SO each non-mixing density is the 2*(2l+1) 4-block matrix; without
+ # SO it is two independent (2l+1) spin blocks (density.f nsp=2 path).
+ nonmix_orbs = [o for o in orbs if not mixing[(o['l'], o['sort'])]]
+ densprint = {}
+
+ if nonmix_orbs:
+ srot_rotreps = {}
+ for o in nonmix_orbs:
+ l, isrt = o['l'], o['sort']
+ for isym in range(nsym):
+ srot_rotreps[(l, isrt, isym)] = _srot_rotrep_nonmixing(
+ ops[isym], l, transmats[(l, isrt)], ifSP, ifSO)
+ dens_raw = {(o['atom'], o['l']): _orbital_densmat_blocks(
+ o, info, spins, ns, below_windows, matn_reps_full)
+ for o in nonmix_orbs}
+ if ifSO:
+ dens_sym = _symmetrize_densmat(nonmix_orbs, info, ops, nsym,
+ srot_rotreps, dens_raw, ns, ifSP, ifSO)
+ for o in nonmix_orbs:
+ key = (o['atom'], o['l'])
+ rotrep_loc, timeinv = rotloc_rotrep[key]
+ densprint[key] = _rotdens_densmat(
+ o, dens_sym[key], rotrep_loc, timeinv)
+ else:
+ dens_sym = _symmetrize_densmat_nonso(nonmix_orbs, info, ops, nsym,
+ srot_rotreps, dens_raw)
+ for o in nonmix_orbs:
+ key = (o['atom'], o['l'])
+ rotrep_loc, _ = rotloc_rotrep[key]
+ densprint[key] = _rotdens_densmat_nonso(
+ dens_sym[key], rotrep_loc)
+
+ mixing_groups = {}
+ for o in orbs:
+ key = (o['l'], o['sort'])
+ if mixing[key]:
+ mixing_groups.setdefault(key, []).append(o)
+ for (l, sort), group in mixing_groups.items():
+ rotreps = []
+ for op in ops:
+ ti = _op_timeinv(op, ifSO)
+ rotreps.append((mixing_rotrep(op, l, transmats[(l, sort)], ti), ti))
+ raw = {o['atom']: _orbital_densmat_mixing(
+ o, info, spins, below_windows, matn_reps_full) for o in group}
+ sym = _symmetrize_densmat_mixing(group, ops, nsym, rotreps, raw)
+ for o in group:
+ key = (o['atom'], l)
+ rotrep_loc, timeinv = rotloc_rotrep[key]
+ densprint[key] = _rotdens_densmat_mixing(sym[o['atom']], rotrep_loc, timeinv)
+
+ # ---- write ----
+ with open(case + '.parproj', 'w') as f:
+ for o in orbs:
+ n = spins[0]['nLO'][(o['l'], o['sort'])] + 2
+ f.write('%6d\n' % n)
+
+ for o in orbs:
+ l = o['l']
+ atom = o['atom']
+ n = spins[0]['nLO'][(l, o['sort'])] + 2
+ ismix = mixing[(l, o['sort'])]
+
+ # (A) Theta projector (outputqmc.f:935-973). Mixing writes the full
+ # 2(2l+1) block from is=1 only; non-mixing the two (2l+1) spin blocks.
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ for ir in range(n):
+ if ismix:
+ P = matn_reps[(atom, l)][ik][:, :, ir]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(ns):
+ P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # (B) density matrix (outputqmc.f:1012-1066). SO writes the single
+ # 2*(2l+1) block; mixing too; non-SO the two (2l+1) spin blocks
+ # (each real then imag, per spin, outputqmc.f:1052-1064).
+ dp = densprint[(atom, l)]
+ if ifSO or ismix:
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, dp[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, dp[m, :].imag)
+ else:
+ for ispin in range(ns):
+ for m in range(2 * l + 1):
+ write_row(f, dp[ispin][m, :].real)
+ for ispin in range(ns):
+ for m in range(2 * l + 1):
+ write_row(f, dp[ispin][m, :].imag)
+
+ # (C) Rloc rotrep (outputqmc.f:1130-1177). SO/mixing 2*(2l+1);
+ # non-SO bare (2l+1). The time-reversal flag follows under ifSP.
+ rotrep_loc, timeinv = rotloc_rotrep[(atom, l)]
+ for m in range(rotrep_loc.shape[0]):
+ write_row(f, rotrep_loc[m, :].real)
+ for m in range(rotrep_loc.shape[0]):
+ write_row(f, rotrep_loc[m, :].imag)
+ if ifSP:
+ f.write('%6d\n' % (1 if timeinv else 0))
diff --git a/python/triqs_dftkit/wien2k/sympar.py b/python/triqs_dftkit/wien2k/sympar.py
new file mode 100644
index 0000000..48708e1
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/sympar.py
@@ -0,0 +1,192 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
+#
+# 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 .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj included-shell symmetry file
+(case.sympar), replacing the corresponding output of the dmftproj Fortran
+executable (outputqmc.f, unit ousympar).
+
+case.sympar is the sibling of case.symqmc. Both write the same per-operation
+spinor symmetry matrix srot(isym)%rotrep(l,isrt)%mat (setsym.f), built
+identically for every shell as transmat . D(R)_lm . transmat^dag. They differ
+only in the shell list: symqmc covers the CORRELATED shells (l_inc==2), sympar
+covers ALL INCLUDED shells (l_inc in {1,2}), iterating orb(1:norb) instead of
+crorb(1:ncrorb). The orb ordering is (sort, l, atom) per dmftproj.f:357-387.
+
+The per-shell basis transmat is per sort: a complex shell uses the exact
+identity, a cubic shell uses the exact double-precision cubic-d harmonics (the
+same full-precision templates dmftproj reads after the precision fix), as
+ctqmcout does. Mixed within one file (Ca complex, Os cubic for CaOs2).
+
+sympar drops three things symqmc/ctqmcout carry: no orbital-description block
+at the top, no ifsplit/irep sub-block selection (always the full matrix), and
+no crorb%ifSOat/correp metadata. It is purely group-theoretic: header, perms,
+optional timeflag line (ifSP), the representation matrices, and (only when
+.not.ifSP) a paramagnetic time-reversal operator per orbital.
+
+This reference case (CaOs2) is SP+SO: indmftpr SO flag=1 sets ifSO, and
+ifSO=>ifSP, so it exercises the non-mixing 2*(2l+1)=10-wide block-diag spinor
+path and writes the timeflag line; the paramagnetic tail is absent.
+"""
+
+import os
+import numpy as np
+
+from ._dmftproj import (dmat, fmt, mixing_rotrep, read_dmftsym, read_fromfile,
+ read_indmftpr, reptrans, timeinv_orbital, tmat,
+ write_row)
+
+
+# --- included shells ---------------------------------------------------------
+
+def _included_shells(info):
+ """One entry per INCLUDED shell (l_inc in {1,2}), one per atom of its sort,
+ in orb order (sort, l, atom). Each shell carries l, basis name, a mixing
+ flag and the transform P: a spin-coupling fromfile basis gives the full
+ 2(2l+1) P (P spinrot P^dag matrix), otherwise the (2l+1) up/up block (the
+ single-precision-cast cubic harmonics dmftproj writes)."""
+ shells = []
+ for isort in range(info['nsort']):
+ s = info['sorts'][isort]
+ for l in s['included_ls']:
+ if s['basis'] == 'fromfile':
+ P, mixing = read_fromfile(s['sourcefile'], l)
+ else:
+ P, mixing = reptrans(s['basis'], l), False
+ for imu in range(1, info['mult'][isort] + 1):
+ atom = sum(info['mult'][:isort]) + imu
+ shells.append(dict(l=l, sort=isort + 1, atom=atom,
+ basis=s['basis'], P=P, mixing=mixing))
+ return shells
+
+
+# --- matrix construction (shared with symqmc) --------------------------------
+
+def _phase(op, ti):
+ a, c = op['a'], op['c']
+ return (c - a) if ti else (a + c) # (g-a) on magnetic ops, else (a+g)
+
+
+def _l0_matrix(op, ti):
+ e = np.exp(1j * _phase(op, ti) / 2)
+ return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)
+
+
+def _nonmixing_matrix(op, shell, ti):
+ """Non-mixing SP+SO whole shell: the up/up block scaled by +-(a+g)/2,
+ block-diagonal over spin, with the orbital time-reversal operator on the
+ magnetic operations (setsym.f, outputqmc.f:1292-1339)."""
+ l, P = shell['l'], shell['P']
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ if ti:
+ rotl = timeinv_orbital(l, rotl)
+ rotrep = P @ rotl @ P.T # antiunitary op: transpose, not dagger
+ else:
+ rotrep = P @ rotl @ np.conj(P.T)
+ e = np.exp(1j * _phase(op, ti) / 2)
+ d = 2 * l + 1
+ mat = np.zeros((2 * d, 2 * d), dtype=complex)
+ mat[:d, :d] = e * rotrep
+ mat[d:, d:] = np.conj(e) * rotrep
+ return mat
+
+
+def _nonmixing_orbital(op, shell):
+ """Non-SO spin-diagonal basis: the bare (2l+1) representation
+ P D(R)_{lm} P^H (setsym.f non-SO branch, outputqmc.f:1342-1349). Under
+ non-SO srot%timeinv is always false, so no time-reversal or phase."""
+ l, P = shell['l'], shell['P']
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ return P @ rotl @ np.conj(P.T)
+
+
+def _shell_matrix(op, shell, ti, ifSO):
+ if shell['mixing'] and not ifSO:
+ raise NotImplementedError('spin-mixing bases without spin-orbit coupling are not implemented')
+ l = shell['l']
+ if l == 0:
+ return _l0_matrix(op, ti)
+ if shell['mixing']:
+ return mixing_rotrep(op, l, shell['P'], bool(ti))
+ if not ifSO:
+ return _nonmixing_orbital(op, shell)
+ return _nonmixing_matrix(op, shell, ti)
+
+
+# --- output ------------------------------------------------------------------
+
+def _write_matrix(f, mat):
+ for m in range(mat.shape[0]):
+ write_row(f, mat[m, :].real)
+ for m in range(mat.shape[0]):
+ write_row(f, mat[m, :].imag)
+
+
+def write_sympar(case):
+ """Write .sympar from .dmftsym, .indmftpr, .struct.
+
+ Detects ifSP/ifSO from the presence of .almblm{up,dn} and the indmftpr
+ SO flag, exactly as oubwin/symqmc/ctqmcout. For SP+SO the file carries the
+ timeflag line and 2*(2l+1)-wide block-diag spinor matrices; the paramagnetic
+ time-reversal tail is written only when .not.ifSP."""
+ info = read_indmftpr(case + '.indmftpr')
+ shells, so = _included_shells(info), info['so']
+ nsym, ops = read_dmftsym(case + '.dmftsym')
+ natom = len(ops[0]['perm'])
+
+ ifSP = os.path.exists(case + '.almblmup') and os.path.exists(case + '.almblmdn')
+ ifSO = bool(so)
+ ifSP = ifSP or ifSO # ifSO => ifSP
+
+ timeflag = []
+ for op in ops:
+ det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+ - op['krotm'][0, 1] * op['krotm'][1, 0])
+ timeflag.append(1 if (ifSO and det2 < 0.0) else 0)
+
+ with open(case + '.sympar', 'w') as f:
+ f.write('%6d %6d\n' % (nsym, natom))
+ for op in ops:
+ f.write(''.join('%6d ' % p for p in op['perm']) + '\n')
+ if ifSP:
+ f.write(''.join('%6d ' % t for t in timeflag) + '\n')
+
+ for isym, op in enumerate(ops):
+ for sh in shells:
+ l = sh['l']
+ ti = bool(timeflag[isym])
+ if l == 0 and not (ifSP and ifSO):
+ f.write(fmt(1.0) + '\n')
+ f.write(fmt(0.0) + '\n')
+ continue
+ _write_matrix(f, _shell_matrix(op, sh, ti, ifSO))
+
+ if not ifSP:
+ for sh in shells:
+ l = sh['l']
+ if l == 0:
+ f.write(fmt(1.0) + '\n')
+ f.write(fmt(0.0) + '\n')
+ continue
+ tm = tmat(l)
+ op = sh['P'] @ tm @ sh['P'].T
+ ident = np.eye(2 * l + 1, dtype=complex)
+ time_op = op @ np.conj(ident)
+ _write_matrix(f, time_op)
diff --git a/python/triqs_dftkit/wien2k/symqmc.py b/python/triqs_dftkit/wien2k/symqmc.py
index bcc85ca..615abed 100644
--- a/python/triqs_dftkit/wien2k/symqmc.py
+++ b/python/triqs_dftkit/wien2k/symqmc.py
@@ -41,195 +41,112 @@
time-reversal operator -i sigma_y (x) T applied to the magnetic operations;
- l = 0 (s) shells: the 2x2 spin phase block.
-The mixing fromfile path is the one dft_tools #148 singles out. dmftproj reads
+The mixing fromfile path is the one dft_tools #148 singles out. Legacy dmftproj reads
both cubic and fromfile bases with single-precision CMPLX casts and a 250-column
line cap (Fortran set_ang_trans.f), and evaluates spin phases in single precision
(outputqmc.f, setsym.f). Its case.symqmc therefore differs at about 1e-7 from
this double-precision generator, whose cubic d basis uses exact coefficients.
"""
-import math
import os
+
import numpy as np
-# --- angular bases (transpose(P) = , m = -l..l) -----------------------
-
-_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
-
-# standard cubic harmonics, Wien2k convention
-_CUBIC = {
- 1: np.array([[0, 1, 0], [-1j, 0, -1j], [1, 0, -1]], dtype=complex) / 1.0,
- 2: np.array([
- [0, 0, 1, 0, 0],
- [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
- [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
- [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
- [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
- ], dtype=complex),
-}
-
-
-def _reptrans(basis, l):
- if basis == 'cubic' and l in _CUBIC:
- return _CUBIC[l]
- return _COMPLEX[l]
-
-
-# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
-
-def _small_d(l, m, n, b):
- f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
- (math.factorial(l + n) * math.factorial(l - n))
- s = 0.0
- for t in range(0, 2 * l + 1):
- if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
- f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
- (math.factorial(l - m - t) * math.factorial(m + n + t) *
- math.factorial(l - n - t) * math.factorial(t))
- f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
- f4 = 1.0 if (2 * t + n + m) == 0 else math.cos(b / 2) ** (2 * t + n + m)
- s += (-1) ** (l - m - t) * f2 * f3 * f4
- return math.sqrt(f1) * s
-
-
-def _dmat(l, a, b, c, det):
- D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
- for m in range(-l, l + 1):
- for n in range(-l, l + 1):
- v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b)
- if det < -0.5:
- v *= (-1) ** l
- D[m + l, n + l] = v
- return D
-
-
-def _tmat(l):
- """Complex-conjugation operator in the spherical-harmonic basis,
- T[m, -m] = (-1)^m (timeinv.f)."""
- T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
- for m in range(-l, l + 1):
- T[-m + l, m + l] = (-1) ** m
- return T
-
-
-def _timeinv_orbital(l, mat):
- return _tmat(l) @ np.conj(mat)
-
-
-# --- fromfile angular basis ---------------------------------------------------
-
-def _read_fromfile(path, l):
- """Parse a dmftproj fromfile basis (each line: the coefficients of a new
- basis vector in {|m,up>, |m,dn>}, m = -l..l, real/imag interleaved; '*'
- marks the end of an irep). Return (P, mixing) where P = is the
- transform matrix (reptrans.transmat), full 2(2l+1) for a spin-mixing basis
- or the (2l+1) up/up block otherwise."""
- n = 2 * (2 * l + 1)
- rows = []
- for line in open(path):
- line = line.rstrip('\n')
- if not line.strip():
- continue
- body = line[1:]
- vals = [float(x) for x in body.split()][:2 * n]
- rows.append([vals[2 * k] + 1j * vals[2 * k + 1] for k in range(n)])
- if len(rows) == n:
- break
- R = np.array(rows) # R[i, :] = |new_i> in old basis
- d = 2 * l + 1
- up_up, up_dn = R[:d, :d], R[:d, d:]
- dn_up, dn_dn = R[d:, :d], R[d:, d:]
- mixing = not (np.allclose(dn_dn, up_up) and
- np.allclose(up_dn, 0) and np.allclose(dn_up, 0))
- P = np.conj(R) if mixing else np.conj(up_up) # = conj()
- return P, mixing
-
-
-# --- input parsing -----------------------------------------------------------
-
-def _read_dmftsym(path):
- lines = open(path).read().split('\n')
- nsym = int(lines[0].split()[0])
- perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
- rest = lines[1 + nsym:]
- starts = [i for i, l in enumerate(rest) if 'Sym. op.' in l]
- ops = []
- for k, s in enumerate(starts):
- ang = rest[s + 1].split()
- a, b, c = (math.radians(float(x)) for x in ang[:3])
- krotm = np.array([[float(x) for x in rest[s + 2 + r].split()] for r in range(3)])
- ops.append(dict(perm=perms[k], a=a, b=b, c=c, krotm=krotm))
- return nsym, ops
+from ._dmftproj import (dmat, mixing_rotrep, read_dmftsym, read_fromfile,
+ read_indmftpr, reptrans, timeinv_orbital, tmat)
-def _read_correlated_shells(indmftpr, struct):
- """Return the list of correlated shells, one entry per correlated atom, plus
- the SO flag, from case.indmftpr and case.struct multiplicities. Each shell
- carries l, the basis name, its transform matrix P = and a mixing
- flag (True for a spin-coupling fromfile basis)."""
- raw = [l.split('!')[0].strip() for l in open(indmftpr)]
- raw = [l for l in raw if l != '']
- nsort = int(raw[0].split()[0])
- mult = [int(x) for x in raw[1].split()][:nsort]
- i = 3
- so = 0
+# --- correlated shells -------------------------------------------------------
+
+def _correlated_shells(info):
+ """One entry per correlated atom (l_inc==2), in sort order. Each shell
+ carries l, basis name, its transform matrix P = and a mixing flag
+ (True for a spin-coupling fromfile basis). symqmc builds its symmetry
+ matrices from the EXACT cubic harmonics (no single-precision cast)."""
shells = []
- for isort in range(nsort):
- basis = raw[i].split()[0]
- i += 1
- sourcefile = None
- if basis == 'fromfile':
- sourcefile = os.path.join(os.path.dirname(indmftpr), raw[i])
- i += 1
- l_inc = [int(x) for x in raw[i].split()]
- i += 1
- ireps = [int(x) for x in raw[i].split()]
- i += 1
- correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
- if any(n > 0 for n in ireps):
- i += 1 # skip the correps line
- if correlated_ls:
- so = int(raw[i].split()[0]) # SO flag follows a correlated sort
- i += 1
- for l in correlated_ls:
- if basis == 'fromfile':
- P, mixing = _read_fromfile(sourcefile, l)
- else:
- P, mixing = _reptrans(basis, l), False
- for _ in range(mult[isort]):
- shells.append(dict(l=l, basis=basis, P=P, mixing=mixing))
- return shells, so
+ for isort in range(info['nsort']):
+ s = info['sorts'][isort]
+ basis = s['basis']
+ for l in s['correlated_ls']:
+ if basis == 'fromfile':
+ P, mixing = read_fromfile(s['sourcefile'], l)
+ else:
+ P, mixing = reptrans(basis, l, cast=False), False
+ for _ in range(info['mult'][isort]):
+ shells.append(dict(l=l, basis=basis, P=P, mixing=mixing))
+ return shells
+
+
+# The symqmc test reaches into these two names directly; keep them as the
+# module's parsing entry points.
+_read_dmftsym = read_dmftsym
+
+
+def _read_correlated_shells(indmftpr, struct):
+ info = read_indmftpr(indmftpr)
+ return _correlated_shells(info), info['so']
def write_symqmc(case):
- """Write .symqmc from .dmftsym, .indmftpr, .struct."""
- nsym, ops = _read_dmftsym(case + '.dmftsym')
- shells, so = _read_correlated_shells(case + '.indmftpr', case + '.struct')
+ """Write .symqmc from .dmftsym, .indmftpr, .struct.
+
+ ifSO comes from the indmftpr SO flag; ifSP from the presence of
+ .almblm{up,dn} (ifSO => ifSP). With SO the matrices are the 2*(2l+1)
+ block-diag spinor representation and a timeflag line is written; without SO
+ they are the bare (2l+1) representation, with a timeflag line (all zero)
+ only when ifSP. The paramagnetic time-reversal tail follows when .not.ifSP."""
+ nsym, ops = read_dmftsym(case + '.dmftsym')
+ info = read_indmftpr(case + '.indmftpr')
+ shells, so = _correlated_shells(info), info['so']
natom = len(ops[0]['perm'])
+ ifSO = bool(so)
+ ifSP = (os.path.exists(case + '.almblmup')
+ and os.path.exists(case + '.almblmdn')) or ifSO
+
timeinv = []
for op in ops:
det2 = op['krotm'][0, 0] * op['krotm'][1, 1] - op['krotm'][0, 1] * op['krotm'][1, 0]
- timeinv.append(1 if (so and det2 < 0.0) else 0)
+ timeinv.append(1 if (ifSO and det2 < 0.0) else 0)
with open(case + '.symqmc', 'w') as f:
f.write('%6d %6d\n' % (nsym, natom))
for op in ops:
f.write(''.join('%6d ' % p for p in op['perm']) + '\n')
- if so:
+ if ifSP:
f.write(''.join('%6d ' % t for t in timeinv) + '\n')
for isym, op in enumerate(ops):
for sh in shells:
- f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym])))
+ f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym], ifSO)))
+ if not ifSP:
+ for sh in shells:
+ f.write(_format_matrix(_time_op_matrix(sh)))
-def _shell_matrix(op, shell, ti):
- """Spinor symmetry matrix for one correlated shell under one operation."""
+def _time_op_matrix(shell):
+ """Paramagnetic time-reversal operator for one correlated shell (.not.ifSP,
+ outputqmc.f:835-887): the (2l+1) operator P T P^T in the new basis; the s
+ shell reduces to the 1x1 identity."""
l = shell['l']
if l == 0:
- return _l0_matrix(op, ti)
+ return np.array([[1.0 + 0j]], dtype=complex)
+ P = shell['P']
+ return (P @ tmat(l) @ P.T) @ np.eye(2 * l + 1, dtype=complex)
+
+
+def _shell_matrix(op, shell, ti, ifSO):
+ """Symmetry matrix for one correlated shell under one operation: the
+ 2*(2l+1) block-diag spinor representation under SO, the bare (2l+1)
+ representation otherwise."""
+ if shell['mixing'] and not ifSO:
+ raise NotImplementedError('spin-mixing bases without spin-orbit coupling are not implemented')
+ l = shell['l']
+ if l == 0:
+ return _l0_matrix(op, ti) if ifSO else np.array([[1.0 + 0j]], dtype=complex)
if shell['mixing']:
return _mixing_matrix(op, shell, ti)
+ if not ifSO:
+ return _nonmixing_orbital(op, shell)
return _nonmixing_matrix(op, shell, ti)
@@ -248,10 +165,12 @@ def _nonmixing_matrix(op, shell, ti):
"""Spin-diagonal basis: the up/up block scaled by +-(a+g)/2, with the
orbital time-reversal operator on the magnetic operations."""
l, P = shell['l'], shell['P']
- rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
if ti:
- rotl = _timeinv_orbital(l, rotl)
- rotrep = P @ rotl @ np.conj(P.T)
+ rotl = timeinv_orbital(l, rotl)
+ rotrep = P @ rotl @ P.T # antiunitary op: transpose, not dagger
+ else:
+ rotrep = P @ rotl @ np.conj(P.T)
e = np.exp(1j * _phase(op, ti) / 2)
d = 2 * l + 1
mat = np.zeros((2 * d, 2 * d), dtype=complex)
@@ -260,30 +179,19 @@ def _nonmixing_matrix(op, shell, ti):
return mat
-def _mixing_matrix(op, shell, ti):
- """Spin-coupling basis: full 2(2l+1) spinor representation. The orbital
- rotation enters block-diagonal (beta=0) or block-antidiagonal (beta=pi, the
- magnetic operations); the spinor time-reversal -i sigma_y (x) T is then
- applied to the magnetic operations (setsym.f, timeinv.f)."""
+def _nonmixing_orbital(op, shell):
+ """Non-SO spin-diagonal basis: the bare (2l+1) representation
+ P D(R)_{lm} P^H (setsym.f non-SO branch). Under non-SO srot%timeinv is
+ always false, so no orbital time-reversal or spin phase enters."""
l, P = shell['l'], shell['P']
- rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
- e = np.exp(1j * _phase(op, ti) / 2)
- d = 2 * l + 1
- spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
- if ti: # beta = pi, block-antidiagonal
- spinrot[:d, d:] = e * rotl
- spinrot[d:, :d] = -np.conj(e) * rotl
- else: # beta = 0, block-diagonal
- spinrot[:d, :d] = e * rotl
- spinrot[d:, d:] = np.conj(e) * rotl
- rotrep = P @ spinrot @ np.conj(P.T)
- if ti:
- tm = _tmat(l)
- tinv = np.zeros((2 * d, 2 * d), dtype=complex)
- tinv[:d, d:] = -tm
- tinv[d:, :d] = tm
- rotrep = (P @ tinv @ P.T) @ np.conj(rotrep)
- return rotrep
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ return P @ rotl @ np.conj(P.T)
+
+
+def _mixing_matrix(op, shell, ti):
+ """Spin-coupling basis: full 2(2l+1) spinor representation (setsym.f,
+ timeinv.f), shared with sympar via _dmftproj.mixing_rotrep."""
+ return mixing_rotrep(op, shell['l'], shell['P'], bool(ti))
def _format_matrix(mat):
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index c43cbed..71df147 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
@@ -42,8 +46,88 @@ add_test(NAME Py_wien2k_symqmc_mixing
set_property(TEST Py_wien2k_symqmc_mixing APPEND PROPERTY ENVIRONMENT
PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+# Generator on a non-centrosymmetric spin-orbit cell (tellurium, p shell), the
+# dft_tools #148 regime the centrosymmetric d-shell fixtures never reached. The
+# reference is the dmftproj matrix data (single precision, ~10 kB).
+file(COPY Te.dmftsym Te.indmftpr Te.struct wien2k_symqmc_te.ref.npy
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+add_test(NAME Py_wien2k_symqmc_noncentrosym
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_noncentrosym.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_symqmc_noncentrosym APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Pure-Python case.oubwin band-window generator vs the dmftproj Fortran output.
+# Reads the existing compressed almblm inputs from the source fixtures and
+# compares both spin outputs against CaOs2.oubwinup in a scratch directory.
+add_test(NAME Py_wien2k_oubwin_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_oubwin_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_oubwin_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Pure-Python correlated-shell projector generator vs the native fixture.
+add_test(NAME Py_wien2k_ctqmcout_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_ctqmcout_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_ctqmcout_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Included-shell symmetry generator, including an uncorrelated Ca d shell.
+add_test(NAME Py_wien2k_sympar_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_sympar_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_sympar_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+add_test(NAME Py_wien2k_parproj_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_parproj_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_parproj_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
add_test(NAME Py_wien2k_kpath_labels
COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_kpath_labels.py
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set_property(TEST Py_wien2k_kpath_labels APPEND PROPERTY ENVIRONMENT
PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+add_test(NAME Py_wien2k_dmftproj_common_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_dmftproj_common_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_dmftproj_common_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+add_test(NAME Py_wien2k_outband_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_outband_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_outband_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Pure-Python generators on the non-spin-orbit path (ifSO=0, the common DMFT
+# case): dmftproj writes (2l+1) matrices, per-spin orthonormalization. Validated
+# against the precision-fixed dmftproj run without -so. Reuses CaOs2.struct/dmftsym.
+file(COPY CaOs2_noso.indmftpr CaOs2_noso.almblmup.gz CaOs2_noso.almblmdn.gz
+ CaOs2_noso.ctqmcout.gz CaOs2_noso.symqmc.gz
+ CaOs2_noso_partial.indmftpr CaOs2_noso_partial.sympar.gz CaOs2_noso_partial.parproj.gz
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+add_test(NAME Py_wien2k_noso_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_noso_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_noso_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Odd-l shells (p l=1, f l=3): the improper-rotation parity (-1)^l is invisible
+# for the even-l (d) fixtures, so these pin dmat's use of iprop. p uses the
+# standard almblm (wide window, full rank); f uses a high-band-cutoff almblm so
+# the high-energy window carries real l=3 character and the Loewdin overlap is
+# full rank. Both reproduce the precision-fixed dmftproj to machine precision.
+file(COPY CaOs2_pshell.indmftpr CaOs2_pshell.ctqmcout.gz CaOs2_pshell.symqmc.gz
+ CaOs2_fshell.indmftpr CaOs2_fshell.almblmup.gz CaOs2_fshell.almblmdn.gz
+ CaOs2_fshell.ctqmcout.gz CaOs2_fshell.symqmc.gz
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+add_test(NAME Py_wien2k_oddl_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_oddl_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_oddl_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
diff --git a/test/python/wien2k/CaOs2_band.almblmdn.gz b/test/python/wien2k/CaOs2_band.almblmdn.gz
new file mode 100644
index 0000000..0329e74
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_band.almblmup.gz b/test/python/wien2k/CaOs2_band.almblmup.gz
new file mode 100644
index 0000000..3eccf62
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_band.indmftpr b/test/python/wien2k/CaOs2_band.indmftpr
new file mode 100644
index 0000000..6aa6f67
--- /dev/null
+++ b/test/python/wien2k/CaOs2_band.indmftpr
@@ -0,0 +1,13 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+fromfile
+jbasis_d.dat
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
+0.6529001225
diff --git a/test/python/wien2k/CaOs2_band.klist_band b/test/python/wien2k/CaOs2_band.klist_band
new file mode 100644
index 0000000..366e653
--- /dev/null
+++ b/test/python/wien2k/CaOs2_band.klist_band
@@ -0,0 +1,2 @@
+GAMMA 0 0 0 8 2.0
+END
diff --git a/test/python/wien2k/CaOs2_band.outband.gz b/test/python/wien2k/CaOs2_band.outband.gz
new file mode 100644
index 0000000..063276c
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.outband.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.almblmdn.gz b/test/python/wien2k/CaOs2_fshell.almblmdn.gz
new file mode 100644
index 0000000..ff9c44c
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.almblmup.gz b/test/python/wien2k/CaOs2_fshell.almblmup.gz
new file mode 100644
index 0000000..0fbf89f
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.ctqmcout.gz b/test/python/wien2k/CaOs2_fshell.ctqmcout.gz
new file mode 100644
index 0000000..2561edd
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.indmftpr b/test/python/wien2k/CaOs2_fshell.indmftpr
new file mode 100644
index 0000000..9624e85
--- /dev/null
+++ b/test/python/wien2k/CaOs2_fshell.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+complex
+0 0 0 2
+0 0 0 0
+1
+3.0 10.0
diff --git a/test/python/wien2k/CaOs2_fshell.symqmc.gz b/test/python/wien2k/CaOs2_fshell.symqmc.gz
new file mode 100644
index 0000000..541d819
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_full.ctqmcout.gz b/test/python/wien2k/CaOs2_full.ctqmcout.gz
new file mode 100644
index 0000000..3cca2cc
Binary files /dev/null and b/test/python/wien2k/CaOs2_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_full.indmftpr b/test/python/wien2k/CaOs2_full.indmftpr
new file mode 100644
index 0000000..227d144
--- /dev/null
+++ b/test/python/wien2k/CaOs2_full.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz b/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz
new file mode 100644
index 0000000..c3c7a7c
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_jbasis_full.indmftpr b/test/python/wien2k/CaOs2_jbasis_full.indmftpr
new file mode 100644
index 0000000..dab1b92
--- /dev/null
+++ b/test/python/wien2k/CaOs2_jbasis_full.indmftpr
@@ -0,0 +1,12 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+fromfile
+jbasis_d.dat
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_jbasis_full.parproj.gz b/test/python/wien2k/CaOs2_jbasis_full.parproj.gz
new file mode 100644
index 0000000..9a60325
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_jbasis_full.sympar.gz b/test/python/wien2k/CaOs2_jbasis_full.sympar.gz
new file mode 100644
index 0000000..a9abae8
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_mode1.indmftpr b/test/python/wien2k/CaOs2_mode1.indmftpr
new file mode 100644
index 0000000..1cabe50
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+-0.15 0.30 1
diff --git a/test/python/wien2k/CaOs2_mode1.oubwindn b/test/python/wien2k/CaOs2_mode1.oubwindn
new file mode 100644
index 0000000..778773d
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.oubwindn
@@ -0,0 +1,5 @@
+ 1
+ 1
+ 1
+ 60 76
+ 1.0000000000000000
diff --git a/test/python/wien2k/CaOs2_mode1.oubwinup b/test/python/wien2k/CaOs2_mode1.oubwinup
new file mode 100644
index 0000000..778773d
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.oubwinup
@@ -0,0 +1,5 @@
+ 1
+ 1
+ 1
+ 60 76
+ 1.0000000000000000
diff --git a/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz b/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz
new file mode 100644
index 0000000..38b1ece
Binary files /dev/null and b/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_mode1_full.indmftpr b/test/python/wien2k/CaOs2_mode1_full.indmftpr
new file mode 100644
index 0000000..13e6507
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1_full.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0 1
diff --git a/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz b/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz
new file mode 100644
index 0000000..38b1ece
Binary files /dev/null and b/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_mode2_full.indmftpr b/test/python/wien2k/CaOs2_mode2_full.indmftpr
new file mode 100644
index 0000000..7732eb7
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode2_full.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+43 92 2
diff --git a/test/python/wien2k/CaOs2_noso.almblmdn.gz b/test/python/wien2k/CaOs2_noso.almblmdn.gz
new file mode 100644
index 0000000..10bba32
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.almblmup.gz b/test/python/wien2k/CaOs2_noso.almblmup.gz
new file mode 100644
index 0000000..8080ae5
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.ctqmcout.gz b/test/python/wien2k/CaOs2_noso.ctqmcout.gz
new file mode 100644
index 0000000..7c8b27d
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.indmftpr b/test/python/wien2k/CaOs2_noso.indmftpr
new file mode 100644
index 0000000..b54bf2f
--- /dev/null
+++ b/test/python/wien2k/CaOs2_noso.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+0
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_noso.symqmc.gz b/test/python/wien2k/CaOs2_noso.symqmc.gz
new file mode 100644
index 0000000..d18da41
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial.indmftpr b/test/python/wien2k/CaOs2_noso_partial.indmftpr
new file mode 100644
index 0000000..19b63bb
--- /dev/null
+++ b/test/python/wien2k/CaOs2_noso_partial.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 1 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+0
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_noso_partial.parproj.gz b/test/python/wien2k/CaOs2_noso_partial.parproj.gz
new file mode 100644
index 0000000..10156b7
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial.sympar.gz b/test/python/wien2k/CaOs2_noso_partial.sympar.gz
new file mode 100644
index 0000000..0547daa
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_partial.indmftpr b/test/python/wien2k/CaOs2_partial.indmftpr
new file mode 100644
index 0000000..afcca02
--- /dev/null
+++ b/test/python/wien2k/CaOs2_partial.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 1 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+-0.15 0.30
diff --git a/test/python/wien2k/CaOs2_partial.parproj.gz b/test/python/wien2k/CaOs2_partial.parproj.gz
new file mode 100644
index 0000000..6507818
Binary files /dev/null and b/test/python/wien2k/CaOs2_partial.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_partial.sympar.gz b/test/python/wien2k/CaOs2_partial.sympar.gz
new file mode 100644
index 0000000..e4e2f53
Binary files /dev/null and b/test/python/wien2k/CaOs2_partial.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_pshell.ctqmcout.gz b/test/python/wien2k/CaOs2_pshell.ctqmcout.gz
new file mode 100644
index 0000000..5669218
Binary files /dev/null and b/test/python/wien2k/CaOs2_pshell.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_pshell.indmftpr b/test/python/wien2k/CaOs2_pshell.indmftpr
new file mode 100644
index 0000000..db81ca3
--- /dev/null
+++ b/test/python/wien2k/CaOs2_pshell.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+complex
+0 2 0 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_pshell.symqmc.gz b/test/python/wien2k/CaOs2_pshell.symqmc.gz
new file mode 100644
index 0000000..29e2b7e
Binary files /dev/null and b/test/python/wien2k/CaOs2_pshell.symqmc.gz differ
diff --git a/test/python/wien2k/Te.dmftsym b/test/python/wien2k/Te.dmftsym
new file mode 100644
index 0000000..a58ddd0
--- /dev/null
+++ b/test/python/wien2k/Te.dmftsym
@@ -0,0 +1,49 @@
+ 6
+ 1 2 3
+ 2 3 1
+ 3 1 2
+ 1 3 2
+ 2 1 3
+ 3 2 1
+
+Sym. op. : 1
+ 0.0 0.0 0.0 1 - euler angles: a,b,c; iprop
+ 1.0000000000000000 0.0000000000000000 0.0000000000000000
+ 0.0000000000000000 0.99999999999999989 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 0.99999999999999989
+
+Sym. op. : 2
+ 240.0 0.0 0.0 1 - euler angles: a,b,c; iprop
+ -0.50000000000000000 -0.86602540378443849 0.0000000000000000
+ 0.86602540378443882 -0.49999999999999994 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 0.99999999999999989
+
+Sym. op. : 3
+ 120.0 0.0 0.0 1 - euler angles: a,b,c; iprop
+ -0.50000000000000000 0.86602540378443849 0.0000000000000000
+ -0.86602540378443882 -0.49999999999999994 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 0.99999999999999989
+
+Sym. op. : 4
+ 0.0 180.0 120.0 1 - euler angles: a,b,c; iprop
+ 0.50000000000000000 0.86602540378443849 0.0000000000000000
+ 0.86602540378443882 -0.49999999999999994 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 -0.99999999999999989
+
+Sym. op. : 5
+ 0.0 180.0 0.0 1 - euler angles: a,b,c; iprop
+ -1.0000000000000000 0.0000000000000000 0.0000000000000000
+ 0.0000000000000000 0.99999999999999989 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 -0.99999999999999989
+
+Sym. op. : 6
+ 0.0 180.0 240.0 1 - euler angles: a,b,c; iprop
+ 0.50000000000000000 -0.86602540378443849 0.0000000000000000
+ -0.86602540378443882 -0.49999999999999994 0.0000000000000000
+ 0.0000000000000000 0.0000000000000000 -0.99999999999999989
+ Global->local coordinates rotation matrices
+ 1
+ -0.0000000000000000 0.50000000000000000 0.86602540378443860
+ 0.0000000000000000 -0.86602540378443860 0.50000000000000000
+ 1.0000000000000000 0.0000000000000000 0.0000000000000000
+ 0.0 90.0 150.0 1 - euler angles: a,b,c; iprop
diff --git a/test/python/wien2k/Te.indmftpr b/test/python/wien2k/Te.indmftpr
new file mode 100644
index 0000000..f2fa44d
--- /dev/null
+++ b/test/python/wien2k/Te.indmftpr
@@ -0,0 +1,8 @@
+1
+3
+3
+cubic
+0 2 0 0
+0 0 0 0
+1
+-0.50 0.50
diff --git a/test/python/wien2k/Te.struct b/test/python/wien2k/Te.struct
new file mode 100644
index 0000000..89c1903
--- /dev/null
+++ b/test/python/wien2k/Te.struct
@@ -0,0 +1,37 @@
+ Tellur (P3121) s-o calc. M|| 0.00 0.00 1.00
+H 1 NEW
+ RELA
+ 9.410000 9.410000 11.199000 0.000000 0.000000 0.000000
+ATOM -1: X=0.76468000 Y=0.76468000 Z=0.00000000
+ MULT= 3 ISPLIT= 8
+ -1: X=0.23532000 Y=0.00000000 Z=0.33333330
+ -1: X=0.00000000 Y=0.23532000 Z=0.66666670
+Te NPT= 381 R0=.000005000 RMT= 2.10000 Z: 52.00000
+LOCAL ROT MATRIX: -0.0000000 0.5000000 0.8660254
+ 0.0000000-0.8660254 0.5000000
+ 1.0000000 0.0000000 0.0000000
+ 6 NUMBER OF SYMMETRY OPERATIONS
+ 1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0 1 0.00000000
+ 1 A 1 so. oper. type orig. index
+ 0-1 0 0.00000000
+ 1-1 0 0.00000000
+ 0 0 1 0.33333330
+ 2 A 3
+-1 1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0 1 0.66666670
+ 3 A 5
+ 0 1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0-1 0.00000000
+ 4 B 2
+-1 0 0 0.00000000
+-1 1 0 0.00000000
+ 0 0-1 0.33333330
+ 5 B 4
+ 1-1 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0-1 0.66666670
+ 6 B 6
diff --git a/test/python/wien2k/wien2k_ctqmcout_python.py b/test/python/wien2k/wien2k_ctqmcout_python.py
new file mode 100644
index 0000000..5c129f7
--- /dev/null
+++ b/test/python/wien2k/wien2k_ctqmcout_python.py
@@ -0,0 +1,100 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.ctqmcout generator (triqs_dftkit.wien2k.ctqmcout)
+# reproduces the dmftproj Fortran correlated-shell projectors to machine
+# precision: regenerate case.ctqmcout from the almblm/struct/dmftsym + a wide
+# energy window and compare to the precision-fixed dmftproj reference.
+#
+# The window (-2.0, 3.0) keeps 50 bands, more than the 20 correlated spin-orbitals
+# of the two SOC Os d shells, so the Loewdin overlap O is full rank.
+# Reference built with the precision-fixed dmftproj (#14) and full-precision
+# templates. The original narrow window is not a valid Loewdin oracle.
+#
+# CaOs2_jbasis_full exercises the spin-mixing fromfile (|j, m_j>) path: the Os d
+# shell is a 2(2l+1)=10 spinor basis (jbasis_d.dat), so the projector, Rloc
+# rotrep and complex-harmonics transform blocks are the full mixing matrices the
+# dft_tools #148 path singles out. Same wide window for the full-rank Loewdin.
+#
+# CaOs2_mode1_full and CaOs2_mode2_full exercise the band-index projection modes
+# (set_projections.f:70-88). Mode 1 ("-2.0 3.0 1") scans all spins/k for the
+# global band-index window; mode 2 ("43 92 2") takes the indices straight from
+# the window line. Both resolve to bands 43..92 (50 bands > 20 correlated
+# spin-orbitals), so the Loewdin overlap stays full rank and the projectors
+# match the precision-fixed dmftproj to machine precision.
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_lines):
+ got = open(got_path).read().split()
+ ref = ''.join(ref_lines).split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ floats = []
+ for a, b in zip(got, ref):
+ if b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ else:
+ floats.append((float(a), float(b)))
+ a, b = np.array(floats).T
+ np.testing.assert_allclose(a, b, atol=1e-11, rtol=0, equal_nan=False)
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ if os.path.exists(os.path.join(HERE, 'jbasis_d.dat')):
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'),
+ os.path.join(tmp, 'jbasis_d.dat'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'CaOs2.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ ctqmcout.write_ctqmcout(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.ctqmcout.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.ctqmcout'), fh.readlines())
+ finally:
+ shutil.rmtree(tmp)
+
+
+_check('CaOs2_full')
+_check('CaOs2_jbasis_full')
+_check('CaOs2_mode1_full') # proj_mode 1: band-index window 43..92
+_check('CaOs2_mode2_full') # proj_mode 2: explicit band indices 43 92
+
+print('wien2k_ctqmcout_python: ok')
+
+# A singular overlap must fail clearly rather than writing NaN projectors.
+try:
+ ctqmcout._sqrt_inv(np.ones((2, 2), dtype=complex))
+except ValueError as error:
+ assert 'rank deficient' in str(error)
+else:
+ raise AssertionError('singular overlap was accepted')
+
+# Rank is invariant under a uniform projector scaling. Small, well-conditioned
+# overlaps must remain valid and normalize to the identity.
+for scale in (1., 1e-20):
+ overlap = np.diag([scale, 2 * scale]).astype(complex)
+ inverse_root = ctqmcout._sqrt_inv(overlap)
+ np.testing.assert_allclose(inverse_root @ overlap @ inverse_root.conj().T,
+ np.eye(2), atol=1e-14, rtol=0)
diff --git a/test/python/wien2k/wien2k_dmftproj_common_python.py b/test/python/wien2k/wien2k_dmftproj_common_python.py
new file mode 100644
index 0000000..cf1097d
--- /dev/null
+++ b/test/python/wien2k/wien2k_dmftproj_common_python.py
@@ -0,0 +1,186 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Unit tests for the shared dmftproj machinery (triqs_dftkit.wien2k._dmftproj),
+# the module the five case.* generators are built on. The generators' golden
+# tests cover the assembled output; these exercise the isolated primitives
+# directly so a regression in a shared unit is caught at its source.
+
+import os
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import _dmftproj as dp
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+# --- list-directed numeric parsing -------------------------------------------
+
+assert dp.to_float('1.0D-3') == 1.0e-3
+assert dp.to_float('4.57E-002') == 4.57e-2
+assert dp.to_float('-2.5d0') == -2.5
+
+assert dp.to_complex('(1.0,-2.0)') == complex(1.0, -2.0)
+
+
+def _reader(text):
+ fh = tempfile.NamedTemporaryFile('w', suffix='.tmp', delete=False)
+ fh.write(text)
+ fh.close()
+ return dp.Reader(fh.name)
+
+
+r = _reader('(1.0,2.0)\n(3.0,-4.0) (5.0,6.0)\n')
+assert dp.read_complex(r) == complex(1.0, 2.0)
+a, b = dp.read_two_complex(r)
+assert a == complex(3.0, -4.0) and b == complex(5.0, 6.0)
+
+
+# --- angular basis: cubic transform, with and without the float32 cast -------
+
+exact = dp.reptrans('cubic', 2, cast=False)
+cast = dp.reptrans('cubic', 2, cast=True)
+assert exact.shape == (5, 5)
+# unitary
+assert np.max(np.abs(exact @ np.conj(exact.T) - np.eye(5))) < 1e-12
+# the l = 1 (p) table is unitary too: it was un-normalized and mis-ordered,
+# and only the d path was covered before (the non-centrosymmetric p bug).
+p1 = dp.reptrans('cubic', 1, cast=False)
+assert p1.shape == (3, 3)
+assert np.max(np.abs(p1 @ np.conj(p1.T) - np.eye(3))) < 1e-12
+# the cast truncates 1/sqrt2 to single precision (legacy constructor rounding)
+assert abs(abs(exact[1, 0]) - 2 ** -0.5) < 1e-15
+assert 0 < abs(abs(cast[1, 0]) - 2 ** -0.5) < 1e-6
+assert dp.reptrans('complex', 2).shape == (5, 5)
+assert np.allclose(dp.reptrans('complex', 2), np.eye(5))
+
+# The tellurium representative has a nontrivial local frame. Its composed
+# orbital rotation must stay unitary and must include more than the atom-mapping
+# symmetry operation alone.
+_, te_ops, te_refs = dp.read_dmftsym(os.path.join(HERE, 'Te.dmftsym'),
+ rotloc=True)
+te_rot = dp.rotloc_orbital(1, te_refs[0], te_ops, 1, 1)
+np.testing.assert_allclose(te_rot @ te_rot.conj().T, np.eye(3),
+ atol=1e-12, rtol=0)
+mapping = next(op for op in te_ops if op['perm'][0] == 1)
+mapping_only = dp.dmat(1, mapping['a'], mapping['b'], mapping['g'],
+ float(mapping['iprop']))
+assert np.max(np.abs(te_rot - mapping_only)) > 0.1
+
+
+# --- Wigner D and orbital time reversal --------------------------------------
+
+for l in (1, 2):
+ D0 = dp.dmat(l, 0.0, 0.0, 0.0, 1.0)
+ assert np.max(np.abs(D0 - np.eye(2 * l + 1))) < 1e-12 # zero rotation
+ D = dp.dmat(l, 0.3, 0.7, 1.1, 1.0)
+ assert np.max(np.abs(D @ np.conj(D.T) - np.eye(2 * l + 1))) < 1e-12 # unitary
+ T = dp.tmat(l)
+ for m in range(-l, l + 1):
+ assert T[-m + l, m + l] == (-1) ** m
+
+
+# --- select_window: contiguous band range in (e1, e2] ------------------------
+
+eband = np.array([-1.0, -0.1, 0.05, 0.2, 0.5]) # bands 10..14, window (-0.15, 0.30]
+incl, lo, hi = dp.select_window(10, 14, eband, -0.15, 0.30)
+assert incl and lo == 11 and hi == 13
+
+
+# --- band-index window (proj_mode 1/2) ---------------------------------------
+
+# set_projections.f:70-88: every k included, indices clamped to nbmin/nbmax.
+incl, lo, hi = dp.select_band_window(1, 92, 60, 76)
+assert incl and lo == 60 and hi == 76
+incl, lo, hi = dp.select_band_window(50, 80, 43, 92) # clamp both ends
+assert incl and lo == 50 and hi == 80
+
+# proj_mode 2 takes b_bot/b_top straight from the (rounded) window line.
+m2 = {'proj_mode': 2, 'e_bot': 43.0, 'e_top': 92.0}
+assert dp.band_index_window(m2, []) == (43, 92)
+
+# proj_mode 1 scans the (e_bot, e_top] energies for the global band-index range.
+m1 = {'proj_mode': 1, 'e_bot': -0.15, 'e_top': 0.30}
+spins = [{'kp': [{'nbmin': 1, 'nbmax': 5,
+ 'eband': np.array([-1.0, -0.1, 0.05, 0.2, 0.5])}]}]
+assert dp.band_index_window(m1, spins) == (2, 4)
+
+
+# --- case.indmftpr / case.dmftsym structured parse ---------------------------
+
+info = dp.read_indmftpr(os.path.join(HERE, 'CaOs2.indmftpr'))
+assert info['nsort'] == 2 and info['mult'] == [1, 2] and info['lmax'] == 3
+assert info['so'] == 1
+assert info['sorts'][1]['correlated_ls'] == [2] # Os d correlated
+assert info['sorts'][0]['correlated_ls'] == [] # Ca nothing
+
+part = dp.read_indmftpr(os.path.join(HERE, 'CaOs2_partial.indmftpr'))
+assert part['sorts'][0]['included_ls'] == [2] # Ca d now an included shell
+assert part['sorts'][0]['correlated_ls'] == []
+assert part['sorts'][1]['correlated_ls'] == [2]
+
+nsym, ops = dp.read_dmftsym(os.path.join(HERE, 'CaOs2.dmftsym'))
+assert nsym == 16 and len(ops) == 16
+assert all(o['krotm'].shape == (3, 3) and len(o['perm']) == 3 for o in ops)
+
+# Atomic sorts and symmetry operations are independent counts (P1 can have
+# many inequivalent atoms but only the identity operation).
+with tempfile.TemporaryDirectory() as tmp:
+ path = os.path.join(tmp, 'P1.dmftsym')
+ with open(path, 'w') as fh:
+ fh.write('1\n1 2\nSym. op. 1\n0 0 0 1\n'
+ '1 0 0\n0 1 0\n0 0 1\nGlobal->local\n'
+ '1\n1 0 0\n0 1 0\n0 0 1\n0 0 0 1\n'
+ '2\n1 0 0\n0 1 0\n0 0 1\n90 0 0 1\n')
+ _, _, refs = dp.read_dmftsym(path, rotloc=True)
+ assert len(refs) == 2
+ np.testing.assert_allclose(refs[1]['a'], np.pi / 2, atol=1e-14, rtol=0)
+
+ # A tiny but nonzero spin coupling must retain the full spinor basis.
+ theta = 1e-9
+ basis = np.eye(6)
+ basis[0, 0] = basis[3, 3] = np.cos(theta)
+ basis[0, 3], basis[3, 0] = np.sin(theta), -np.sin(theta)
+ path = os.path.join(tmp, 'basis.dat')
+ with open(path, 'w') as fh:
+ for row in basis:
+ fh.write(' ' + ' '.join(f'{value:.17E} 0' for value in row).replace('E', 'D') + '\n')
+ transform, mixing = dp.read_fromfile(path, 1)
+ assert mixing and transform.shape == (6, 6)
+ np.testing.assert_array_equal(transform, basis.conjugate())
+
+ # Each correlated l with an explicit irep owns its own correps record.
+ path = os.path.join(tmp, 'multi.indmftpr')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n1\n-2 3\n')
+ parsed = dp.read_indmftpr(path)
+ assert parsed['so'] == 1 and parsed['sorts'][0]['correlated_ls'] == [1, 2]
+
+ # Selecting individual irreps is not implemented: do not silently export
+ # the entire shell in place of the requested correlated subspace.
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncubic\n0 0 2\n0 0 2\n10\n1\n-2 3\n')
+ try:
+ dp.read_indmftpr(path)
+ except NotImplementedError as error:
+ assert 'irreducible' in str(error)
+ else:
+ raise AssertionError('split correlated shell was silently accepted')
+
+ # The input lacks the cross overlaps needed for two local orbitals;
+ # the native implementation rejects this case as well.
+ path = os.path.join(tmp, 'multi.almblm')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\n0\n1\n2\n')
+ try:
+ dp.read_almblm(path, dict(nsort=1, lmax=0, mult=[1]))
+ except NotImplementedError as error:
+ assert 'local orbital' in str(error)
+ else:
+ raise AssertionError('multiple local orbitals were silently accepted')
+
+print('wien2k_dmftproj_common: ok')
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')
diff --git a/test/python/wien2k/wien2k_noso_python.py b/test/python/wien2k/wien2k_noso_python.py
new file mode 100644
index 0000000..50f17e6
--- /dev/null
+++ b/test/python/wien2k/wien2k_noso_python.py
@@ -0,0 +1,127 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.* generators on the non-spin-orbit path (ifSO=0),
+# the common DMFT case (e.g. SrVO3). dmftproj writes (2l+1) matrices there, not
+# the 2*(2l+1) spinor matrices of the SO path: the orthonormalization is
+# per-spin (orthogonal_wannier), the symmetry and Rloc matrices are the bare
+# (2l+1) representation, the density matrix splits into two independent spin
+# blocks, and srot%timeinv is always false (setsym.f sets it only under SP+SO).
+#
+# CaOs2_noso is the spin-polarized, non-SO CaOs2 (indmftpr SO flag 0, wide
+# window -2.0 3.0): two equivalent Os d shells, cubic harmonics. The wide window
+# keeps the Loewdin overlap full rank. ctqmcout and symqmc are checked against
+# the precision-fixed dmftproj reference for this case.
+#
+# CaOs2_noso_partial adds an uncorrelated Ca d shell, so sympar and parproj have
+# an included-but-not-correlated shell to symmetrize. Same struct/dmftsym/almblm
+# as the correlated-only case. References built with the precision-fixed
+# dmftproj run without -so.
+
+import gzip
+import os
+import shutil
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout, symqmc, sympar, parproj
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_lines):
+ got = open(got_path).read().split()
+ ref = ''.join(ref_lines).split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ floats = []
+ for a, b in zip(got, ref):
+ if b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ else:
+ floats.append((float(a), float(b)))
+ actual, expected = np.array(floats).T
+ np.testing.assert_allclose(actual, expected, atol=1e-11, rtol=0,
+ equal_nan=False)
+
+
+def _stage(case):
+ tmp = tempfile.mkdtemp()
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'CaOs2_noso.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ return tmp
+
+
+def _check_output(case, ext, writer):
+ tmp = _stage(case)
+ try:
+ writer(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.{ext}.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.{ext}'), fh.readlines())
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Correlated-only: the projector and correlated-shell symmetry files.
+_check_output('CaOs2_noso', 'ctqmcout', ctqmcout.write_ctqmcout)
+_check_output('CaOs2_noso', 'symqmc', symqmc.write_symqmc)
+
+# Partial: the included-shell symmetry and partial-projector files.
+_check_output('CaOs2_noso_partial', 'sympar', sympar.write_sympar)
+_check_output('CaOs2_noso_partial', 'parproj', parproj.write_parproj)
+
+# Shells with different l on one atom remain independent in the non-SO path.
+orbs = [dict(atom=1, sort=1, l=l) for l in (0, 1)]
+raw = {(1, l): [np.eye(2 * l + 1) * value for value in (1., 2.)]
+ for l in (0, 1)}
+representations = {(l, 1, 0): (np.eye(2 * l + 1), False, 0.) for l in (0, 1)}
+result = parproj._symmetrize_densmat_nonso(
+ orbs, {}, [dict(perm=[1])], 1, representations, raw)
+for shell in raw:
+ for actual, expected in zip(result[shell], raw[shell]):
+ np.testing.assert_allclose(actual, expected, atol=1e-14, rtol=0)
+
+# SO metadata requires both spin channels, including when a single-channel
+# almblm happens to be present. Never export a half-sized SO projector.
+tmp = _stage('CaOs2_noso')
+try:
+ case = os.path.join(tmp, 'CaOs2_noso')
+ os.rename(case + '.almblmup', case + '.almblm')
+ os.unlink(case + '.almblmdn')
+ with open(case + '.indmftpr') as fh:
+ lines = fh.read().splitlines()
+ lines[-2] = '1'
+ with open(case + '.indmftpr', 'w') as fh:
+ fh.write('\n'.join(lines) + '\n')
+ for writer in (ctqmcout.write_ctqmcout, parproj.write_parproj):
+ try:
+ writer(case)
+ except ValueError as error:
+ assert 'both up and down' in str(error)
+ else:
+ raise AssertionError('single-channel SO projectors were accepted')
+finally:
+ shutil.rmtree(tmp)
+
+for module in (symqmc, sympar):
+ try:
+ module._shell_matrix({}, dict(l=1, mixing=True, P=np.eye(6)), 0, False)
+ except NotImplementedError as error:
+ assert 'spin-mixing' in str(error)
+ else:
+ raise AssertionError('non-SO symmetry accepted a spinor basis')
+
+print('wien2k_noso_python: ok')
diff --git a/test/python/wien2k/wien2k_oddl_python.py b/test/python/wien2k/wien2k_oddl_python.py
new file mode 100644
index 0000000..74e65e1
--- /dev/null
+++ b/test/python/wien2k/wien2k_oddl_python.py
@@ -0,0 +1,83 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the generators on the odd-l shells (p, l=1; f, l=3). The d-only fixtures
+# (l=2) cannot catch the improper-rotation parity: the Wigner-D parity factor is
+# (-1)^l, which is +1 for even l and only flips sign for odd l. case.dmftsym
+# stores the proper part of each operation in krotm (det(krotm)=+1 even for
+# improper ops) and the true parity in iprop; dmat must use iprop. The p- and
+# f-shell symqmc here pin that sign.
+#
+# Both are the spin-orbit CaOs2 cell (two equivalent Os atoms) with the Os shell
+# promoted to a complex p or f shell. The Loewdin overlap must be full rank:
+# p (l=1): the standard almblm, wide window -2..3 (6 spin-orbitals/atom).
+# f (l=3): Os carries no f weight near E_F, so the projector is rank-deficient
+# on the standard bands. CaOs2_fshell.almblm is recomputed with a high band
+# cutoff (lapw1/lapwso EMAX raised on the converged density); the window
+# 3..10 then sits on high-energy bands with real l=3 character, full rank
+# (overlap condition number ~1.3).
+
+import gzip
+import os
+import shutil
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout, symqmc
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_gz):
+ got = open(got_path).read().split()
+ with gzip.open(ref_gz, 'rt') as fh:
+ ref = fh.read().split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ floats = []
+ for a, b in zip(got, ref):
+ if b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ else:
+ floats.append((float(a), float(b)))
+ actual, expected = np.array(floats).T
+ np.testing.assert_allclose(actual, expected, atol=1e-11, rtol=0,
+ equal_nan=False)
+
+
+def _check(case, almblm_case, exts):
+ tmp = tempfile.mkdtemp()
+ try:
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for spin in ('up', 'dn'):
+ _gunzip(os.path.join(HERE, f'{almblm_case}.almblm{spin}.gz'),
+ os.path.join(tmp, f'{case}.almblm{spin}'))
+ case_path = os.path.join(tmp, case)
+ for ext, writer in exts:
+ writer(case_path)
+ _compare(case_path + '.' + ext,
+ os.path.join(HERE, f'{case}.{ext}.gz'))
+ finally:
+ shutil.rmtree(tmp)
+
+
+# p shell: standard SO almblm, full rank on the wide window.
+_check('CaOs2_pshell', 'CaOs2',
+ [('ctqmcout', ctqmcout.write_ctqmcout), ('symqmc', symqmc.write_symqmc)])
+
+# f shell: high-band-cutoff almblm, full rank on the high-energy window.
+_check('CaOs2_fshell', 'CaOs2_fshell',
+ [('ctqmcout', ctqmcout.write_ctqmcout), ('symqmc', symqmc.write_symqmc)])
+
+print('wien2k_oddl_python: ok')
diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py
new file mode 100644
index 0000000..048f070
--- /dev/null
+++ b/test/python/wien2k/wien2k_oubwin_python.py
@@ -0,0 +1,84 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.oubwin generator (triqs_dftkit.wien2k.oubwin)
+# reproduces the dmftproj Fortran output: regenerate case.oubwinup / case.oubwindn
+# from case.almblm{up,dn} + case.indmftpr and compare to the native reference,
+# ignoring only trailing whitespace removed from the shipped fixture.
+#
+# CaOs2 is spin-orbit + spin-polarized, so both outputs match CaOs2.oubwinup.
+# Reuse the compressed almblm inputs shipped with the SOC fixture.
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import oubwin
+from triqs_dftkit.wien2k import _dmftproj as dp
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ for ext in ('almblmup', 'almblmdn'):
+ with gzip.open(os.path.join(HERE, f'{case}.{ext}.gz'), 'rb') as src, \
+ open(os.path.join(tmp, f'{case}.{ext}'), 'wb') as dst:
+ shutil.copyfileobj(src, dst)
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'), tmp)
+ oubwin.write_oubwin(os.path.join(tmp, case))
+ with open(os.path.join(HERE, f'{case}.oubwinup')) as f:
+ ref = [line.rstrip() for line in f]
+ for ext in ('oubwinup', 'oubwindn'):
+ with open(os.path.join(tmp, f'{case}.{ext}')) as f:
+ got = [line.rstrip() for line in f]
+ assert got == ref, f'{case}.{ext} differs from dmftproj reference'
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Independent expected band ranges for set_projections.f's branch ordering.
+for nbmin, energies, expected in (
+ (10, [-3., -2., -1., 0.5], (False, 0, 0)), # opens only at nbmax
+ (5, [0.5], (False, 0, 0)), # nbmin == nbmax
+ (10, [0., 0.5, 1., 2.], (True, 11, 12)), # strict lower, inclusive upper
+ (10, [0.5, 1.], (True, 10, 11)), # open window reaches nbmax
+ (10, [-1., 0., 2.], (False, 0, 0))): # no in-window band
+ got = dp.select_window(nbmin, nbmin + len(energies) - 1,
+ np.array(energies), 0., 1.)
+ assert got == expected, (energies, got, expected)
+
+# Three k-points with distinct weights, a Fermi shift of 10 Ry, and one
+# excluded point. The lowest-band weight is used even if that band is outside
+# the window; an excluded point emits no band range or weight record.
+with tempfile.TemporaryDirectory() as tmp:
+ almblm = os.path.join(tmp, 'case.almblm')
+ with open(almblm, 'w') as f:
+ f.write('0\n3\n0\n10.0\n1.0\n0\n')
+ for ik, nbmin, energies, weights in (
+ (1, 10, [9., 10.25, 10.75, 12.], [.25, .80, .60, .25]),
+ (2, 20, [8., 9.], [.10, .10]),
+ (3, 30, [10.25, 10.75], [.75, .50])):
+ f.write(f'IK = {ik}\n0 0 0\n{ik} {nbmin} {nbmin + len(energies) - 1}\n')
+ for energy, weight in zip(energies, weights):
+ f.write(f'{weight} {energy}\n')
+ f.write('ATOM = 1\n1\n')
+ f.write('(0.,0.) (0.,0.)\n' * len(energies))
+ windows = oubwin._windows_from_spin(dp.read_almblm(almblm, dict(nsort=1, lmax=0, mult=[1])), dict(proj_mode=0, e_bot=0., e_top=1.), None)
+ assert windows == [(True, 11, 12, 0.25), (False, 0, 0, 0.10),
+ (True, 30, 31, 0.75)], windows
+ output = os.path.join(tmp, 'case.oubwin')
+ oubwin._write_oubwin_file(output, 0, windows)
+ with open(output) as f:
+ assert [line.split() for line in f] == [
+ ['3'], ['0'], ['1'], ['11', '12'], ['0.2500000000000000'],
+ ['0'], ['1'], ['30', '31'], ['0.7500000000000000']]
+
+_check('CaOs2')
+
+print('wien2k_oubwin_python: ok')
diff --git a/test/python/wien2k/wien2k_outband_python.py b/test/python/wien2k/wien2k_outband_python.py
new file mode 100644
index 0000000..a05fa66
--- /dev/null
+++ b/test/python/wien2k/wien2k_outband_python.py
@@ -0,0 +1,98 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.outband generator (triqs_dftkit.wien2k.outband)
+# reproduces the dmftproj Fortran band-structure projectors to machine
+# precision: regenerate case.outband from the band-mode almblm/struct/dmftsym,
+# the k-path (klist_band) and the band indmftpr, and compare to the
+# precision-fixed dmftproj -band reference.
+#
+# case.outband feeds convert_bands_input. It is the band analog of
+# case.ctqmcout: the same correlated-shell projector construction (raw Alm/Clm
+# projector -> rot_projectmat -> fromfile/cubic transmat -> Loewdin
+# orthonormalization) plus the raw Theta projectors of the parproj path, written
+# in the outband layout. The CaOs2_band fixture is SP+SO with the spin-mixing
+# fromfile (|j, m_j>) Os d shell (jbasis_d.dat); the wide window (-2.0, 3.0)
+# keeps 50 bands > 20 correlated spin-orbitals, so the Loewdin overlap is full
+# rank and the projectors match to machine precision.
+#
+# Band-mode specifics (dmftproj.f:557-581): nkband, the number of k-points along
+# the plotted path, is read from CaOs2_band.klist_band; the Fermi energy is the
+# last line of CaOs2_band.indmftpr (the band almblm keeps a placeholder eferm
+# record that the reader skips).
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import outband
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_lines):
+ got = open(got_path).read().split()
+ ref = ''.join(ref_lines).split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ floats = []
+ for a, b in zip(got, ref):
+ if a.lstrip('-').isdigit() and b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ elif _isfloat(a) and _isfloat(b):
+ floats.append((float(a), float(b)))
+ else:
+ assert a == b, (a, b) # k-label tokens
+ a, b = np.array(floats).T
+ np.testing.assert_allclose(a, b, atol=1e-11, rtol=0, equal_nan=False)
+
+
+def _isfloat(tok):
+ try:
+ float(tok)
+ return True
+ except ValueError:
+ return False
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ for ext in ('indmftpr', 'klist_band'):
+ shutil.copy(os.path.join(HERE, f'{case}.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'),
+ os.path.join(tmp, 'jbasis_d.dat'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'{case}.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ outband.write_outband(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.outband.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.outband'), fh.readlines())
+ with open(os.path.join(tmp, f'{case}.klist_band'), 'w') as fh:
+ fh.write('GAMMA 0 0 0 8 2.0\nX 4 0 0 8 2.0\nEND\n')
+ assert outband.read_k_list(os.path.join(tmp, f'{case}.klist_band')) == (2, [(1, 'GAMMA'), (2, 'X')])
+ try:
+ outband.write_outband(os.path.join(tmp, case))
+ except ValueError as exc:
+ assert 'k-point counts differ' in str(exc)
+ else:
+ raise AssertionError('inconsistent band-path dimensions must be rejected')
+ finally:
+ shutil.rmtree(tmp)
+
+
+_check('CaOs2_band')
+
+print('wien2k_outband_python: ok')
diff --git a/test/python/wien2k/wien2k_parproj_python.py b/test/python/wien2k/wien2k_parproj_python.py
new file mode 100644
index 0000000..a759fae
--- /dev/null
+++ b/test/python/wien2k/wien2k_parproj_python.py
@@ -0,0 +1,96 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.parproj generator (triqs_dftkit.wien2k.parproj)
+# reproduces the dmftproj Fortran output. parproj holds the partial projectors
+# and density matrices for the INCLUDED orbitals; the CaOs2_partial fixture adds
+# an uncorrelated Ca d-shell. Inputs are the shared CaOs2 almblm/struct/dmftsym
+# renamed to the CaOs2_partial case.
+#
+# Floats compared at 1e-7 (dmftproj single-precision basis-transform floor).
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import parproj
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_lines):
+ got = open(got_path).read().split()
+ ref = ''.join(ref_lines).split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ floats = []
+ for a, b in zip(got, ref):
+ if b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ else:
+ floats.append((float(a), float(b)))
+ a, b = np.array(floats).T
+ np.testing.assert_allclose(a, b, atol=1e-11, rtol=0, equal_nan=False)
+
+
+def _check(case, mode=0):
+ tmp = tempfile.mkdtemp()
+ try:
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ if mode:
+ path = os.path.join(tmp, f'{case}.indmftpr')
+ with open(path) as fh:
+ lines = fh.read().splitlines()
+ lines[-1] = lines[-1] + ' 1' if mode == 1 else '60 76 2'
+ with open(path, 'w') as fh:
+ fh.write('\n'.join(lines) + '\n')
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'), os.path.join(tmp, 'jbasis_d.dat'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'CaOs2.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ parproj.write_parproj(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.parproj.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.parproj'), fh.readlines())
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Identity symmetry must preserve distinct s and p density blocks on one atom.
+orbs = [dict(atom=1, sort=1, l=l) for l in (0, 1)]
+raw = {(1, l): [np.eye(2 * l + 1, dtype=complex) * value
+ for value in (1., 2., 0.25, 0.25)] for l in (0, 1)}
+representations = {(l, 1, 0): (np.eye(2 * l + 1), False, 0.) for l in (0, 1)}
+result = parproj._symmetrize_densmat(orbs, {}, [dict(perm=[1])], 1,
+ representations, raw, 2, True, True)
+for shell in raw:
+ for actual, expected in zip(result[shell], raw[shell]):
+ np.testing.assert_allclose(actual, expected, atol=1e-14, rtol=0)
+
+_check('CaOs2_partial')
+_check('CaOs2_jbasis_full')
+# All definitions select bands 60..76 here; the partial projectors and
+# below-window density must match the native energy-window reference.
+_check('CaOs2_partial', mode=1)
+_check('CaOs2_partial', mode=2)
+
+print('wien2k_parproj_python: ok')
+
+# Equivalent atoms must share the average after identity and swap operations.
+raw_mixing = {1: np.diag([1., 3.]), 2: np.diag([5., 7.])}
+actual = parproj._symmetrize_densmat_mixing(
+ [dict(atom=1), dict(atom=2)], [dict(perm=[1, 2]), dict(perm=[2, 1])], 2,
+ [(np.eye(2), False)] * 2, raw_mixing)
+for matrix in actual.values():
+ np.testing.assert_allclose(matrix, np.diag([3., 5.]), atol=1e-14, rtol=0)
diff --git a/test/python/wien2k/wien2k_sympar_python.py b/test/python/wien2k/wien2k_sympar_python.py
new file mode 100644
index 0000000..a2e01b8
--- /dev/null
+++ b/test/python/wien2k/wien2k_sympar_python.py
@@ -0,0 +1,81 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.sympar generator (triqs_dftkit.wien2k.sympar)
+# reproduces the dmftproj Fortran output. sympar holds the symmetry matrices for
+# all INCLUDED orbitals (correlated + partial), so the CaOs2_partial fixture adds
+# an uncorrelated Ca d-shell to the correlated Os d-shells. Inputs are the shared
+# CaOs2 almblm/struct/dmftsym renamed to the CaOs2_partial case.
+#
+# Floats compared at 1e-7: dmftproj stores the cubic basis transform through a
+# single-precision CMPLX cast, so its matrices carry ~1e-7 noise the double-
+# precision generator reproduces to the same floor.
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import sympar
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _gunzip(src, dst):
+ with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+ shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_lines):
+ got = open(got_path).read().split()
+ ref = ''.join(ref_lines).split()
+ assert len(got) == len(ref), (len(got), len(ref))
+ # 16 operations, three atoms, time flags, then three 10x10 shell matrices.
+ header = 2 + 16 * 3 + 16
+ np.testing.assert_array_equal(np.array(got[:header], int),
+ np.array(ref[:header], int))
+ data = np.array(got[header:], float).reshape(16, 3, 2, 10, 10)
+ reference = np.array(ref[header:], float).reshape(data.shape)
+ matrices = data[:, :, 0] + 1j * data[:, :, 1]
+ expected = reference[:, :, 0] + 1j * reference[:, :, 1]
+ np.testing.assert_allclose(matrices, expected, atol=1e-11, rtol=0,
+ equal_nan=False)
+ for operation in matrices:
+ for matrix in operation:
+ np.testing.assert_allclose(matrix @ matrix.conj().T, np.eye(10),
+ atol=1e-11, rtol=0)
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ for ext in ('struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'CaOs2.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ sympar.write_sympar(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.sympar.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.sympar'), fh.readlines())
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Inversion has proper part I but parity -1: odd-l matrices must be -I.
+inversion = dict(a=0., b=0., c=0., iprop=-1, krotm=np.eye(3))
+for l in (1, 3):
+ shell = dict(l=l, mixing=False, P=np.eye(2 * l + 1, dtype=complex))
+ np.testing.assert_allclose(sympar._shell_matrix(inversion, shell, False, True),
+ -np.eye(2 * (2 * l + 1)), atol=1e-14, rtol=0)
+p = sympar.reptrans('cubic', 1)
+np.testing.assert_allclose(p @ p.conj().T, np.eye(3), atol=1e-11, rtol=0)
+
+_check('CaOs2_partial')
+
+print('wien2k_sympar_python: ok')
diff --git a/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy
index b7e2e01..4cce82c 100644
Binary files a/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy and b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_l0.ref.npy b/test/python/wien2k/wien2k_symqmc_l0.ref.npy
index 679d0e0..ae517ee 100644
Binary files a/test/python/wien2k/wien2k_symqmc_l0.ref.npy and b/test/python/wien2k/wien2k_symqmc_l0.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_mixing.py b/test/python/wien2k/wien2k_symqmc_mixing.py
index eb729a8..dc10a30 100644
--- a/test/python/wien2k/wien2k_symqmc_mixing.py
+++ b/test/python/wien2k/wien2k_symqmc_mixing.py
@@ -65,10 +65,10 @@ def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
data, shells, ops, so = _run(indmftpr, extra)
ref = np.load(os.path.join(HERE, ref_npy))
assert data.shape == ref.shape, (data.shape, ref.shape)
- assert np.max(np.abs(data - ref)) < 1e-6, np.max(np.abs(data - ref))
+ assert np.max(np.abs(data - ref)) < 1e-11, np.max(np.abs(data - ref))
for op in ops:
for sh in shells:
- mat = symqmc._shell_matrix(op, sh, _timeinv(op, so))
+ mat = symqmc._shell_matrix(op, sh, _timeinv(op, so), bool(so))
dev = np.max(np.abs(mat @ np.conj(mat.T) - np.eye(mat.shape[0])))
assert dev < unit_tol, ('not unitary', dev)
diff --git a/test/python/wien2k/wien2k_symqmc_noncentrosym.py b/test/python/wien2k/wien2k_symqmc_noncentrosym.py
new file mode 100644
index 0000000..c80bb0c
--- /dev/null
+++ b/test/python/wien2k/wien2k_symqmc_noncentrosym.py
@@ -0,0 +1,78 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.symqmc generator on a non-centrosymmetric
+# spin-orbit cell (tellurium, P3_1 21, point group 32, a correlated p shell),
+# the regime dft_tools #148 is about and the centrosymmetric d-shell fixtures
+# (CaOs2, SrVO3) never exercised. Two convention bugs hid there:
+#
+# - the cubic l = 1 (p) transform must be unitary and px,py,pz ordered;
+# - the magnetic (time-inversion) operations need the antiunitary basis
+# change P D P^T, not P D P^dag, which differs only for a complex transform
+# (the p_y row), so the real cubic-d path never showed it.
+#
+# Reference: the dmftproj Fortran case.symqmc for the same symmetry input. The
+# Fortran reads the cubic template at single precision, so its floor is ~1e-7;
+# the generator is double precision. We compare the matrix data at 1e-6 and
+# separately assert every spinor matrix is unitary.
+
+import os
+import shutil
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import symqmc
+from triqs_dftkit.wien2k._dmftproj import reptrans
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+# The cubic l = 1 transform itself must be unitary (P P^dag = I). The d (l = 2)
+# table was correct; only this one carried the un-normalized, mis-ordered rows.
+P1 = reptrans('cubic', 1)
+assert np.max(np.abs(P1 @ np.conj(P1.T) - np.eye(3))) < 1e-12, 'cubic l=1 not unitary'
+
+
+def _matrix_data(path, nsym, natom):
+ """Flat array of the matrix blocks, skipping header/perm/timeflag lines."""
+ toks = iter(open(path).read().split())
+ assert int(next(toks)) == nsym and int(next(toks)) == natom
+ for _ in range(nsym * natom + nsym): # perm rows + the time-reversal flags
+ next(toks)
+ return np.array([float(x) for x in toks])
+
+
+tmp = tempfile.mkdtemp()
+try:
+ for src, dst in (('Te.dmftsym', 'case.dmftsym'),
+ ('Te.struct', 'case.struct'),
+ ('Te.indmftpr', 'case.indmftpr')):
+ shutil.copy(os.path.join(HERE, src), os.path.join(tmp, dst))
+ case = os.path.join(tmp, 'case')
+ symqmc.write_symqmc(case)
+ nsym, ops = symqmc._read_dmftsym(case + '.dmftsym')
+ shells, so = symqmc._read_correlated_shells(case + '.indmftpr', case + '.struct')
+ natom = len(ops[0]['perm'])
+ data = _matrix_data(case + '.symqmc', nsym, natom)
+finally:
+ shutil.rmtree(tmp)
+
+ref = np.load(os.path.join(HERE, 'wien2k_symqmc_te.ref.npy'))
+assert data.shape == ref.shape, (data.shape, ref.shape)
+np.testing.assert_allclose(data, ref, atol=1e-7, rtol=0, equal_nan=False)
+
+
+def _timeinv(op, so):
+ k = op['krotm']
+ return 1 if (so and k[0, 0] * k[1, 1] - k[0, 1] * k[1, 0] < 0.0) else 0
+
+
+for op in ops:
+ for sh in shells:
+ mat = symqmc._shell_matrix(op, sh, _timeinv(op, so), bool(so))
+ dev = np.max(np.abs(mat @ np.conj(mat.T) - np.eye(mat.shape[0])))
+ assert dev < 1e-7, ('not unitary', dev)
+
+print('wien2k_symqmc_noncentrosym: ok')
diff --git a/test/python/wien2k/wien2k_symqmc_te.ref.npy b/test/python/wien2k/wien2k_symqmc_te.ref.npy
new file mode 100644
index 0000000..4d006d4
Binary files /dev/null and b/test/python/wien2k/wien2k_symqmc_te.ref.npy differ