broyden.F90 6.98 KB
 Markus Betzinger committed Apr 26, 2016 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  Gregor Michalicek committed Nov 08, 2017 14 15  SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,& mmap,nmaph,mapmt,mapvac2,nmap,fm,mit,sm,lpot)  Markus Betzinger committed Apr 26, 2016 16 17  #include"cpp_double.h"  Gregor Michalicek committed Nov 08, 2017 18   Markus Betzinger committed Apr 26, 2016 19 20  USE m_metric USE m_types  Gregor Michalicek committed Nov 08, 2017 21   Markus Betzinger committed Apr 26, 2016 22  IMPLICIT NONE  Gregor Michalicek committed Nov 08, 2017 23   Markus Betzinger committed Apr 26, 2016 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  Gregor Michalicek committed Nov 08, 2017 33 34 35 36 37  ! Scalar Arguments INTEGER, INTENT (IN) :: mmap,nmap INTEGER, INTENT (IN) :: mapmt,mapvac2 INTEGER, INTENT (INOUT) :: mit  Markus Betzinger committed Apr 26, 2016 38  LOGICAL,OPTIONAL,INTENT(IN) :: lpot  Gregor Michalicek committed Nov 08, 2017 39 40 41  ! Array Arguments REAL, INTENT (IN) :: fm(nmap)  Markus Betzinger committed Apr 26, 2016 42  REAL, INTENT (INOUT) :: sm(nmap)  Gregor Michalicek committed Nov 08, 2017 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  Markus Betzinger committed Apr 26, 2016 51 52  REAL, ALLOCATABLE :: am(:) REAL, ALLOCATABLE :: fm1(:),sm1(:),ui(:),um(:),vi(:),vm(:)  Gregor Michalicek committed Nov 08, 2017 53 54  ! External Functions  Markus Betzinger committed Apr 26, 2016 55 56  REAL CPP_BLAS_sdot EXTERNAL CPP_BLAS_sdot  Gregor Michalicek committed Nov 08, 2017 57 58  ! External Subroutines  Markus Betzinger committed Apr 26, 2016 59  EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal  Gregor Michalicek committed Nov 08, 2017 60   Markus Betzinger committed Apr 26, 2016 61 62  dfivi = zero  Gregor Michalicek committed Nov 08, 2017 63 64  l_pot = .FALSE. IF (PRESENT(lpot)) l_pot = lpot  Markus Betzinger committed Apr 26, 2016 65 66  ALLOCATE (fm1(mmap),sm1(mmap),ui(mmap),um(mmap),vi(mmap),vm(mmap))  Gregor Michalicek committed Nov 08, 2017 67  ALLOCATE (am(input%maxiter+1))  Gregor Michalicek committed May 03, 2016 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  Markus Betzinger committed Apr 26, 2016 77  IF (mit.NE.1) THEN  Gregor Michalicek committed Nov 08, 2017 78 79  ! load input charge density (sm1) and difference of ! in and out charge densities (fm1) from previous iteration (m-1)  Markus Betzinger committed Apr 26, 2016 80 81  REWIND 57  Gregor Michalicek committed Nov 08, 2017 82 83  READ (57) mit, alphan, (fm1(i),i=1,nmap), (sm1(i),i=1,nmap) IF (ABS(input%alpha-alphan) > 0.0001) THEN  Markus Betzinger committed Apr 26, 2016 84  WRITE (6,*) 'mixing parameter has been changed; reset'  Daniel Wortmann committed Apr 29, 2016 85  WRITE (6,*) 'broyden algorithm or set alpha to',alphan  Gregor Michalicek committed Nov 08, 2017 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)  Markus Betzinger committed Apr 26, 2016 93  END IF  Gregor Michalicek committed Nov 08, 2017 94 95  ! save F_m and rho_m for next iteration  Markus Betzinger committed Apr 26, 2016 96 97 98 99  REWIND 57 nit = mit +1 IF (nit > input%maxiter+1) nit = 1 WRITE (57) nit,input%alpha,fm,sm  Gregor Michalicek committed Nov 08, 2017 100   Markus Betzinger committed Apr 26, 2016 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 ui  Gregor Michalicek committed Nov 08, 2017 108 109  DO k = 1, nmap um(k) = input%alpha * fm1(k) + sm1(k)  Markus Betzinger committed Apr 26, 2016 110  END DO  Gregor Michalicek committed Nov 08, 2017 111 112  iread = MIN(mit-1,input%maxiter+1) DO it = 2, iread  Daniel Wortmann committed Jul 11, 2017 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  Markus Betzinger committed Apr 26, 2016 115  am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)  Gregor Michalicek committed Nov 10, 2017 116  ! calculate um(:) = -am(it)*ui(:) + um  Markus Betzinger committed Apr 26, 2016 117 118 119  CALL CPP_BLAS_saxpy(nmap,-am(it),ui,1,um,1) WRITE(6,FMT='(5x," for it",i2,5x,f10.6)')it,am(it) END DO  Gregor Michalicek committed Nov 08, 2017 120   Markus Betzinger committed Apr 26, 2016 121  IF (input%imix.EQ.3) THEN  Gregor Michalicek committed Nov 08, 2017 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 :  Markus Betzinger committed Apr 26, 2016 131  smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1)  Gregor Michalicek committed Nov 08, 2017 132 133 134  ! generate vm = alpha*sm1 - \sum vi vm(:) = input%alpha * fm1(:)  Markus Betzinger committed Apr 26, 2016 135  DO it = 2,iread  Gregor Michalicek committed Nov 08, 2017 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  Markus Betzinger committed Apr 26, 2016 138  bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)  Gregor Michalicek committed Nov 10, 2017 139  ! calculate vm(:) = -bm*vi(:) + vm  Markus Betzinger committed Apr 26, 2016 140 141 142  CALL CPP_BLAS_saxpy(nmap,-bm,vi,1,vm,1) !write(6,FMT='(5x," for it",i2,5x,f10.6)') it, bm END DO  Gregor Michalicek committed Nov 08, 2017 143 144 145  ! complete evaluation of vm ! vmnorm = -  Markus Betzinger committed Apr 26, 2016 146  vmnorm = CPP_BLAS_sdot(nmap,fm1,1,um,1) - smnorm  Gregor Michalicek committed Nov 08, 2017 147 148  ! if (vmnorm.lt.tol_10) stop  Markus Betzinger committed Apr 26, 2016 149  CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)  Gregor Michalicek committed Nov 08, 2017 150   Markus Betzinger committed Apr 26, 2016 151  ELSE IF (input%imix.EQ.5) THEN  Gregor Michalicek committed Nov 08, 2017 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 / vmnorm = one / CPP_BLAS_sdot(nmap,fm1,1,vm,1)  Markus Betzinger committed Apr 26, 2016 162 163 164  CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1) ELSE IF (input%imix.EQ.7) THEN  Gregor Michalicek committed Nov 08, 2017 165 166 167 168 169 170 171 172 173  !**************************************** ! generalized anderson method !**************************************** ! calculate vm = alpha*wfm1 -\sum  Markus Betzinger committed Apr 26, 2016 201  fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1)  Gregor Michalicek committed Nov 10, 2017 202  ! calculate sm(:) = (1.0-fmvm)*ui(:) + sm  Markus Betzinger committed Apr 26, 2016 203 204  CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1) END IF  Gregor Michalicek committed Nov 08, 2017 205   Markus Betzinger committed Apr 26, 2016 206  mit = mit + 1  Gregor Michalicek committed Nov 08, 2017 207  DEALLOCATE (fm1,sm1,ui,um,vi,vm,am)  Markus Betzinger committed Apr 26, 2016 208 209 210  END SUBROUTINE broyden END MODULE m_broyden