broyden.F90 7.07 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
MODULE m_broyden
  USE m_juDFT
  !################################################################
  !     IMIX = 3 : BROYDEN'S FIRST METHOD
  !     IMIX = 5 : BROYDEN'S SECOND METHOD
  !     IMIX = 7 : GENERALIZED ANDERSEN METHOD
  !     sm   : input charge density of iteration m
  !            afterwards update rho(m+1)
  !     fm   : output minus input charge density of iteration m
  !     sm1  : input charge density of iteration m-1
  !     fm1   : output minus inputcharge density of iteration m-1
  !################################################################
CONTAINS
14
  SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
15
                     hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,lpot)
16 17

#include"cpp_double.h"
18

19 20
    USE m_metric
    USE m_types
21
    USE m_broyd_io
22

23
    IMPLICIT NONE
24

25 26 27 28 29 30 31 32 33
    TYPE(t_oneD),INTENT(IN)    :: oneD
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_noco),INTENT(IN)    :: noco
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_stars),INTENT(IN)   :: stars
    TYPE(t_cell),INTENT(IN)    :: cell
    TYPE(t_sphhar),INTENT(IN)  :: sphhar
    TYPE(t_atoms),INTENT(IN)   :: atoms
34
    TYPE(t_hybrid),INTENT(IN)  :: hybrid
35 36 37

    ! Scalar Arguments
    INTEGER, INTENT (IN)        :: mmap,nmap
38
    INTEGER, INTENT (IN)        :: mapmt,mapvac2
39
    LOGICAL,OPTIONAL,INTENT(IN) :: lpot
40 41 42

    ! Array Arguments
    REAL,    INTENT (IN)    :: fm(nmap) 
43
    REAL,    INTENT (INOUT) :: sm(nmap)
44 45

    ! Local Scalars
46
    INTEGER         :: i,it,k,nit,iread,nmaph, mit
47
    REAL            :: bm,dfivi,fmvm,smnorm,vmnorm,alphan
48
    LOGICAL         :: l_pot, l_exist
49 50 51
    REAL, PARAMETER :: one=1.0, zero=0.0

    ! Local Arrays
52 53
    REAL, ALLOCATABLE :: am(:)
    REAL, ALLOCATABLE :: fm1(:),sm1(:),ui(:),um(:),vi(:),vm(:)
54 55

    ! External Functions
56 57
    REAL CPP_BLAS_sdot
    EXTERNAL CPP_BLAS_sdot
58 59

    ! External Subroutines
60
    EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal
61

62 63
    dfivi = zero

64 65
    l_pot = .FALSE.
    IF (PRESENT(lpot)) l_pot = lpot
66 67

    ALLOCATE (fm1(mmap),sm1(mmap),ui(mmap),um(mmap),vi(mmap),vm(mmap))
68
    ALLOCATE (am(input%maxiter+1))
69 70 71 72 73 74 75 76 77

    fm1 = 0.0
    sm1 = 0.0
    ui  = 0.0
    um  = 0.0
    vi  = 0.0
    vm  = 0.0
    am  = 0.0

78 79 80 81
    mit = 0
    l_exist = initBroydenHistory(input,hybrid,nmap) ! returns true if there already exists a Broyden history
    IF(.NOT.l_exist) mit = 1

82
    IF (mit.NE.1) THEN
83 84
       ! load input charge density (sm1) and difference of 
       ! in and out charge densities (fm1) from previous iteration (m-1)
85

86
       CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))
87
       IF (ABS(input%alpha-alphan) > 0.0001) THEN
88
          WRITE (6,*) 'mixing parameter has been changed; reset'
Daniel Wortmann's avatar
Daniel Wortmann committed
89
          WRITE (6,*) 'broyden algorithm or set alpha to',alphan
90 91 92 93 94 95 96
          CALL juDFT_error("mixing parameter (input) changed", calledby ="broyden")
       END IF

       ! generate F_m   - F_(m-1)  ... sm1
       !      and rho_m - rho_(m-1) .. fm1
       sm1(1:nmap) = sm(1:nmap) - sm1(1:nmap)
       fm1(1:nmap) = fm(1:nmap) - fm1(1:nmap)
97
    END IF
98 99

    ! save F_m and rho_m for next iteration
100 101
    nit = mit +1
    IF (nit > input%maxiter+1) nit = 1
102
    CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm)
