LCOV - code coverage report
Current view: top level - eigen - od_vacfun.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 145 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.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_od_vacfun
       8             : CONTAINS
       9           0 :   SUBROUTINE od_vacfun(&
      10             :        m_cyl,cell,vacuum,DIMENSION,stars,&
      11           0 :        jsp,input,noco,ipot,oneD, n2d_1, ivac,evac,bkpt,MM,vM,&
      12           0 :        vxy,vz,kvac3,nv2, tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv)
      13             :     !*********************************************************************
      14             :     !     determines the necessary values and derivatives on the cylindrical
      15             :     !     vacuum boundary for energy parameter evac. also sets up the
      16             :     !     hamiltonian matrices necessary to update the full hamiltonian
      17             :     !     matrix.
      18             :     !     Y. Mokrousov, June 2002
      19             :     !*********************************************************************
      20             :     USE m_vacuz
      21             :     USE m_vacudz
      22             :     USE m_intgr, ONLY : intgz0
      23             :     USE m_types
      24             :     IMPLICIT NONE
      25             :     TYPE(t_dimension),INTENT(IN):: DIMENSION
      26             :     TYPE(t_oneD),INTENT(IN)     :: oneD
      27             :     TYPE(t_input),INTENT(IN)    :: input
      28             :     TYPE(t_vacuum),INTENT(IN)   :: vacuum
      29             :     TYPE(t_noco),INTENT(IN)     :: noco
      30             :     TYPE(t_stars),INTENT(IN)    :: stars
      31             :     TYPE(t_cell),INTENT(IN)     :: cell
      32             :     !     ..
      33             :     !     .. Scalar Arguments ..
      34             : 
      35             :     INTEGER, INTENT (IN) :: jsp   ,ipot,MM ,vM
      36             :     INTEGER, INTENT (IN) :: ivac,n2d_1,m_cyl
      37             :     !     ..
      38             :     !     .. Array Arguments ..
      39             :     INTEGER, INTENT (IN) :: nv2(input%jspins)
      40             :     INTEGER, INTENT (IN) :: kvac3(DIMENSION%nv2d,input%jspins)
      41             :     COMPLEX, INTENT (IN) :: vxy(vacuum%nmzxyd,n2d_1-1)
      42             :     COMPLEX, INTENT (OUT):: tddv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d)
      43             :     COMPLEX, INTENT (OUT):: tduv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d)
      44             :     COMPLEX, INTENT (OUT):: tudv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d)
      45             :     COMPLEX, INTENT (OUT):: tuuv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d)
      46             :     REAL,    INTENT (IN) :: vz(vacuum%nmzd,2,4) ,evac(2,input%jspins)
      47             :     REAL,    INTENT (IN) :: bkpt(3) 
      48             :     REAL,    INTENT (OUT):: udz(-vM:vM,DIMENSION%nv2d,input%jspins),uz(-vM:vM,DIMENSION%nv2d,input%jspins)
      49             :     REAL,    INTENT (OUT):: dudz(-vM:vM,DIMENSION%nv2d,input%jspins)
      50             :     REAL,    INTENT (OUT):: duz(-vM:vM,DIMENSION%nv2d,input%jspins)
      51             :     REAL,    INTENT (OUT):: ddnv(-vM:vM,DIMENSION%nv2d,input%jspins)
      52             :     !     ..
      53             :     !     .. Local Scalars ..
      54             :     REAL ev,scale,xv,yv,vzero,v1,wronk
      55             :     INTEGER i,ik,jk,np1,jspin,jsp1,jsp2 ,l,m
      56             :     INTEGER i1,i2,i3,ind1,ind3
      57             :     LOGICAL tail
      58             :     !     ..
      59             :     !     .. Local Arrays ..
      60           0 :     REAL wdz(-vM:vM,DIMENSION%nv2d,input%jspins),wz(-vM:vM,DIMENSION%nv2d,input%jspins)
      61           0 :     REAL dwdz(-vM:vM,DIMENSION%nv2d,input%jspins),dwz(-vM:vM,DIMENSION%nv2d,input%jspins)
      62           0 :     REAL u(vacuum%nmzd,-vM:vM,DIMENSION%nv2d,input%jspins),ud(vacuum%nmzd,-vM:vM,DIMENSION%nv2d,input%jspins)
      63           0 :     REAL v(3),x(vacuum%nmzd)
      64           0 :     REAL vr0(vacuum%nmzd,2,4)
      65           0 :     REAL w(vacuum%nmzd,-vM:vM,DIMENSION%nv2d,input%jspins),wd(vacuum%nmzd,-vM:vM,DIMENSION%nv2d,input%jspins)
      66             :     REAL qssbti(2)
      67             :     !     ..
      68             : 
      69             : 
      70           0 :     tail = .TRUE.
      71           0 :     np1 = vacuum%nmzxy + 1
      72             : 
      73             :     !     wronksian for the schrodinger equation given by an identity
      74             : 
      75           0 :     wronk = 2.0
      76             : 
      77           0 :     qssbti(1) = - noco%qss(3)/2.
      78           0 :     qssbti(2) = + noco%qss(3)/2.
      79             : 
      80           0 :     tuuv(:,:,:,:) = CMPLX(0.,0.)
      81           0 :     tudv(:,:,:,:) = CMPLX(0.,0.)
      82           0 :     tduv(:,:,:,:) = CMPLX(0.,0.)
      83           0 :     tddv(:,:,:,:) = CMPLX(0.,0.)
      84             : 
      85             :     !     generate basis functions for each 1-d k_z+g_z and m
      86             : 
      87           0 :     DO jspin = 1,input%jspins
      88           0 :        DO  ik = 1,nv2(jspin)
      89           0 :           DO  m = 0,vM
      90           0 :              v(1) = 0.0
      91           0 :              v(2) = 0.0
      92           0 :              v(3) = bkpt(3) + kvac3(ik,jspin) + qssbti(jspin)
      93           0 :              ev = evac(ivac,jspin) - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
      94             :              !     constructing of the 'pseudopotential'
      95           0 :              DO  i=1,vacuum%nmzd
      96             :                 v1 = 1./(8.*((cell%z1+(i-1)*vacuum%delz)**2))&
      97           0 :                      -(m*m)/(2.*((cell%z1+(i-1)*vacuum%delz)**2))
      98           0 :                 vr0(i,ivac,jspin) = vz(i,ivac,jspin)-v1
      99             :              ENDDO
     100           0 :              vzero = vr0(vacuum%nmzd,ivac,jspin)
     101             :              !     obtaining solutions with the 'pseudopotential'
     102             : 
     103             :              CALL vacuz(ev,vr0(:,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
     104           0 :                   wz(m,ik,jspin),dwz(m,ik,jspin),w(1,m,ik,jspin))
     105             :              CALL vacudz(ev,vr0(1,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
     106             :                   wdz(m,ik,jspin),dwdz(m,ik,jspin),ddnv(m,ik,jspin),&
     107           0 :                   wd(1,m,ik,jspin),dwz(m,ik,jspin),w(1,m,ik,jspin))
     108             :              !     make sure the solutions satisfy the wronksian
     109             :              scale = wronk/ (wdz(m,ik,jspin)*dwz(m,ik,jspin)-&
     110           0 :                   dwdz(m,ik,jspin)*wz(m,ik,jspin))
     111           0 :              wdz(m,ik,jspin) = scale*wdz(m,ik,jspin)
     112           0 :              dwdz(m,ik,jspin) = scale*dwdz(m,ik,jspin)
     113           0 :              ddnv(m,ik,jspin) = scale*ddnv(m,ik,jspin)
     114           0 :              IF (m.GT.0) THEN
     115           0 :                 wdz(-m,ik,jspin) = wdz(m,ik,jspin)
     116           0 :                 dwdz(-m,ik,jspin) = dwdz(m,ik,jspin)
     117           0 :                 ddnv(-m,ik,jspin) = ddnv(m,ik,jspin)
     118             :              END IF
     119           0 :              DO  i = 1,vacuum%nmz
     120           0 :                 wd(i,m,ik,jspin) = scale*wd(i,m,ik,jspin)
     121           0 :                 w(i,m,ik,jspin) = scale*w(i,m,ik,jspin)
     122           0 :                 IF (m.GT.0) THEN
     123           0 :                    wd(i,-m,ik,jspin) = wd(i,m,ik,jspin)
     124           0 :                    w(i,-m,ik,jspin) = w(i,m,ik,jspin)
     125             :                 END IF
     126             :              ENDDO
     127             :              !     constructing 'real' solutions
     128           0 :              DO  i=1,vacuum%nmz
     129           0 :                 u(i,m,ik,jspin)=w(i,m,ik,jspin)/SQRT(cell%z1+(i-1)*vacuum%delz)
     130           0 :                 ud(i,m,ik,jspin)=wd(i,m,ik,jspin)/SQRT(cell%z1+(i-1)*vacuum%delz)
     131           0 :                 IF (m.GT.0) THEN
     132           0 :                    u(i,-m,ik,jspin) = u(i,m,ik,jspin)
     133           0 :                    ud(i,-m,ik,jspin) = ud(i,m,ik,jspin)
     134             :                 END IF
     135             :              ENDDO
     136             :              duz(m,ik,jspin)=(-dwz(m,ik,jspin))/SQRT(cell%z1)-&
     137           0 :                   wz(m,ik,jspin)/(2.0*((cell%z1)**(1.5)))
     138           0 :              uz(m,ik,jspin)=wz(m,ik,jspin)/SQRT(cell%z1)
     139             :              dudz(m,ik,jspin)=(-dwdz(m,ik,jspin))/SQRT(cell%z1)-&
     140           0 :                   wdz(m,ik,jspin)/(2.0*((cell%z1)**(1.5))) 
     141           0 :              udz(m,ik,jspin)=wdz(m,ik,jspin)/SQRT(cell%z1)
     142           0 :              IF (m.GT.0) THEN
     143           0 :                 duz(-m,ik,jspin) = duz(m,ik,jspin)
     144           0 :                 uz(-m,ik,jspin) = uz(m,ik,jspin)
     145           0 :                 dudz(-m,ik,jspin) = dudz(m,ik,jspin)
     146           0 :                 udz(-m,ik,jspin) = udz(m,ik,jspin)
     147             :              END IF
     148             :           ENDDO
     149             :        ENDDO
     150             :     ENDDO
     151             : 
     152             :     !     set up the tuuv, etc. matrices
     153             : 
     154           0 :     IF (noco%l_noco) THEN
     155           0 :        IF (ipot.EQ.1) THEN
     156             :           jsp1 = 1
     157             :           jsp2 = 1
     158           0 :        ELSEIF (ipot.EQ.2) THEN
     159             :           jsp1 = 2
     160             :           jsp2 = 2
     161           0 :        ELSEIF (ipot.EQ.3) THEN
     162           0 :           jsp1 = 2
     163           0 :           jsp2 = 1
     164             :        ENDIF
     165             :     ELSE
     166           0 :        jsp1 = jsp
     167           0 :        jsp2 = jsp
     168             :     ENDIF
     169             : 
     170           0 :     DO  ik = 1,nv2(jsp1)
     171           0 :        DO  jk = 1,nv2(jsp2)
     172           0 :           i1 = 0
     173           0 :           i2 = 0
     174           0 :           i3 = kvac3(ik,jsp1) - kvac3(jk,jsp2)
     175           0 :           ind3 = stars%ig(i1,i2,i3)
     176           0 :           IF (ind3.EQ.0) CYCLE
     177           0 :           DO  m = -vM,vM
     178           0 :              DO  l = -vM,vM
     179           0 :                 IF (l.EQ.m .OR. (iabs(m).LE.m_cyl&
     180           0 :                      &                             .AND. iabs(l).LE.m_cyl)) THEN
     181             :                    !     determine the warping component of the potential
     182           0 :                    ind1 = oneD%ig1(i3,m-l)
     183           0 :                    IF (ind1.NE.0) THEN
     184           0 :                       IF(ind1.NE.1) THEN
     185           0 :                          ind1 = ind1 - 1
     186             :                          !     only the warping part
     187             :                          !--->             tuuv
     188           0 :                          DO i = 1,vacuum%nmzxy
     189           0 :                             x(np1-i) = w(i,m,ik,jsp1)*w(i,l,jk,jsp2) *REAL(vxy(i,ind1))
     190             :                          ENDDO
     191           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     192           0 :                          DO i = 1,vacuum%nmzxy
     193           0 :                             x(np1-i) = w(i,m,ik,jsp1)*w(i,l,jk,jsp2) *AIMAG(vxy(i,ind1))
     194             :                          ENDDO
     195           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     196           0 :                          tuuv(m,l,ik,jk) = CMPLX(xv,yv)
     197             :                          !--->             tddv
     198           0 :                          DO i = 1,vacuum%nmzxy
     199           0 :                             x(np1-i) = wd(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *REAL(vxy(i,ind1))
     200             :                          ENDDO
     201           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     202           0 :                          DO i = 1,vacuum%nmzxy
     203           0 :                             x(np1-i) = wd(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *AIMAG(vxy(i,ind1))
     204             :                          ENDDO
     205           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     206           0 :                          tddv(m,l,ik,jk) = CMPLX(xv,yv)
     207             :                          !--->             tudv
     208           0 :                          DO i = 1,vacuum%nmzxy
     209           0 :                             x(np1-i) = w(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *REAL(vxy(i,ind1))
     210             :                          ENDDO
     211           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     212           0 :                          DO i = 1,vacuum%nmzxy
     213           0 :                             x(np1-i) = w(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *AIMAG(vxy(i,ind1))
     214             :                          ENDDO
     215           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     216           0 :                          tudv(m,l,ik,jk) = CMPLX(xv,yv)
     217             :                          !--->             tduv
     218           0 :                          DO i = 1,vacuum%nmzxy
     219           0 :                             x(np1-i) = wd(i,m,ik,jsp1)*w(i,l,jk,jsp2) *REAL(vxy(i,ind1))
     220             :                          ENDDO
     221           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     222           0 :                          DO i = 1,vacuum%nmzxy
     223           0 :                             x(np1-i) = wd(i,m,ik,jsp1)*w(i,l,jk,jsp2) *AIMAG(vxy(i,ind1))
     224             :                          ENDDO
     225           0 :                          CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     226           0 :                          tduv(m,l,ik,jk) = CMPLX(xv,yv)
     227             :                       ELSE
     228             :                          !--->          diagonal terms
     229           0 :                          IF ((ipot.EQ.1).OR.(ipot.EQ.2)) THEN
     230           0 :                             tuuv(m,m,ik,ik) = CMPLX(evac(ivac,jsp1),0.0)
     231             :                             tddv(m,m,ik,ik) = CMPLX(evac(ivac,jsp1)*&
     232           0 :                                  &                       ddnv(m,ik,jsp1),0.0)
     233           0 :                             tudv(m,m,ik,ik) = CMPLX(0.5,0.0)
     234           0 :                             tduv(m,m,ik,ik) = CMPLX(0.5,0.0)
     235             :                          ELSE
     236             :                             !--->             tuuv
     237           0 :                             DO i = 1,vacuum%nmz
     238           0 :                                x(vacuum%nmz+1-i) = w(i,m,ik,jsp1)*w(i,l,jk,jsp2) *vz(i,ivac,3)
     239             :                             ENDDO
     240           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     241           0 :                             DO i = 1,vacuum%nmz
     242           0 :                                x(vacuum%nmz+1-i) = w(i,m,ik,jsp1)*w(i,l,jk,jsp2) *vz(i,ivac,4)
     243             :                             ENDDO
     244           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     245           0 :                             tuuv(m,l,ik,jk) = CMPLX(xv,yv)
     246             :                             !--->             tddv
     247           0 :                             DO i = 1,vacuum%nmz
     248           0 :                                x(vacuum%nmz+1-i) = wd(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *vz(i,ivac,3)
     249             :                             ENDDO
     250           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     251           0 :                             DO i = 1,vacuum%nmz
     252           0 :                                x(vacuum%nmz+1-i) = wd(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *vz(i,ivac,4)
     253             :                             ENDDO
     254           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     255           0 :                             tddv(m,l,ik,jk) = CMPLX(xv,yv)
     256             :                             !--->             tudv
     257           0 :                             DO i = 1,vacuum%nmz
     258           0 :                                x(vacuum%nmz+1-i) = w(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *vz(i,ivac,3)
     259             :                             ENDDO
     260           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     261           0 :                             DO i = 1,vacuum%nmz
     262           0 :                                x(vacuum%nmz+1-i) = w(i,m,ik,jsp1)*wd(i,l,jk,jsp2) *vz(i,ivac,4)
     263             :                             ENDDO
     264           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     265           0 :                             tudv(m,l,ik,jk) = CMPLX(xv,yv)
     266             :                             !--->             tduv
     267           0 :                             DO i = 1,vacuum%nmz
     268           0 :                                x(vacuum%nmz+1-i) = wd(i,m,ik,jsp1)*w(i,l,jk,jsp2) *vz(i,ivac,3)
     269             :                             ENDDO
     270           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     271           0 :                             DO i = 1,vacuum%nmz
     272           0 :                                x(vacuum%nmz+1-i) = wd(i,m,ik,jsp1)*w(i,l,jk,jsp2) *vz(i,ivac,4)
     273             :                             ENDDO
     274           0 :                             CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     275           0 :                             tduv(m,l,ik,jk) = CMPLX(xv,yv)
     276             : 
     277             :                          ENDIF !ipot
     278             :                       END IF !ind1 ne 1
     279             :                    ENDIF    ! ind1 ne 0
     280             :                 END IF
     281             :              ENDDO
     282             :           ENDDO
     283             :        ENDDO
     284             :     ENDDO
     285             : 
     286             : 
     287           0 :     RETURN
     288             :   END SUBROUTINE od_vacfun
     289             : END MODULE m_od_vacfun
     290             : 
     291             : 
     292             : 
     293             : 
     294             :       

Generated by: LCOV version 1.13