broyden.F90 7.07 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  SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&  Gregor Michalicek committed Nov 10, 2017 15  hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,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 10, 2017 21  USE m_broyd_io  Gregor Michalicek committed Nov 08, 2017 22   Markus Betzinger committed Apr 26, 2016 23  IMPLICIT NONE  Gregor Michalicek committed Nov 08, 2017 24   Markus Betzinger committed Apr 26, 2016 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  Gregor Michalicek committed Nov 10, 2017 34  TYPE(t_hybrid),INTENT(IN) :: hybrid  Gregor Michalicek committed Nov 08, 2017 35 36 37  ! Scalar Arguments INTEGER, INTENT (IN) :: mmap,nmap  Gregor Michalicek committed Nov 10, 2017 38  INTEGER, INTENT (IN) :: mapmt,mapvac2  Markus Betzinger committed Apr 26, 2016 39  LOGICAL,OPTIONAL,INTENT(IN) :: lpot  Gregor Michalicek committed Nov 08, 2017 40 41 42  ! Array Arguments REAL, INTENT (IN) :: fm(nmap)  Markus Betzinger committed Apr 26, 2016 43  REAL, INTENT (INOUT) :: sm(nmap)  Gregor Michalicek committed Nov 08, 2017 44 45  ! Local Scalars  46  INTEGER :: i,it,k,nit,iread,nmaph, mit  Gregor Michalicek committed Nov 08, 2017 47  REAL :: bm,dfivi,fmvm,smnorm,vmnorm,alphan  Gregor Michalicek committed Nov 10, 2017 48  LOGICAL :: l_pot, l_exist  Gregor Michalicek committed Nov 08, 2017 49 50 51  REAL, PARAMETER :: one=1.0, zero=0.0 ! Local Arrays  Markus Betzinger committed Apr 26, 2016 52 53  REAL, ALLOCATABLE :: am(:) REAL, ALLOCATABLE :: fm1(:),sm1(:),ui(:),um(:),vi(:),vm(:)  Gregor Michalicek committed Nov 08, 2017 54 55  ! External Functions  Markus Betzinger committed Apr 26, 2016 56 57  REAL CPP_BLAS_sdot EXTERNAL CPP_BLAS_sdot  Gregor Michalicek committed Nov 08, 2017 58 59  ! External Subroutines  Markus Betzinger committed Apr 26, 2016 60  EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal  Gregor Michalicek committed Nov 08, 2017 61   Markus Betzinger committed Apr 26, 2016 62 63  dfivi = zero  Gregor Michalicek committed Nov 08, 2017 64 65  l_pot = .FALSE. IF (PRESENT(lpot)) l_pot = lpot  Markus Betzinger committed Apr 26, 2016 66 67  ALLOCATE (fm1(mmap),sm1(mmap),ui(mmap),um(mmap),vi(mmap),vm(mmap))  Gregor Michalicek committed Nov 08, 2017 68  ALLOCATE (am(input%maxiter+1))  Gregor Michalicek committed May 03, 2016 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  Gregor Michalicek committed Nov 10, 2017 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  Markus Betzinger committed Apr 26, 2016 82  IF (mit.NE.1) THEN  Gregor Michalicek committed Nov 08, 2017 83 84  ! 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 85   Gregor Michalicek committed Nov 10, 2017 86  CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))  Gregor Michalicek committed Nov 08, 2017 87  IF (ABS(input%alpha-alphan) > 0.0001) THEN  Markus Betzinger committed Apr 26, 2016 88  WRITE (6,*) 'mixing parameter has been changed; reset'  Daniel Wortmann committed Apr 29, 2016 89  WRITE (6,*) 'broyden algorithm or set alpha to',alphan  Gregor Michalicek committed Nov 08, 2017 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)  Markus Betzinger committed Apr 26, 2016 97  END IF  Gregor Michalicek committed Nov 08, 2017 98 99  ! save F_m and rho_m for next iteration  Markus Betzinger committed Apr 26, 2016 100 101  nit = mit +1 IF (nit > input%maxiter+1) nit = 1  Gregor Michalicek committed Nov 10, 2017 102  CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm)  Gregor Michalicek committed Nov 08, 2017 103   Markus Betzinger committed Apr 26, 2016 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 ui  Gregor Michalicek committed Nov 14, 2017 111  um(:nmap) = input%alpha * fm1(:nmap) + sm1(:nmap)  112   Gregor Michalicek committed Nov 08, 2017 113 114  iread = MIN(mit-1,input%maxiter+1) DO it = 2, iread  Gregor Michalicek committed Nov 10, 2017 115 116 117  CALL readUVec(input,hybrid,nmap,it-mit,mit,ui) CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)  Markus Betzinger committed Apr 26, 2016 118  am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)  119  ! calculate um(:) = -am(it)*ui(:) + um(:)  Markus Betzinger committed Apr 26, 2016 120 121 122  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 123   Markus Betzinger committed Apr 26, 2016 124  IF (input%imix.EQ.3) THEN  Gregor Michalicek committed Nov 08, 2017 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 :  Markus Betzinger committed Apr 26, 2016 134  smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1)  Gregor Michalicek committed Nov 08, 2017 135 136 137  ! generate vm = alpha*sm1 - \sum vi vm(:) = input%alpha * fm1(:)  138   Markus Betzinger committed Apr 26, 2016 139  DO it = 2,iread  Gregor Michalicek committed Nov 10, 2017 140 141  CALL readUVec(input,hybrid,nmap,it-mit,mit,ui) CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)  Markus Betzinger committed Apr 26, 2016 142  bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)  Gregor Michalicek committed Nov 10, 2017 143  ! calculate vm(:) = -bm*vi(:) + vm  Markus Betzinger committed Apr 26, 2016 144 145 146  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 147 148 149  ! complete evaluation of vm ! vmnorm = -  Markus Betzinger committed Apr 26, 2016 150  vmnorm = CPP_BLAS_sdot(nmap,fm1,1,um,1) - smnorm  Gregor Michalicek committed Nov 08, 2017 151 152  ! if (vmnorm.lt.tol_10) stop  Markus Betzinger committed Apr 26, 2016 153  CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)  Gregor Michalicek committed Nov 08, 2017 154   Markus Betzinger committed Apr 26, 2016 155  ELSE IF (input%imix.EQ.5) THEN  Gregor Michalicek committed Nov 08, 2017 156 157 158 159  !**************************************** ! broyden's second method !****************************************  160  ! multiply fm1 with metric matrix and store in vm: w |fm1>  Gregor Michalicek committed Nov 08, 2017 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 / vmnorm = one / CPP_BLAS_sdot(nmap,fm1,1,vm,1)  Markus Betzinger committed Apr 26, 2016 166 167 168  CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1) ELSE IF (input%imix.EQ.7) THEN  Gregor Michalicek committed Nov 08, 2017 169 170 171 172 173 174 175 176 177  !**************************************** ! generalized anderson method !**************************************** ! calculate vm = alpha*wfm1 -\sum  Markus Betzinger committed Apr 26, 2016 202  fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1)  203  ! calculate sm(:) = (1.0-fmvm)*um(:) + sm  Markus Betzinger committed Apr 26, 2016 204 205  CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1) END IF  Gregor Michalicek committed Nov 08, 2017 206 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