LCOV - code coverage report
Current view: top level - eigen - vacfun.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 73 103 70.9 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_vacfun
       2             :   use m_juDFT
       3             : CONTAINS
       4          64 :   SUBROUTINE vacfun(&
       5             :        vacuum,stars,input,noco,jspin1,jspin2,&
       6         128 :        sym, cell,ivac,evac,bkpt, vxy,vz,kvac,nv2,&
       7         128 :        tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
       8             :     !*********************************************************************
       9             :     !     determines the necessary values and derivatives on the vacuum
      10             :     !     boundary (ivac=1 upper vacuum; ivac=2, lower) for energy
      11             :     !     parameter evac.  also sets up the 2d hamiltonian matrices
      12             :     !     necessary to update the full hamiltonian matrix.
      13             :     !               m. weinert
      14             :     !*********************************************************************
      15             : 
      16             :     USE m_intgr, ONLY : intgz0
      17             :     USE m_vacuz
      18             :     USE m_vacudz
      19             :     USE m_types
      20             :     IMPLICIT NONE
      21             : 
      22             :     TYPE(t_input),INTENT(IN)       :: input
      23             :     TYPE(t_vacuum),INTENT(IN)      :: vacuum
      24             :     TYPE(t_noco),INTENT(IN)        :: noco
      25             :     TYPE(t_sym),INTENT(IN)         :: sym
      26             :     TYPE(t_stars),INTENT(IN)       :: stars
      27             :     TYPE(t_cell),INTENT(IN)        :: cell
      28             :     !     ..
      29             :     !     .. Scalar Arguments ..
      30             :     INTEGER, INTENT (IN) :: ivac,jspin1,jspin2
      31             :     REAL,    INTENT (OUT) :: wronk
      32             :     !     ..
      33             :     !     .. Array Arguments ..
      34             :     INTEGER, INTENT (IN) :: nv2(:)!(input%jspins)
      35             :     INTEGER, INTENT (IN) :: kvac(:,:,:)!(2,dimension%nv2d,input%jspins)
      36             :     COMPLEX, INTENT (IN) :: vxy(:,:,:,:) !(vacuum%nmzxyd,stars%ng2-1,nvac,:)
      37             :     COMPLEX, INTENT (OUT):: tddv(:,:),tduv(:,:)!(dimension%nv2d,dimension%nv2d)
      38             :     COMPLEX, INTENT (OUT):: tudv(:,:),tuuv(:,:)!(dimension%nv2d,dimension%nv2d)
      39             :     REAL,ALLOCATABLE,INTENT (IN) :: vz(:,:,:) !(vacuum%nmzd,2,4) ,
      40             :     REAL,    INTENT (IN) :: evac(:,:)!(2,input%jspins)
      41             :     REAL,    INTENT (IN) :: bkpt(3) 
      42             :     REAL,    INTENT (OUT):: udz(:,:),uz(:,:)!(dimension%nv2d,input%jspins)
      43             :     REAL,    INTENT (OUT):: dudz(:,:),duz(:,:)!(dimension%nv2d,input%jspins)
      44             :     REAL,    INTENT (OUT):: ddnv(:,:)!(dimension%nv2d,input%jspins)
      45             :     !     ..
      46             :     !     .. Local Scalars ..
      47             :     REAL ev,scale,xv,yv,vzero,fac
      48             :     COMPLEX phase
      49             :     INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,ipot
      50             :     LOGICAL tail
      51             :     !     ..
      52             :     !     .. Local Arrays ..
      53         128 :     REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins)
      54         128 :     REAL v(3),x(vacuum%nmzd), qssbti(2,2)
      55             :     !     ..
      56          64 :     fac=MERGE(1.0,-1.0,jspin1>=jspin2)
      57          64 :     ipot=MERGE(jspin1,3,jspin1==jspin2)
      58             : 
      59        2784 :     tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0
      60         168 :     udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0.
      61          64 :     tail = .true.
      62          64 :     np1 = vacuum%nmzxy + 1
      63             :     !--->    wronksian for the schrodinger equation given by an identity
      64          64 :     wronk = 2.0
      65             :     !---> setup the spin-spiral q-vector
      66         192 :     qssbti(1:2,1) = - noco%qss(1:2)/2
      67         192 :     qssbti(1:2,2) = + noco%qss(1:2)/2
      68             :     !--->    generate basis functions for each 2-d k+g
      69         128 :     DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
      70        2728 :        DO  ik = 1,nv2(jspin)
      71        2600 :           v(1:2) = bkpt(1:2) + kvac(:,ik,jspin) + qssbti(1:2,jspin)
      72        2600 :           v(3) = 0.0
      73        2600 :           ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat))
      74        2600 :           vzero = vz(vacuum%nmzd,ivac,jspin)
      75             :           CALL vacuz(ev,vz(1,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
      76        2600 :                uz(ik,jspin),duz(ik,jspin),u(1,ik,jspin))
      77             :           CALL vacudz(ev,vz(1,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
      78             :                udz(ik,jspin),dudz(ik,jspin),ddnv(ik,jspin),&
      79        2600 :                ud(1,ik,jspin),duz(ik,jspin),u(1,ik,jspin))
      80             :           !--->       make sure the solutions satisfy the wronksian
      81             :           scale = wronk/ (udz(ik,jspin)*duz(ik,jspin)-&
      82        2600 :                &                         dudz(ik,jspin)*uz(ik,jspin))
      83        2600 :           udz(ik,jspin)  = scale*udz(ik,jspin)
      84        2600 :           dudz(ik,jspin) = scale*dudz(ik,jspin)
      85        2600 :           ddnv(ik,jspin) = scale*ddnv(ik,jspin)
      86        2664 :           ud(:,ik,jspin) = scale*ud(:,ik,jspin)
      87             :        enddo
      88             :     ENDDO
      89             :     !--->    set up the tuuv, etc. matrices
      90        2664 :     DO  ik = 1,nv2(jspin1)
      91      111968 :        DO  jk = 1,nv2(jspin2)
      92             : 
      93             :           !--->     determine the warping component of the potential
      94      109304 :           i1 = kvac(1,ik,jspin1) - kvac(1,jk,jspin2)
      95      109304 :           i2 = kvac(2,ik,jspin1) - kvac(2,jk,jspin2)
      96      109304 :           i3 = 0
      97      109304 :           ind3 = stars%ig(i1,i2,i3)
      98      109304 :           IF (ind3.EQ.0) CYCLE
      99      109304 :           phase = stars%rgphs(i1,i2,i3)
     100      109304 :           ind2 = stars%ig2(ind3)
     101      109304 :           IF (ind2.EQ.0) THEN
     102           0 :              WRITE (6,FMT=8000) ik,jk
     103             : 8000         FORMAT (' **** error in map2 for 2-d stars',2i5)
     104           0 :              CALL juDFT_error("error in map2 for 2-d stars",calledby ="vacfun")
     105             :           END IF
     106             :           !--->     get the proper warping index (vxy starts with the 2nd star)
     107      109304 :           ind2 = ind2 - 1
     108      111904 :           IF (ind2.NE.0) THEN
     109             :              !--->       only the warping part, 1st star (G=0) is done later
     110             : 
     111             :              !--->       obtain the warping matrix elements
     112             :              !--->       note that the tail correction (tail=.true.) is included for
     113             :              !--->       the integrals, i.e. the integrand is from infinity inward
     114             : 
     115             :              !--->       tuuv
     116    21447504 :              DO  i = 1,vacuum%nmzxy
     117    10777104 :                 x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
     118             :              enddo
     119      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     120    10777104 :              DO  i = 1,vacuum%nmzxy
     121    10777104 :                 x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
     122             :              enddo
     123      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     124      106704 :              tuuv(ik,jk) = phase*cmplx(xv,yv)
     125             : 
     126             :              !--->       tddv
     127    10777104 :              DO  i = 1,vacuum%nmzxy
     128    10777104 :                 x(np1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
     129             :              enddo
     130      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     131    10777104 :              DO  i = 1,vacuum%nmzxy
     132    10777104 :                 x(np1-i) =ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
     133             :              enddo
     134      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     135      106704 :              tddv(ik,jk) = phase*cmplx(xv,yv)
     136             : 
     137             :              !--->       tudv
     138    10777104 :              DO  i = 1,vacuum%nmzxy
     139    10777104 :                 x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
     140             :              enddo
     141      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     142    10777104 :              DO  i = 1,vacuum%nmzxy
     143    10777104 :                 x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
     144             :              enddo
     145      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     146      106704 :              tudv(ik,jk) = phase*cmplx(xv,yv)
     147             : 
     148             :              !--->       tduv
     149    10777104 :              DO  i = 1,vacuum%nmzxy
     150    10777104 :                 x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
     151             :              enddo
     152      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
     153    10777104 :              DO  i = 1,vacuum%nmzxy
     154    10777104 :                 x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
     155             :              enddo
     156      106704 :              CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
     157      106704 :              tduv(ik,jk) = phase*cmplx(xv,yv)
     158             : 
     159             :           ELSE
     160             :              !--->       diagonal (film muffin-tin) terms
     161        2600 :              IF (jspin1==jspin2) THEN
     162        2600 :                 tuuv(ik,ik) = cmplx(evac(ivac,jspin1),0.0)
     163        2600 :                 tddv(ik,ik) = cmplx(evac(ivac,jspin1)*ddnv(ik,jspin1),0.0)
     164        2600 :                 tudv(ik,ik) = cmplx(0.5,0.0)
     165        2600 :                 tduv(ik,ik) = cmplx(0.5,0.0)
     166             :              ELSE
     167             :                 !--->          tuuv
     168           0 :                 DO i = 1,vacuum%nmz
     169           0 :                    x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
     170             :                 ENDDO
     171           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     172           0 :                 DO i = 1,vacuum%nmz
     173           0 :                    x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
     174             :                 ENDDO
     175           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     176           0 :                 tuuv(ik,jk) = cmplx(xv,yv)
     177             : 
     178             :                 !--->          tddv
     179           0 :                 DO i = 1,vacuum%nmz
     180           0 :                    x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
     181             :                 ENDDO
     182           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     183           0 :                 DO i = 1,vacuum%nmz
     184           0 :                    x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
     185             :                 ENDDO
     186           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     187           0 :                 tddv(ik,jk) = cmplx(xv,yv)
     188             : 
     189             :                 !--->          tudv
     190           0 :                 DO i = 1,vacuum%nmz
     191           0 :                    x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
     192             :                 ENDDO
     193           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     194           0 :                 DO i = 1,vacuum%nmz
     195           0 :                    x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
     196             :                 ENDDO
     197           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     198           0 :                 tudv(ik,jk) = cmplx(xv,yv)
     199             : 
     200             :                 !--->          tduv
     201           0 :                 DO i = 1,vacuum%nmz
     202           0 :                    x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
     203             :                 ENDDO
     204           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
     205           0 :                 DO i = 1,vacuum%nmz
     206           0 :                    x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
     207             :                 ENDDO
     208           0 :                 CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
     209           0 :                 tduv(ik,jk) = cmplx(xv,yv)
     210             :              ENDIF
     211             : 
     212             :           ENDIF
     213             :        enddo
     214             :     enddo
     215             : 
     216          64 :   END SUBROUTINE vacfun
     217             : END MODULE m_vacfun
     218             : 

Generated by: LCOV version 1.13