TABLE OF CONTENTS


ABINIT/m_ddk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddk

FUNCTION

  Objects and methods to extract data from DDK files.
  The DDK files are binary (soon also netcdf) files with Hamiltonian derivatives
  wrt k, and the corresponding wave functions

COPYRIGHT

 Copyright (C) 2016-2022 ABINIT group (MJV, HM, MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 MODULE m_ddk
26 
27  use defs_basis
28  use m_abicore
29  use m_errors
30  use m_xmpi
31  use m_nctk
32  use m_hdr
33  use m_dtset
34  use m_krank
35  use m_crystal
36  use m_mpinfo
37  use m_cgtools
38  use m_hamiltonian
39  use m_initylmg
40  use m_ebands
41  use m_pawcprj
42  use m_getgh1c
43  use netcdf
44 
45  use m_fstrings,      only : strcat, sjoin, itoa, ktoa
46  use m_io_tools,      only : iomode_from_fname
47  use m_time,          only : cwtime, cwtime_report
48  use defs_abitypes,   only : MPI_type
49  use defs_datatypes,  only : ebands_t, pseudopotential_type
50  use m_vkbr,          only : vkbr_t, nc_ihr_comm, vkbr_init, vkbr_free
51  use m_pawtab,        only : pawtab_type
52  use m_wfk,           only : wfk_read_ebands !, wfk_read_h1mat
53  use m_wfd,           only : wfd_t, wfd_init, wave_t
54 
55  implicit none
56 
57  private

m_ddk/ddk_red2car [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_red2car

FUNCTION

  Convert ddk matrix element from reduced coordinates to cartesian coordinates.

INPUTS

OUTPUT

SOURCE

655 pure subroutine ddk_red2car(rprimd, vred, vcar)
656 
657 !Arguments -------------------------------------
658  real(dp),intent(in) :: rprimd(3,3)
659  real(dp),intent(in) :: vred(2,3)
660  real(dp),intent(out) :: vcar(2,3)
661 
662 !Local variables -------------------------------
663  real(dp) :: vtmp(2,3)
664 
665 !************************************************************************
666 
667  ! Go to Cartesian coordinates (same as pmat2cart routine)
668  ! V_cart = 1/(2pi) * Rprimd x V_red
669  ! where V_red is the derivative computed in the DFPT routines (derivative wrt reduced component).
670  vtmp(1,:) = rprimd(:,1)*vred(1,1) &
671             +rprimd(:,2)*vred(1,2) &
672             +rprimd(:,3)*vred(1,3)
673  vtmp(2,:) = rprimd(:,1)*vred(2,1) &
674             +rprimd(:,2)*vred(2,2) &
675             +rprimd(:,3)*vred(2,3)
676  vcar = vtmp / two_pi
677 
678 end subroutine ddk_red2car

m_ddk/ddkop_apply [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_apply

FUNCTION

  Apply velocity operator dH/dk to wavefunction in G-space. Store results in object.

INPUTS

  eig0nk: Eigenvalue associated to the wavefunction.
  npw_k: Number of planewaves.
  nspinor: Number of spinor components.
  cwave(2,npw_k*nspinor)=input wavefunction in reciprocal space
  cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives)
     if not allocated or size=0, they are locally computed (and not sorted)!!

SIDE EFFECTS

 Stores:
  gh1c(2,npw1*nspinor)= <G|H^(1)|C> or  <G|H^(1)-lambda.S^(1)|C> on the k+q sphere
                        (only kinetic+non-local parts if optlocal=0)

SOURCE

875 subroutine ddkop_apply(self, eig0nk, npw_k, nspinor, cwave, cwaveprj)
876 
877 !Arguments ------------------------------------
878 !scalars
879  class(ddkop_t),intent(inout) :: self
880  integer,intent(in) :: npw_k, nspinor
881  real(dp),intent(in) :: eig0nk
882 !arrays
883  real(dp),intent(inout) :: cwave(2,npw_k*nspinor)
884  type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
885 
886 !Local variables-------------------------------
887 !scalars
888  integer,parameter :: berryopt0 = 0, optlocal0 = 0, tim_getgh1c = 1, usevnl0 = 0, opt_gvnlx1 = 0
889  integer :: idir, sij_opt, ispinor, ipws, ipw, optnl
890  real(dp) :: eshift
891 !arrays
892  real(dp) :: grad_berry(2,(berryopt0/4)), gvnlx1(2,usevnl0)
893  real(dp),pointer :: dkinpw(:),kinpw1(:)
894 
895 !************************************************************************
896 
897  self%eig0nk = eig0nk
898 
899  if (self%inclvkb /= 0) then
900    ! optlocal0 = 0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
901    ! optnl = 2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
902    ! opt_gvnlx1 = option controlling the use of gvnlx1 array:
903    optnl = 2 !; if (self%inclvkb == 0) optnl = 0
904 
905    eshift = self%eig0nk - self%dfpt_sciss
906    do idir=1,3
907      sij_opt = self%gs_hamkq(idir)%usepaw
908      call getgh1c(berryopt0, cwave, cwaveprj, self%gh1c(:,:,idir), &
909        grad_berry, self%gs1c(:,:,idir), self%gs_hamkq(idir), gvnlx1, idir, self%ipert, eshift, self%mpi_enreg, optlocal0, &
910        optnl, opt_gvnlx1, self%rf_hamkq(idir), sij_opt, tim_getgh1c, usevnl0)
911    end do
912 
913  else
914    ! optnl 0 with DDK does not work as expected.
915    ! So I treat the kinetic term explicitly without calling getgh1c.
916    do idir=1,3
917      kinpw1 => self%gs_hamkq(idir)%kinpw_kp
918      dkinpw => self%rf_hamkq(idir)%dkinpw_k
919      do ispinor=1,nspinor
920        do ipw=1,npw_k
921          ipws = ipw + npw_k*(ispinor-1)
922          if (kinpw1(ipw) < huge(zero)*1.d-11) then
923            self%gh1c(:,ipws,idir) = dkinpw(ipw) * cwave(:,ipws)
924          else
925            self%gh1c(:,ipws,idir) = zero
926          end if
927        end do
928      end do
929    end do
930  end if
931 
932 end subroutine ddkop_apply

m_ddk/ddkop_free [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_free

FUNCTION

  Free memory

INPUTS

SOURCE

1103 subroutine ddkop_free(self)
1104 
1105 !Arguments ------------------------------------
1106 !scalars
1107  class(ddkop_t),intent(inout) :: self
1108 
1109 !Local variables-------------------------------
1110 !scalars
1111  integer :: idir
1112 
1113 !************************************************************************
1114 
1115  ABI_SFREE(self%gh1c)
1116  ABI_SFREE(self%gs1c)
1117 
1118  do idir=1,3
1119    call self%gs_hamkq(idir)%free()
1120    call self%htg(idir)%free()
1121    call self%rf_hamkq(idir)%free()
1122  end do
1123 
1124  self%mpi_enreg => null()
1125 
1126 end subroutine ddkop_free

m_ddk/ddkop_get_braket [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_get_braket

FUNCTION

  Compute diagonal matrix element in Cartesian coordinates.

INPUTS

  eig0mk: Eigenvalue associated to the "bra" wavefunction
  istwkf_k: defines storage of wavefunctions for this k-point
  npw_k: Number of planewaves.
  nspinor: Number of spinor components.
  brag(2,npw_k*nspinor)=input wavefunction in reciprocal space

SOURCE

 953 function ddkop_get_braket(self, eig0mk, istwf_k, npw_k, nspinor, brag, mode) result(vk)
 954 
 955 !Arguments ------------------------------------
 956 !scalars
 957  class(ddkop_t),intent(in) :: self
 958  integer,intent(in) :: istwf_k, npw_k, nspinor
 959  real(dp),intent(in) :: eig0mk
 960  character(len=*),optional,intent(in) :: mode
 961 !arrays
 962  real(dp),intent(in) :: brag(2*npw_k*nspinor)
 963  real(dp) :: vk(2,3)
 964 
 965 !Local variables-------------------------------
 966 !scalars
 967  integer :: idir
 968  real(dp) :: doti
 969 !arrays
 970  real(dp) :: dotarr(2), vk_red(2, 3)
 971  character(len=50) :: my_mode
 972 
 973 !************************************************************************
 974 
 975  if (self%usepaw == 0) then
 976    ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)>  (NC psps)
 977    do idir=1,3
 978      dotarr = cg_zdotc(npw_k * nspinor, brag, self%gh1c(:,:,idir))
 979      if (istwf_k > 1) then
 980        !dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term)
 981        doti = two * dotarr(2)
 982        if (istwf_k == 2 .and. self%mpi_enreg%me_g0 == 1) then
 983          ! nspinor always 1
 984          ! TODO: Recheck this part but it should be ok.
 985          doti = doti - (brag(1) * self%gh1c(2,1,idir) - brag(2) * self%gh1c(1,1,idir))
 986        end if
 987        dotarr(2) = doti; dotarr(1) = zero
 988      end if
 989      vk(:, idir) = dotarr
 990    end do
 991  else
 992    ABI_ERROR("PAW Not Implemented")
 993    ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(jband,k)^(0)> (PAW)
 994    ! eshiftkq = half * (eig0mk - self%eig0nk)
 995    ABI_UNUSED(eig0mk)
 996  end if
 997 
 998  my_mode = "cart"; if (present(mode)) my_mode = mode
 999  select case (mode)
1000  case ("cart")
1001    vk_red = vk
1002    call ddk_red2car(self%rprimd, vk_red, vk)
1003  case ("reduced")
1004    continue
1005  case default
1006    ABI_ERROR(sjoin("Invalid vaue for mode:", mode))
1007  end select
1008 
1009 end function ddkop_get_braket

m_ddk/ddkop_get_vdiag [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_get_vdiag

FUNCTION

  Simplified interface to compute the diagonal matrix element of the velocity operator in cartesian coords.

INPUTS

SOURCE

1025 function ddkop_get_vdiag(self, eig0nk, istwf_k, npw_k, nspinor, cwave, cwaveprj, mode) result(vk)
1026 
1027 !Arguments ------------------------------------
1028 !scalars
1029  class(ddkop_t),intent(inout) :: self
1030  integer,intent(in) :: istwf_k, npw_k, nspinor
1031  real(dp),intent(in) :: eig0nk
1032  character(len=*),optional,intent(in) :: mode
1033 !arrays
1034  real(dp),intent(inout) :: cwave(2,npw_k*nspinor)
1035  type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
1036  real(dp) :: vk(3)
1037 
1038 !Local variables-------------------------------
1039  character(len=50) :: my_mode
1040 !arrays
1041  real(dp) :: cvk(2, 3)
1042 
1043 !************************************************************************
1044 
1045  my_mode = "cart"; if (present(mode)) my_mode = mode
1046  call self%apply(eig0nk, npw_k, nspinor, cwave, cwaveprj)
1047  cvk = self%get_braket(eig0nk, istwf_k, npw_k, nspinor, cwave, mode=my_mode)
1048  vk = cvk(1, :)
1049 
1050 end function ddkop_get_vdiag

m_ddk/ddkop_get_vnondiag [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_get_vdiag

FUNCTION

  Simplified interface to compute the off-diagonal matrix elemente of the velocity operator in cartesian coords.

INPUTS

SOURCE

1064 function ddkop_get_vnondiag(self, eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, cwave_ket, cwaveprj, mode) result(cvk)
1065 
1066 !Arguments ------------------------------------
1067 !scalars
1068  class(ddkop_t),intent(inout) :: self
1069  integer,intent(in) :: istwf_k, npw_k, nspinor
1070  real(dp),intent(in) :: eig0nk_bra
1071  character(len=*),optional,intent(in) :: mode
1072 !arrays
1073  real(dp),intent(inout) :: cwave_bra(2,npw_k*nspinor),cwave_ket(2,npw_k*nspinor)
1074  type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
1075  real(dp) :: cvk(2,3)
1076 
1077 !Local variables-------------------------------
1078  character(len=50) :: my_mode
1079 !arrays
1080 
1081 !************************************************************************
1082 
1083  my_mode = "cart"; if (present(mode)) my_mode = mode
1084  call self%apply(eig0nk_bra, npw_k, nspinor, cwave_ket, cwaveprj)
1085  cvk = self%get_braket(eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, mode=my_mode)
1086 
1087 end function ddkop_get_vnondiag

m_ddk/ddkop_new [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_new

FUNCTION

  Build new object. Use dtset%inclvkb to determine whether non-local part should be included.

INPUTS

 dtset<dataset_type>=All input variables for this dataset.
 cryst<crystal_t>=Crystal structure.
 pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
 psps<pseudopotential_type>=Variables related to pseudopotentials.
 mpi_enreg=information about MPI parallelization
 mpw=Maximum number of plane-waves over k-points.
 ngfft(18)=contain all needed information about 3D FFT

OUTPUT

SOURCE

703 type(ddkop_t) function ddkop_new(dtset, cryst, pawtab, psps, mpi_enreg, mpw, ngfft) result(new)
704 
705 !Arguments ------------------------------------
706 !scalars
707  type(dataset_type),intent(in) :: dtset
708  type(crystal_t),intent(in) :: cryst
709  type(pseudopotential_type),intent(in) :: psps
710  type(MPI_type),target,intent(in) :: mpi_enreg
711  integer,intent(in) :: mpw
712 !arrays
713  integer,intent(in) :: ngfft(18)
714  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
715 
716 !Local variables-------------------------------
717 !scalars
718  integer,parameter :: cplex1 = 1
719  integer :: nfft, mgfft, idir
720 
721 ! *************************************************************************
722 
723  ABI_CHECK(dtset%usepaw == 0, "PAW not tested/implemented!")
724 
725  new%inclvkb = dtset%inclvkb
726  new%usepaw = dtset%usepaw
727  new%ipert = cryst%natom + 1
728  new%dfpt_sciss = dtset%dfpt_sciss
729  new%mpw = mpw
730 
731  new%rprimd = cryst%rprimd
732  new%mpi_enreg => mpi_enreg
733 
734  ! Not used because vlocal1 is not applied.
735  nfft = product(ngfft(1:3))
736  mgfft = maxval(ngfft(1:3))
737 
738  ABI_MALLOC(new%gh1c, (2, new%mpw*dtset%nspinor, 3))
739  ABI_MALLOC(new%gs1c, (2, new%mpw*dtset%nspinor, 3))
740 
741  do idir=1,3
742    ! ==== Initialize most of the Hamiltonian (and derivative) ====
743    ! 1) Allocate all arrays and initialize quantities that do not depend on k and spin.
744    ! 2) Perform the setup needed for the non-local factors:
745    ! * Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
746    ! * PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
747    call init_hamiltonian(new%gs_hamkq(idir), psps, pawtab, dtset%nspinor, dtset%nsppol, dtset%nspden, cryst%natom,&
748      cryst%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg)
749      !paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
750      !usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
751 
752    ! Prepare application of the NL part.
753    call init_rf_hamiltonian(cplex1, new%gs_hamkq(idir), new%ipert, new%rf_hamkq(idir), has_e1kbsc=.true.)
754      !&paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
755      !&mpi_spintab=mpi_enreg%my_isppoltab)
756  end do
757 
758 end function ddkop_new

m_ddk/ddkop_setup_spin_kpoint [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkop_setup_spin_kpoint

FUNCTION

  Prepare internal tables that depend on k-point/spin

INPUTS

  dtset<dataset_type>=All input variables for this dataset.
  cryst<crystal_t>=Crystal structure.
  psps<pseudopotential_type>=Variables related to pseudopotentials.
  spin: spin index
  kpoint(3): K-point in reduced coordinates.
  istwkf_k: defines storage of wavefunctions for this k-point
  npw_k: Number of planewaves.
  kg_k(3,npw_k)=reduced planewave coordinates.

SOURCE

782 subroutine ddkop_setup_spin_kpoint(self, dtset, cryst, psps, spin, kpoint, istwf_k, npw_k, kg_k)
783 
784 !Arguments ------------------------------------
785 !scalars
786  integer,intent(in) :: spin, npw_k, istwf_k
787  class(ddkop_t),intent(inout) :: self
788  type(crystal_t) :: cryst
789  type(dataset_type),intent(in) :: dtset
790  type(pseudopotential_type),intent(in) :: psps
791 !arrays
792  integer,intent(in) :: kg_k(3,npw_k)
793  real(dp),intent(in) :: kpoint(3)
794 
795 !Local variables-------------------------------
796 !scalars
797  integer,parameter :: nkpt1=1, nsppol1=1
798  type(mpi_type) :: mpienreg_seq
799 !arrays
800  integer :: npwarr(nkpt1), dummy_nband(nkpt1*nsppol1)
801  integer :: idir, nkpg, nkpg1, useylmgr1, optder !, nylmgr1
802  real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:)
803 
804 !************************************************************************
805 
806  ABI_CHECK(npw_k <= self%mpw, "npw_k > mpw!")
807  self%kpoint = kpoint
808 
809  ! Set up the spherical harmonics (Ylm) at k+q if useylm = 1
810  useylmgr1 = 0; optder = 0
811  if (psps%useylm == 1) then
812    useylmgr1 = 1; optder = 1
813  end if
814 
815  ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm))
816  ABI_MALLOC(ylmgr1_k, (npw_k, 3+6*(optder/2), psps%mpsang**2*psps%useylm*useylmgr1))
817 
818  if (psps%useylm == 1) then
819    ! Fake MPI_type for sequential part. dummy_nband and nsppol1 are not used in sequential mode.
820    call initmpi_seq(mpienreg_seq)
821    dummy_nband = 0; npwarr = npw_k
822    call initylmg(cryst%gprimd, kg_k, kpoint, nkpt1, mpienreg_seq, psps%mpsang, npw_k, dummy_nband, nkpt1, &
823       npwarr, nsppol1, optder, cryst%rprimd, ylm_k, ylmgr1_k)
824    call destroy_mpi_enreg(mpienreg_seq)
825  end if
826 
827  do idir=1,3
828    call self%htg(idir)%free()
829 
830    ! Continue to initialize the Hamiltonian
831    call self%gs_hamkq(idir)%load_spin(spin, with_nonlocal=.true.)
832    call self%rf_hamkq(idir)%load_spin(spin, with_nonlocal=.true.)
833 
834    ! We need ffnl1 and dkinpw for 3 dirs. Note that the Hamiltonian objects use pointers to keep a reference
835    ! to the output results of this routine.
836    ! This is the reason why we need to store the targets in self%htg
837    call getgh1c_setup(self%gs_hamkq(idir), self%rf_hamkq(idir), dtset, psps, kpoint, kpoint, idir, self%ipert, & ! In
838      cryst%natom, cryst%rmet, cryst%gprimd, cryst%gmet, istwf_k, npw_k, npw_k, &            ! In
839      useylmgr1, kg_k, ylm_k, kg_k, ylm_k, ylmgr1_k, &                                       ! In
840      self%htg(idir)%dkinpw, nkpg, nkpg1, self%htg(idir)%kpg_k, self%htg(idir)%kpg1_k, &     ! Out
841      self%htg(idir)%kinpw1, self%htg(idir)%ffnlk, self%htg(idir)%ffnl1, &                   ! Out
842      self%htg(idir)%ph3d, self%htg(idir)%ph3d1)                                             ! Out
843  end do
844 
845  ABI_FREE(ylm_k)
846  ABI_FREE(ylmgr1_k)
847 
848 end subroutine ddkop_setup_spin_kpoint

