cdnsp.f90 7.13 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17
      MODULE m_cdnsp
      USE m_juDFT
!     *******************************************************
!     sets up the starting density for the spin-polarized
!     calculation from a paramagnetic density
!     changed to suit both ferromagnetic and antiferro-
!     magnetic case. changes only in mt-part - r.pentcheva Jan'96
!     *******************************************************
      CONTAINS
        SUBROUTINE cdnsp(&
             &                 atoms,input,vacuum,sphhar,&
18
             &                 stars,sym,noco,oneD,cell,DIMENSION)
19 20

          USE m_intgr, ONLY : intgr3
21
          USE m_constants
22
          USE m_cdn_io
23 24 25
          USE m_types
          IMPLICIT NONE
          !     ..
26 27 28 29
          TYPE(t_stars),INTENT(IN)     :: stars
          TYPE(t_vacuum),INTENT(IN)    :: vacuum
          TYPE(t_atoms),INTENT(IN)     :: atoms
          TYPE(t_sphhar),INTENT(IN)    :: sphhar
30
          TYPE(t_input),INTENT(INOUT)  :: input
31
          TYPE(t_sym),INTENT(IN)       :: sym
32
          TYPE(t_noco),INTENT(IN)      :: noco
33 34 35
          TYPE(t_oneD),INTENT(IN)      :: oneD
          TYPE(t_cell),INTENT(IN)      :: cell
          TYPE(t_dimension),INTENT(IN) :: DIMENSION
36 37 38 39

          ! local type instances
          TYPE(t_potden)               :: den

40
          !     .. Local Scalars ..
41
          REAL dummy,p,pp,qtot1,qtot2,spmtot,qval,sfp,fermiEnergyTemp
42
          INTEGER i,ivac,j,k,lh,n,na,jsp_new
43
          INTEGER ios 
44
          LOGICAL n_exist,l_qfix
45 46
          !     ..
          !     .. Local Arrays ..
47 48
          REAL rhoc(atoms%jmtd,atoms%ntype,dimension%jspd)
          REAL tec(atoms%ntype,dimension%jspd),qintc(atoms%ntype,dimension%jspd)
49 50 51 52 53 54 55 56 57
          CHARACTER(len=140), ALLOCATABLE :: clines(:)
          CHARACTER(len=140)              :: lineread
          !      ..
          sfp = 2 * SQRT(pi_const)
          !sphhar%nlhd = MAXVAL(sphhar%nlh(:))

          IF (input%jspins/=2)  CALL juDFT_error&
               &     ("cdnsp: set jspins = 2 and remove fl7para!",calledby&
               &     ="cdnsp")
58

59
          CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
60 61 62
          ALLOCATE (den%cdom(1),den%cdomvz(1,1),den%cdomvxy(1,1,1))
          ALLOCATE (den%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,atoms%n_u),input%jspins))
          den%mmpMat = CMPLX(0.0,0.0)
63 64

          input%jspins=1
65
          CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
66
          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
67
                           CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
68
          input%jspins=2
69

70 71 72 73 74 75 76 77
          qval = 0.
          na = 1
          !
          !     ---> set jspins=2
          jsp_new = 2
          !
          DO n = 1,atoms%ntype
             DO j = 1,atoms%jri(n)
78
                den%mt(j,0,n,1) = den%mt(j,0,n,1) - rhoc(j,n,1)/sfp
79
             ENDDO
80 81
             !         WRITE (16,FMT='(8f10.4)') (den%mt(i,0,n,1),i=1,16)
             CALL intgr3(den%mt(1,0,n,1),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qval)
82 83 84
             p = (atoms%bmu(n)+sfp*qval)/ (2.*sfp*qval)
             pp = 1. - p
             DO j = 1,atoms%jri(n)
85 86
                den%mt(j,0,n,jsp_new) = pp*den%mt(j,0,n,1) + rhoc(j,n,1)/ (2.*sfp)
                den%mt(j,0,n,1)       =  p*den%mt(j,0,n,1) + rhoc(j,n,1)/ (2.*sfp)
