scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Simple output file for noncollinear case was not correct

This commit is contained in:
giannozz 2021-12-16 19:58:35 +00:00
parent bd6f745af2
commit 6551b7ad43
1 changed files with 203 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

322
PP/src/pw2gt.f90
View File

@ -9,33 +9,41 @@
PROGRAM pw2gt
!-----------------------------------------------------------------------
! Convert QE data files to the file format used by the GreenT code.
! The output file is a single, formatted and self-contained file
! (BEWARE! may be quite large!) that can be useful for anybody else
! not willing to figure out how read QE data files.
! IMPORTANT NOTICE: must be run on a single processor ONLY.
! Uncomment next line to test the file: the code will re-read the file
! just written, fill and diagonalize the hamiltonian matrix
!#define __DEBUG
! The output is a single, simple, formatted and self-contained file.
! Can be useful for anybody wanting to access and process QE data
! with little hassle.
!
! IMPORTANT NOTICE 1: must be run on a single processor ONLY
! IMPORTANT NOTICE 2: the file may easily become VERY large!
!
! Input: namelist &inputpp ... /
! outdir as in QE input (default values as in QE).
! prefix " " " " "
! check .false. (default) / .true.
! if .true. do not write the file but read it,
! fill and diagonalize the hamiltonian matrix
! BEWARE! do not attempt to do that for systems
! exceeding a few thousands plane waves
!
! Input: namelist &inputpp [outdir=...] [prefix=...] / as in QE input
! (default values as in QE).
! Output file written in "outdir"/"prefix".save/output.dat
!
! Written by Paolo Giannozzi, with help from Marco Pala
!
USE io_global, ONLY : ionode, ionode_id
USE io_files, ONLY : tmp_dir, prefix
USE mp_global, ONLY : mp_startup
USE io_files, ONLY : tmp_dir, prefix, restart_dir
USE mp_global, ONLY : mp_startup, mp_global_end
USE mp_images, ONLY : intra_image_comm, nproc_image
USE mp, ONLY : mp_bcast
USE environment,ONLY : environment_start, environment_end
!
IMPLICIT NONE
!
LOGICAL :: needwf = .true.
LOGICAL :: needwf = .true., check
INTEGER :: ios
CHARACTER(LEN=256) :: outdir
CHARACTER(LEN=256) :: outdir, fileout
CHARACTER(LEN=256), EXTERNAL :: trimcheck
!
NAMELIST / inputpp / outdir, prefix
NAMELIST / inputpp / outdir, prefix, check
!
! initialise environment
!
@ -49,6 +57,7 @@ PROGRAM pw2gt
prefix = 'pwscf'
CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir )
IF ( trim( outdir ) == ' ' ) outdir = './'
check = .false.
!
IF ( ionode ) THEN
!
@ -71,23 +80,27 @@ PROGRAM pw2gt
IF ( nproc_image > 1 ) CALL errore ('pw2gt', &
'must be run on a single processor only', nproc_image )
!
! Read xml file, allocate and initialize general variables
!
CALL read_file_new ( needwf )
!
CALL simple_output ()
!
#if defined (__DEBUG)
CALL simple_diag ()
#endif
fileout = TRIM ( restart_dir() ) // 'output.dat'
IF ( check ) THEN
!
CALL simple_diag ( fileout )
!
ELSE
!
! Read xml file, allocate and initialize general variables
!
CALL read_file_new ( needwf )
CALL simple_output ( fileout )
!
END IF
!
CALL environment_end ( 'PW2GT' )
!
CALL stop_pp()
CALL mp_global_end ()
!
END PROGRAM pw2gt
!----------------------------------------------------------------------------
SUBROUTINE simple_output ( )
SUBROUTINE simple_output ( fileout )
!----------------------------------------------------------------------------
!
! Not-so-smart but easy-to-read output file for simple cases (e.g. Si)
@ -101,26 +114,25 @@ SUBROUTINE simple_output ( )
USE gvecs, ONLY : doublegrid
USE scf, ONLY : vrs, vltot, v, kedtau
USE fft_rho, ONLY : rho_r2g
USE lsda_mod, ONLY : nspin
USE lsda_mod, ONLY : nspin, isk
USE klist, ONLY : nks, xk, ngk, igk_k
USE uspp, ONLY : okvan, nkb, vkb, dvan, dvan_so
USE uspp_param, ONLY : nh
USE uspp_init, ONLY : init_us_2
USE fft_base, ONLY : dfftp, dffts
USE wvfct, ONLY : nbnd, et
USE wvfct, ONLY : nbnd, et, npwx
USE wavefunctions, ONLY : evc
USE noncollin_module, ONLY : lspinorb, domag
USE pw_restart_new, ONLY : read_collected_wfc
!
IMPLICIT NONE
!
CHARACTER(LEN=*), intent(in) :: fileout
COMPLEX(dp), ALLOCATABLE :: vaux(:,:)
CHARACTER(LEN=256) :: fileout
INTEGER :: iun, ig, is, ik, ikb, ibnd, na, nt, npw
!
fileout = TRIM ( restart_dir() ) // 'output.dat'
WRITE( UNIT = stdout, FMT = '(/,5X,"Writing simple output data file ",A)' ) &
fileout
trim(fileout)
OPEN ( NEWUNIT = iun, FORM = 'formatted', STATUS = 'unknown', &
FILE = fileout )
WRITE(iun,'("# Primitive lattice vectors a_1, a_2, a_3 (a.u.)")')
@ -164,13 +176,13 @@ SUBROUTINE simple_output ( )
DO nt = 1, nsp
WRITE(iun,'("# Atom type, number of projectors")')
WRITE(iun,*) nt, nh(nt)
IF ( nspin /= 4 ) THEN
WRITE(iun,*) dvan(1:nh(nt),1:nh(nt),nt)
ELSE
IF ( lspinorb ) THEN
DO is = 1, nspin
WRITE(iun,'("# spin component n.",i4)') is
WRITE(iun,*) dvan_so(1:nh(nt),1:nh(nt),is,nt)
END DO
ELSE
WRITE(iun,*) dvan(1:nh(nt),1:nh(nt),nt)
END IF
END DO
WRITE(iun,'("# number of beta functions")')
@ -181,7 +193,11 @@ SUBROUTINE simple_output ( )
WRITE(iun,'("# number of plane waves for all k-points")')
WRITE(iun,*) ngk(1:nks)
DO ik=1,nks
WRITE(iun,'("# k-point n.",i4)') ik
IF( nspin == 2 ) THEN
WRITE(iun,'("# k-point n.",i4," spin ",i1)') ik,isk(ik)
ELSE
WRITE(iun,'("# k-point n.",i4)') ik
ENDIF
WRITE(iun,*) tpiba*xk(:,ik)
WRITE(iun,'("# number of plane waves")')
npw = ngk(ik)
@ -200,6 +216,8 @@ SUBROUTINE simple_output ( )
WRITE(iun,'("# band n.",i4)') ibnd
WRITE(iun,'("# eigenvalue (eV):",e25.15)') et(ibnd,ik)*13.6058
WRITE(iun,'(2e25.15)') evc(1:npw,ibnd)
IF ( nspin == 4 ) &
WRITE(iun,'(2e25.15)') evc(1+npwx:npw+npwx,ibnd)
END DO
END DO
WRITE(iun,'("# end of file")')
@ -209,68 +227,66 @@ SUBROUTINE simple_output ( )
END SUBROUTINE simple_output
!
!----------------------------------------------------------------------------
SUBROUTINE simple_diag ( )
SUBROUTINE simple_diag ( fileout )
!----------------------------------------------------------------------------
!
! Check: read the simple data file, re-diagonalize the matrix
!
USE kinds, ONLY : dp
USE io_global, ONLY : stdout
USE io_files, ONLY : restart_dir
USE cell_base, ONLY : at, bg, alat, tpiba
USE ions_base, ONLy : ityp, atm, tau
USE gvect, ONLY : mill
USE klist, ONLY : xk, igk_k
USE uspp, ONLY : vkb, dvan, dvan_so
USE uspp_param, ONLY : nh
USE wavefunctions,ONLY : evc
!
IMPLICIT NONE
!
CHARACTER(LEN=*), intent(in) :: fileout
!
! These variables are read from file
INTEGER :: npw, nbnd, nspin, ngm, nat, nkb, nsp, npol, nks
INTEGER :: npw, nbnd, nspin, current_spin, ngm, nat, nkb, nsp, npol, nks
INTEGER, ALLOCATABLE :: ngk(:), ityp(:), nh(:), mill(:,:), igk(:)
LOGICAL :: lspinorb, domag, okvan
REAL(dp), ALLOCATABLE :: et(:)
COMPLEX(dp), ALLOCATABLE :: vaux(:,:)
COMPLEX(dp) :: gfortran_merda
REAL(dp) :: at(3,3), bg(3,3), xk(3)
REAL(dp), ALLOCATABLE :: et(:), dvan(:,:,:)
COMPLEX(dp), ALLOCATABLE :: vaux(:,:), evc(:,:), vkb(:,:), dvan_so(:,:,:,:)
!
CHARACTER(LEN=256) :: fileout
CHARACTER(LEN=80) :: line
COMPLEX(dp), ALLOCATABLE :: h(:,:), v(:,:)
INTEGER :: iun, ig, is, ik, ikb, ibnd, na, nt, nt_, i, j, ii, jj, ij, &
ih, jh, i1, i2, i3, ipol
nhm, ih, jh, i1, i2, i3, ipol
INTEGER, ALLOCATABLE :: limm(:,:,:)
LOGICAL :: debug = .false., skip_diag = .false.
REAL(dp) :: g(3)
LOGICAL :: debug = .true.
REAL(dp), ALLOCATABLE :: e(:)
COMPLEX(dp), ALLOCATABLE :: h(:,:), v(:,:)
!
fileout = TRIM ( restart_dir() ) // 'output.dat'
WRITE( UNIT = stdout, FMT = '(/,5X,"Checking simple output data file ",A)' ) &
fileout
trim(fileout)
OPEN ( NEWUNIT = iun, FORM = 'formatted', STATUS = 'old', FILE = fileout )
READ(iun,'(a)') line
READ(iun,*) at(:,1), at(:,2), at(:,3)
IF ( debug) WRITE(stdout,*) trim(line), at(:,1), at(:,2), at(:,3)
WRITE(stdout,*) trim(line), at(:,1), at(:,2), at(:,3)
READ(iun,'(a)') line
READ(iun,*) bg(:,1), bg(:,2), bg(:,3)
IF ( debug) WRITE(stdout,*) trim(line), bg(:,1), bg(:,2), bg(:,3)
WRITE(stdout,*) trim(line), bg(:,1), bg(:,2), bg(:,3)
READ(iun,'(a)') line
READ(iun,*) nsp
IF ( debug) WRITE(stdout,*) trim(line), nsp
WRITE(stdout,*) trim(line), nsp
READ(iun,'(a)') line
READ(iun,*) nat
IF ( debug) WRITE(stdout,*) trim(line), nat
WRITE(stdout,*) trim(line), nat
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
WRITE(stdout,'(a)') line
ALLOCATE (ityp(nat))
DO na =1, nat
READ(iun,*) nt
ityp(na) = nt
READ(iun,'(a)') line
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
! READ(iun,'(a,3e25.15)') atm(nt), tau(:,na)
IF (debug) WRITE(stdout,'(a)') line
END DO
READ(iun,'(a)') line
READ(iun,*) ngm
IF ( debug) WRITE(stdout,*) trim(line), ngm
WRITE(stdout,*) trim(line), ngm
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
WRITE(stdout,'(a)') line
ALLOCATE ( mill(3,ngm) )
READ(iun,'(3i8)') (mill(:,ig), ig=1,ngm)
! if i1=mill(1,ig), i2=mill(2,ig), i3=mill(3,ig), then:
! limm(i1,i2,i3) = ig
@ -285,95 +301,119 @@ SUBROUTINE simple_diag ( )
!
READ(iun,'(a)') line
READ(iun,*) nspin
IF ( debug) WRITE(stdout,*) trim(line), nspin
WRITE(stdout,*) trim(line), nspin
IF ( nspin <= 2 ) THEN
npol = 1
domag=.FALSE.
lspinorb=.FALSE.
ELSE
READ(iun,'(a)') line
READ(iun,*) domag, lspinorb
IF ( domag .OR. .NOT.lspinorb ) &
CALL errore ('simple_diag','spin-orbit with no magnetization only',1)
IF ( .NOT.lspinorb ) skip_diag = .true.
npol = 2
ENDIF
ALLOCATE (vaux(ngm,nspin))
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
WRITE(stdout,'(a)') line
DO is=1,nspin
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
IF (debug) WRITE(stdout,'(a)') line
READ(iun,'(2e25.15)') (vaux(ig,is), ig=1,ngm)
! should be READ(iun,*) (vaux(ig,is), ig=1,ngm)
END DO
READ(iun,'(a)') line
READ(iun,*) okvan
IF ( okvan ) CALL errore ('simple_diag','US PP not implemented',1)
IF ( okvan ) skip_diag = .true.
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
WRITE(stdout,'(a)',advance='no') trim(line)
ALLOCATE (nh(nsp))
READ(iun,*) nh(:)
WRITE(stdout,*) nh(:)
nhm = maxval(nh(:))
IF ( lspinorb ) THEN
ALLOCATE (dvan_so(nhm,nhm,nspin,nsp) )
ELSE
ALLOCATE (dvan(nhm,nhm,nsp) )
END IF
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
WRITE(stdout,'(a)') line
DO nt = 1, nsp
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
IF (debug) WRITE(stdout,'(a)') line
READ(iun,*) nt_, nh(nt)
IF ( nspin /= 4 ) THEN
READ(iun,*) dvan(1:nh(nt),1:nh(nt),nt)
ELSE
IF ( lspinorb ) THEN
DO is=1,nspin
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
IF (debug) WRITE(stdout,'(a)') line
READ(iun,*) dvan_so(1:nh(nt),1:nh(nt),is,nt)
END DO
ELSE
READ(iun,*) dvan(1:nh(nt),1:nh(nt),nt)
END IF
END DO
READ(iun,'(a)') line
READ(iun,*) nkb
IF ( debug) WRITE(stdout,*) trim(line), nkb
WRITE(stdout,*) trim(line), nkb
!
READ(iun,'(a)') line
READ(iun,*) nks
IF ( debug) WRITE(stdout,*) trim(line), nks
READ(iun,'(a)') line
WRITE(stdout,*) trim(line), nks
READ(iun,'(a)') line
ALLOCATE ( ngk(nks) )
READ(iun,*) ngk(1:nks)
DO ik=1,nks
READ(iun,'(a)') line
READ(iun,*) xk(:,ik)
IF ( debug) WRITE(stdout,*) trim(line), xk(:,ik)
! For LSDA, read spin index for this k-point
IF ( nspin == 2 ) THEN
READ(line(21:),*) current_spin
IF (current_spin < 1 .or. current_spin > 2 ) &
CALL errore ('simple_diag','mismatch in spin index',ik)
ELSE
current_spin = 1
END IF
READ(iun,*) xk(:)
WRITE(stdout,*) trim(line), xk(:)
READ(iun,'(a)') line
READ(iun,*) npw
IF ( debug) WRITE(stdout,*) trim(line), npw
WRITE(stdout,*) trim(line), npw
IF ( npw /= ngk(ik) ) CALL errore ('simple_diag','mismatch in Npw',ik)
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
READ(iun,'(i8)') (igk_k(ig,ik), ig=1,npw)
WRITE(stdout,'(a)') line
ALLOCATE ( igk(npw) )
READ(iun,'(i8)') (igk(ig), ig=1,npw)
ALLOCATE ( vkb(npw,nkb) )
DO ikb=1,nkb
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
IF (debug) WRITE(stdout,'(a)') line
READ(iun,'(2e25.15)') vkb(1:npw,ikb)
END DO
READ(iun,'(a)') line
READ(iun,*) nbnd
IF ( debug) WRITE(stdout,*) trim(line), nbnd
WRITE(stdout,*) trim(line), nbnd
ALLOCATE ( et(nbnd), evc(npol*npw,nbnd) )
DO ibnd=1,nbnd
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
IF (debug) WRITE(stdout,'(a)') line
READ(iun,'(a)') line
IF ( debug) WRITE(stdout,'(a)') line
READ(line(19:),*) et(ibnd)
READ(iun,'(2e25.15)') evc(1:npw,ibnd)
IF ( npol == 2 ) &
READ(iun,'(2e25.15)') evc(1+npw:2*npw,ibnd)
END DO
WRITE(stdout,'("# Data read for k-point #",i4,", diagonalizing...")') ik
ALLOCATE ( h(npol*npw,npol*npw), v(npol*npw,npol*npw), et(npol*npw) )
WRITE(stdout,'("Data read for k-point #",i4,", eigenvalues :")') ik
WRITE(stdout,'(6f12.6)') et(:)
ALLOCATE ( h(npol*npw,npol*npw), v(npol*npw,npol*npw), e(npol*npw) )
h(:,:) = (0.0_dp,0.0_dp)
DO j=1,npw
!
! kinetic energy
!
g(:) = mill(1,igk_k(j,ik))*bg(:,1) + &
mill(2,igk_k(j,ik))*bg(:,2) + &
mill(3,igk_k(j,ik))*bg(:,3)
h(j,j)= ( xk(1,ik)+g(1) )**2 + &
( xk(2,ik)+g(2) )**2 + &
( xk(3,ik)+g(3) )**2
g(:) = mill(1,igk(j))*bg(:,1) + &
mill(2,igk(j))*bg(:,2) + &
mill(3,igk(j))*bg(:,3)
h(j,j)= ( xk(1)+g(1) )**2 + &
( xk(2)+g(2) )**2 + &
( xk(3)+g(3) )**2
IF ( npol == 2 ) h(npw+j,npw+j) = h(j,j)
!
! nonlocal potential
@ -384,13 +424,11 @@ SUBROUTINE simple_diag ( )
IF ( nt == ityp(na) ) THEN
DO ih=1,nh(nt)
DO jh=1,nh(nt)
IF ( nspin /= 4 ) THEN
DO i=j,npw
h(i,j) = h(i,j) + vkb(i,ikb+ih) * &
dvan(ih,jh,nt) *DCONJG(vkb(j,ikb+jh))
END DO
ELSE
DO i=j,npw
IF ( lspinorb ) THEN
!
! noncollinear spin-orbit
!
DO i=1,npw
h(i,j) = h(i,j) + vkb(i,ikb+ih) * &
dvan_so(ih,jh,1,nt) *DCONJG(vkb(j,ikb+jh))
h(i+npw,j+npw) = h(i+npw,j+npw) + vkb(i,ikb+ih)* &
@ -402,6 +440,24 @@ SUBROUTINE simple_diag ( )
h(i+npw,j) = h(i+npw,j) + vkb(i,ikb+ih) * &
dvan_so(ih,jh,3,nt) *DCONJG(vkb(j,ikb+jh))
END DO
ELSE IF ( npol == 2 ) THEN
!
! noncollinear but no spin-orbit
!
DO i=1,npw
h(i,j) = h(i,j) + vkb(i,ikb+ih) * &
dvan(ih,jh,nt) *DCONJG(vkb(j,ikb+jh))
h(i+npw,j+npw) = h(i+npw,j+npw) + vkb(i,ikb+ih)* &
dvan(ih,jh,nt) *DCONJG(vkb(j,ikb+jh))
END DO
ELSE
!
! collinear (LSDA) or unpolarized
!
DO i=1,npw
h(i,j) = h(i,j) + vkb(i,ikb+ih) * &
dvan(ih,jh,nt) *DCONJG(vkb(j,ikb+jh))
END DO
END IF
END DO
END DO
@ -409,13 +465,13 @@ SUBROUTINE simple_diag ( )
END IF
END DO
END DO
!
!
! local potential
!
DO i=j,npw
i1 = mill(1,igk_k(i,ik)) - mill(1,igk_k(j,ik))
i2 = mill(2,igk_k(i,ik)) - mill(2,igk_k(j,ik))
i3 = mill(3,igk_k(i,ik)) - mill(3,igk_k(j,ik))
DO i=1,npw
i1 = mill(1,igk(i)) - mill(1,igk(j))
i2 = mill(2,igk(i)) - mill(2,igk(j))
i3 = mill(3,igk(i)) - mill(3,igk(j))
IF ( ABS(i1) > SIZE(limm,1) .OR. &
ABS(i2) > SIZE(limm,2) .OR. &
ABS(i3) > SIZE(limm,3) ) &
@ -423,20 +479,48 @@ SUBROUTINE simple_diag ( )
ij = limm ( i1,i2,i3 )
IF ( ij <= 0 .OR. ij > ngm ) &
CALL errore ('simple_diag','internal error (2)',i)
DO ipol = 1, npol
ii = (ipol-1)*npw + i
jj = (ipol-1)*npw + j
h(ii,jj) = h(ii,jj) + vaux(ij,1)
IF ( i > j ) h(jj,ii) =DCONJG(h(ii,jj))
END DO
IF (npol == 2) THEN
ii = npw + i
jj = npw + j
IF (domag) THEN
!
! noncollinear magnetic
!
h( i, j) = h( i, j) + vaux(ij,1) + vaux(ij,4)
h( i,jj) = h( i,jj) + vaux(ij,2) - (0.0_dp,1.0_dp)*vaux(ij,3)
h(ii, j) = h( i,jj) + vaux(ij,2) + (0.0_dp,1.0_dp)*vaux(ij,3)
h(ii,jj) = h(ii,jj) + vaux(ij,1) - vaux(ij,4)
ELSE
!
! noncollinear (spin-orbit) not magnetic
!
h( i, j) = h( i, j) + vaux(ij,1)
h(ii,jj) = h(ii,jj) + vaux(ij,1)
END IF
ELSE
!
! collinear (LSDA) or unpolarized
!
h( i, j) = h( i, j) + vaux(ij,current_spin)
END IF
END DO
END DO
!
CALL cdiagh ( npol*npw, h, npol*npw, et, v)
WRITE(stdout,'(4f12.6)') (et(ibnd)*13.6058, ibnd=1,nbnd)
DEALLOCATE ( et, v, h )
IF ( .not. skip_diag ) THEN
CALL cdiagh ( npol*npw, h, npol*npw, e, v)
WRITE(stdout,'("Recomputed eigenvalues:")')
WRITE(stdout,'(6f12.6)') e(1:nbnd)*13.6058
ELSE
CALL infomsg ('simple_diag','not implemented: diagonalization skipped')
END IF
DEALLOCATE ( e, v, h, evc, et, vkb, igk )
END DO
DEALLOCATE (limm, vaux)
DEALLOCATE ( ngk, nh, vaux, limm, mill, ityp )
IF ( lspinorb ) THEN
DEALLOCATE ( dvan_so )
ELSE
DEALLOCATE ( dvan )
END IF
READ(iun,'(a)') line
WRITE(stdout,'(a)') line
!