m_ddk/ddkop_t [ Types ]

[ Top ] [ m_ddk ] [ Types ]

NAME

  ddkop_t

FUNCTION

  This object provides a simplified interface to compute matrix elements of the
  velocity operator with the DFPT routines.

SOURCE

 84  type,public :: ddkop_t
 85 
 86   integer :: ipert
 87   ! Perturbation type: natom + 1
 88 
 89   integer :: inclvkb
 90   ! Option for calculating the matrix elements of [Vnl,r].
 91   ! 0 to exclude commutator, 2 to include it
 92 
 93   integer :: usepaw
 94   ! 0 for NC, 1 for PAW
 95 
 96   integer :: mpw
 97   ! Maximum number of plane-waves over k-points (used to dimension arrays)
 98 
 99   real(dp) :: kpoint(3)
100   ! K-point (set in setup_spin_kpoint)
101 
102   real(dp) :: eig0nk
103 
104   real(dp) :: dfpt_sciss = zero
105 
106   real(dp) :: rprimd(3,3)
107 
108   type(MPI_type),pointer :: mpi_enreg => null()
109 
110   type(gs_hamiltonian_type) :: gs_hamkq(3)
111 
112   type(rf_hamiltonian_type) :: rf_hamkq(3)
113 
114   type(ham_targets_t), private :: htg(3)
115   ! Store arrays targetted by the hamiltonians.
116 
117   real(dp), allocatable :: gh1c(:,:,:)
118    !gh1c, (2, mpw*nspinor, 3))
119 
120   real(dp), allocatable :: gs1c(:,:,:)
121    ! gs1c, (2, mpw*nspinor, 3*psps%usepaw))
122 
123  contains
124 
125    procedure :: setup_spin_kpoint => ddkop_setup_spin_kpoint
126     ! Prepare application of dH/dk for given spin, k-point.
127 
128    procedure :: apply => ddkop_apply
129     ! Apply dH/dk to input wavefunction.
130 
131    procedure :: get_braket => ddkop_get_braket
132     ! Compute matrix element (complex results) in cartesian coords.
133 
134    procedure :: get_vdiag => ddkop_get_vdiag
135     ! Compute diagonal matrix element (real) in cartesian coords.
136 
137    procedure :: get_vnondiag => ddkop_get_vnondiag
138     ! Compute off diagonal matrix elements in cartesian coords.
139 
140    procedure :: free => ddkop_free
141     ! Free memory.
142 
143  end type ddkop_t
144 
145  public :: ddkop_new   ! Build object
146 
147  !interface ddkop_t
148  !  procedure ddkop_new
149  !end interface ddkop_t

