hs_int.F90 3.58 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
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 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_hs_int
CONTAINS
  !Subroutine to construct the interstitial Hamiltonian and overlap matrix
  SUBROUTINE hs_int(input,noco,stars,lapw,mpi,cell,isp,vpw,&
       smat,hmat)
    USE m_types
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)      :: input
    TYPE(t_noco),INTENT(IN)       :: noco
    TYPE(t_stars),INTENT(IN)      :: stars
    TYPE(t_cell),INTENT(IN)       :: cell
    TYPE(t_lapw),INTENT(IN)       :: lapw
    TYPE(t_mpi),INTENT(IN)        :: mpi
    INTEGER,INTENT(IN)            :: isp
    COMPLEX,INTENT(IN)            :: vpw(:,:)
    CLASS(t_mat),INTENT(INOUT)     :: smat(:,:),hmat(:,:)


    INTEGER :: ispin,jspin,vpw_spin !spin indices
    INTEGER :: i,j,ii(3),iispin,jjspin,i0
    INTEGER :: in
    COMPLEX :: th,ts,phase
    REAL    :: b1(3),b2(3),r2

    IF (noco%l_noco.AND.isp==2) RETURN !was done already
    DO ispin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
       iispin=MIN(ispin,SIZE(smat,1))
       DO jspin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
          jjspin=MIN(jspin,SIZE(smat,1))
          IF (jspin==ispin) THEN
             vpw_spin=ispin
          ELSE
             vpw_spin=3
          ENDIF
          !$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
          !$OMP SHARED(mpi,lapw,stars,input,cell,vpw) &
          !$OMP SHARED(jjspin,iispin,ispin,jspin,vpw_spin)&
          !$OMP SHARED(hmat,smat)&
          !$OMP PRIVATE(ii,i0,i,j,in,phase,b1,b2,r2,th,ts)
          DO  i = mpi%n_rank+1,lapw%nv(ispin),mpi%n_size
             i0=(i-1)/mpi%n_size+1
             !--->    loop over (k+g)
             DO  j = 1,i  
                !-->     determine index and phase factor
                ii = lapw%gvec(:,i,ispin) - lapw%gvec(:,j,jspin)
                in = stars%ig(ii(1),ii(2),ii(3))
                IF (in.EQ.0) CYCLE
                phase = stars%rgphs(ii(1),ii(2),ii(3))
                !+APW_LO
                ts = phase*stars%ustep(in)
                IF (input%l_useapw) THEN
                   b1=lapw%bkpt+lapw%gvec(:,i,ispin)
                   b2=lapw%bkpt+lapw%gvec(:,j,jspin)
                   r2 = DOT_PRODUCT(MATMUL(b2,cell%bbmat),b1)   
                   th = phase*(0.5*r2*stars%ustep(in)+vpw(in,vpw_spin))
                ELSE
                   IF (vpw_spin==3.AND.jspin==2) THEN !The off-diagonal part is only due
                      th = vpw(in,vpw_spin)           !to potential
                      ts=0.0
                   ELSEIF(vpw_spin==3) THEN
                      th = CONJG(vpw(in,vpw_spin))
                      ts=0.0
                   ELSE
                      th = phase* (0.25* (lapw%rk(i,ispin)**2+lapw%rk(j,jspin)**2)*stars%ustep(in) + vpw(in,vpw_spin))
                   ENDIF
                ENDIF
                !-APW_LO
                !--->    determine matrix element and store
                IF (hmat(1,1)%l_real) THEN
                   hmat(jjspin,iispin)%data_r(j,i0) = REAL(th)
                   smat(jjspin,iispin)%data_r(j,i0) = REAL(ts)
                else
                   hmat(jjspin,iispin)%data_c(j,i0) = th
                   smat(jjspin,iispin)%data_c(j,i0) = ts
                endif
             ENDDO
          ENDDO
          !$OMP END PARALLEL DO
       ENDDO
    ENDDO
  END SUBROUTINE hs_int
END MODULE m_hs_int