coredr.F90 5.72 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
MODULE m_coredr
CONTAINS
  SUBROUTINE coredr(input,atoms,seig, rho,DIMENSION,sphhar, vrs, qints,rhc)
    !     *******************************************************
    !     *****   set up the core densities for compounds   *****
    !     *****   for relativistic core                     *****
    !     *******************************************************

    USE m_etabinit
    USE m_spratm
    USE m_ccdnup
12
    USE m_cdn_io
13 14 15 16 17 18 19 20 21 22 23
    USE m_types
    IMPLICIT NONE
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_atoms),INTENT(IN)     :: atoms
    !
    !     .. Scalar Arguments ..
    REAL seig
    !     ..
    !     .. Array Arguments ..
24 25 26
    REAL   , INTENT (IN) :: vrs(atoms%jmtd,atoms%ntype,input%jspins)
    REAL,    INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
    REAL,    INTENT (OUT) :: rhc(DIMENSION%msh,atoms%ntype,input%jspins),qints(atoms%ntype,input%jspins)
27 28 29 30 31 32 33
    !     ..
    !     .. Local Scalars ..
    REAL dxx,rnot,sume,t2,t2b,z,t1,rr,d,v1,v2
    INTEGER i,j,jatom,jspin,k,n,ncmsh
    LOGICAL exetab
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
34
    REAL br(atoms%jmtd,atoms%ntype),brd(DIMENSION%msh),etab(100,atoms%ntype),&
35 36
         rhcs(atoms%jmtd,atoms%ntype,input%jspins),rhochr(DIMENSION%msh),rhospn(DIMENSION%msh),&
         tecs(atoms%ntype,input%jspins),vr(atoms%jmtd,atoms%ntype),vrd(DIMENSION%msh)
Daniel Wortmann's avatar
Daniel Wortmann committed
37
    INTEGER nkmust(atoms%ntype),ntab(100,atoms%ntype),ltab(100,atoms%ntype)
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

    !     ..
    ntab(:,:) = -1 ; ltab(:,:) = -1 ; etab(:,:) = 0.0
    !
    ! setup potential and field
    !
    IF (input%jspins.EQ.1) THEN
       DO n = 1,atoms%ntype
          DO j = 1,atoms%jmtd
             vr(j,n) = vrs(j,n,1)
             br(j,n) = 0.0
          END DO
       END DO
    ELSE
       DO n = 1,atoms%ntype
          DO j = 1,atoms%jmtd
             vr(j,n) = (vrs(j,n,1)+vrs(j,n,input%jspins))/2.
             br(j,n) = (vrs(j,n,input%jspins)-vrs(j,n,1))/2.
          END DO
       END DO
    END IF
    !
    ! setup eigenvalues
    exetab = .FALSE.
    INQUIRE (file='core.dat',exist=exetab)
    IF (exetab) THEN
       OPEN (58,file='core.dat',form='formatted',status='old')
       REWIND 58
       DO n = 1,atoms%ntype
          READ (58,FMT=*) nkmust(n)
          DO k = 1,nkmust(n)
             READ (58,FMT='(f12.6,2i3)') etab(k,n),ntab(k,n),&
                  &                                               ltab(k,n)

          END DO
       END DO
    ELSE
       OPEN (58,file='core.dat',form='formatted',status='new')
76
       CALL etabinit(atoms,DIMENSION,input, vr, etab,ntab,ltab,nkmust)
77 78 79 80 81 82 83 84 85 86 87 88
    END IF
    !
    ncmsh = DIMENSION%msh
    seig = 0.
    ! ---> set up densities
    DO jatom = 1,atoms%ntype
       !
       DO j = 1,atoms%jri(jatom)
          vrd(j) = vr(j,jatom)
          brd(j) = br(j,jatom)
       END DO

