LCOV - code coverage report
Current view: top level - eigen_soc - sorad.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 79 110 71.8 %
Date: 2019-09-08 04:53:50 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         192 :   SUBROUTINE sorad(atoms,input,ntyp,vr,enpara,spav,rsoc,usdus)
      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             :     !     ..
      28             :     !     .. Scalar Arguments ..
      29             :     INTEGER, INTENT (IN) :: ntyp
      30             :     LOGICAL, INTENT (IN) :: spav ! if T, spin-averaged pot is used
      31             :     !     ..
      32             :     !     .. Array Arguments ..
      33             :     REAL,    INTENT (IN) :: vr(:,:)!(atoms%jmtd,input%jspins),
      34             :     !     ..
      35             :     !     .. Local Scalars ..
      36             :     REAL ddn1,e ,ulops,dulops,duds1
      37             :     INTEGER i,j,ir,jspin,l,noded,nodeu,ilo,ilop
      38             :     !     ..
      39             :     !     .. Local Arrays ..
      40         384 :     REAL, ALLOCATABLE :: p(:,:),pd(:,:),q(:,:),qd(:,:),plo(:,:)
      41         384 :     REAL, ALLOCATABLE :: plop(:,:),glo(:,:),fint(:),pqlo(:,:)
      42         192 :     REAL, ALLOCATABLE :: filo(:,:)
      43         192 :     REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:)
      44             :     !     ..
      45             :     
      46         192 :     IF (atoms%jri(ntyp)>atoms%jmtd)  CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
      47             :     ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
      48         192 :          &   qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2) )
      49             : 
      50         576 :     p = 0.0 ; pd = 0.0 ; q = 0.0 ; plo = 0.0 ; fint = 0.0
      51         576 :     qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0
      52             :     
      53             : 
      54        1924 :     DO l = 0,atoms%lmax(ntyp) 
      55             : 
      56        5052 :        DO jspin = 1,input%jspins
      57             :           !
      58             :           !--->    calculate normalized function at e: p and q 
      59             :           !
      60        3320 :           e = enpara%el0(l,ntyp,jspin)
      61             :           CALL radsra(&
      62             :                e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
      63             :                usdus%us(l,ntyp,jspin),usdus%dus(l,ntyp,jspin),&
      64        3320 :                nodeu,p(:,jspin),q(:,jspin))
      65             :           !                     
      66             :           !--->    calculate orthogonal energy derivative at e : pd and qd
      67             :           !
      68             :           CALL radsrd(&
      69             :                e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
      70             :                usdus%uds(l,ntyp,jspin),usdus%duds(l,ntyp,jspin),&
      71             :                usdus%ddn(l,ntyp,jspin),noded,pd(:,jspin),qd(:,jspin),&
      72        5052 :                p(:,jspin),q(:,jspin),usdus%dus(l,ntyp,jspin))
      73             : 
      74             :        END DO     ! end of spin loop
      75             :        !
      76             :        !---> in case of jspins=1
      77             :        !
      78             : 
      79        1732 :        IF (input%jspins.EQ.1) THEN
      80       79520 :           DO i = 1,atoms%jri(ntyp)
      81       79376 :              p(i,2) =  p(i,1)
      82       79520 :              pd(i,2) = pd(i,1)
      83             :           ENDDO
      84             :        ENDIF
      85             :        !
      86             :        !---> common spin-orbit integrant V   (average spin directions)
      87             :        !                                  SO
      88     1043580 :        v0(:) = 0.0
      89        1732 :        IF (input%jspins.EQ.1) THEN
      90         144 :           v0(1:atoms%jri(ntyp)) = vr(1:atoms%jri(ntyp),1)
      91         144 :           e = enpara%el0(l,ntyp,1)
      92             :        ELSE
      93      920604 :           DO i = 1,atoms%jri(ntyp)
      94      920604 :              v0(i) = (vr(i,1)+vr(i,input%jspins))/2.
      95             :           END DO
      96        1588 :           e = (enpara%el0(l,ntyp,1)+enpara%el0(l,ntyp,input%jspins))/2.
      97             :        END IF
      98             : 
      99        1732 :        CALL sointg(ntyp,e,vr,v0,atoms,input,vso)
     100        1732 :        IF (spav) THEN
     101           0 :           DO i= 1,atoms%jmtd
     102           0 :              vso(i,1)= (vso(i,1)+vso(i,2))/2.
     103           0 :              vso(i,2)= vso(i,1)
     104             :           ENDDO
     105             :        ENDIF
     106             : 
     107             :        !                        s       s'            .s       s'
     108             :        !-->  radial integrals <u  |V  |u  > = rsopp, <u  |V  |u  > = rsopdp etc.
     109             :        !                            SO                     SO
     110             : 
     111        1732 :        IF (l.GT.0) THEN ! there is no spin-orbit for s-states
     112        4620 :           DO i = 1, 2
     113       16940 :              DO j = 1, 2
     114        6160 :                 rsoc%rsopp(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     115        6160 :                 rsoc%rsopdp(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     116        6160 :                 rsoc%rsoppd(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     117        9240 :                 rsoc%rsopdpd(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     118             :              ENDDO
     119             :           ENDDO
     120             :        ENDIF ! l>0
     121             :        !
     122             :        !--->  Check for local orbitals with same l
     123             :        !
     124        2030 :        DO ilo = 1, atoms%nlo(ntyp)
     125        1838 :           IF (atoms%llo(ilo,ntyp).EQ.l) THEN
     126             : 
     127          22 :              DO jspin = 1,input%jspins
     128          12 :                 e = enpara%ello0(ilo,ntyp,jspin)
     129             :                 CALL radsra(&
     130             :                      e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     131             :                      usdus%ulos(ilo,ntyp,jspin),usdus%dulos(ilo,ntyp,jspin),&
     132          12 :                      nodeu,plo(:,jspin),qlo(:,jspin))
     133             : 
     134             :                 !+apw+lo
     135          22 :                 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
     136           0 :                    ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
     137           0 :                    glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
     138           0 :                    pqlo(1:atoms%jri(ntyp),1)=plo(1:atoms%jri(ntyp),jspin)
     139           0 :                    pqlo(1:atoms%jri(ntyp),2)=qlo(1:atoms%jri(ntyp),jspin)
     140           0 :                    i = atoms%ulo_der(ilo,ntyp)
     141           0 :                    IF(atoms%l_dulo(ilo,ntyp)) i=1
     142             :                    CALL radsrdn(&
     143             :                         e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     144             :                         usdus%ulos(ilo,ntyp,jspin),duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
     145           0 :                         pqlo,usdus%dulos(ilo,ntyp,jspin),i)
     146           0 :                    ddn1 = SQRT(ddn1)
     147           0 :                    IF(atoms%l_dulo(ilo,ntyp)) ddn1=1.0
     148           0 :                    plo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)/ddn1
     149           0 :                    qlo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),2)/ddn1
     150           0 :                    usdus%dulos(ilo,ntyp,jspin) = duds1/ddn1
     151           0 :                    DEALLOCATE (glo,pqlo,filo)
     152             :                 ENDIF
     153             :                 !-apw+lo
     154             :              ENDDO
     155             : 
     156          10 :              IF (input%jspins.EQ.1) THEN
     157           8 :                 plo(1:atoms%jri(ntyp),2) = plo(1:atoms%jri(ntyp),1)
     158           8 :                 e = (enpara%ello0(ilo,ntyp,1) + enpara%el0(l,ntyp,1) )/2
     159             :              ELSE
     160             :                 e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +&
     161           2 :                      enpara%el0(l,ntyp,1) + enpara%el0(l,ntyp,input%jspins) )/4
     162             :              END IF
     163          10 :              CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
     164          10 :              IF (spav) THEN
     165           0 :                 DO i= 1,atoms%jmtd
     166           0 :                    vso(i,1)= (vso(i,1)+vso(i,2))/2.
     167           0 :                    vso(i,2)= vso(i,1)
     168             :                 ENDDO
     169             :              ENDIF
     170             : 
     171          30 :              DO i = 1, 2
     172         110 :                 DO j = 1, 2
     173          40 :                    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))
     174          40 :                    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))
     175          40 :                    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))
     176          60 :                    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))
     177             :                 ENDDO
     178             :              ENDDO
     179             : 
     180          22 :              DO i = 1,input%jspins
     181          12 :                 fint(:) = plo(:,i) *  p(:,i) + qlo(:,i) *  q(:,i)
     182          12 :                 CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,i))
     183        9096 :                 fint(:) = plo(:,i) * pd(:,i) + qlo(:,i) * qd(:,i)
     184          22 :                 CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,i))
     185             :              ENDDO
     186             : 
     187          20 :              DO ilop = 1, atoms%nlo(ntyp)
     188          20 :                 IF (atoms%llo(ilop,ntyp).EQ.l) THEN
     189             : 
     190          22 :                    DO jspin = 1,input%jspins
     191          12 :                       e = enpara%ello0(ilop,ntyp,jspin)
     192             :                       CALL radsra(&
     193             :                            e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     194          12 :                            ulops,dulops,nodeu,plop(:,jspin),q(:,1))
     195             :                       !+apw+lo
     196          22 :                       IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate orthogonal energy derivative at e
     197           0 :                          ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
     198           0 :                          glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
     199           0 :                          pqlo(1:atoms%jri(ntyp),1)=plop(1:atoms%jri(ntyp),jspin)
     200           0 :                          pqlo(1:atoms%jri(ntyp),2)=q(1:atoms%jri(ntyp),1)
     201           0 :                          i = atoms%ulo_der(ilo,ntyp)
     202           0 :                          IF(atoms%l_dulo(ilo,ntyp)) i=1
     203             :                          CALL radsrdn(&
     204             :                               e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
     205             :                               ulops,duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
     206           0 :                               pqlo,dulops,i)
     207           0 :                          plop(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)
     208           0 :                          DEALLOCATE (glo,pqlo,filo)
     209             :                       ENDIF
     210             :                       !-apw+lo
     211             :                    ENDDO
     212             : 
     213          10 :                    IF (input%jspins.EQ.1) THEN
     214           8 :                       plop(1:atoms%jri(ntyp),2) = plop(1:atoms%jri(ntyp),1)
     215           8 :                       e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilop,ntyp,1) )/2
     216             :                    ELSE
     217             :                       e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +  &
     218           2 :                            enpara%ello0(ilop,ntyp,1) + enpara%ello0(ilop,ntyp,input%jspins) )/4
     219             :                    END IF
     220          10 :                    CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
     221          10 :                    IF (spav) THEN
     222           0 :                       DO i= 1,atoms%jmtd
     223           0 :                          vso(i,1)= (vso(i,1)+vso(i,2))/2.
     224           0 :                          vso(i,2)= vso(i,1)
     225             :                       ENDDO
     226             :                    ENDIF
     227             : 
     228          30 :                    DO i = 1, 2
     229         110 :                       DO j = 1, 2
     230             :                          rsoc%rsoploplop(ntyp,ilo,ilop,i,j) =&
     231          60 :                               radso(plo(:atoms%jri(ntyp),i),plop(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
     232             :                       ENDDO
     233             :                    ENDDO
     234             : 
     235             :                 ENDIF
     236             :              ENDDO
     237             : 
     238             :           ENDIF
     239             :        ENDDO ! end of lo-loop
     240             : 
     241             :     ENDDO ! end of l-loop
     242             : 
     243         192 :     DEALLOCATE ( p,pd,q,qd,plo,plop,qlo,fint,v0,vso )
     244             :     !      rsoplop (:,:,:,:) = 0.0
     245             :     !      rsoplopd(:,:,:,:) = 0.0
     246             :     !      rsopplo (:,:,:,:) = 0.0
     247             :     !      rsopdplo(:,:,:,:) = 0.0
     248             :     !      rsoploplop(:,:,:,:,:) = 0.0
     249             : 
     250         192 :   END SUBROUTINE sorad
     251             :   !--------------------------------------------------------------------
     252       24840 :   REAL FUNCTION radso(a,b,vso,dx,r0)
     253             :     !
     254             :     !     compute radial spin-orbit integrals
     255             :     !
     256             :     USE m_intgr, ONLY : intgr0
     257             :     IMPLICIT NONE
     258             :     !
     259             :     !     .. Scalar Arguments ..
     260             :     REAL,    INTENT (IN) :: r0,dx
     261             :     !     ..
     262             :     !     .. Array Arguments ..
     263             :     REAL,    INTENT (IN) :: a(:),b(:),vso(:)
     264             :     !     ..
     265             :     !     .. Local Arrays ..
     266       49680 :     REAL q(size(a))
     267             :     !     ..
     268    14376752 :     q = a*b*vso
     269       24840 :     CALL intgr0(q,r0,dx,size(a),radso)
     270             : 
     271       24840 :     RETURN
     272             :   END FUNCTION radso
     273             :   !--------------------------------------------------------------------
     274             : END MODULE m_sorad

Generated by: LCOV version 1.13