hsmt.F90 9.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
MODULE m_hsmt
  USE m_juDFT
  IMPLICIT NONE
CONTAINS
  SUBROUTINE hsmt(DIMENSION,atoms,sphhar,sym,enpara,SUB_COMM,n_size,n_rank,&
       jsp,input,matsize,mpi,lmaxb,gwc,noco,cell,&
       lapw, bkpt, vr,vs_mmp, oneD,usdus, kveclo,aa,bb,tlmplm)
    !=====================================================================
    !     redesign of old hssphn into hsmt
    !           splitted into several parts:
    !                ->hsmt_init_soc
    !                ->hsmt_sph (spherical part)
    !                ->hsmt_extra (LOs and LDA+U)
    !                ->hsmt_nonsp (non-spherical part)
    !
    !                     Daniel Wortmann 2014
    !=====================================================================
    !*********************************************************************
    !     updates the hamiltonian and overlap matrices with the
    !     contributions from the spheres, both spherical and non-
    !     spherical.
    !                m. weinert  1986
    !     modified for inversion symmetry (real version), rp.p nov. 92
    !
    !     Since this non-collinear version allows only one definite
    !     direction of magnetic field in each muffin-tin, no matrix
    !     potential has to be used inside the MT's. However, the boundary
    !     conditions change, because the local (spin-) frame is different
    !     from the global frame. As a result, the usual spin-up and -down
    !     Hamitonian- and overlapp-matrix elements have to be
    !     premultiplied by a phasefactor, depending on the direction of
    !     the magentic field of each MT, to form the elements of the full
    !     complex matrices.
    !
    !     Philipp Kurz 98/01/27
    !*********************************************************************
    !------------------------------------------------------------------+
    ! Note for ev-parallelization:                                     |
    !                             unlike in hsint and hsvac, we do not |
    ! move through the H-matrix block by block (i.e. first the up/up,  |
    ! then the down/down, and finally the up/down spin blocks) but     |
    ! go  through one block only and update all blocks simultaniously. |
    ! Therefore, virtually we move through the whole matrix (up to     |
    ! nv(1)+nv(2)) and project back to the first block.                |
    !                                                   gb-00          |
    !------------------------------------------------------------------+
    !******** ABBREVIATIONS ***********************************************
    !     aa       : hamitonian matrix
    !     bb       : overlapp  matrix
    !     alph,beta: Euler angles of the local magnetic field direction of
    !                each atom (-type).
    !     chi      : Pauli spinors of the local spin-coordinate-frames
    !     chinn    : prefactors to be multiplied with the Hamiltonian- and
    !                overlappmatrix elements
    !**********************************************************************
    !
    !******** Spin-orbit interaction *************************************
    !   When l_soc=true & l_noco=true & l_ss=false, spin-orbit interaction
    !   is added in the first variation
    !   The formula by Youn et al. J.Comp.Phys. 172, 387 (2001) is used
    !   (note the wrong sign in the paper):
    !   <G|Hso|G'>= -i sigma . (G x G') sum_l (2l+1)/4pi Rso P'_l(G.G')
    !
    !   Jussi Enkovaara 2004, Juelich
    !********************************************************************
#include"cpp_double.h"
    USE m_hsmt_socinit
    USE m_hsmt_nonsph
    USE m_hsmt_sph
    USE m_hsmt_extra
    USE m_types
72 73
    USE m_hsmt_fjgj
    USE m_hsmt_overlap
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(IN)       :: mpi
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
    TYPE(t_oneD),INTENT(IN)      :: oneD
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_noco),INTENT(IN)      :: noco
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_enpara),INTENT(IN)    :: enpara
    TYPE(t_lapw),INTENT(INOUT)   :: lapw !nlotot,nv_tot&nmat can be updated
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: SUB_COMM,n_size,n_rank 
    INTEGER, INTENT (IN) :: jsp  ,matsize 
    INTEGER, INTENT (IN) :: lmaxb
    INTEGER, INTENT (IN) :: gwc

    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: bkpt(3) 
    REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,DIMENSION%jspd)
    COMPLEX,INTENT(IN):: vs_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u,input%jspins)
    TYPE(t_usdus),INTENT(INOUT)  :: usdus

    INTEGER, INTENT (OUT) :: kveclo(atoms%nlotot)
#ifdef CPP_INVERSION
    REAL,    INTENT (INOUT) :: aa(matsize),bb(matsize)
#else
    COMPLEX, INTENT (INOUT) :: aa(matsize),bb(matsize)
#endif

    TYPE(t_tlmplm),INTENT(INOUT) :: tlmplm


    !     ..
    !     .. Local Scalars ..
    INTEGER :: k,i,hlpmsize,isp,jsp_start,jsp_end,nintsp,iintsp,nc,ab_dim
    LOGICAL :: l_socfirst
    !     ..
    !     .. Local Arrays ..
    REAL                 :: v(3)
    COMPLEX              :: isigma(2,2,3)
    REAL, ALLOCATABLE    :: fj(:,:,:,:),gj(:,:,:,:)
    REAL, ALLOCATABLE    :: gk(:,:,:),vk(:,:,:)
    REAL, PARAMETER      :: eps = 1.0e-30

    TYPE(t_rsoc):: rsoc

124 125 126
#if 1==2
    REAL, ALLOCATABLE    :: fj_test(:,:,:,:),gj_test(:,:,:,:)
    COMPLEX,ALLOCATABLE  :: bb_test(:)
127
#endif
128

