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.
7 ! calculate new occupation numbers using derivatives of total energy
11 integer, parameter :: maxit
=10000
13 real(8), parameter :: eps
=1.d
-12
14 real(8) tau
,sum
,gs
,gsp
,dgs
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
25 sum
=sum
+wkpt(ik
)*occsv(ist
,ik
)
28 t1
=(chgval
-sum
)/dble(nstsv
)
29 occsv(:,:)=occsv(:,:)+t1
30 ! redistribute charge so that occupancies are in the interval [0,occmax]
34 if (occsv(ist
,ik
).gt
.occmax
) then
35 sum
=sum
+wkpt(ik
)*(occsv(ist
,ik
)-occmax
)
38 if (occsv(ist
,ik
).lt
.0.d0
) then
39 sum
=sum
+wkpt(ik
)*occsv(ist
,ik
)
47 t1
=wkpt(ik
)*(occmax
-occsv(ist
,ik
))
49 occsv(ist
,ik
)=occsv(ist
,ik
)+t1
/wkpt(ik
)
52 t1
=wkpt(ik
)*occsv(ist
,ik
)
54 occsv(ist
,ik
)=occsv(ist
,ik
)-t1
/wkpt(ik
)
61 ! find suitable value of kapa such that sum of gamma is 0
72 gamma(ist
,ik
)=t1
*(occmax
-occsv(ist
,ik
))
74 gamma(ist
,ik
)=t1
*occsv(ist
,ik
)
76 gs
=gs
+wkpt(ik
)*gamma(ist
,ik
)
77 sum
=sum
+wkpt(ik
)*gamma(ist
,ik
)**2
83 if (t1
.lt
.eps
) goto 10
86 if (gs
*dgs
.gt
.0.d0
) dkapa
=-dkapa
87 if (gs
*gsp
.lt
.0.d0
) then
97 write(*,'("Error(rdmdedn): could not find offset")')
101 ! normalize gamma if sum of squares is greater than 1
105 sum
=sum
+wkpt(ik
)*gamma(ist
,ik
)**2
108 if (sum
.gt
.1.d0
) then
110 gamma(:,:)=t1
*gamma(:,:)
112 ! find step size which keeps occupancies in the interval [0,occmax]
115 if (abs(tau
).lt
.eps
) goto 30
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
125 if (gamma(ist
,ik
).lt
.0.d0
) then
137 occsv(ist
,ik
)=occsv(ist
,ik
)+tau
*gamma(ist
,ik
)
140 ! write derivatives and occupancies to a file
141 call rdmwritededn(dedn
)
142 deallocate(dedn
,gamma
)