coredr.F90 5.73 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 24

    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 ..
Daniel Wortmann's avatar
Daniel Wortmann committed
25 26 27
    REAL   , INTENT (IN) :: vrs(atoms%jmtd,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT (OUT) :: rhc(DIMENSION%msh,atoms%ntype,DIMENSION%jspd),qints(atoms%ntype,DIMENSION%jspd)
28 29 30 31 32 33 34
    !     ..
    !     .. 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
35 36 37 38
    REAL br(atoms%jmtd,atoms%ntype),brd(DIMENSION%msh),etab(100,atoms%ntype),&
         rhcs(atoms%jmtd,atoms%ntype,DIMENSION%jspd),rhochr(DIMENSION%msh),rhospn(DIMENSION%msh),&
         tecs(atoms%ntype,DIMENSION%jspd),vr(atoms%jmtd,atoms%ntype),vrd(DIMENSION%msh)
    INTEGER nkmust(atoms%ntype),ntab(100,atoms%ntype),ltab(100,atoms%ntype)
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

    !     ..
    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')
77
       CALL etabinit(atoms,DIMENSION,input, vr, etab,ntab,ltab,nkmust)
78 79 80 81 82 83 84 85 86 87 88 89
    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

90 91 92 93 94 95 96 97 98 99 100 101 102
       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
103 104
       IF (atoms%jri(jatom).LT.ncmsh) THEN
          DO i = atoms%jri(jatom) + 1,ncmsh
105 106 107 108 109 110 111 112 113 114
             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
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
          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)
    !
169 170
    !----> store core charge densities
    CALL writeCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
171 172
  END SUBROUTINE coredr
END MODULE m_coredr