subroutine plot_write_vdw_data(num_wann, wannier_data, real_lattice, u_matrix, u_matrix_opt, &
have_disentangled, w90_system, error, comm, stdout, seedname)
!================================================!
!
! Write a file with Wannier centres, spreads and occupations for
! post-processing computation of vdW C6 coeffients.
!
! Based on code written by Lampros Andrinopoulos.
!================================================!
use w90_constants, only: cmplx_0
use w90_constants, only: dp
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_input
use w90_io, only: io_date
use w90_types, only: wannier_data_type, w90_system_type
use w90_utility, only: utility_translate_home
implicit none
type(wannier_data_type), intent(in) :: wannier_data
type(w90_system_type), intent(in) :: w90_system
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), intent(in) :: u_matrix(:, :, :)
complex(kind=dp), intent(in) :: u_matrix_opt(:, :, :)
integer, intent(in) :: stdout
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
integer :: iw, vdw_unit, r, s, k, m, ierr, ndim
real(kind=dp) :: wc(3, num_wann)
real(kind=dp) :: ws(num_wann)
complex(kind=dp), allocatable :: f_w(:, :), v_matrix(:, :) !f_w2(:,:)
wc = wannier_data%centres
ws = wannier_data%spreads
! translate Wannier centres to the home unit cell
do iw = 1, num_wann
call utility_translate_home(wc(:, iw), real_lattice)
end do
allocate (f_w(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating f_w in plot_write_vdw_data', comm)
return
end if
!~ ! aam: remove f_w2 at end
!~ allocate(f_w2(num_wann, num_wann),stat=ierr)
!~ if (ierr/=0) call io_error('Error in allocating f_w2 in plot_write_vdw_data')
if (have_disentangled) then
! dimension of occupied subspace
if (w90_system%num_valence_bands .le. 0) then
call set_error_input(error, 'Please set num_valence_bands in seedname.win', comm)
return
end if
ndim = w90_system%num_valence_bands
allocate (v_matrix(ndim, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating V_matrix in plot_write_vdw_data', comm)
return
end if
! aam: initialise
f_w(:, :) = cmplx_0
v_matrix(:, :) = cmplx_0
!~ f_w2(:,:) = cmplx_0
! aam: IN THE END ONLY NEED DIAGONAL PART, SO COULD SIMPLIFY...
! aam: calculate V = U_opt . U
do s = 1, num_wann
do k = 1, ndim
do m = 1, num_wann
v_matrix(k, s) = v_matrix(k, s) + u_matrix_opt(k, m, 1)*u_matrix(m, s, 1)
end do
end do
end do
! aam: calculate f = V^dagger . V
do r = 1, num_wann
do s = 1, num_wann
do k = 1, ndim
f_w(r, s) = f_w(r, s) + v_matrix(k, s)*conjg(v_matrix(k, r))
end do
end do
end do
!~ ! original formulation
!~ do r=1,num_wann
!~ do s=1,num_wann
!~ do nkp=1,num_kpts
!~ do k=1,ndimfroz(nkp)
!~ do m=1,num_wann
!~ do l=1,num_wann
!~ f_w2(r,s) = f_w2(r,s) + &
!~ u_matrix_opt(k,m,nkp) * u_matrix(m,s,nkp) * &
!~ conjg(u_matrix_opt(k,l,nkp)) * conjg(u_matrix(l,r,nkp))
!~ end do
!~ end do
!~ end do
!~ end do
!~ end do
!~ end do
!~ ! test equivalence
!~ do r=1,num_wann
!~ do s=1,num_wann
!~ if (abs(real(f_w(r,s),dp)-real(f_w2(r,s),dp)).gt.eps6) then
!~ write(*,'(i6,i6,f16.10,f16.10)') r,s,real(f_w(r,s),dp),real(f_w2(r,s),dp)
!~ end if
!~ if (abs(aimag(f_w(r,s))-aimag(f_w2(r,s))).gt.eps6) then
!~ write(*,'(a,i6,i6,f16.10,f16.10)') 'Im: ',r,s,aimag(f_w(r,s)),aimag(f_w2(r,s))
!~ end if
!~ end do
!~ end do
else
! for valence only, all occupancies are unity
f_w(:, :) = 1.0_dp
end if
! aam: write the seedname.vdw file directly here
open (newunit=vdw_unit, file=trim(seedname)//'.vdw', action='write')
if (have_disentangled) then
write (vdw_unit, '(a)') 'disentangle T'
else
write (vdw_unit, '(a)') 'disentangle F'
end if
write (vdw_unit, '(a)') 'amalgamate F'
write (vdw_unit, '(a,i3)') 'degeneracy', w90_system%num_elec_per_state
write (vdw_unit, '(a)') 'num_frag 2'
write (vdw_unit, '(a)') 'num_wann'
write (vdw_unit, '(i3,1x,i3)') num_wann/2, num_wann/2
write (vdw_unit, '(a)') 'tol_occ 0.9'
write (vdw_unit, '(a)') 'pxyz'
write (vdw_unit, '(a)') 'F F F'
write (vdw_unit, '(a)') 'F F F'
write (vdw_unit, '(a)') 'tol_dist 0.05'
write (vdw_unit, '(a)') 'centres_spreads_occ'
write (vdw_unit, '(a)') 'ang'
do iw = 1, num_wann
write (vdw_unit, '(4(f13.10,1x),1x,f11.8)') wc(1:3, iw), ws(iw), real(f_w(iw, iw))
end do
close (vdw_unit)
write (stdout, '(/a/)') ' vdW data written to file '//trim(seedname)//'.vdw'
if (have_disentangled) then
deallocate (v_matrix, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating v_matrix in plot_write_vdw_data', comm)
return
end if
end if
deallocate (f_w, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating f_w in plot_write_vdw_data', comm)
return
end if
return
end subroutine plot_write_vdw_data
