Commit 1368dfc3 authored by Matthias Redies's avatar Matthias Redies

fix calc_cmt

parent e6189aff
......@@ -30,7 +30,7 @@ contains
complex, intent(inout) :: cmt_out(:,:,:)
complex, allocatable :: acof(:,:,:), bcof(:,:,:), ccof(:,:,:,:)
complex, allocatable :: cmt(:,:,:), cmt_out2(:,:,:)
complex, allocatable :: cmt(:,:,:)
integer :: ikp, nbands, ok(4) ! index of parent k-point
integer :: iatom, itype, ieq, indx, i, j, idum, iop, l, ll, lm, m
......@@ -39,7 +39,6 @@ contains
complex :: cdum
type(t_lapw) :: lapw_ik, lapw_ikp
allocate(cmt_out2, source=cmt_out)
ikp = kpts%bkp(ik)
nbands = hybdat%nbands(ikp)
......@@ -105,19 +104,16 @@ contains
END DO
END DO
! write cmt at irreducible k-points in direct-access file cmt
if(ik <= kpts%nkpt) then
cmt_out = cmt
else
iop = kpts%bksym(ik)
call waveftrafo_gen_cmt(cmt, c_phase, zmat_ikp%l_real, ikp, iop, atoms, &
mpdata, hybinp, kpts, sym, nbands, cmt_out)
endif
call read_cmt(cmt_out2, ik)
if(maxval(abs(cmt_out - cmt_out2))> 1e-10) then
write (*,*) "ik = ", ik
write (*,*) "max diff = ", maxval(abs(cmt_out - cmt_out2))
call juDFT_error("too big off a diff")
endif
end subroutine calc_cmt
end module m_calc_cmt
......@@ -88,8 +88,9 @@
if(allocated(c_phase)) deallocate(c_phase)
allocate(c_phase(hybdat%nbands(ikpt)))
if(ikpt /= kpts%bkp(ikpt)) call juDFT_error("We should be reading the parent z-mat here!")
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, ikpt, &
jsp, z(ikpt), c_phase)
jsp, z(ikpt), c_phase=c_phase)
call calc_cmt(atoms, cell, input, noco, hybinp, hybdat, mpdata, kpts, &
sym, oneD, z(kpts%bkp(ikpt)), jsp, ikpt, c_phase, &
cmt(:hybdat%nbands(ikpt),:,:,ikpt))
......
......@@ -83,7 +83,8 @@ CONTAINS
! read in mt wavefunction coefficients from file cmt
nbasfcn = calc_number_of_basis_functions(lapw, atoms, noco)
CALL zmat%init(sym%invs, nbasfcn, hybdat%nbands(kpts%bkp(nk)))
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, kpts%bkp(nk), jsp, zmat, c_phase)
if(nk /= kpts%bkp(nk)) call juDFT_error("We should be reading the parent z-mat here!")
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, kpts%bkp(nk), jsp, zmat, c_phase=c_phase)
call calc_cmt(atoms, cell, input, noco, hybinp, hybdat, mpdata, kpts, &
sym, oneD, zmat, jsp, nk, c_phase, cmt)
call zmat%free()
......
......@@ -19,7 +19,6 @@ CONTAINS
hybdat, noco, oneD, mpi, input, jsp, zmat)
! nkpt :: number of all k-points
USE m_radfun
USE m_radflo
USE m_abcof
......@@ -184,12 +183,12 @@ CONTAINS
! abcof calculates the wavefunction coefficients
! stored in acof,bcof,ccof
lapw(ikpt0)%nmat = lapw(ikpt0)%nv(jsp) + atoms%nlotot
CALL abcof(input, atoms, sym, cell, lapw(ikpt0), hybdat%nbands(ikpt0), hybdat%usdus, noco, jsp, &!hybdat%kveclo_eig(:,ikpt0),&
oneD, acof(:hybdat%nbands(ikpt0), :, :), bcof(:hybdat%nbands(ikpt0), :, :), &
ccof(:, :hybdat%nbands(ikpt0), :, :), zmat(ikpt0))
! MT wavefunction coefficients are calculated in a local coordinate system rotate them in the global one
CALL hyb_abcrot(hybinp, atoms, hybdat%nbands(ikpt0), sym, acof(:hybdat%nbands(ikpt0), :, :), &
bcof(:hybdat%nbands(ikpt0), :, :), ccof(:, :hybdat%nbands(ikpt0), :, :))
......@@ -239,6 +238,7 @@ CONTAINS
END DO
END DO
! write cmt at irreducible k-points in direct-access file cmt
CALL write_cmt(cmt, ikpt0)
CALL zhlp%alloc(zmat(ikpt0)%l_real, zmat(ikpt0)%matsize1, zmat(ikpt0)%matsize2)
......@@ -263,6 +263,7 @@ CONTAINS
call waveftrafo_gen_zmat(zmat(ikpt0), ikpt0, iop, &
kpts, sym, jsp, input, hybdat%nbands(ikpt0), &
lapw(ikpt0), lapw(ikpt), zhlp, c_phase)
call waveftrafo_gen_cmt(cmt, c_phase, zmat(ikpt0)%l_real, ikpt0, iop, atoms, &
mpdata, hybinp, kpts, sym, hybdat%nbands(ikpt0), cmthlp)
CALL write_cmt(cmthlp, ikpt)
......
......@@ -151,7 +151,9 @@ CONTAINS
call z%init(olap%l_real, nbasfcn, hybdat%nbands(nk))
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, nk, jsp, z, c_phase)
if(nk /= kpts%bkp(nk)) call juDFT_error("We should be reading the parent z-mat here!")
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, nk, jsp, z, c_phase=c_phase)
CALL timestart("symm_hf")
CALL symm_hf_init(sym, kpts, nk, nsymop, rrot, psym)
......
......@@ -837,7 +837,9 @@ CONTAINS
! read in cmt coefficients from direct access file cmt at kpoint nk
momentum = cmplx_0
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, nk, jsp, z, c_phase)
if(nk /= kpts%bkp(nk)) call juDFT_error("We should be reading the parent z-mat here!")
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, nk, jsp, z, c_phase=c_phase)
call calc_cmt(atoms, cell, input, noco, hybinp, hybdat, mpdata, kpts, &
sym, oneD, z, jsp, nk, c_phase, cmt)
......
This diff is collapsed.
......@@ -179,7 +179,7 @@ contains
CALL write_matrix(mat, rec, id_olap)
END subroutine write_olap
subroutine read_z(atoms, cell, hybdat, kpts, sym, noco, input, ik, jsp, z_out, c_phase)
subroutine read_z(atoms, cell, hybdat, kpts, sym, noco, input, ik, jsp, z_out, parent_z, c_phase)
USE m_eig66_io
use m_types
use m_trafo
......@@ -194,30 +194,39 @@ contains
integer, intent(in) :: ik, jsp
TYPE(t_mat), INTENT(INOUT) :: z_out
complex, intent(inout), optional :: c_phase(:)
type(t_mat), intent(inout), target, optional :: parent_z
complex, intent(inout), optional :: c_phase(:)
INTEGER :: ikp, iop
type(t_mat) :: tmp_mat
complex :: cmt(input%neig,hybdat%maxlmindx,atoms%nat)
complex :: cmthlp(input%neig,hybdat%maxlmindx,atoms%nat)
type(t_lapw) :: lapw_ik, lapw_ikp
INTEGER :: ikp, iop
type(t_mat), pointer :: ptr_mat
type(t_mat), target :: tmp_mat
complex :: cmt(input%neig,hybdat%maxlmindx,atoms%nat)
complex :: cmthlp(input%neig,hybdat%maxlmindx,atoms%nat)
type(t_lapw) :: lapw_ik, lapw_ikp
cmt=0;cmthlp=0
call timestart("read_z")
if(ik <= kpts%nkpt) then
call read_eig(hybdat%eig_id,ik,jsp,zmat=z_out)
if(present(parent_z)) call parent_z%copy(z_out,1,1)
else
if(present(parent_z)) then
ptr_mat => parent_z
else
call tmp_mat%init(z_out)
ptr_mat => tmp_mat
endif
ikp = kpts%bkp(ik) ! parrent k-point
iop = kpts%bksym(ik) ! connecting symm
call tmp_mat%init(z_out)
call read_eig(hybdat%eig_id,ikp, jsp,zmat=tmp_mat)
call read_eig(hybdat%eig_id,ikp, jsp,zmat=ptr_mat)
CALL lapw_ik%init(input, noco, kpts, atoms, sym, ik, cell, sym%zrfs)
CALL lapw_ikp%init(input, noco, kpts, atoms, sym, ikp, cell, sym%zrfs)
call waveftrafo_gen_zmat(tmp_mat, ikp, iop, kpts, sym, jsp, input, &
call waveftrafo_gen_zmat(ptr_mat, ikp, iop, kpts, sym, jsp, input, &
hybdat%nbands(ikp), lapw_ikp, lapw_ik, z_out, &
c_phase)
endif
......
......@@ -42,7 +42,8 @@ check_arguments.F90
check_para.F90
xmlOutput.F90
round.f90
#npy.f90
npy.f90
endian_swap.f90
round.f90
)
if (JUDFT_USE_MPI)
......
module endian_swap
implicit none
PRIVATE
PUBLIC :: Big_Endian
PUBLIC :: Swap_Endian
INTERFACE Swap_Endian
module procedure SWAP_I1
module procedure SWAP_I2
module procedure SWAP_I4
module procedure SWAP_I8
module procedure SWAP_F4
module procedure SWAP_F8
module procedure SWAP_F16
module procedure SWAP_C4
module procedure SWAP_C8
END INTERFACE Swap_Endian
CONTAINS
FUNCTION Big_Endian()
LOGICAL :: Big_Endian
Big_Endian = ichar(transfer(1,'a')) == 0
END FUNCTION Big_Endian
function SWAP_I4(input) result(output)
implicit none
integer(4), parameter :: b_sz = 4
integer(b_sz), intent(in) :: input
integer(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_I4
function SWAP_I2(input) result(output)
implicit none
integer(4), parameter :: b_sz = 2
integer(b_sz), intent(in) :: input
integer(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_I2
function SWAP_I1(input) result(output)
implicit none
integer(4), parameter :: b_sz = 1
integer(b_sz), intent(in) :: input
integer(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_I1
function SWAP_I8(input) result(output)
implicit none
integer(4), parameter :: b_sz = 8
integer(b_sz), intent(in) :: input
integer(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_I8
function SWAP_F4(input) result(output)
implicit none
integer(4), parameter :: b_sz = 4
real(b_sz), intent(in) :: input
real(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_F4
function SWAP_F8(input) result(output)
implicit none
integer(4), parameter :: b_sz = 8
real(b_sz), intent(in) :: input
real(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_F8
function SWAP_F16(input) result(output)
implicit none
integer(4), parameter :: b_sz = 16
real(b_sz), intent(in) :: input
real(b_sz) :: output
integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
integer(1) :: i
byte_arr_tmp = transfer(input, byte_arr_tmp)
do i = 1,b_sz
byte_arr(i) = byte_arr_tmp(1 + b_sz - i)
enddo
output = transfer(byte_arr, output)
end function SWAP_F16
function SWAP_C8(input) result(output)
implicit none
integer(4), parameter :: b_sz = 8
complex(b_sz), intent(in) :: input
complex(b_sz) :: output
real(b_sz) :: re, im
re = Swap_Endian(real(input))
im = Swap_Endian(aimag(input))
output = cmplx(re, im)
end function swap_c8
function SWAP_C4(input) result(output)
implicit none
integer(4), parameter :: b_sz = 4
complex(b_sz), intent(in) :: input
complex(b_sz) :: output
real(b_sz) :: re, im
re = Swap_Endian(real(input))
im = Swap_Endian(aimag(input))
output = cmplx(re, im)
end function swap_c4
END module endian_swap
......@@ -14,5 +14,6 @@ MODULE m_juDFT
USE m_judft_usage
use m_judft_para
use m_juDFT_round
use m_npy
! use m_npy
END MODULE m_juDFT
This diff is collapsed.
......@@ -532,7 +532,9 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
parent = 0
CALL zMat%init(olap%l_real,nbasfcn,input%neig)
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, ikpt, jsp, zMat, c_phase)
if(ikpt /= kpts%bkp(ikpt)) call juDFT_error("We should be reading the parent z-mat here!")
call read_z(atoms, cell, hybdat, kpts, sym, noco, input, ikpt, jsp, zMat, c_phase=c_phase)
CALL symm_hf_init(sym,kpts,ikpt,nsymop,rrot,psym)
CALL symm_hf(kpts,ikpt,sym,hybdat,eig_irr,input,atoms,mpdata,hybinp,cell,lapw,&
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment