u_mix.f90 6.71 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
MODULE m_umix
  USE m_juDFT
  !
  ! mix the old and new density matrix for the lda+U method
  !                                                 gb.2001
Daniel Wortmann's avatar
Daniel Wortmann committed
12 13
  ! --------------------------------------------------------
  ! Extension to multiple U per atom type by G.M. 2017
14
CONTAINS
15
  SUBROUTINE u_mix(input,atoms,n_mmp_in,n_mmp_out)
16

Daniel Wortmann's avatar
Daniel Wortmann committed
17
    USE m_types
18
    USE m_cdn_io
19
    USE m_nmat_rot
Daniel Wortmann's avatar
Daniel Wortmann committed
20 21
    USE m_xmlOutput

22
    ! ... Arguments
Daniel Wortmann's avatar
Daniel Wortmann committed
23

24
    IMPLICIT NONE
25
    TYPE(t_input),INTENT(IN)   :: input
26
    TYPE(t_atoms),INTENT(IN)   :: atoms
27 28
    COMPLEX, INTENT (INOUT)    :: n_mmp_out(-3:3,-3:3,atoms%n_u,input%jspins)
    COMPLEX, INTENT (INOUT)    :: n_mmp_in (-3:3,-3:3,atoms%n_u,input%jspins)
29 30
    !
    ! ... Locals ...
Daniel Wortmann's avatar
Daniel Wortmann committed
31 32 33
    INTEGER j,k,iofl,l,itype,ios,i_u,jsp,lty(atoms%n_u)
    REAL alpha,spinf,gam,del,sum1,sum2,mix_u, uParam, jParam
    REAL    theta(atoms%n_u),phi(atoms%n_u),zero(atoms%n_u)
34
    LOGICAL n_exist
Daniel Wortmann's avatar
Daniel Wortmann committed
35
    CHARACTER(LEN=20)   :: attributes(6)
36
    COMPLEX,ALLOCATABLE :: n_mmp(:,:,:,:)
37 38 39 40 41 42
    !
    ! check for possible rotation of n_mmp
    !
    INQUIRE (file='n_mmp_rot',exist=n_exist)
    IF (n_exist) THEN
       OPEN (68,file='n_mmp_rot',status='old',form='formatted')
Daniel Wortmann's avatar
Daniel Wortmann committed
43 44 45 46 47 48 49 50 51 52 53
       DO i_u = 1, atoms%n_u
          l = atoms%lda_u(i_u)%l
          READ(68,*,iostat=ios) theta(i_u),phi(i_u)
          IF (ios == 0) THEN
             lty(i_u) = l
          ELSE
             IF (i_u == 1)  CALL juDFT_error("ERROR reading n_mmp_rot", calledby ="u_mix")
             theta(i_u) = theta(i_u-1) ; phi(i_u) = phi(i_u-1)
             lty(i_u) = lty(i_u-1)
          END IF
       END DO
54 55
       CLOSE (68)
       zero = 0.0
56
       CALL nmat_rot(zero,-theta,-phi,3,atoms%n_u,input%jspins,lty,n_mmp_out)
Daniel Wortmann's avatar
Daniel Wortmann committed
57 58
    END IF

59
    ! Write out n_mmp_out to out.xml file
Daniel Wortmann's avatar
Daniel Wortmann committed
60 61

    CALL openXMLElementNoAttributes('ldaUDensityMatrix')
62
    DO jsp = 1, input%jspins
Daniel Wortmann's avatar
Daniel Wortmann committed
63 64 65 66 67 68 69 70 71 72 73 74 75 76
       DO i_u = 1, atoms%n_u
          l = atoms%lda_u(i_u)%l
          itype = atoms%lda_u(i_u)%atomType
          uParam = atoms%lda_u(i_u)%u
          jParam = atoms%lda_u(i_u)%j
          attributes = ''
          WRITE(attributes(1),'(i0)') jsp
          WRITE(attributes(2),'(i0)') itype
          WRITE(attributes(3),'(i0)') i_u
          WRITE(attributes(4),'(i0)') l
          WRITE(attributes(5),'(f15.8)') uParam
          WRITE(attributes(6),'(f15.8)') jParam
          CALL writeXMLElementMatrixPoly('densityMatrixFor',&
                                         (/'spin    ','atomType','uIndex  ','l       ','U       ','J       '/),&
77
                                         attributes,n_mmp_out(-l:l,-l:l,i_u,jsp))
Daniel Wortmann's avatar
Daniel Wortmann committed
78 79 80 81
       END DO
    END DO
    CALL closeXMLElement('ldaUDensityMatrix')

82 83
    ! Check whether density matrix is present and open density-matrix - file
    n_exist = isDensityMatrixPresent()
84 85 86
    OPEN (69,file='n_mmp_mat',status='unknown',form='formatted')

    IF (n_exist) THEN
87 88 89 90 91 92 93 94 95
       IF (input%ldauLinMix) THEN

          ! mix here straight with given mixing factors

          ALLOCATE (n_mmp(-3:3,-3:3,atoms%n_u,input%jspins))

          alpha = input%ldauMixParam
          spinf = input%ldauSpinf

