hsmt.F90 10 KB
Newer Older
1 2 3 4 5
MODULE m_hsmt
  USE m_juDFT
  IMPLICIT NONE
CONTAINS
  SUBROUTINE hsmt(DIMENSION,atoms,sphhar,sym,enpara,SUB_COMM,n_size,n_rank,&
Daniel Wortmann's avatar
Daniel Wortmann committed
6
       jsp,input,mpi,lmaxb,noco,cell,&
7
       lapw, bkpt, vr,vs_mmp, oneD,usdus, kveclo,tlmplm,l_real,hamOvlp)
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
    !=====================================================================
    !     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
69 70 71
#ifdef CPP_GPU    
    USE m_hsmt_nonsph_gpu
#endif
72 73 74
    USE m_hsmt_sph
    USE m_hsmt_extra
    USE m_types
75 76
    USE m_hsmt_fjgj
    USE m_hsmt_overlap
77
    IMPLICIT NONE
78 79 80 81 82 83 84 85 86 87 88 89
    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
    TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp
90 91 92
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: SUB_COMM,n_size,n_rank 
93
    INTEGER, INTENT (IN) :: jsp  
94 95 96 97 98
    INTEGER, INTENT (IN) :: lmaxb

    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: bkpt(3) 
Daniel Wortmann's avatar
Daniel Wortmann committed
99
    REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
100 101 102 103
    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)
104 105

    LOGICAL,INTENT(IN)    :: l_real
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123

    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
    !
    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
Daniel Wortmann's avatar
Daniel Wortmann committed
184
       ALLOCATE ( fj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntype,2),gj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntype,2))
185
    ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
186 187
       ALLOCATE(fj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntype,ab_dim))
       ALLOCATE(gj(DIMENSION%nvd,0:atoms%lmaxd,atoms%ntype,ab_dim))
188 189 190 191
    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
       CALL timestart("hsmt spherical")
200
       CALL hsmt_sph(sym,DIMENSION,atoms,SUB_COMM,n_size,n_rank,sphhar,isp,ab_dim,&
201
            input,hlpmsize,noco,l_socfirst,cell,nintsp,lapw,enpara%el0,usdus,&
202
            vr,gk,rsoc,isigma,fj,gj,l_real,hamOvlp)
203
       CALL timestop("hsmt spherical")
204 205 206 207 208 209 210 211 212 213 214
#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)
215 216
       bb=0.0
       bb_test=0.0
217
       CALL timestart("hsmt_overlap")
218
       CALL  hsmt_overlap(input,atoms,n_size,n_rank,isp,l_socfirst,hlpmsize,ab_dim,&
219 220
       noco,cell,nintsp, lapw,usdus,gk,fj,gj,bb_test)
       CALL timestop("hsmt_overlap")
221 222 223 224
       CALL timestart("hsmt_overlap_zherk")
       CALL  hsmt_overlap_zherk(input,sym,atoms,n_size,n_rank,isp,l_socfirst,hlpmsize,ab_dim,&
       noco,cell,nintsp, lapw,usdus,gk,vk,fj,gj,bb)
       CALL timestop("hsmt_overlap_zherk")
225
       print *,"bb_diff",maxval(abs(bb-bb_test))
226
       
227 228
       deallocate(bb_test)
#endif
229 230
       IF (.NOT.input%secvar) THEN
          CALL timestart("hsmt extra")
Daniel Wortmann's avatar
Daniel Wortmann committed
231
          IF (ANY(atoms%nlo>0).OR.(atoms%n_u.GT.0)) &
232
               CALL hsmt_extra(DIMENSION,atoms,sym,isp,n_size,n_rank,input,nintsp,sub_comm,&
Daniel Wortmann's avatar
Daniel Wortmann committed
233
               hlpmsize,lmaxb,noco,l_socfirst,lapw,cell,enpara%el0,&
234
               fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in
235
               kveclo,l_real,hamOvlp%a_r,hamOvlp%b_r,hamOvlp%a_c,hamOvlp%b_c) !out/in
236 237
          CALL timestop("hsmt extra")
          CALL timestart("hsmt non-spherical")
238
#ifndef CPP_GPU
239
          CALL hsmt_nonsph(DIMENSION,atoms,sym,SUB_COMM,n_size,n_rank,input,isp,nintsp,&
240
               hlpmsize,noco,l_socfirst,lapw,cell,tlmplm,fj,gj,gk,vk,oneD,l_real,hamOvlp%a_r,hamOvlp%a_c)
241 242

          CALL timestop("hsmt non-spherical")
243 244 245 246 247 248 249
#else
          CALL timestart("hsmt non-spherical-GPU")
          CALL hsmt_nonsph_gpu(DIMENSION,atoms,sym,SUB_COMM,n_size,n_rank,input,isp,nintsp,&
               hlpmsize,noco,l_socfirst,lapw,cell,tlmplm,fj,gj,gk,vk,oneD,l_real,hamOvlp%a_r,hamOvlp%a_c)

          CALL timestop("hsmt non-spherical-GPU")
#endif
250 251
       ENDIF
    ENDDO
252 253 254 255
#if 1==2
       print *,"hsmt checkmode"
       call judft_end("test",0)
#endif
256 257 258
    RETURN
  END SUBROUTINE hsmt
END MODULE m_hsmt