src/rdmvaryn.f90

   1
   2 ! Copyright (C) 2007-2008 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
   3 ! This file is distributed under the terms of the GNU Lesser General Public
   4 ! License. See the file COPYING for license details.
   5
   6 subroutine rdmvaryn
   7 ! calculate new occupation numbers using derivatives of total energy
   8 use modmain
   9 implicit none
  10 ! local variables
  11 integer, parameter :: maxit=10000
  12 integer it,ik,ist
  13 real(8), parameter :: eps=1.d-12
  14 real(8) tau,sum,gs,gsp,dgs
  15 real(8) kapa,dkapa,t1
  16 ! allocatable arrays
  17 real(8), allocatable :: dedn(:,:)
  18 real(8), allocatable :: gamma(:,:)
  19 allocate(dedn(nstsv,nkpt))
  20 allocate(gamma(nstsv,nkpt))
  21 ! add constant to occupancies for charge conservation
  22 sum=0.d0
  23 do ik=1,nkpt
  24   do ist=1,nstsv
  25     sum=sum+wkpt(ik)*occsv(ist,ik)
  26   end do
  27 end do
  28 t1=(chgval-sum)/dble(nstsv)
  29 occsv(:,:)=occsv(:,:)+t1
  30 ! redistribute charge so that occupancies are in the interval [0,occmax]
  31 sum=0.d0
  32 do ik=1,nkpt
  33   do ist=1,nstsv
  34     if (occsv(ist,ik).gt.occmax) then
  35       sum=sum+wkpt(ik)*(occsv(ist,ik)-occmax)
  36       occsv(ist,ik)=occmax
  37     end if
  38     if (occsv(ist,ik).lt.0.d0) then
  39       sum=sum+wkpt(ik)*occsv(ist,ik)
  40       occsv(ist,ik)=0.d0
  41     end if
  42   end do
  43 end do
  44 do ist=1,nstsv
  45   do ik=1,nkpt
  46     if (sum.gt.0.d0) then
  47       t1=wkpt(ik)*(occmax-occsv(ist,ik))
  48       t1=min(t1,sum)
  49       occsv(ist,ik)=occsv(ist,ik)+t1/wkpt(ik)
  50       sum=sum-t1
  51     else
  52       t1=wkpt(ik)*occsv(ist,ik)
  53       t1=min(t1,-sum)
  54       occsv(ist,ik)=occsv(ist,ik)-t1/wkpt(ik)
  55       sum=sum+t1
  56     end if
  57   end do
  58 end do
  59 ! get the derivatives
  60 call rdmdedn(dedn)
  61 ! find suitable value of kapa such that sum of gamma is 0
  62 gsp=0.d0
  63 kapa=0.d0
  64 dkapa=0.1d0
  65 do it=1,maxit
  66   gs=0.d0
  67   sum=0.d0
  68   do ik=1,nkpt
  69     do ist=1,nstsv
  70       t1=dedn(ist,ik)-kapa
  71       if (t1.gt.0.d0) then
  72         gamma(ist,ik)=t1*(occmax-occsv(ist,ik))
  73       else
  74         gamma(ist,ik)=t1*occsv(ist,ik)
  75       end if
  76       gs=gs+wkpt(ik)*gamma(ist,ik)
  77       sum=sum+wkpt(ik)*gamma(ist,ik)**2
  78     end do
  79   end do
  80   sum=sqrt(sum)
  81   sum=max(sum,1.d0)
  82   t1=abs(gs)/sum
  83   if (t1.lt.eps) goto 10
  84   if (it.ge.2) then
  85     dgs=gs-gsp
  86     if (gs*dgs.gt.0.d0) dkapa=-dkapa
  87     if (gs*gsp.lt.0.d0) then
  88       dkapa=0.5d0*dkapa
  89     else
  90       dkapa=1.1d0*dkapa
  91     end if
  92   end if
  93   gsp=gs
  94   kapa=kapa+dkapa
  95 end do
  96 write(*,*)
  97 write(*,'("Error(rdmdedn): could not find offset")')
  98 write(*,*)
  99 stop
 100 10 continue
 101 ! normalize gamma if sum of squares is greater than 1
 102 sum=0.d0
 103 do ik=1,nkpt
 104   do ist=1,nstsv
 105     sum=sum+wkpt(ik)*gamma(ist,ik)**2
 106   end do
 107 end do
 108 if (sum.gt.1.d0) then
 109   t1=1.d0/sqrt(sum)
 110   gamma(:,:)=t1*gamma(:,:)
 111 end if
 112 ! find step size which keeps occupancies in the interval [0,occmax]
 113 tau=taurdmn
 114 20 continue
 115 if (abs(tau).lt.eps) goto 30
 116 do ik=1,nkpt
 117   do ist=1,nstsv
 118     t1=occsv(ist,ik)+tau*gamma(ist,ik)
 119     if (gamma(ist,ik).gt.0.d0) then
 120       if (t1.gt.occmax+eps) then
 121         tau=0.75d0*tau
 122         goto 20
 123       end if
 124     end if
 125     if (gamma(ist,ik).lt.0.d0) then
 126       if (t1.lt.-eps) then
 127         tau=0.75d0*tau
 128         goto 20
 129       end if
 130     end if
 131   end do
 132 end do
 133 30 continue
 134 ! update occupancies
 135 do ik=1,nkpt
 136   do ist=1,nstsv
 137     occsv(ist,ik)=occsv(ist,ik)+tau*gamma(ist,ik)
 138   end do
 139 end do
 140 ! write derivatives and occupancies to a file
 141 call rdmwritededn(dedn)
 142 deallocate(dedn,gamma)
 143 return
 144 end subroutine
 145