LCOV - code coverage report
Current view: top level - eigen_soc - sorad.f90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 98 123 79.7 %
Date: 2024-03-29 04:21:46 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             : MODULE m_sorad
       7             :   USE m_juDFT
       8             :   !*********************************************************************
       9             :   !     generates radial spin-orbit matrix elements
      10             :   !*********************************************************************
      11             : CONTAINS
      12         204 :   SUBROUTINE sorad(atoms,input,ntyp,vr,enpara,spav,rsoc,usdus,hub1data)
      13             : 
      14             :     USE m_constants, ONLY : c_light
      15             :     USE m_intgr,     ONLY : intgr0
      16             :     USE m_sointg
      17             :     USE m_radsra
      18             :     USE m_radsrd
      19             :     USE m_radsrdn
      20             :     USE m_types
      21             :     IMPLICIT NONE
      22             :     TYPE(t_enpara),INTENT(IN)   :: enpara
      23             :     TYPE(t_input),INTENT(IN)    :: input
      24             :     TYPE(t_atoms),INTENT(IN)    :: atoms
      25             :     TYPE(t_usdus),INTENT(INOUT) :: usdus
      26             :     TYPE(t_rsoc),INTENT(INOUT)  :: rsoc
      27             :     TYPE(t_hub1data),OPTIONAL,INTENT(IN)  :: hub1data
      28             :     !     ..
      29             :     !     .. Scalar Arguments ..
      30             :     INTEGER, INTENT (IN) :: ntyp
      31             :     LOGICAL, INTENT (IN) :: spav ! if T, spin-averaged pot is used
      32             :     !     ..
      33             :     !     .. Array Arguments ..
      34             :     REAL,    INTENT (IN) :: vr(:,:)!(atoms%jmtd,input%jspins),
      35             :     !     ..
      36             :     !     .. Local Scalars ..
      37             :     REAL ddn1,e ,ulops,dulops,duds1
      38             :     INTEGER i,j,ir,jspin,l,noded,nodeu,ilo,ilop
      39             :     LOGICAL l_hia
      40             :     !     ..
      41             :     !     .. Local Arrays ..
      42         204 :     REAL, ALLOCATABLE :: p(:,:),pd(:,:),q(:,:),qd(:,:),plo(:,:)
      43         204 :     REAL, ALLOCATABLE :: plop(:,:),glo(:,:),fint(:),pqlo(:,:)
      44         204 :     REAL, ALLOCATABLE :: filo(:,:)
      45         204 :     REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:),vrTmp(:)
      46             :     !     ..
      47             :     
      48         204 :     IF (atoms%jri(ntyp)>atoms%jmtd)  CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
      49             :     ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
      50        2856 :          &   qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2),vrTmp(atoms%jmtd) )
      51             : 
      52     1322712 :     p = 0.0 ; pd = 0.0 ; q = 0.0 ; plo = 0.0 ; fint = 0.0
      53     1469476 :     qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0; vrTmp = 0.0
      54             :     
      55             : 
      56        2144 :     DO l = 0,atoms%lmax(ntyp) 
      57        1940 :        l_hia=.FALSE.
      58        1940 :        DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
      59        1940 :           IF(atoms%lda_u(i)%atomType.EQ.ntyp.AND.atoms%lda_u(i)%l.EQ.l) THEN
      60           0 :              l_hia=.TRUE.
      61             :           ENDIF
      62             :        ENDDO
      63        5374 :        DO jspin = 1,input%jspins
      64             :           !TODO: here genMTBasis should be used
      65     2432390 :          vrTmp = vr(:,jspin)
      66      221738 :          if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
      67        3434 :           IF(l_hia.AND.input%jspins.EQ.2) THEN
      68           0 :              IF(PRESENT(hub1data)) THEN
      69           0 :                 IF(hub1data%l_performSpinavg) vrTmp = (vr(:,1)+vr(:,2))/2.0
      70             :              ENDIF
      71             :           ENDIF
      72             :           !
      73             :           !--->    calculate normalized function at e: p and q 
      74             :           !
      75        3434 :           e = enpara%el0(l,ntyp,jspin)
      76             :           CALL radsra(&
      77             :                e,l,vrTmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
      78             :                usdus%us(l,ntyp,jspin),usdus%dus(l,ntyp,jspin),&
      79        3434 :                nodeu,p(:,jspin),q(:,jspin))
      80             :           !                     
      81             :           !--->    calculate orthogonal energy derivative at e : pd and qd
      82             :           !
      83             :           CALL radsrd(&
      84             :                e,l,vrTmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
      85             :                usdus%uds(l,ntyp,jspin),usdus%duds(l,ntyp,jspin),&
      86             :                usdus%ddn(l,ntyp,jspin),noded,pd(:,jspin),qd(:,jspin),&
      87        5374 :                p(:,jspin),q(:,jspin),usdus%dus(l,ntyp,jspin))
      88             : 
      89             :        END DO     ! end of spin loop
      90             :        !
      91             :        !---> in case of jspins=1
      92             :        !
      93             : 
      94        1940 :        IF (input%jspins.EQ.1) THEN
      95      382036 :           DO i = 1,atoms%jri(ntyp)
      96      381590 :              p(i,2) =  p(i,1)
      97      382036 :              pd(i,2) = pd(i,1)
      98             :           ENDDO
      99             :        ENDIF
     100             :        !
     101             :        !---> common spin-orbit integrant V   (average spin directions)
     102             :        !                                  SO
     103     1416360 :        v0(:) = 0.0
     104        1940 :        IF (input%jspins.EQ.1) THEN
     105      382036 :           v0(1:atoms%jri(ntyp)) = vr(1:atoms%jri(ntyp),1)
     106         446 :           e = enpara%el0(l,ntyp,1)
     107             :        ELSE
     108     1012596 :           DO i = 1,atoms%jri(ntyp)
     109     1012596 :              v0(i) = (vr(i,1)+vr(i,input%jspins))/2.
     110             :           END DO
     111        1494 :           e = (enpara%el0(l,ntyp,1)+enpara%el0(l,ntyp,input%jspins))/2.
     112             :        END IF
     113             : 
     114     2682112 :        CALL sointg(ntyp,e,vr,v0,atoms,input,vso)
     115        1940 :        IF (spav) THEN
     116      109152 :           DO i= 1,atoms%jmtd
     117      109008 :              vso(i,1)= (vso(i,1)+vso(i,2))/2.
     118      109152 :              vso(i,2)= vso(i,1)
     119             :           ENDDO
     120             :        ENDIF
     121             : 
     122             :        !                        s       s'            .s       s'
     123             :        !-->  radial integrals <u  |V  |u  > = rsopp, <u  |V  |u  > = rsopdp etc.
     124             :        !                            SO                     SO
     125             : 
     126        1940 :        IF (l.GT.0) THEN ! there is no spin-orbit for s-states
     127        5208 :           DO i = 1, 2
     128       12152 :              DO j = 1, 2
     129     5009248 :                 rsoc%rsopp(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     130     5009248 :                 rsoc%rsopdp(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     131     5009248 :                 rsoc%rsoppd(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     132     5012720 :                 rsoc%rsopdpd(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     133             :              ENDDO
     134             :           ENDDO
     135             :        ENDIF ! l>0
     136             :        !
     137             :        !--->  Check for local orbitals with same l
     138             :        !
     139        4044 :        DO ilo = 1, atoms%nlo(ntyp)
     140        3840 :           IF (atoms%llo(ilo,ntyp).EQ.l) THEN
     141             : 
     142         550 :              DO jspin = 1,input%jspins
     143      281150 :                vrTmp = vr(:,jspin)
     144       48866 :                if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
     145         354 :                 e = enpara%ello0(ilo,ntyp,jspin)
     146             :                 CALL radsra(&
     147             :                      e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     148             :                      usdus%ulos(ilo,ntyp,jspin),usdus%dulos(ilo,ntyp,jspin),&
     149         354 :                      nodeu,plo(:,jspin),qlo(:,jspin))
     150             : 
     151             :                 !+apw+lo
     152         550 :                 IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN !  calculate energy derivative (of order atoms%ulo_der) at e
     153           0 :                    ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
     154           0 :                    glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
     155           0 :                    pqlo(1:atoms%jri(ntyp),1)=plo(1:atoms%jri(ntyp),jspin)
     156           0 :                    pqlo(1:atoms%jri(ntyp),2)=qlo(1:atoms%jri(ntyp),jspin)
     157           0 :                    i = atoms%ulo_der(ilo,ntyp)
     158           0 :                    IF(atoms%l_dulo(ilo,ntyp)) i=1
     159             :                    CALL radsrdn(&
     160             :                         e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     161             :                         usdus%ulos(ilo,ntyp,jspin),duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
     162           0 :                         pqlo,usdus%dulos(ilo,ntyp,jspin),i)
     163           0 :                    ddn1 = SQRT(ddn1)
     164           0 :                    IF(atoms%l_dulo(ilo,ntyp)) ddn1=1.0
     165           0 :                    plo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)/ddn1
     166           0 :                    qlo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),2)/ddn1
     167           0 :                    usdus%dulos(ilo,ntyp,jspin) = duds1/ddn1
     168           0 :                    DEALLOCATE (glo,pqlo,filo)
     169             :                 ENDIF
     170             :                 !-apw+lo
     171             :              ENDDO
     172             : 
     173         196 :              IF (input%jspins.EQ.1) THEN
     174       33860 :                 plo(1:atoms%jri(ntyp),2) = plo(1:atoms%jri(ntyp),1)
     175          38 :                 e = (enpara%ello0(ilo,ntyp,1) + enpara%el0(l,ntyp,1) )/2
     176             :              ELSE
     177             :                 e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +&
     178         158 :                      enpara%el0(l,ntyp,1) + enpara%el0(l,ntyp,input%jspins) )/4
     179             :              END IF
     180      335488 :              CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
     181         196 :              IF (spav) THEN
     182       24256 :                 DO i= 1,atoms%jmtd
     183       24224 :                    vso(i,1)= (vso(i,1)+vso(i,2))/2.
     184       24256 :                    vso(i,2)= vso(i,1)
     185             :                 ENDDO
     186             :              ENDIF
     187             : 
     188         588 :              DO i = 1, 2
     189        1372 :                 DO j = 1, 2
     190         784 :                    rsoc%rsoplop (ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),p (:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     191         784 :                    rsoc%rsoplopd(ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
     192         784 :                    rsoc%rsopplo (ntyp,ilo,i,j) = radso(p (:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
     193        1176 :                    rsoc%rsopdplo(ntyp,ilo,i,j) = radso(pd(:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
     194             :                 ENDDO
     195             :              ENDDO
     196             : 
     197         550 :              DO i = 1,input%jspins
     198      280796 :                 fint(:) = plo(:,i) *  p(:,i) + qlo(:,i) *  q(:,i)
     199         354 :                 CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,i))
     200      280796 :                 fint(:) = plo(:,i) * pd(:,i) + qlo(:,i) * qd(:,i)
     201         550 :                 CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,i))
     202             :              ENDDO
     203             : 
     204         548 :              DO ilop = 1, atoms%nlo(ntyp)
     205         548 :                 IF (atoms%llo(ilop,ntyp).EQ.l) THEN
     206             : 
     207         550 :                    DO jspin = 1,input%jspins
     208      281150 :                      vrTmp = vr(:,jspin)
     209       48866 :                      if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
     210         354 :                       e = enpara%ello0(ilop,ntyp,jspin)
     211             :                       CALL radsra(&
     212             :                            e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     213         354 :                            ulops,dulops,nodeu,plop(:,jspin),q(:,1))
     214             :                       !+apw+lo
     215         550 :                       IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate orthogonal energy derivative at e
     216           0 :                          ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
     217           0 :                          glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
     218           0 :                          pqlo(1:atoms%jri(ntyp),1)=plop(1:atoms%jri(ntyp),jspin)
     219           0 :                          pqlo(1:atoms%jri(ntyp),2)=q(1:atoms%jri(ntyp),1)
     220           0 :                          i = atoms%ulo_der(ilo,ntyp)
     221           0 :                          IF(atoms%l_dulo(ilo,ntyp)) i=1
     222             :                          CALL radsrdn(&
     223             :                               e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     224             :                               ulops,duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
     225           0 :                               pqlo,dulops,i)
     226           0 :                          plop(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)
     227           0 :                          DEALLOCATE (glo,pqlo,filo)
     228             :                       ENDIF
     229             :                       !-apw+lo
     230             :                    ENDDO
     231             : 
     232         196 :                    IF (input%jspins.EQ.1) THEN
     233       33860 :                       plop(1:atoms%jri(ntyp),2) = plop(1:atoms%jri(ntyp),1)
     234          38 :                       e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilop,ntyp,1) )/2
     235             :                    ELSE
     236             :                       e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +  &
     237         158 :                            enpara%ello0(ilop,ntyp,1) + enpara%ello0(ilop,ntyp,input%jspins) )/4
     238             :                    END IF
     239      335488 :                    CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
     240         196 :                    IF (spav) THEN
     241       24256 :                       DO i= 1,atoms%jmtd
     242       24224 :                          vso(i,1)= (vso(i,1)+vso(i,2))/2.
     243       24256 :                          vso(i,2)= vso(i,1)
     244             :                       ENDDO
     245             :                    ENDIF
     246             : 
     247         588 :                    DO i = 1, 2
     248        1372 :                       DO j = 1, 2
     249             :                          rsoc%rsoploplop(ntyp,ilo,ilop,i,j) =&
     250        1176 :                               radso(plo(:atoms%jri(ntyp),i),plop(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     251             :                       ENDDO
     252             :                    ENDDO
     253             : 
     254             :                 ENDIF
     255             :              ENDDO
     256             : 
     257             :           ENDIF
     258             :        ENDDO ! end of lo-loop
     259             : 
     260             :     ENDDO ! end of l-loop
     261             : 
     262         204 :     DEALLOCATE ( p,pd,q,qd,plo,plop,qlo,fint,v0,vso )
     263             :     !      rsoplop (:,:,:,:) = 0.0
     264             :     !      rsoplopd(:,:,:,:) = 0.0
     265             :     !      rsopplo (:,:,:,:) = 0.0
     266             :     !      rsopdplo(:,:,:,:) = 0.0
     267             :     !      rsoploplop(:,:,:,:,:) = 0.0
     268             : 
     269         204 :   END SUBROUTINE sorad
     270             :   !--------------------------------------------------------------------
     271       31696 :   REAL FUNCTION radso(a,b,vso,dx,r0)
     272             :     !
     273             :     !     compute radial spin-orbit integrals
     274             :     !
     275             :     USE m_intgr, ONLY : intgr0
     276             :     IMPLICIT NONE
     277             :     !
     278             :     !     .. Scalar Arguments ..
     279             :     REAL,    INTENT (IN) :: r0,dx
     280             :     !     ..
     281             :     !     .. Array Arguments ..
     282             :     REAL,    INTENT (IN) :: a(:),b(:),vso(:)
     283             :     !     ..
     284             :     !     .. Local Arrays ..
     285       31696 :     REAL q(size(a))
     286             :     !     ..
     287    23149056 :     q = a*b*vso
     288       31696 :     CALL intgr0(q,r0,dx,size(a),radso)
     289             : 
     290       31696 :     RETURN
     291             :   END FUNCTION radso
     292             :   !--------------------------------------------------------------------
     293             : END MODULE m_sorad

Generated by: LCOV version 1.14