coredr.F90 5.68 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
    REAL   , INTENT (IN) :: vrs(atoms%jmtd,atoms%ntype,input%jspins)
    REAL,    INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
26
    REAL,    INTENT (OUT) :: rhc(atoms%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 ..
34 35 36
    REAL br(atoms%jmtd,atoms%ntype),brd(atoms%msh),etab(100,atoms%ntype),&
         rhcs(atoms%jmtd,atoms%ntype,input%jspins),rhochr(atoms%msh),rhospn(atoms%msh),&
         tecs(atoms%ntype,input%jspins),vr(atoms%jmtd,atoms%ntype),vrd(atoms%msh)
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
    END IF
    !
79
    ncmsh = atoms%msh
80 81 82 83 84 85 86 87 88
    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
          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)

128
       CALL spratm(atoms%msh,vrd,brd,z,rnot,dxx,ncmsh,&
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
            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
147
          DO j = 1,atoms%msh
148 149 150 151
             rhc(j,jatom,input%jspins) = (rhochr(j)+rhospn(j))*0.5
             rhc(j,jatom,1) = (rhochr(j)-rhospn(j))*0.5
          ENDDO
       ELSE
152
          DO j = 1,atoms%msh
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
             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