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 SUBROUTINE broyden(& & cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,& & mmap,nmaph,mapmt,mapvac2,nmap,fm,mit,sm,lpot) #include"cpp_double.h" USE m_metric USE m_types IMPLICIT NONE 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 ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: mmap,nmap INTEGER, INTENT (IN) :: mapmt,mapvac2 INTEGER, INTENT (INOUT) :: mit LOGICAL,OPTIONAL,INTENT(IN) :: lpot ! .. ! .. Array Arguments .. REAL, INTENT (IN) :: fm(nmap) REAL, INTENT (INOUT) :: sm(nmap) !-odim !+odim ! .. ! .. 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 .. REAL, ALLOCATABLE :: am(:) REAL, ALLOCATABLE :: fm1(:),sm1(:),ui(:),um(:),vi(:),vm(:) ! .. ! .. External Functions .. REAL CPP_BLAS_sdot EXTERNAL CPP_BLAS_sdot ! .. ! .. External Subroutines .. EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal ! .. dfivi = zero IF (PRESENT(lpot)) THEN l_pot=lpot ELSE l_pot=.FALSE. ENDIF ALLOCATE (fm1(mmap),sm1(mmap),ui(mmap),um(mmap),vi(mmap),vm(mmap)) ALLOCATE ( am(input%maxiter+1) ) fm1 = 0.0 sm1 = 0.0 ui = 0.0 um = 0.0 vi = 0.0 vm = 0.0 am = 0.0 ! IF (mit.NE.1) THEN ! ! load input charge density (sm1) and difference of ! in and out charge densities (fm1) from previous iteration (m-1) REWIND 57 READ (57) mit,alphan,(fm1(i),i=1,nmap),(sm1(i),i=1,nmap) IF ( ABS(input%alpha-alphan) > 0.0001 ) THEN WRITE (6,*) 'mixing parameter has been changed; reset' WRITE (6,*) 'broyden algorithm or set alpha to',alphan CALL juDFT_error("mixing parameter (input) changed", calledby& & ="broyden") ENDIF ! ! loop to 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) END IF ! ! save F_m and rho_m for next iteration ! REWIND 57 nit = mit +1 IF (nit > input%maxiter+1) nit = 1 WRITE (57) nit,input%alpha,fm,sm ! 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 ! DO k = 1,nmap um(k) = input%alpha*fm1(k) + sm1(k) END DO iread=MIN(mit-1,input%maxiter+1) DO it = 2,iread READ (59,rec=(it-1)*2-1) (ui(i),i=1,nmap) READ (59,rec=(it-1)*2) (vi(i),i=1,nmap),dfivi am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1) 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 ! ! **************************************** ! broyden's first method ! **************************************** ! IF (input%imix.EQ.3) THEN ! ! 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 : smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1) ! ! loop to generate vm = alpha*sm1 - \sum vi ! vm(:) = input%alpha*fm1(:) DO it = 2,iread READ (59,rec=2*(it-1)-1)(ui(i),i=1,nmap) READ (59,rec=2*(it-1))(vi(i),i=1,nmap),dfivi bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1) CALL CPP_BLAS_saxpy(nmap,-bm,vi,1,vm,1) !write(6,FMT='(5x," for it",i2,5x,f10.6)') it, bm END DO ! ! complete evaluation of vm ! vmnorm = - ! vmnorm = CPP_BLAS_sdot(nmap,fm1,1,um,1) - smnorm ! * if (vmnorm.lt.tol_10) stop CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1) ! write bm(it) ! ELSE IF (input%imix.EQ.5) THEN ! ! **************************************** ! 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) CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1) ELSE IF (input%imix.EQ.7) THEN ! ! **************************************** ! generalized anderson method ! **************************************** ! ! calculate vm = alpha*wfm1 -\sum ! fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1) CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1) END IF mit = mit + 1 DEALLOCATE ( fm1,sm1,ui,um,vi,vm,am ) END SUBROUTINE broyden END MODULE m_broyden