prg_pulaymixer_mod.F90

Go to the documentation of this file.

 
module prg_pulaymixer_mod
 
  use bml
  use prg_kernelparser_mod
  use prg_extras_mod
 
  implicit none
 
  private
 
  integer, parameter :: dp = kind(1.0d0)
 
  type, public :: mx_type
 
    character(20) :: mixertype
 
    integer       :: verbose
 
    integer       :: mpulay
 
    real(dp)      :: mixcoeff
 
    logical       :: mixeron
 
  end type mx_type
 
  public :: prg_parse_mixer, prg_qmixer, prg_linearmixer
 
contains
 
  subroutine prg_parse_mixer(input,filename)
    implicit none
    type(mx_type), intent(inout) :: input
    integer, parameter :: nkey_char = 1, nkey_int = 2, nkey_re = 1, nkey_log = 1
    character(len=*) :: filename
 
    !Library of keywords with the respective defaults.
    character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
         'MixerType=']
    character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
         'Linear']
 
    character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
         'Verbose=','MPulay=']
    integer :: valvector_int(nkey_int) = (/ &
         0,5/)
 
    character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
         'MixCoeff=']
    real(dp) :: valvector_re(nkey_re) = (/&
         0.25 /)
 
    character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
         'MixerON=']
    logical :: valvector_log(nkey_log) = (/&
         .true./)
 
    !Start and stop characters
    character(len=50), parameter :: startstop(2) = [character(len=50) :: &
         'MIXER{', '}']
 
    call prg_parsing_kernel(keyvector_char,valvector_char&
         ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
         keyvector_log,valvector_log,trim(filename),startstop)
 
    !Characters
    input%mixertype = valvector_char(1)
 
    !Integers
    input%verbose = valvector_int(1)
    input%mpulay = valvector_int(2)
 
    !Logicals
    input%mixeron = valvector_log(1)
 
    !Reals
    input%mixcoeff = valvector_re(1)
 
  end subroutine prg_parse_mixer
 
 
  subroutine prg_qmixer(charges,oldcharges,dqin,dqout,scferror,piter,pulaycoef,mpulay,verbose)
 
    implicit none
    integer :: i,j,info,s,k,n
    real(dp) :: alpha,coeff
    real(dp), intent(in) :: pulaycoef
    real(dp), intent(inout) :: scferror
    real(dp) :: error, errop
    integer :: piter
    real(dp) :: dt, atenuacion
    integer :: nint, na, nb
    integer, intent(in) :: mpulay,verbose
    real(dp), allocatable :: d(:),dl(:),dnewin(:)
    real(dp), allocatable :: dnewout(:)
    real(dp), intent(inout) :: charges(:)
    real(dp), allocatable, intent(inout) :: oldcharges(:)
    real(dp), allocatable, intent(inout) :: dqin(:,:),dqout(:,:)
    real(dp), allocatable :: coef(:,:),b(:),ipiv(:)
 
    if(verbose > 0)write(*,*)"Performing Pulay mixing scheme ..."
    n=size(charges)
 
    alpha = pulaycoef !the coefficient for mixing
 
    if(allocated(oldcharges).eqv..false.)then
      allocate(oldcharges(n),dqin(n,mpulay),dqout(n,mpulay))
    endif
 
    if(allocated(dqin).eqv..false.)then
      allocate(dqin(n,mpulay),dqout(n,mpulay))
    endif
 
    s=min(piter-1,mpulay) !mpulay is the iteration number
 
    if(piter.eq.1) then
      charges=(1.0_dp-alpha)*oldcharges + alpha*charges
      scferror = prg_norm2(charges(:)-oldcharges(:))
      if(verbose.ge.1)then
        write(*,*)"SCF error =", scferror
      endif
      oldcharges=charges
    else
 
      allocate(d(n),dnewin(n),dnewout(n))
 
      d=charges
 
      allocate(coef(s+1,s+1)) !Allocating the coeffs matrix
      allocate(b(s+1))
      allocate(ipiv(s+1))
 
      if(piter.le.mpulay+1)then  !If piter=6 => mpulay=5
        dqin(:,piter-1)=oldcharges(:)
        dqout(:,piter-1)=d(:)
      endif
 
      if(piter.gt.mpulay+1)then
 
        do j=1,s-1
          dqin(:,j)=dqin(:,j+1)
          dqout(:,j)=dqout(:,j+1)
        enddo
 
        dqin(:,s)=oldcharges(:)
        dqout(:,s)=d(:)
 
      endif
 
      coef=0.0_dp
 
      do i=1,s+1
        coef(s+1,i)=-1.0d0
        coef(i,s+1)=-1.0d0
        b(i)=0
      enddo
      b(s+1)=-1.0d0
      coef(s+1,s+1)=0.0_dp
 
      do i=1,s
        do j=1,s
          do k=1,n
            coef(i,j)=coef(i,j)+(dqout(k,i)-dqin(k,i))*(dqout(k,j)-dqin(k,j))
          enddo
        enddo
      enddo
 
      if(verbose.ge.1)then
        write(*,*)"coefs"
        do i=1,s+1
          write(*,'(10f12.5)')(coef(i,j),j=1,s+1)
        enddo
        write(*,*)"dqin"
        write(*,'(10f12.5)')(dqin(n,j),j=1,s)
      endif
 
      call dgesv(s+1,1,coef,s+1,ipiv,b,s+1,info)
 
      if(info.ne.0)then 
        write(*,*)'WARNING: (At prg_qmixer) Singular matrix in pulay. Doing linear mixing'
        !STOP 
        charges=(1.0_dp-alpha*0.1_dp)*oldcharges + alpha*0.1_dp*charges
        scferror = prg_norm2(charges(:)-oldcharges(:))
        oldcharges = charges
        !piter = 1
      else
 
        dnewin=0.0_dp
        dnewout=0.0_dp
 
        if(verbose.ge.1)then
          write(*,*)"eigen coefs"
          write(*,'(6f10.5)')(b(j),j=1,s)
        endif
 
        do j=1,s
          dnewin(:)=dnewin(:)+b(j)*dqin(:,j)
          dnewout(:)= dnewout(:)+b(j)*dqout(:,j)
        enddo
 
        d=(1.0_dp-alpha)*dnewin + alpha*dnewout
 
        scferror = prg_norm2(charges(:)-oldcharges(:))
 
        if(verbose.ge.1)then
          write(*,*)"SCF error =", scferror
        endif
 
        charges=d
 
        oldcharges=d
      endif
 
    endif
 
  end subroutine prg_qmixer
 
  subroutine prg_linearmixer(charges,oldcharges,scferror,linmixcoef,verbose)
    implicit none
    real(dp), intent(in) :: linmixcoef
    real(dp), intent(inout) :: scferror
    integer, intent(in) :: verbose
    real(dp), allocatable, intent(inout) :: charges(:),oldcharges(:)
 
    scferror = prg_norm2(charges(:)-oldcharges(:))
 
    if(verbose.ge.1)then
      write(*,*)"SCF error =", scferror
    endif
 
    charges = (1.0_dp - linmixcoef)*oldcharges + linmixcoef*charges
    oldcharges = charges
 
  end subroutine prg_linearmixer
 
end module prg_pulaymixer_mod