89 90 91 92 93 94 95 96 97 98 99 100 101
       IF (input%l_core_confpot) THEN
          !--->    linear extension of the potential with slope t1 / a.u.
          rr = atoms%rmt(jatom)
          d = EXP(atoms%dx(jatom))
          t1=0.125
          !         t2  = vrd(jri(jatom))/rr - rr*t1
          !         t2b = brd(jri(jatom))/rr - rr*t1
          t2  = vrs(atoms%jri(jatom),jatom,1)     /rr - rr*t1
          t2b = vrs(atoms%jri(jatom),jatom,input%jspins)/rr - rr*t1
       ELSE
          t2 = vrd(atoms%jri(jatom))/ (atoms%jri(jatom)-ncmsh)
          t2b = brd(atoms%jri(jatom))/ (atoms%jri(jatom)-ncmsh)
       ENDIF
102 103
       IF (atoms%jri(jatom).LT.ncmsh) THEN
          DO i = atoms%jri(jatom) + 1,ncmsh
104 105 106 107 108 109 110 111 112 113
             IF (input%l_core_confpot) THEN
                rr = d*rr
                v1 = rr*( t2  + rr*t1 )
                v2 = rr*( t2b + rr*t1 )
                vrd(i) = 0.5*(v2 + v1)
                brd(i) = 0.5*(v2 - v1)
             ELSE
                vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
                brd(i) = brd(atoms%jri(jatom)) + t2b* (i-atoms%jri(jatom))
             ENDIF
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
          END DO
       END IF

       !        rr = rmsh(1,jatom)
       !        do i =1, ncmsh
       !          rr = d*rr
       !         write(*,'(3f20.10)') rr,vrd(i),brd(i)
       !        enddo

       !
       rnot = atoms%rmsh(1,jatom)
       z = atoms%zatom(jatom)
       dxx = atoms%dx(jatom)

       CALL spratm(DIMENSION%msh,vrd,brd,z,rnot,dxx,ncmsh,&
            etab(1,jatom),ntab(1,jatom),ltab(1,jatom), sume,rhochr,rhospn)

       seig = seig + atoms%neq(jatom)*sume
       !
       !     rho_up=2(ir) = (rhochr(ir)  + rhospn(ir))*0.5
       !     rho_dw=1(ir) = (rhochr(ir)  - rhospn(ir))*0.5
       !
       IF (input%jspins.EQ.2) THEN
          DO j = 1,atoms%jri(jatom)
             rhcs(j,jatom,input%jspins) = (rhochr(j)+rhospn(j))*0.5
             rhcs(j,jatom,1) = (rhochr(j)-rhospn(j))*0.5
          END DO
       ELSE
          DO j = 1,atoms%jri(jatom)
             rhcs(j,jatom,1) = rhochr(j)
          END DO
       END IF
       IF (input%jspins.EQ.2) THEN
          DO j = 1,DIMENSION%msh
             rhc(j,jatom,input%jspins) = (rhochr(j)+rhospn(j))*0.5
             rhc(j,jatom,1) = (rhochr(j)-rhospn(j))*0.5
          ENDDO
       ELSE
          DO j = 1,DIMENSION%msh
             rhc(j,jatom,1) = rhochr(j)
          END DO
       END IF
       !
       ! store atomic eigenvalues to file.58
       IF (jatom.EQ.1) REWIND 58
       WRITE (58,FMT=*) nkmust(jatom)
       DO k = 1,nkmust(jatom)
          WRITE (58,FMT='(f12.6,2i3)') etab(k,jatom),ntab(k,jatom), ltab(k,jatom)
       END DO
       !---->update spherical charge density rho with the core density.
       CALL ccdnup(atoms,sphhar,input,jatom, rho, sume,vrs,rhochr,rhospn, tecs,qints)

    END DO ! loop over atoms (jatom)
    !
168 169
    !----> store core charge densities
    CALL writeCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
170 171
  END SUBROUTINE coredr
END MODULE m_coredr