129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
    !
    CALL timestart("hsmt init")
    l_socfirst= noco%l_soc .AND. noco%l_noco .AND. (.NOT. noco%l_ss)
    IF (l_socfirst) CALL hsmt_socinit(mpi,atoms,sphhar,enpara,input,vr,noco,& !in
         rsoc,usdus,isigma) !out

#ifdef CPP_MPI
    !
    ! determine size of close-packed matrix
    !
    nc = 0
    k = 0
    DO  i=1+n_rank, lapw%nv(1)+atoms%nlotot, n_size
       nc = nc + 1
       k = k + n_size*(nc-1) + n_rank + 1
    ENDDO
    hlpmsize = k
#else
    hlpmsize = (DIMENSION%nvd+atoms%nlotot)*(DIMENSION%nvd+atoms%nlotot+1)/2
#endif


    !!---> pk non-collinear
    !---> loop over the local spins and interstitial spins
    nintsp = 1
    IF (noco%l_noco) THEN
       jsp_start = 1
       jsp_end   = 2
       !--->    in a spin-spiral calculations the a- and b-coeff and other
       !--->    quantities depend on the interstitial spin. therefore, an
       !--->    additional loop over the interstitial spin is needed in some
       !--->    places.
       IF (noco%l_ss) nintsp = 2
    ELSE
       jsp_start = jsp
       jsp_end   = jsp
    ENDIF
    !Set up the k+G+qss vectors
    ab_dim = 1
    IF (noco%l_ss) ab_dim = 2
    ALLOCATE(vk(DIMENSION%nvd,3,ab_dim),gk(DIMENSION%nvd,3,ab_dim))

    DO iintsp = 1,nintsp
       DO k = 1,lapw%nv(iintsp)
          IF (iintsp ==1) THEN
             v=bkpt+(/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)-noco%qss/2
          ELSE
             v=bkpt+(/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)+noco%qss/2
          ENDIF
          vk(k,:,iintsp) = v
          gk(k,:,iintsp) = MATMUL(TRANSPOSE(cell%bmat),v)/MAX (lapw%rk(k,iintsp),eps)
       END DO
    ENDDO
    !! the fj and gj arrays are constructed in hssphn_sph and used later
    IF (noco%l_constr.OR.l_socfirst) THEN
       ALLOCATE ( fj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntypd,2),gj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntypd,2))
    ELSE
       ALLOCATE(fj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntypd,ab_dim))
       ALLOCATE(gj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntypd,ab_dim))
    ENDIF
    CALL timestop("hsmt init")

    DO isp = jsp_start,jsp_end
192 193 194 195 196 197 198
#if 1==2
       ALLOCATE(bb_test(size(bb)))
       bb_test=0.0
       bb_test=bb
       fj=0.0
       gj=0.0
#endif
199 200 201 202 203
       CALL timestart("hsmt spherical")
       CALL hsmt_sph(DIMENSION,atoms,SUB_COMM,n_size,n_rank,sphhar,isp,ab_dim,&
            input,hlpmsize,noco,l_socfirst,cell,nintsp,lapw,enpara%el0,usdus,&
            vr,gk,rsoc,isigma, aa,bb,fj,gj)
       CALL timestop("hsmt spherical")
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
#if 1==2
       ALLOCATE(fj_test(size(fj,1),0:size(fj,2)-1,size(fj,3),size(fj,4)))
       ALLOCATE(gj_test(size(fj,1),0:size(fj,2)-1,size(fj,3),size(fj,4)))
       fj_test=0.0
       gj_test=0.0
       CALL timestart("hsmt_fjgj")
       CALL hsmt_fjgj(atoms,isp,noco,l_socfirst,cell,nintsp, lapw,usdus,fj_test,gj_test)
       CALL timestop("hsmt_fjgj")
       print *,"fj diff",maxval(abs(fj-fj_test))
       print *,"gj diff",maxval(abs(gj-gj_test))
       deallocate(fj_test,gj_test)
       CALL timestart("hsmt_overlap")
        CALL  hsmt_overlap(input,atoms,n_size,n_rank,isp,l_socfirst,hlpmsize,ab_dim,&
       noco,cell,nintsp, lapw,usdus,gk,fj,gj,bb_test)
       CALL timestop("hsmt_overlap")
       print *,"bb_diff",maxval(abs(bb-bb_test))
       deallocate(bb_test)
#endif
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
       IF (.NOT.input%secvar) THEN
          CALL timestart("hsmt extra")
          IF (ANY(atoms%nlo>0).OR.ANY(atoms%lda_u%l.GE.0)) &
               CALL hsmt_extra(DIMENSION,atoms,sym,isp,n_size,n_rank,input,nintsp,sub_comm,&
               hlpmsize,lmaxb,gwc,noco,l_socfirst,lapw,cell,enpara%el0,&
               fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in
               kveclo,aa,bb) !out/in
          CALL timestop("hsmt extra")
          CALL timestart("hsmt non-spherical")
          CALL hsmt_nonsph(DIMENSION,atoms,sym,SUB_COMM,n_size,n_rank,input,isp,nintsp,&
               hlpmsize,noco,l_socfirst,lapw,cell,tlmplm,fj,gj,gk,vk,oneD,aa)

          CALL timestop("hsmt non-spherical")
       ENDIF
    ENDDO
237 238 239 240
#if 1==2
       print *,"hsmt checkmode"
       call judft_end("test",0)
#endif
241 242 243
    RETURN
  END SUBROUTINE hsmt
END MODULE m_hsmt