broyden.F90 6.98 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 15
  SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                     mmap,nmaph,mapmt,mapvac2,nmap,fm,mit,sm,lpot)
16 17

#include"cpp_double.h"
18

19 20
    USE m_metric
    USE m_types
21

22
    IMPLICIT NONE
23

24 25 26 27 28 29 30 31 32
    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
33 34 35 36 37

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

    ! Array Arguments
    REAL,    INTENT (IN)    :: fm(nmap) 
42
    REAL,    INTENT (INOUT) :: sm(nmap)
43 44 45 46 47 48 49 50

    ! Local Scalars
    INTEGER         :: i,it,k,nit,npos,iread,nmaph
    REAL            :: bm,dfivi,fmvm,smnorm,vmnorm,alphan
    LOGICAL         :: l_pot
    REAL, PARAMETER :: one=1.0, zero=0.0

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

    ! External Functions
55 56
    REAL CPP_BLAS_sdot
    EXTERNAL CPP_BLAS_sdot
57 58

    ! External Subroutines
59
    EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal
60

61 62
    dfivi = zero

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

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

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

77
    IF (mit.NE.1) THEN
78 79
       ! load input charge density (sm1) and difference of 
       ! in and out charge densities (fm1) from previous iteration (m-1)
80 81

       REWIND 57
82 83
       READ (57) mit, alphan, (fm1(i),i=1,nmap), (sm1(i),i=1,nmap)
       IF (ABS(input%alpha-alphan) > 0.0001) THEN
84
          WRITE (6,*) 'mixing parameter has been changed; reset'
Daniel Wortmann's avatar
Daniel Wortmann committed
85
          WRITE (6,*) 'broyden algorithm or set alpha to',alphan
86 87 88 89 90 91 92
          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)
93
    END IF
94 95

    ! save F_m and rho_m for next iteration
96 97 98 99
    REWIND 57
    nit = mit +1
    IF (nit > input%maxiter+1) nit = 1
    WRITE (57) nit,input%alpha,fm,sm
100

101 102 103 104 105 106 107
    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
108 109
       DO k = 1, nmap
          um(k) = input%alpha * fm1(k) + sm1(k)
110
       END DO
111 112
       iread = MIN(mit-1,input%maxiter+1)
       DO it = 2, iread
Daniel Wortmann's avatar
Daniel Wortmann committed
113 114
          READ (59,rec=(it-1)*2-1) (ui(i),i=1,nmap)
          READ (59,rec=(it-1)*2)   (vi(i),i=1,nmap),dfivi
115
          am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)
116
          ! calculate um(:) = -am(it)*ui(:) + um
117 118 119
          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
120

121
       IF (input%imix.EQ.3) THEN
122 123 124 125 126 127 128 129 130
          !****************************************
          !        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>
131
          smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1)
132 133 134

          ! generate vm = alpha*sm1  - \sum <ui|w|sm1> vi
          vm(:) = input%alpha * fm1(:)
135
          DO it = 2,iread
136 137
             READ (59,rec=2*(it-1)-1) (ui(i),i=1,nmap)
             READ (59,rec=2*(it-1))   (vi(i),i=1,nmap), dfivi
138
             bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)
139
             ! calculate vm(:) = -bm*vi(:) + vm
140 141 142
             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
143 144 145

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

149
          CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
150

151
       ELSE IF (input%imix.EQ.5) THEN
152 153 154 155 156 157 158 159 160 161
          !****************************************
          !      broyden's second method
          !****************************************

          ! multiply fm1 with metric matrix and store in vm:  w |fm1>  
          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)
162 163 164
          CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1)

       ELSE IF (input%imix.EQ.7) THEN
165 166 167 168 169 170 171 172 173
          !****************************************
          !      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)

174
          DO it = 2,iread
175 176
             READ (59,rec=2*(it-1)) (vi(i),i=1,nmap), dfivi
             ! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
177 178
             CALL CPP_BLAS_saxpy(nmap,-am(it)*dfivi,vi,1,vm,1)
          END DO
179

180
          vmnorm = CPP_BLAS_sdot(nmap,fm1,1,vm,1)
181
          ! if (vmnorm.lt.tol_10) stop
182

183
          ! calculate vm(:) = (1.0/vmnorm)*vm(:)
184
          CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
185 186

          ! save dfivi(mit) for next iteration
187
          dfivi = vmnorm
188

189
       END IF
190 191 192

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

193
       npos=mit-1
194 195
       IF (mit.GT.input%maxiter+1) npos = MOD(mit-2,input%maxiter)+1

Daniel Wortmann's avatar
Daniel Wortmann committed
196
       WRITE (59,rec=2*npos-1) (um(i),i=1,nmap)
197 198 199 200
       WRITE (59,rec=2*npos)   (vm(i),i=1,nmap), dfivi

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

206
    mit = mit + 1
207
    DEALLOCATE (fm1,sm1,ui,um,vi,vm,am)
208 209 210

  END SUBROUTINE broyden
END MODULE m_broyden