96
          sum1 = 0.0
97
          IF (input%jspins.EQ.1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
98
             DO i_u = 1, atoms%n_u
99 100
                DO j = -3,3
                   DO k = -3,3
101 102
                      sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
                      n_mmp(k,j,i_u,1) = alpha * n_mmp_out(k,j,i_u,1) + (1.0-alpha) * n_mmp_in(k,j,i_u,1)
Daniel Wortmann's avatar
Daniel Wortmann committed
103 104 105
                   END DO
                END DO
             END DO
106 107 108
             WRITE (6,'(a16,f12.6)') 'n_mmp distance =',sum1
          ELSE
             sum2 = 0.0
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110 111
             gam = 0.5 * alpha * (1.0 + spinf)
             del = 0.5 * alpha * (1.0 - spinf)
             DO i_u = 1,atoms%n_u
112 113
                DO j = -3,3
                   DO k = -3,3
114 115 116 117 118 119 120 121 122 123 124 125
                      sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
                      sum2 = sum2 + ABS(n_mmp_out(k,j,i_u,2) - n_mmp_in(k,j,i_u,2))

                      n_mmp(k,j,i_u,1) =       gam * n_mmp_out(k,j,i_u,1) + &
                                         (1.0-gam) * n_mmp_in (k,j,i_u,1) + &
                                               del * n_mmp_out(k,j,i_u,2) - &
                                               del * n_mmp_in (k,j,i_u,2)

                      n_mmp(k,j,i_u,2) =       gam * n_mmp_out(k,j,i_u,2) + &
                                         (1.0-gam) * n_mmp_in (k,j,i_u,2) + &
                                               del * n_mmp_out(k,j,i_u,1) - &
                                               del * n_mmp_in (k,j,i_u,1)
Daniel Wortmann's avatar
Daniel Wortmann committed
126 127 128
                   END DO
                END DO
             END DO
129 130 131 132 133
             WRITE (6,'(a23,f12.6)') 'n_mmp distance spin 1 =',sum1
             WRITE (6,'(a23,f12.6)') 'n_mmp distance spin 2 =',sum2
          ENDIF
          WRITE (69,9000) n_mmp
          WRITE (69,'(2(a6,f5.3))') 'alpha=',alpha,'spinf=',spinf
134
          n_mmp_in = n_mmp
135 136
          DEALLOCATE (n_mmp)
       ELSE ! input%ldauLinMix
137 138

          ! calculate distance and write new n_mmp to mix in broyden.F
139

140
          sum1 = 0.0
Daniel Wortmann's avatar
Daniel Wortmann committed
141
          DO i_u = 1, atoms%n_u
142 143
             DO j = -3,3
                DO k = -3,3
144
                   sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
Daniel Wortmann's avatar
Daniel Wortmann committed
145 146 147
                END DO
             END DO
          END DO
148
          IF (input%jspins.EQ.1) THEN
149 150 151 152
             WRITE (6,'(a16,f12.6)') 'n_mmp distance =',sum1
          ELSE
             sum2 = 0.0
             WRITE (6,'(a23,f12.6)') 'n_mmp distance spin 1 =',sum1
Daniel Wortmann's avatar
Daniel Wortmann committed
153
             DO i_u = 1, atoms%n_u
154 155
                DO j = -3,3
                   DO k = -3,3
156
                      sum2 = sum2 + ABS(n_mmp_out(k,j,i_u,2) - n_mmp_in(k,j,i_u,2))
Daniel Wortmann's avatar
Daniel Wortmann committed
157 158 159
                   END DO
                END DO
             END DO
160
             DO j=-3,3
161
                WRITE(6,'(14f12.6)') (n_mmp_in(k,j,1,2),k=-3,3)
Daniel Wortmann's avatar
Daniel Wortmann committed
162
             END DO
163 164
             WRITE (6,'(a23,f12.6)') 'n_mmp distance spin 2 =',sum2
             DO j=-3,3
165
                WRITE(6,'(14f12.6)') (n_mmp_out(k,j,1,2),k=-3,3)
Daniel Wortmann's avatar
Daniel Wortmann committed
166 167
             END DO
          END IF
168 169
          WRITE (69,9000) n_mmp_in
          WRITE (69,9000) n_mmp_out
170 171
       END IF ! input%ldauLinMix
    ELSE ! isDensityMatrixPresent()
172
       ! first time with lda+u; write new n_mmp  
173
       WRITE (69,9000) n_mmp_out
174
       WRITE (69,'(2(a6,f5.3))') 'alpha=',input%ldauMixParam,'spinf=',input%ldauSpinf
175
       n_mmp_in = n_mmp_out
176
    END IF ! isDensityMatrixPresent()
177 178 179 180 181 182

9000 FORMAT(7f20.13)

    CLOSE (69)
  END SUBROUTINE u_mix
END MODULE m_umix