LCOV - code coverage report
Current view: top level - eigen - slomat.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 58 68 85.3 %
Date: 2024-03-28 04:22:06 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_slomat
       8             :   !***********************************************************************
       9             :   ! updates the overlap matrix with the contributions from the local
      10             :   ! orbitals.
      11             :   !                                                p.kurz sept. 1996
      12             :   !***********************************************************************
      13             : CONTAINS
      14       11286 :    SUBROUTINE slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,ntyp,na,&
      15       11286 :                      isp,ud, alo1,blo1,clo1,fjgj,&
      16             :                      igSpinPr,igSpin,chi,smat,l_fullj,lapwq,fjgjq)
      17             :     !***********************************************************************
      18             :     ! locol stores the number of columns already processed; on parallel
      19             :     !       computers this decides, whether the LO-contribution is
      20             :     !       done on this node                                          gb00
      21             :     !
      22             :     ! function legpol() at end of module
      23             :     !***********************************************************************
      24             : 
      25             :       USE m_constants,ONLY: fpi_const
      26             :       !USE m_types_mpimat
      27             :       USE m_types
      28             :       USE m_hsmt_fjgj
      29             : 
      30             :       IMPLICIT NONE
      31             : 
      32             :       TYPE(t_input),INTENT(IN)  :: input
      33             :       TYPE(t_atoms),INTENT(IN)  :: atoms
      34             :       TYPE(t_sym),INTENT(IN)    :: sym
      35             :       TYPE(t_lapw),INTENT(IN),Target   :: lapw
      36             :       TYPE(t_mpi),INTENT(IN)    :: fmpi
      37             :       TYPE(t_cell),INTENT(IN)   :: cell
      38             :       TYPE(t_nococonv),INTENT(IN)   :: nococonv
      39             :       TYPE(t_fjgj),INTENT(IN),target   :: fjgj
      40             : 
      41             :       ! Scalar Arguments
      42             :       INTEGER, INTENT (IN)      :: na,ntyp
      43             :       INTEGER, INTENT (IN)      :: igSpin,igSpinPr
      44             :       COMPLEX, INTENT (IN)      :: chi
      45             :       INTEGER, INTENT(IN)       :: isp
      46             : 
      47             :       ! Array Arguments
      48             :       REAL,   INTENT (IN)       :: alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
      49             :       TYPE(t_usdus),INTENT(IN)  :: ud
      50             :       CLASS(t_mat),INTENT(INOUT) :: smat
      51             :       LOGICAL, INTENT(IN) :: l_fullj
      52             : 
      53             :       TYPE(t_lapw), OPTIONAL, TARGET, INTENT(IN) :: lapwq
      54             :       TYPE(t_fjgj), OPTIONAL, Target, INTENT(IN) :: fjgjq
      55             : 
      56             :       ! Local Scalars
      57             :       REAL    :: con,dotp,fact1,fact2,fact3,fl2p1
      58             :       INTEGER :: invsfct,k ,l,lo,lop,lp,nkvec,nkvecp,kp,i
      59             :       INTEGER :: locol,lorow
      60             :       LOGICAL :: l_samelapw
      61             : 
      62       11286 :       COMPLEX,   ALLOCATABLE  :: cphPr(:), cph(:)
      63             : 
      64             :       TYPE(t_lapw),pointer :: lapwPr
      65             :       TYPE(t_fjgj),pointer :: fjgjPr
      66             : 
      67       11286 :       l_samelapw = .FALSE.
      68       11286 :       IF (.NOT.PRESENT(lapwq)) l_samelapw = .TRUE.
      69             :       IF (.NOT.l_samelapw) THEN
      70             :          lapwPr => lapwq
      71             :          fjgjPr => fjgjq
      72             :       ELSE
      73             :          lapwPr => lapw
      74             :          fjgjPr => fjgj
      75             :       END IF
      76             :       
      77       56430 :       ALLOCATE(cph(MAXVAL(lapw%nv)))
      78       56430 :       ALLOCATE(cphPr(MAXVAL(lapwPr%nv)))
      79             : 
      80             :       ! TODO: Introduce the logic for different lapw and the full rectangular
      81             :       !       instead of triangular construction...
      82             : 
      83       11286 :       CALL lapwPr%phase_factors(igSpinPr,atoms%taual(:,na),nococonv%qss,cphPr)
      84       11286 :       IF (l_samelapw.and.igspin==igSpinPr) THEN
      85     1444464 :          cph = cphPr
      86             :       ELSE
      87           0 :          CALL lapw%phase_factors(igSpin,atoms%taual(:,na),nococonv%qss,cph)
      88             :       END IF
      89             : 
      90       11286 :       IF ((sym%invsat(na) == 0) .OR. (sym%invsat(na) == 1)) THEN
      91             :          ! If this atom is the first of two atoms related by inversion, the
      92             :          ! contributions to the overlap matrix of both atoms are added simultaneously.
      93             :          ! Where it is made use of the fact, that the sum of these contributions
      94             :          ! is twice the real part of the contribution of each atom. Note, that in
      95             :          ! this case there are twice as many (2*(2*l+1)) k-vectors.
      96             :          ! (compare abccoflo and comments there).
      97             :          IF (sym%invsat(na) == 0) invsfct = 1
      98       11286 :          IF (sym%invsat(na) == 1) invsfct = 2
      99             : 
     100       11286 :          con = fpi_const/SQRT(cell%omtil)* ((atoms%rmt(ntyp))**2)/2.0
     101             : 
     102             :          !$acc kernels present(smat,smat%data_c,smat%data_r)&
     103             :          !$acc & copyin(fjgj,fjgj%fj,fjgj%gj,fjgjPr,fjgjPr%fj,fjgjPr%gj,l,lapw,lapw%kvec(:,:,na),lapwPr,lapwPr%kvec(:,:,na),ud,clo1(:),dotp,cph(:),cphPr(:),atoms,lapw%index_lo(:,na),lapw%gk(:,:,:)) &
     104             :          !$acc & copyin(lapwPr%index_lo(:,na),lapwPr%gk(:,:,:),ud%dulon(:,ntyp,isp),ud%ddn(:,ntyp,isp),ud%uulon(:,ntyp,isp),ud%uloulopn(:,:,ntyp,isp),blo1(:)) &
     105             :          !$acc & copyin(atoms,atoms%nlo(ntyp),lapw%nv(:),lapwPr%nv(:),atoms%llo(:,ntyp),alo1(:), fmpi, fmpi%n_size, fmpi%n_rank)&
     106             :          !$acc & default(none)
     107             :          
     108             :          !$acc loop gang private(fact1,fl2p1,l) independent
     109       28724 :          DO lo = 1,atoms%nlo(ntyp) !loop over all LOs for this atom
     110       17438 :             l = atoms%llo(lo,ntyp)
     111       17438 :             fl2p1 = (2*l+1)/fpi_const
     112             :             fact1 = (con**2) * fl2p1 * ( &
     113             :                   & alo1(lo)*(alo1(lo) &
     114             :                         & + 2*clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
     115             :                   & blo1(lo)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
     116             :                         & + 2*clo1(lo)*ud%dulon(lo,ntyp,isp)) + &
     117       17438 :                   & clo1(lo)*clo1(lo) )
     118             :             
     119       70144 :             DO nkvec = 1,invsfct* (2*l+1) !Each LO can have several functions
     120             :                ! lop-kG part (l_fullj). |====|
     121       41420 :                IF (l_fullj) THEN
     122           0 :                   lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvec
     123           0 :                   kp = lapwPr%kvec(nkvec,lo,na)
     124             :                   !$acc loop vector private(fact2,dotp) independent
     125           0 :                   DO k = fmpi%n_rank + 1, lapw%nv(igSpin), fmpi%n_size
     126             :                      fact2 = con * fl2p1 * ( &
     127             :                            & fjgj%fj(k,l,isp,igSpin)*(alo1(lo) &
     128             :                                                   & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
     129             :                            & fjgj%gj(k,l,isp,igSpin)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
     130           0 :                                                   & + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
     131           0 :                      dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
     132             : 
     133           0 :                      IF (smat%l_real) THEN
     134             :                         smat%data_r(lorow,(k-1) / fmpi%n_size + 1) = smat%data_r(lorow,(k-1) / fmpi%n_size + 1) &
     135             :                                             & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
     136           0 :                                             & * CONJG(cphPr(kp)) * cph(k)
     137             :                      ELSE
     138             :                         smat%data_c(lorow,(k-1) / fmpi%n_size + 1) = smat%data_c(lorow,(k-1) / fmpi%n_size + 1) &
     139             :                                             & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
     140           0 :                                             & * CONJG(cphPr(kp)) * cph(k)
     141             :                      END IF
     142             :                   END DO
     143             :                   !$acc end loop
     144             :                END IF
     145       41420 :                locol = lapw%nv(igSpin)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
     146       58858 :                IF (MOD(locol-1,fmpi%n_size) == fmpi%n_rank) THEN
     147       26146 :                   locol=(locol-1)/fmpi%n_size+1 !this is the column in local storage!
     148       26146 :                   k = lapw%kvec(nkvec,lo,na)
     149             :                   ! Calculate the overlap matrix elements with the regular Flapw
     150             :                   ! basis functions
     151             :                   !$acc loop vector private(fact2,dotp,kp) independent
     152     3109227 :                   DO kp = 1,lapwPr%nv(igSpinPr)
     153             :                      fact2 = con * fl2p1 * ( &
     154             :                            & fjgjPr%fj(kp,l,isp,igSpinPr)*(alo1(lo) &
     155             :                                                      & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
     156             :                            & fjgjPr%gj(kp,l,isp,igSpinPr)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
     157     3083081 :                                                      & + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
     158    12332324 :                      dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
     159             : 
     160     3109227 :                      IF (smat%l_real) THEN
     161             :                         smat%data_r(kp,locol) = smat%data_r(kp,locol) &
     162             :                                             & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
     163     1122558 :                                             & * CONJG(cphPr(kp)) * cph(k)
     164             :                      ELSE
     165             :                         smat%data_c(kp,locol) = smat%data_c(kp,locol) &
     166             :                                             & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
     167     1960523 :                                             & * CONJG(cphPr(kp)) * cph(k)
     168             :                      END IF
     169             :                   END DO
     170             :                   !$acc end loop
     171             :                   ! Calculate the overlap matrix elements with other local orbitals
     172             :                   ! of the same atom, if they have the same l
     173             :                   !$acc loop vector private(lp,fact3,nkvecp,kp,lorow,dotp) independent
     174       39764 :                   DO lop = 1, MERGE(lo-1,atoms%nlo(ntyp),igSpinPr==igSpin.AND..NOT.l_fullj)
     175       13618 :                      IF (lop==lo) CYCLE !Do later
     176       13618 :                      lp = atoms%llo(lop,ntyp)
     177       39764 :                      IF (l == lp) THEN
     178             :                         fact3 = con**2 * fl2p1 * ( &
     179             :                               & alo1(lop)*(alo1(lo) &
     180             :                                        & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
     181             :                               & blo1(lop)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
     182             :                                        & + clo1(lo)*ud%dulon(lo,ntyp,isp)) + &
     183             :                               & clo1(lop)*(alo1(lo)*ud%uulon(lop,ntyp,isp) &
     184             :                                        & + blo1(lo)*ud%dulon(lop,ntyp,isp) &
     185          60 :                                        & + clo1(lo)*ud%uloulopn(lop,lo,ntyp,isp)) )
     186         240 :                         DO nkvecp = 1,invsfct* (2*lp+1)
     187         180 :                            kp = lapwPr%kvec(nkvecp,lop,na)
     188         180 :                            lorow=lapwPr%nv(igSpinPr)+lapwPr%index_lo(lop,na)+nkvecp
     189         720 :                            dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
     190             : 
     191         240 :                            IF (smat%l_real) THEN
     192             :                               smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
     193             :                                                      & + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
     194         180 :                                                      & * CONJG(cphPr(kp)) * cph(k)
     195             :                            ELSE
     196             :                               smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
     197             :                                                      & + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
     198           0 :                                                      & * CONJG(cphPr(kp)) * cph(k)
     199             :                            END IF
     200             :                         END DO
     201             :                      END IF
     202             :                   END DO
     203             :                   ! Calculate the overlap matrix elements of one local
     204             :                   ! orbital with itself
     205       26146 :                   lop=lo
     206             :                   !$acc loop vector private(kp,lorow,dotp) independent
     207       75926 :                   DO nkvecp = 1,MERGE(nkvec,invsfct* (2*l+1),igSpinPr==igSpin.AND..NOT.l_fullj)
     208       49780 :                      kp = lapwPr%kvec(nkvecp,lo,na)
     209       49780 :                      lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvecp
     210      199120 :                      dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
     211             : 
     212       75926 :                      IF (smat%l_real) THEN
     213             :                         smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
     214             :                                                & + chi * invsfct * fact1 * legpol(l,dotp) &
     215       18047 :                                                & * CONJG(cphPr(kp)) * cph(k)
     216             :                      ELSE
     217             :                         smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
     218             :                                                & + chi * invsfct * fact1 * legpol(l,dotp) &
     219       31733 :                                                & * CONJG(cphPr(kp)) * cph(k)
     220             :                      END IF
     221             :                   END DO
     222             :                   !$acc end loop
     223             :                END IF ! mod(locol-1,n_size) = nrank
     224             :             END DO
     225             :             !$acc end loop
     226             :          END DO
     227             :          !$acc end loop
     228             :          !$acc end kernels
     229             :       END IF
     230             :     
     231       11286 :    END SUBROUTINE slomat
     232             : 
     233     3133041 :    PURE REAL FUNCTION legpol(l,arg)
     234             :       !$acc routine seq
     235             : 
     236             :       IMPLICIT NONE
     237             : 
     238             :       REAL,INTENT(IN)   :: arg
     239             :       INTEGER,INTENT(IN):: l
     240             : 
     241             :       INTEGER :: lp
     242             : 
     243     3133041 :       REAL :: plegend(0:l)
     244             : 
     245     3133041 :       plegend(0) = 1.0
     246     3133041 :       IF (l.GE.1) THEN
     247     2654253 :          plegend(1) = arg
     248     2870268 :          DO lp = 1,l - 1
     249     2870268 :             plegend(lp+1) = (lp+lp+1)*arg*plegend(lp)/ (lp+1) -lp*plegend(lp-1)/ (lp+1)
     250             :          END DO
     251             :       END IF
     252     3133041 :       legpol = plegend(l)
     253     3133041 :    END FUNCTION legpol
     254             : 
     255             : END MODULE m_slomat

Generated by: LCOV version 1.14