87 88 89
             ENDDO
             DO lh = 1,sphhar%nlh(atoms%ntypsy(na))
                DO j = 1,atoms%jri(n)
90 91
                   den%mt(j,lh,n,jsp_new) = pp*den%mt(j,lh,n,1)
                   den%mt(j,lh,n,1)       =  p*den%mt(j,lh,n,1)
92 93 94 95 96
                ENDDO
             ENDDO
             na = na + atoms%neq(n)
          ENDDO
          DO k = 1,stars%ng3
97 98
             den%pw(k,jsp_new) = 0.5 * den%pw(k,1)
             den%pw(k,1)       = den%pw(k,jsp_new)
99 100 101 102
          ENDDO
          IF (input%film) THEN
             DO ivac = 1,vacuum%nvac
                DO j = 1, vacuum%nmz
103 104
                   den%vacz(j,ivac,jsp_new) = 0.5 * den%vacz(j,ivac,1)
                   den%vacz(j,ivac,1)       = den%vacz(j,ivac,jsp_new)
105 106 107
                ENDDO
                DO k = 2, stars%ng2
                   DO j = 1,vacuum%nmzxy
108 109
                      den%vacxy(j,k-1,ivac,jsp_new) = 0.5 * den%vacxy(j,k-1,ivac,1)
                      den%vacxy(j,k-1,ivac,1)       = den%vacxy(j,k-1,ivac,jsp_new)
110 111 112 113
                   ENDDO
                ENDDO
             ENDDO
          ENDIF
114
          !     ----> write the spin-polarized density
115
          CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
116
                            CDN_INPUT_DEN_const,0,-1.0,0.0,.FALSE.,den)
117 118 119 120 121 122 123
          !
          !     -----> This part is only used for testing th e magnetic moment in 
          !     ----->   each sphere
          !
          DO n = 1,atoms%ntype
             qtot1=0.00
             qtot2=0.00
124 125
             CALL intgr3(den%mt(1,0,n,1),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qtot1)
             CALL intgr3(den%mt(1,0,n,jsp_new),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qtot2)
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
             spmtot=sfp*(qtot1-qtot2)
             WRITE (6,'('' moment in sphere '',2x,'':'',f8.4)') spmtot
          ENDDO

          !--->   read enpara and then double it
          INQUIRE(file='enpara',exist=n_exist)
          IF (n_exist) THEN
             OPEN(40,file ='enpara',status='old',form='formatted')
             REWIND 40
             n = 0
             DO
                READ (40,'(a)',iostat = ios) lineread
                IF (ios/=0) EXIT
                n          = n+1
             ENDDO

             ALLOCATE (clines(n))

             REWIND 40
             DO i = 1,n
                READ (40,'(a)') clines(i)
             ENDDO

             REWIND 40
             DO i = 1,n
                WRITE (40,'(a)') TRIM(clines(i))
             ENDDO
             DO i = 1,n
                WRITE (40,'(a)') TRIM(clines(i))
             ENDDO

157
             DEALLOCATE (clines)
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
             CLOSE(40)
          ENDIF
          !
          ! for lda+U: flip n-matrix
          !
          IF (atoms%n_u.GT.0) THEN
             INQUIRE (file='n_mmp_mat',exist=n_exist)
             IF (n_exist) THEN
                OPEN (69,file='n_mmp_mat',status='old',form='formatted')
                REWIND 69

                n=0
                DO
                   READ (69,'(a)',iostat=ios) lineread
                   IF (ios.NE.0) EXIT
                   n = n+1
                ENDDO
                ALLOCATE (clines(n))
                REWIND 69
                DO i=1,n
                   WRITE (69,'(a)') TRIM(clines(i))
                ENDDO
                DO i=1,n
                   WRITE (69,'(a)') TRIM(clines(i))
                ENDDO
                DEALLOCATE (clines)

                CLOSE(69)
             ENDIF
          ENDIF


        END SUBROUTINE cdnsp
      END MODULE m_cdnsp