m_ddk/ddkstore_compute_ddk [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkstore_compute_ddk

FUNCTION

  Calculate the DDK matrix elements using the commutator formulation.

INPUTS

  prefix: Prefix for output EVK file. Empty if output files are not wanted

SOURCE

211 subroutine ddkstore_compute_ddk(ds, wfk_path, prefix, dtset, psps, pawtab, ngfftc, comm)
212 
213 !Arguments ------------------------------------
214 !scalars
215  class(ddkstore_t),intent(inout) :: ds
216  character(len=*),intent(in) :: wfk_path, prefix
217  integer,intent(in) :: comm
218  type(dataset_type),intent(in) :: dtset
219  type(pseudopotential_type),intent(in) :: psps
220  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
221 !arrays
222  integer,intent(in) :: ngfftc(18)
223 
224 !Local variables ------------------------------
225 !scalars
226  integer,parameter :: master = 0
227  integer :: mband, nbcalc, nsppol, ib_v, ib_c, mpw, spin, nspinor, nkpt, nband_k, npw_k
228  integer :: ii, ik, bmin, bmax, istwf_k, idir, my_rank, nproc, ierr, bstop
229  real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all
230  integer :: ncerr, ncid
231  character(len=500) :: msg
232  character(len=fnlen) :: fname
233  logical :: write_ncfile
234  type(wfd_t) :: wfd
235  type(vkbr_t) :: vkbr
236  type(ebands_t) :: ebands
237  type(crystal_t) :: cryst
238  type(hdr_type) :: tmp_hdr, hdr
239  type(ddkop_t) :: ddkop
240  type(wave_t),pointer :: wave_v, wave_c
241 !arrays
242  integer,allocatable :: distrib_mat(:,:,:,:), distrib_diago(:,:,:), nband(:,:), kg_k(:,:)
243  logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
244  real(dp) :: kpt(3), vv(2, 3)
245  real(dp),allocatable :: cg_c(:,:), cg_v(:,:)
246  complex(dpc) :: vg(3), vr(3)
247  complex(gwpc),allocatable :: ihrc(:,:), ug_c(:), ug_v(:)
248  type(pawcprj_type),allocatable :: cwaveprj(:,:)
249 
250 !************************************************************************
251 
252  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
253 
254  if (my_rank == master) call wrtout([std_out, ab_out], " Computation of velocity matrix elements (DDK)", newlines=1)
255 
256  ABI_CHECK(psps%usepaw == 0, "PAW not implemented")
257 
258  ! Get ebands and hdr from WFK file.
259  ebands = wfk_read_ebands(wfk_path, comm, out_hdr=hdr)
260  cryst = hdr%get_crystal()
261 
262  ! Extract important dimensions from hdr%
263  nkpt    = hdr%nkpt
264  nsppol  = hdr%nsppol
265  nspinor = hdr%nspinor
266  mband   = hdr%mband
267 
268  ! Define band range
269  ! TODO: Perhaps one should allocate output arrays using bmin:bmax
270  ! and allow for nc output only if bmin == 1 and bmax == mband
271  if (ds%bmax == -1) ds%bmax = mband
272 
273  if (ds%bmin < 1 .or. ds%bmin > mband .or. ds%bmax > mband .or. ds%bmin > ds%bmax) then
274    ABI_ERROR(sjoin("Invalid value for bmin, bmax", itoa(ds%bmin), itoa(ds%bmax), "with mband:", itoa(mband)))
275  end if
276 
277  bmin = ds%bmin; bmax = ds%bmax
278  nbcalc  = bmax - bmin + 1
279  write_ncfile = len_trim(prefix) > 0
280  if (write_ncfile .and. .not. (bmin == 1 .and. bmax == mband) ) then
281    write_ncfile = .False.
282    ABI_WARNING("Cannot write ncfile if .not. (bmin == 1 .and. bmax == mband)")
283  end if
284 
285  if (my_rank == master) then
286    write(ab_out, "(a)")" Parameters extracted from the Abinit header:"
287    write(ab_out, "(a, f5.1)") '    ecut:    ', hdr%ecut
288    write(ab_out, "(a, i0)")   '    nkpt:    ', nkpt
289    write(ab_out, "(a, i0)")   '    mband:   ', mband
290    write(ab_out, "(a, i0)")   '    nsppol:  ', nsppol
291    write(ab_out, "(a, i0)")   '    nspinor: ', nspinor
292    write(ab_out, "(a, i0)")   '    useylm:  ', dtset%useylm
293    write(ab_out, "(a, i0)")   '    inclvkb: ', dtset%inclvkb
294    write(ab_out, "(2(a, i0))")'    bmin: ', bmin, ", bmax: ", bmax
295    if (ds%only_diago) then
296      write(ab_out, "(a)")'    Computing diagonal matrix elements only'
297    else
298      write(ab_out, "(a)")'    Computing diagonal and off-diagonal matrix elements'
299    end if
300    write(ab_out, "(2(a, i0))")'    Between band index bmin: ', bmin, ", bmax: ", bmax
301    write(ab_out, "(a)")""
302  end if
303 
304  ! Create distribution of the wavefunctions mask.
305  ABI_MALLOC(nband, (nkpt, nsppol))
306  ABI_MALLOC(keep_ur, (mband, nkpt, nsppol))
307  ABI_MALLOC(bks_mask, (mband, nkpt, nsppol))
308  keep_ur = .false.; bks_mask = .false.; nband = mband
309 
310  if (ds%only_diago) then
311    ! Distribute k-points, spin and (b, b) diagonal over MPI processors.
312    ABI_MALLOC(distrib_diago, (bmin:bmax, nkpt, nsppol))
313    distrib_diago = -1
314 
315    ! Create bks_mask to load the wavefunctions.
316    ii = 0
317    do spin=1,nsppol
318      do ik=1,nkpt
319        do ib_v=bmin,bmax
320           ii = ii + 1; if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism.
321           distrib_diago(ib_v, ik, spin) = my_rank
322           bks_mask(ib_v, ik, spin) = .true.
323        end do
324      end do
325    end do
326    call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_diago == my_rank))))
327 
328  else
329    ! Distribute k-points, spin and (b, b') pairs over the processors
330    ABI_MALLOC(distrib_mat, (bmin:bmax, bmin:bmax, nkpt, nsppol))
331    call xmpi_distab(nproc, distrib_mat)
332 
333    ! Create bks_mask to load the wavefunctions
334    do spin=1,nsppol
335      do ik=1,nkpt
336        ! Loop over v bands
337        do ib_v=bmin,bmax
338         ! Loop over c bands
339          do ib_c=bmin,bmax
340            if (distrib_mat(ib_c, ib_v, ik, spin) == my_rank) then
341              bks_mask(ib_v, ik, spin) = .true.
342              bks_mask(ib_c, ik, spin) = .true.
343            end if
344          end do
345        end do
346      end do
347    end do
348 
349    call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_mat == my_rank))))
350  end if
351 
352  ! Initialize distributed wavefunctions object
353  call wfd_init(wfd, cryst, pawtab, psps, keep_ur, mband, nband, nkpt, nsppol,&
354    bks_mask, dtset%nspden, nspinor, hdr%ecut, dtset%ecutsm, dtset%dilatmx, ebands%istwfk, ebands%kptns,&
355    ngfftc, dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
356 
357  ABI_FREE(bks_mask)
358  ABI_FREE(keep_ur)
359  ABI_FREE(nband)
360 
361  call wfd%print(header="Wavefunctions on the k-points grid")
362 
363  ! Read wavefunctions from WFK file.
364  call wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path))
365 
366  ! Allocate workspace arrays
367  mpw = maxval(wfd%npwarr)
368  ABI_MALLOC(kg_k, (3, mpw))
369  ABI_MALLOC(ug_c, (mpw*nspinor))
370  ABI_MALLOC(ug_v, (mpw*nspinor))
371  if (dtset%useria /= 666) then
372    ABI_MALLOC(cg_c, (2, mpw*nspinor))
373    ABI_MALLOC(cg_v, (2, mpw*nspinor))
374  end if
375 
376  ABI_MALLOC(cwaveprj, (0, 0))
377  ABI_CALLOC(ds%dipoles, (3, 2, bmin:bmax, bmin:bmax, nkpt, nsppol))
378  ABI_MALLOC(ihrc, (3, nspinor**2))
379 
380  if (ds%only_diago) then
381    ABI_CALLOC(ds%vdiago, (3, bmin:bmax, nkpt, nsppol))
382  else
383    ABI_CALLOC(ds%vmat, (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol))
384  end if
385 
386  if (dtset%useria /= 666) then
387    ddkop = ddkop_new(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft)
388    !if (my_rank == master) call ddkop%print(ab_out)
389  end if
390 
391  call cwtime(cpu_all, wall_all, gflops_all, "start")
392 
393  do spin=1,nsppol
394    do ik=1,nkpt
395 
396      ! Only do a subset a k-points
397      if (ds%only_diago) then
398        if (all(distrib_diago(:, ik, spin) /= my_rank)) cycle
399      else
400        if (all(distrib_mat(bmin:bmax, bmin:bmax, ik, spin) /= my_rank)) cycle
401      end if
402      call cwtime(cpu, wall, gflops, "start")
403 
404      nband_k  = wfd%nband(ik, spin)
405      istwf_k  = wfd%istwfk(ik)
406      npw_k    = wfd%npwarr(ik)
407      kpt      = wfd%kibz(:,ik)
408      kg_k(:,1:npw_k) = wfd%kdata(ik)%kg_k
409 
410      if (dtset%useria /= 666) then
411        call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kpt, istwf_k, npw_k, kg_k)
412      else
413        ! Allocate KB form factors
414        ! Prepare term i <n,k|[Vnl,r]|n"k>
415        if (dtset%inclvkb /= 0) call vkbr_init(vkbr, cryst, psps, dtset%inclvkb, istwf_k, npw_k, kpt, kg_k)
416      end if
417 
418      ! Loop over bands
419      do ib_v=bmin,bmax
420        if (ds%only_diago) then
421          if (distrib_diago(ib_v,ik,spin) /= my_rank) cycle
422        else
423          if (all(distrib_mat(:,ib_v,ik,spin) /= my_rank)) cycle
424        end if
425 
426        if (dtset%useria /= 666) then
427          call wfd%copy_cg(ib_v, ik, spin, cg_v)
428          call ddkop%apply(ebands%eig(ib_v, ik, spin), npw_k, wfd%nspinor, cg_v, cwaveprj)
429        else
430          ABI_CHECK(wfd%get_wave_ptr(ib_v, ik, spin, wave_v, msg) == 0, msg)
431          ug_v(1:npw_k*nspinor) = wave_v%ug
432        end if
433 
434        ! Loop over bands
435        bstop = bmax; if (ds%only_diago) bstop = ib_v
436        do ib_c=ib_v,bstop
437          if (.not. ds%only_diago) then
438            if (distrib_mat(ib_c, ib_v, ik, spin) /= my_rank) cycle
439          end if
440 
441          if (dtset%useria /= 666) then
442            call wfd%copy_cg(ib_c, ik, spin, cg_c)
443            vv = ddkop%get_braket(ebands%eig(ib_c, ik, spin), istwf_k, npw_k, nspinor, cg_c, mode=ds%mode)
444            !if (ib_v == ib_c) vv(2, :) = zero
445 
446            if (ds%only_diago) then
447              ds%vdiago(:,ib_c,ik,spin) = vv(1, :)
448            else
449              ds%vmat(:,:,ib_c,ib_v,ik,spin) = vv
450              ! Hermitian conjugate
451              if (ib_v /= ib_c) then
452                ds%vmat(1,:,ib_v,ib_c,ik,spin) =  vv(1, :)
453                ds%vmat(2,:,ib_v,ib_c,ik,spin) = -vv(2, :)
454              end if
455            end if
456 
457            do idir=1,3
458              ds%dipoles(idir,:,ib_c,ib_v,ik,spin) = vv(:, idir)
459              ! Hermitian conjugate
460              if (ib_v /= ib_c) ds%dipoles(idir,:,ib_v,ib_c,ik,spin) = [vv(1, idir), -vv(2, idir)]
461            end do
462 
463          else
464            ABI_CHECK(wfd%get_wave_ptr(ib_c, ik, spin, wave_c, msg) == 0, msg)
465            ug_c(1:npw_k*nspinor) = wave_c%ug
466 
467            ! Calculate matrix elements of i[H,r] for NC pseudopotentials.
468            ihrc = nc_ihr_comm(vkbr, cryst, psps, npw_k, nspinor, istwf_k, dtset%inclvkb, kpt, ug_c, ug_v, kg_k)
469 
470            ! HM: 24/07/2018
471            ! Transform dipoles to be consistent with results from DFPT
472            ! Perturbations with DFPT are along the reciprocal lattice vectors
473            ! Perturbations with Commutator are along real space lattice vectors
474            ! dot(A, DFPT) = X
475            ! dot(B, COMM) = X
476            ! B = 2 pi (A^{-1})^T =>
477            ! dot(B^T B,COMM) = 2 pi DFPT
478            vr = (2*pi)*(2*pi)*sum(ihrc(:,:),dim=2)
479            vg(1) = dot_product(Cryst%gmet(1,:), vr)
480            vg(2) = dot_product(Cryst%gmet(2,:), vr)
481            vg(3) = dot_product(Cryst%gmet(3,:), vr)
482 
483            ! Save matrix elements of i*r in the IBZ
484            ds%dipoles(:,1,ib_c,ib_v,ik,spin) = real(vg, kind=dp)
485            ds%dipoles(:,1,ib_v,ib_c,ik,spin) = real(vg, kind=dp) ! Hermitian conjugate
486            if (ib_v == ib_c) then
487              ds%dipoles(:,2,ib_c,ib_v,ik,spin) = zero
488              ds%dipoles(:,2,ib_v,ib_c,ik,spin) = zero
489            else
490              ds%dipoles(:,2,ib_c,ib_v,ik,spin) =  aimag(vg)
491              ds%dipoles(:,2,ib_v,ib_c,ik,spin) = -aimag(vg) ! Hermitian conjugate
492            end if
493          end if
494 
495        end do
496      end do
497 
498      ! Free KB form factors
499      call vkbr_free(vkbr)
500 
501      if (nkpt < 1000 .or. (nkpt > 1000 .and. mod(ik, 200) == 0) .or. ik <= nproc) then
502        write(msg,'(2(a,i0),a)')" k-point [", ik, "/", nkpt, "]"
503        call cwtime_report(msg, cpu, wall, gflops)
504      end if
505 
506    end do ! k-points
507  end do ! spin
508 
509  call cwtime_report(msg, cpu_all, wall_all, gflops_all)
510 
511  ABI_FREE(ug_c)
512  ABI_FREE(ug_v)
513  ABI_FREE(kg_k)
514  ABI_FREE(ihrc)
515  ABI_FREE(cwaveprj)
516  ABI_SFREE(distrib_mat)
517  ABI_SFREE(distrib_diago)
518 
519  if (dtset%useria /= 666) then
520    ABI_FREE(cg_c)
521    ABI_FREE(cg_v)
522    call ddkop%free()
523  end if
524 
525  ! Gather the k-points computed by all processes
526  call xmpi_sum_master(ds%dipoles, master, comm, ierr)
527 
528  if (ds%only_diago) then
529    call xmpi_sum_master(ds%vdiago, master, comm, ierr)
530  else
531    call xmpi_sum_master(ds%vmat, master, comm, ierr)
532  end if
533 
534  ! Write matrix elements to disk.
535 
536  ! Output EVK file in netcdf format.
537  if (my_rank == master .and. write_ncfile) then
538    ! Have to build hdr on k-grid with info about perturbation.
539    call hdr_copy(hdr, tmp_hdr)
540    tmp_hdr%qptn = zero
541 
542    !fname = strcat(prefix, "NEW_EVK.nc")
543    !call wrtout(ab_out, sjoin("- Writing file: ", fname))
544    !NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
545    !tmp_hdr%pertcase = 0
546    !NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.))
547    !NCF_CHECK(cryst%ncwrite(ncid))
548    !NCF_CHECK(ebands_ncwrite(ebands, ncid))
549    !if (ds%only_diago) then
550    !  ncerr = nctk_def_arrays(ncid, [ &
551    !    nctkarr_t('vred_diagonal', "dp", "three, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
552    !else
553    !  ncerr = nctk_def_arrays(ncid, [ nctkarr_t('vred_matrix', "dp", &
554    !      "two, three, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
555    !end if
556    !NCF_CHECK(ncerr)
557    !NCF_CHECK(nctk_set_datamode(ncid))
558    !if (ds%only_diago) then
559    !  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_diagonal"), ds%vdiago))
560    !else
561    !  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_matrix"), ds%vmat))
562    !end if
563    !NCF_CHECK(nf90_close(ncid))
564 
565    do ii=1,3
566      fname = strcat(prefix, '_', itoa(ii), "_EVK.nc")
567      call wrtout(ab_out, sjoin("- Writing EVK file: ", fname, "for reduced direction:", itoa(ii)))
568      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
569      tmp_hdr%pertcase = 3 * cryst%natom + ii
570      NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.))
571      NCF_CHECK(cryst%ncwrite(ncid))
572      NCF_CHECK(ebands_ncwrite(ebands, ncid))
573      ncerr = nctk_def_arrays(ncid, [ &
574        nctkarr_t('h1_matrix_elements', "dp", &
575         "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
576      NCF_CHECK(ncerr)
577      NCF_CHECK(nctk_set_datamode(ncid))
578      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), ds%dipoles(ii,:,:,:,:,:)))
579      NCF_CHECK(nf90_close(ncid))
580    end do
581    call tmp_hdr%free()
582  end if
583 
584  if (my_rank == master .and. dtset%prtvol > 0) then
585    write(ab_out, "(2a)")ch10,"Writing velocity matrix elements (only diagonal terms, real part) for testing purpose:"
586    do spin=1,nsppol
587      do ik=1,min(nkpt, 4)
588        write(ab_out, "(2(a,1x,i0),2x,2a)")"For spin: ", spin, ", ikbz: ", ik, ", kpt: ", trim(ktoa(wfd%kibz(:,ik)))
589        do ib_c=bmin,min(bmin+8, bmax)
590          write(ab_out, "(3(es16.6,2x))") ds%dipoles(:,1,ib_c,ib_c,ik,spin)
591        end do
592        write(ab_out,*)""
593        !do ib_c=bmin,min(bmin+8, bmax)
594        !  write(ab_out, "(a, 6(es16.6,2x))")"Sum_k: ", sum(ds%dipoles(:,:,ib_c,ib_c,:,spin), dim=3) / nkpt
595        !end do
596      end do
597    end do
598  end if
599 
600  ! Free memory
601  call wfd%free()
602  call ebands_free(ebands)
603  call cryst%free()
604  call hdr%free()
605 
606  ! Block all procs here so that we know output files are available when code returns.
607  call xmpi_barrier(comm)
608 
609 end subroutine ddkstore_compute_ddk

m_ddk/ddkstore_free [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddkstore_free

FUNCTION

  Free memory

INPUTS

SOURCE

625 subroutine ddkstore_free(self)
626 
627 !Arguments ------------------------------------
628 !scalars
629  class(ddkstore_t),intent(inout) :: self
630 
631 !************************************************************************
632 
633  ABI_SFREE(self%vdiago)
634  ABI_SFREE(self%vmat)
635  ABI_SFREE(self%dipoles)
636 
637 end subroutine ddkstore_free

m_ddk/ddkstore_t [ Types ]

[ Top ] [ m_ddk ] [ Types ]

NAME

  ddkstore_t

FUNCTION

  This object stores the matrix elements of the velocity operator computed with the DFPT routines.

SOURCE

161  type,public :: ddkstore_t
162 
163    integer :: bmin = 1, bmax = -1
164    ! Min and max band index
165 
166    character(len=50) :: mode = "reduced"
167     ! "cart" or "reduced"
168 
169    logical :: only_diago = .False.
170    ! True if we are computing only the diagonal elements
171 
172    real(dp),allocatable :: dipoles(:,:,:,:,:,:)
173     ! (3, 2, mband, mband, nkpt, nsppol))
174 
175    real(dp),allocatable :: vdiago(:,:,:,:)
176     ! (3, bmin:bmax, nkpt, nsppol)
177 
178    real(dp),allocatable :: vmat(:,:,:,:,:,:)
179     ! (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol))
180 
181  contains
182 
183    procedure :: compute_ddk => ddkstore_compute_ddk
184     ! Calculate DDK matrix elements (diago or full b,b' matrix).
185     ! Return results in datatype. Optionally, save results to disk in EVK format.
186 
187    procedure :: free => ddkstore_free
188     ! Free memory.
189 
190  end type ddkstore_t

m_ddk/ham_targets_free [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

FUNCTION

INPUTS

SOURCE

1140 subroutine ham_targets_free(self)
1141 
1142 !Arguments ------------------------------------
1143 !scalars
1144  class(ham_targets_t),intent(inout) :: self
1145 
1146 !************************************************************************
1147 
1148  ABI_SFREE(self%ffnlk)
1149  ABI_SFREE(self%ffnl1)
1150  ABI_SFREE(self%kpg_k)
1151  ABI_SFREE(self%kpg1_k)
1152  ABI_SFREE(self%dkinpw)
1153  ABI_SFREE(self%kinpw1)
1154  ABI_SFREE(self%ph3d)
1155  ABI_SFREE(self%ph3d1)
1156 
1157 end subroutine ham_targets_free