src/mixpulay.f90

   1
   2 ! Copyright (C) 2008 S. Suehara. This file is distributed under the terms of the
   3 ! GNU Lesser General Public License. See the file COPYING for license details.
   4
   5 !BOP
   6 ! !ROUTINE: mixpulay
   7 ! !INTERFACE:
   8 subroutine mixpulay(iscl,n,maxsd,nu,mu,f,d)
   9 ! !INPUT/OUTPUT PARAMETERS:
  10 !   iscl  : self-consistent loop number (in,integer)
  11 !   n     : vector length (in,integer)
  12 !   maxsd : maximum subspace dimension (in,integer)
  13 !   nu    : current output vector as well as the next input vector of the
  14 !           self-consistent loop (inout,real(n))
  15 !   mu    : used internally (inout,real(n,maxsd))
  16 !   f     : used internally (inout,real(n,maxsd))
  17 !   d     : RMS difference between old and new output vectors (out,real)
  18 ! !DESCRIPTION:
  19 !   Pulay's mixing scheme which uses direct inversion in the iterative subspace
  20 !   (DIIS). See {\it Chem. Phys. Lett.} {\bf 73}, 393 (1980).
  21 !
  22 ! !REVISION HISTORY:
  23 !   Created October 2008 (S. Suehara, NIMS)
  24 !   Modified, October 2008 (JKD)
  25 !EOP
  26 !BOC
  27 implicit none
  28 ! arguments
  29 integer, intent(in) :: iscl
  30 integer, intent(in) :: n
  31 integer, intent(in) :: maxsd
  32 real(8), intent(inout) :: nu(n)
  33 real(8), intent(inout) :: mu(n,maxsd)
  34 real(8), intent(inout) :: f(n,maxsd)
  35 real(8), intent(out) :: d
  36 ! local variables
  37 integer i,j,k,m,jc,jn,info
  38 ! initial mixing parameter
  39 real(8), parameter :: beta=0.1d0
  40 ! allocatable arrays
  41 integer, allocatable :: ipiv(:)
  42 real(8), allocatable :: alpha(:),a(:,:),work(:)
  43 ! external functions
  44 real(8) ddot
  45 external ddot
  46 if (n.lt.1) then
  47   write(*,*)
  48   write(*,'("Error(mixpulay): n < 1 : ",I8)') n
  49   write(*,*)
  50   stop
  51 end if
  52 if (maxsd.lt.2) then
  53   write(*,*)
  54   write(*,'("Error(mixpulay): maxsd < 2 : ",I8)') maxsd
  55   write(*,*)
  56   stop
  57 end if
  58 if (iscl.le.1) then
  59   mu(:,1)=nu(:)
  60   f(:,1)=0.d0
  61   d=1.d0
  62   return
  63 end if
  64 ! current index
  65 jc=mod(iscl-1,maxsd)+1
  66 ! next index
  67 jn=mod(iscl,maxsd)+1
  68 if (iscl.le.2) then
  69   nu(:)=beta*nu(:)+(1.d0-beta)*mu(:,1)
  70   f(:,2)=nu(:)-mu(:,1)
  71   mu(:,2)=nu(:)
  72   if (maxsd.ge.3) mu(:,3)=0.d0
  73   d=0.d0
  74   do k=1,n
  75     d=d+f(k,2)**2
  76   end do
  77   d=sqrt(d/dble(n))
  78   return
  79 end if
  80 ! matrix size
  81 m=min(iscl,maxsd)+1
  82 allocate(ipiv(m),alpha(m),a(m,m),work(m))
  83 ! compute f and RMS difference
  84 d=0.d0
  85 do k=1,n
  86   f(k,jc)=nu(k)-mu(k,jc)
  87   d=d+f(k,jc)**2
  88 end do
  89 d=sqrt(d/dble(n))
  90 ! solve the linear system
  91 a(:,:)=0.d0
  92 do i=1,m-1
  93   do j=i,m-1
  94     a(i,j)=a(i,j)+ddot(n,f(:,i),1,f(:,j),1)
  95   end do
  96   a(i,m)=1.d0
  97 end do
  98 alpha(:)=0.d0
  99 alpha(m)=1.d0
 100 call dsysv('U',m,1,a,m,ipiv,alpha,m,work,m,info)
 101 if (info.ne.0) then
 102   write(*,*)
 103   write(*,'("Error(mixpulay): could not solve linear system")')
 104   write(*,'(" DSYSV returned INFO = ",I8)') info
 105   write(*,*)
 106   stop
 107 end if
 108 nu(:)=0.d0
 109 do i=1,m-1
 110   nu(:)=nu(:)+alpha(i)*(mu(:,i)+f(:,i))
 111 end do
 112 mu(:,jn)=nu(:)
 113 deallocate(ipiv,alpha,a,work)
 114 return
 115 end subroutine
 116 !EOC
 117