mcd_init.f90 6.06 KB
Newer Older
1 2
MODULE m_mcdinit
CONTAINS
3
  SUBROUTINE mcd_init(atoms,input,DIMENSION,vr,g,f,mcd,itype,jspin)
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

    !-----------------------------------------------------------------------
    !
    ! 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
24
    TYPE(t_mcd),INTENT(INOUT)    :: mcd
25 26 27 28 29 30 31

    INTEGER, PARAMETER :: l_max = 3

    ! Arguments ...

    INTEGER, INTENT (IN)  :: itype
    INTEGER, INTENT (IN)  :: jspin
Daniel Wortmann's avatar
Daniel Wortmann committed
32
    REAL,    INTENT (IN)  :: vr(atoms%jmtd,atoms%ntype,input%jspins)
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
    REAL,    INTENT (IN)  :: f(atoms%jmtd,2,0:atoms%lmaxd,jspin:jspin)
    REAL,    INTENT (IN)  :: g(atoms%jmtd,2,0:atoms%lmaxd,jspin:jspin)

    ! 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

54
    mcd%ncore(itype) = 0
55
    bmu = 0.0
56
    CALL setcor(itype,1,atoms,input,bmu, nst,kappa,nprnc,occ)
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

    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")
84
             IF ( (e.LE.mcd%emcd_up).AND.(e.GE.mcd%emcd_lo) ) THEN
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
                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
111
                   mcd%m_mcd(icore,l,iri,i) = CMPLX(0.0,0.0)
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
                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
Daniel Wortmann's avatar
Daniel Wortmann committed
142
                CALL nabla(itype,icore,atoms%jri(itype),atoms%dx(itype),DIMENSION%nstd,atoms%ntype,&
143
                     j_core(icore),l_core(icore),l_max,ms,atoms%rmsh(:,itype),gc(:,icore,ispin),&
144
                     gv(:,0:,ispin,i),dgv(:,0:,ispin,i), mcd%m_mcd(:,:,:,i) )
145 146 147
             ENDDO

             DO i = 1, 2*icore*l_core(icore)
148 149 150
                mcd%ncore(itype) = mcd%ncore(itype) + 1
                IF (mcd%ncore(itype)>DIMENSION%nstd)  CALL juDFT_error("dimension%nstd too small" ,calledby ="mcd_init")
                mcd%e_mcd(itype,ispin,mcd%ncore(itype)) = e_mcd1(icore)
151 152 153 154 155 156 157 158 159 160 161 162
             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
163 164
    !         DO icore = 1, mcd%ncore(itype)
    !           write (*,'(10f10.5)') (mcd%m_mcd(icore,l,iri,i),l=1,9)
165 166 167 168 169
    !         ENDDO
    !       ENDDO
    !      ENDDO
  END SUBROUTINE mcd_init
END MODULE m_mcdinit