mcd_init.f90 6.21 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
MODULE m_mcdinit
CONTAINS
  SUBROUTINE mcd_init(atoms,input,DIMENSION,&
       vr,g,f,emcd_up,emcd_lo,itype,jspin, ncore,e_mcd,m_mcd)

    !-----------------------------------------------------------------------
    !
    ! For a given atom-type 'itype' look, whether a core state is in the
    ! energy interval [emcd_lo,emcd_up] and, if found, calculate the 
    ! MCD-matrix elements 'm_mcd'.
    !          
    !-----------------------------------------------------------------------

    USE m_nabla
    USE m_dr2fdr
    USE m_constants, ONLY : c_light
    USE m_setcor
    USE m_differ
    USE m_types
    IMPLICIT NONE

    TYPE(t_dimension),INTENT(IN) :: DIMENSION
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_atoms),INTENT(IN)     :: atoms

    INTEGER, PARAMETER :: l_max = 3

    ! Arguments ...

    INTEGER, INTENT (IN)  :: itype
    INTEGER, INTENT (IN)  :: jspin
    REAL,    INTENT (IN)  :: emcd_up,emcd_lo
    REAL,    INTENT (IN)  :: vr(atoms%jmtd,atoms%ntypd,input%jspins)
    REAL,    INTENT (IN)  :: f(atoms%jmtd,2,0:atoms%lmaxd,jspin:jspin)
    REAL,    INTENT (IN)  :: g(atoms%jmtd,2,0:atoms%lmaxd,jspin:jspin)
    INTEGER, INTENT (OUT) :: ncore(atoms%ntypd)
    REAL,    INTENT (OUT) :: e_mcd(atoms%ntypd,input%jspins,DIMENSION%nstd)
    COMPLEX, INTENT (OUT) :: m_mcd(:,:,:,:)

    ! Locals ...

    INTEGER kap,mue,iri,l,ispin,i,icore,korb,nst,n_core,ierr
    REAL  c,t2,e,fj,fl,fn ,d,ms,rn ,bmu
    INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd),l_core(DIMENSION%nstd)
    REAL vrd(DIMENSION%msh),occ(DIMENSION%nstd,1),a(DIMENSION%msh),b(DIMENSION%msh),j_core(DIMENSION%nstd),e_mcd1(DIMENSION%nstd)
    REAL gv1(atoms%jmtd)
    REAL, ALLOCATABLE :: gc(:,:,:),fc(:,:,:)
    REAL, ALLOCATABLE :: gv(:,:,:,:),fv(:,:,:,:),dgv(:,:,:,:)

    !-----------------------------------------------------------------------

    c = c_light(1.0)
    ALLOCATE ( gc(atoms%jri(itype),atoms%ncst(itype),input%jspins) )
    ALLOCATE ( fc(atoms%jri(itype),atoms%ncst(itype),input%jspins) )

    ! core setup

    ncore(itype) = 0
    bmu = 0.0
    CALL setcor(itype,1,atoms,bmu, nst,kappa,nprnc,occ)

    DO ispin = jspin, jspin

       ! extend core potential

       DO iri = 1, atoms%jri(itype)
          vrd(iri) = vr(iri,itype,ispin)
       ENDDO
       t2 = vrd(atoms%jri(itype)) / (atoms%jri(itype) - DIMENSION%msh)
       DO iri = atoms%jri(itype) + 1, DIMENSION%msh
          vrd(iri) =  vrd(atoms%jri(itype))  + t2* ( iri-atoms%jri(itype) )
       ENDDO

       ! calculate core

       n_core = 0
       DO korb = 1, atoms%ncst(itype)
          IF (occ(korb,1).GT.0) THEN
             fn = nprnc(korb)
             fj = iabs(kappa(korb)) - .5e0
             fl = fj + (.5e0)*isign(1,kappa(korb))
             e = -2* (atoms%zatom(itype)/ (fn+fl))**2
             d = EXP(atoms%dx(itype))
             rn = atoms%rmsh(1,itype)*( d**(DIMENSION%msh-1) )
             CALL differ(fn,fl,fj,c,atoms%zatom(itype),atoms%dx(itype),atoms%rmsh(1,itype),&
                  rn,d,DIMENSION%msh,vrd, e, a,b,ierr)
             IF (ierr/=0)  CALL juDFT_error("error in core-levels", calledby="mcd_init")
             IF ( (e.LE.emcd_up).AND.(e.GE.emcd_lo) ) THEN
                WRITE(*,*) 'good    ev = ',e
                n_core = n_core + 1
                j_core(n_core) = fj
                l_core(n_core) = NINT( fl )
                e_mcd1(n_core) = e
                DO iri = 1, atoms%jri(itype)
                   gc(iri,n_core,ispin) = a(iri)
                   fc(iri,n_core,ispin) = b(iri)
                ENDDO
             ENDIF
          ENDIF
       ENDDO

    ENDDO

    !-----------------------------------------------------------------------

    IF (n_core.GT.0) THEN

       ALLOCATE ( gv(atoms%jri(itype),0:l_max,input%jspins,2) )
       ALLOCATE (dgv(atoms%jri(itype),0:l_max,input%jspins,2) )
       ALLOCATE ( fv(atoms%jri(itype),0:l_max,input%jspins,2) )
       DO i = 1, 2
          DO iri = 3*(itype-1)+1 , 3*(itype-1)+3
             DO l = 1, (l_max+1)**2
                DO icore = 1, DIMENSION%nstd
                   m_mcd(icore,l,iri,i) = CMPLX(0.0,0.0)
                ENDDO
             ENDDO
          ENDDO
       ENDDO
       !
       ! bring LAPW wavefunctions in a proper form:
       !
       DO ispin = jspin, jspin
          ms = ispin - 1.5
          DO l = 0, l_max
             DO iri = 1, atoms%jri(itype)
                gv(iri,l,ispin,1) = f(iri,1,l,ispin)   ! large component of u
                fv(iri,l,ispin,1) = f(iri,2,l,ispin)   ! small              .
                gv(iri,l,ispin,2) = g(iri,1,l,ispin)   ! large component of u
                fv(iri,l,ispin,2) = g(iri,2,l,ispin)   ! small
             ENDDO
             gv1(:) = atoms%rmsh(:,itype) * gv(:,l,ispin,1)
             CALL dr2fdr(&                                          ! deriative of u (large)&
                  gv1,atoms%rmsh(1,itype),atoms%jri(itype), dgv(1,l,ispin,1) )
             gv1(:) = atoms%rmsh(:,itype) * gv(:,l,ispin,2)              !              .
             CALL dr2fdr(&                                          ! deriative of u (large)&
                  gv1,atoms%rmsh(1,itype),atoms%jri(itype), dgv(1,l,ispin,2) )
          ENDDO
          !
          !
          !
          DO icore = 1, n_core

             DO i = 1, 2
                !              write(*,*) j_core(icore),l_core(icore),l_max,ms
                CALL nabla(itype,icore,atoms%jri(itype),atoms%dx(itype),DIMENSION%nstd,atoms%ntypd,&
                     j_core(icore),l_core(icore),l_max,ms,atoms%rmsh(:,itype),gc(:,icore,ispin),&
                     gv(:,0:,ispin,i),dgv(:,0:,ispin,i), m_mcd(:,:,:,i) )
             ENDDO

             DO i = 1, 2*icore*l_core(icore)
                ncore(itype) = ncore(itype) + 1
                IF (ncore(itype)>DIMENSION%nstd)  CALL juDFT_error("dimension%nstd too small" ,calledby ="mcd_init")
                e_mcd(itype,ispin,ncore(itype)) = e_mcd1(icore)
             ENDDO
          ENDDO
       ENDDO

       DEALLOCATE (gv,fv,dgv)
    ENDIF
    DEALLOCATE (gc,fc)


    !      DO i = 1, 2
    !       DO iri = 3*(itype-1)+1 , 3*(itype-1)+3
    !         write (*,*) iri
    !         DO icore = 1, ncore(itype)
    !           write (*,'(10f10.5)') (m_mcd(icore,l,iri,i),l=1,9)
    !         ENDDO
    !       ENDDO
    !      ENDDO
  END SUBROUTINE mcd_init
END MODULE m_mcdinit