103

104 105 106 107 108 109 110
    IF (mit.EQ.1) THEN 
       !     update for rho for mit=1 is straight mixing
       !     sm = sm + alpha*fm
       CALL CPP_BLAS_saxpy(nmap,input%alpha,fm,1,sm,1)
    ELSE
       !     |vi> = w|vi> 
       !     loop to generate um : um = sm1 + alpha*fm1 - \sum <fm1|w|vi> ui
111
       um(:nmap) = input%alpha * fm1(:nmap) + sm1(:nmap)
112

113 114
       iread = MIN(mit-1,input%maxiter+1)
       DO it = 2, iread
115 116 117
          CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
          CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)

118
          am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)
119
          ! calculate um(:) = -am(it)*ui(:) + um(:)
120 121 122
          CALL CPP_BLAS_saxpy(nmap,-am(it),ui,1,um,1)
          WRITE(6,FMT='(5x,"<vi|w|Fm> for it",i2,5x,f10.6)')it,am(it) 
       END DO
123

124
       IF (input%imix.EQ.3) THEN
125 126 127 128 129 130 131 132 133
          !****************************************
          !        broyden's first method
          !****************************************

          ! convolute drho(m) with the metric: |fm1> = w|sm1>
          CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
                      mmap,nmaph,mapmt,mapvac2,sm1,fm1,l_pot)

          ! calculate the norm of sm1 : <sm1|w|sm1>
134
          smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1)
135 136 137

          ! generate vm = alpha*sm1  - \sum <ui|w|sm1> vi
          vm(:) = input%alpha * fm1(:)
138

139
          DO it = 2,iread
140 141
             CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
             CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
142
             bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)
143
             ! calculate vm(:) = -bm*vi(:) + vm
144 145 146
             CALL CPP_BLAS_saxpy(nmap,-bm,vi,1,vm,1)
             !write(6,FMT='(5x,"<ui|w|Fm> for it",i2,5x,f10.6)') it, bm 
          END DO
147 148 149

          ! complete evaluation of vm
          ! vmnorm = <um|w|sm1>-<sm1|w|sm1>
150
          vmnorm = CPP_BLAS_sdot(nmap,fm1,1,um,1) - smnorm
151 152
          ! if (vmnorm.lt.tol_10) stop

153
          CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
154

155
       ELSE IF (input%imix.EQ.5) THEN
156 157 158 159
          !****************************************
          !      broyden's second method
          !****************************************

160
          ! multiply fm1 with metric matrix and store in vm:  w |fm1>
161 162 163 164 165
          CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
                      mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)

          ! calculate the norm of fm1 and normalize vm it: vm = wfm1 / <fm1|w|fm1>
          vmnorm = one / CPP_BLAS_sdot(nmap,fm1,1,vm,1)
166 167 168
          CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1)

       ELSE IF (input%imix.EQ.7) THEN
169 170 171 172 173 174 175 176 177
          !****************************************
          !      generalized anderson method
          !****************************************

          ! calculate vm = alpha*wfm1 -\sum <fm1|w|vi> <fi1|w|vi><vi|
          ! convolute fm1 with the metrik and store in vm
          CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
                      mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)

178
          DO it = 2,iread
179
             CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
180
             ! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
181 182
             CALL CPP_BLAS_saxpy(nmap,-am(it)*dfivi,vi,1,vm,1)
          END DO
183

184
          vmnorm = CPP_BLAS_sdot(nmap,fm1,1,vm,1)
185
          ! if (vmnorm.lt.tol_10) stop
186

187
          ! calculate vm(:) = (1.0/vmnorm)*vm(:)
188
          CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
189 190

          ! save dfivi(mit) for next iteration
191
          dfivi = vmnorm
192

193
       END IF
194 195 196

       ! write um,vm and dfivi on file broyd.?

197 198
       CALL writeUVec(input,hybrid,nmap,mit,um)
       CALL writeVVec(input,hybrid,nmap,mit,dfivi,vm)
199 200 201

       ! update rho(m+1)
       ! calculate <fm|w|vm>
202
       fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1)
203
       ! calculate sm(:) = (1.0-fmvm)*um(:) + sm
204 205
       CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1)
    END IF
206 207

    DEALLOCATE (fm1,sm1,ui,um,vi,vm,am)
208 209 210

  END SUBROUTINE broyden
END MODULE m_broyden