Skip to content

GitLab

fleur
fleur
Commit addfc05e authored by Daniel Wortmann's avatar Daniel Wortmann
Browse files
Options

Merge branch 'develop' into mixer

parents ebcccf8e 0cf0bbf4
Show whitespace changes
Inline Side-by-side
Showing with 703 additions and 741 deletions
cdn/rot_den_mat.f
View file @ addfc05e
......@@ -15,6 +15,7 @@ c This subroutine rotates the direction of the magnetization of the
c density matrix by multiplying with the unitary 2x2 spin rotation
c matrix. --> U*rho*U^dagger
c Philipp Kurz 2000-02-03
c new method for improved stability (l_new=t) gb'19
c***********************************************************************
use m_constants
......@@ -28,11 +29,29 @@ C .. Scalar Arguments ..
C ..
C .. Local Scalars ..
INTEGER ispin
REAL eps
REAL eps,r11n,r22n
COMPLEX r21n
LOGICAL l_new
C ..
C .. Local Arrays ..
COMPLEX u2(2,2),rho(2,2),rhoh(2,2)
C ..
l_new = .true.
IF (l_new) THEN
r11n = 0.5*(1.0+cos(beta))*rho11 - sin(beta)*real(rho21) +
+ 0.5*(1.0-cos(beta))*rho22
r22n = 0.5*(1.0-cos(beta))*rho11 + sin(beta)*real(rho21) +
+ 0.5*(1.0+cos(beta))*rho22
r21n = CMPLX(cos(alph),-sin(alph))*(sin(beta)*(rho11-rho22) +
+ 2.0*(cos(beta)*real(rho21)-cmplx(0.0,aimag(rho21))))*0.5
rho11 = r11n
rho22 = r22n
rho21 = r21n
ELSE
eps = 1.0e-10
......@@ -79,6 +98,8 @@ c---> are real.
rho22 = real(rho(2,2))
rho21 = rho(2,1)
ENDIF
END SUBROUTINE rot_den_mat
END MODULE m_rotdenmat
eigen/hsvac.F90
View file @ addfc05e
......@@ -40,14 +40,14 @@ CONTAINS
! .. Local Scalars ..
COMPLEX hij,sij,apw_lo,c_1
REAL d2,gz,sign,th,wronk
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ipot,ii0,i0
INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ii0,i0
INTEGER ivac,irec,imz,igvm2,igvm2i,s1,s2
INTEGER jspin1,jspin2,jmax,jsp_start,jsp_end
INTEGER i_start,nc,nc_0
! ..
! .. Local Arrays ..
INTEGER:: nv2(input%jspins)
INTEGER kvac1(DIMENSION%nv2d,input%jspins),kvac2(DIMENSION%nv2d,input%jspins)
INTEGER kvac(2,DIMENSION%nv2d,input%jspins)
INTEGER map2(DIMENSION%nvd,input%jspins)
COMPLEX tddv(DIMENSION%nv2d,DIMENSION%nv2d),tduv(DIMENSION%nv2d,DIMENSION%nv2d)
COMPLEX tudv(DIMENSION%nv2d,DIMENSION%nv2d),tuuv(DIMENSION%nv2d,DIMENSION%nv2d)
......@@ -66,16 +66,14 @@ CONTAINS
nv2(jspin) = 0
k_loop:DO k = 1,lapw%nv(jspin)
DO j = 1,nv2(jspin)
IF (lapw%k1(k,jspin).EQ.kvac1(j,jspin)&
.AND. lapw%k2(k,jspin).EQ.kvac2(j,jspin)) THEN
IF (all(lapw%gvec(1:2,k,jspin)==kvac(1:2,j,jspin))) THEN
map2(k,jspin) = j
CYCLE k_loop
END IF
ENDDO
nv2(jspin) = nv2(jspin) + 1
IF (nv2(jspin)>DIMENSION%nv2d) CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac")
kvac1(nv2(jspin),jspin) = lapw%k1(k,jspin)
kvac2(nv2(jspin),jspin) = lapw%k2(k,jspin)
kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,k,jspin)
map2(k,jspin) = nv2(jspin)
ENDDO k_loop
ENDDO
......@@ -86,20 +84,16 @@ CONTAINS
s1=MIN(SIZE(hmat,1),jspin1) !in colinear case s1=1
DO jspin2=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
s2=MIN(SIZE(hmat,1),jspin2) !in colinear case s2=1
ipot=3
IF (jspin1==jspin2) ipot=jspin1
!---> get the wavefunctions and set up the tuuv, etc matrices
jspin=jsp
CALL vacfun(&
vacuum,DIMENSION,stars,&
jsp,input,noco,jspin1,jspin2,&
sym, cell,ivac,evac(1,1),lapw%bkpt,v%vacxy(:,:,ivac,ipot),v%vacz(:,:,:),kvac1,kvac2,nv2,&
vacuum,stars,&
input,noco,jspin1,jspin2,&
sym, cell,ivac,evac,lapw%bkpt,v%vacxy,v%vacz,kvac,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
!
!---> generate a and b coeffficients
!
IF (noco%l_noco) THEN
DO jspin = 1,input%jspins
DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
DO k = 1,lapw%nv(jspin)
gz = sign*cell%bmat(3,3)*lapw%k3(k,jspin)
i2 = map2(k,jspin)
......@@ -109,16 +103,6 @@ CONTAINS
b(k,jspin) = c_1 * CMPLX(duz(i2,jspin), gz* uz(i2,jspin) )
ENDDO
ENDDO
ELSE
DO k = 1,lapw%nv(jsp)
gz = sign*cell%bmat(3,3)*lapw%gvec(3,k,jsp)
i2 = map2(k,jsp)
th = gz*cell%z1
c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
a(k,jsp) = - c_1 * CMPLX(dudz(i2,jsp), gz*udz(i2,jsp) )
b(k,jsp) = c_1 * CMPLX(duz(i2,jsp), gz* uz(i2,jsp) )
ENDDO
ENDIF
!---> update hamiltonian and overlap matrices
IF (jspin1==jspin2) THEN
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
......@@ -153,30 +137,25 @@ CONTAINS
!Diagonal term of Overlapp matrix, Hamiltonian later
sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2)
IF (hmat(1,1)%l_real) THEN
smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
smat(s2,s1)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
ELSE
smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
smat(s2,s1)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
ENDIF
ENDDO
ENDIF
!---> hamiltonian update
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
DO i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
i0=(i-1)/mpi%n_size+1 !local column index
ik = map2(i,jspin2)
DO j = 1,MIN(i,lapw%nv(jspin1))
jk = map2(j,jspin1)
IF (jspin2>jspin1) THEN
hij = CONJG(CONJG(a(j,jspin1))* (tuuv(jk,ik)*a(i,jspin2) +tudv(jk,ik)*b(i,jspin2))&
+ CONJG(b(j,jspin1))* (tddv(jk,ik)*b(i,jspin2) +tduv(jk,ik)*a(i,jspin2)))
ELSE
hij = CONJG(a(i,jspin2))* (tuuv(ik,jk)*a(j,jspin1) +tudv(ik,jk)*b(j,jspin1))&
+ CONJG(b(i,jspin2))* (tddv(ik,jk)*b(j,jspin1) +tduv(ik,jk)*a(j,jspin1))
ENDIF
ik = map2(i,jspin1)
DO j = 1,MERGE(i,lapw%nv(jspin2),jspin1==jspin2)
jk = map2(j,jspin2)
hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))&
+ CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2))
IF (hmat(1,1)%l_real) THEN
hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij)
hmat(s2,s1)%data_r(j,i0) = hmat(s2,s1)%data_r(j,i0) + REAL(hij)
ELSE
hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + hij
hmat(s2,s1)%data_c(j,i0) = hmat(s2,s1)%data_c(j,i0) + hij
ENDIF
ENDDO
ENDDO
......
eigen/vacfun.f90
View file @ addfc05e
......@@ -2,8 +2,8 @@ MODULE m_vacfun
use m_juDFT
CONTAINS
SUBROUTINE vacfun(&
vacuum,DIMENSION,stars, jsp,input,noco,jsp1,jsp2,&
sym, cell,ivac,evac,bkpt, vxy,vz,kvac1,kvac2,nv2,&
vacuum,stars,input,noco,jspin1,jspin2,&
sym, cell,ivac,evac,bkpt, vxy,vz,kvac,nv2,&
tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
!*********************************************************************
! determines the necessary values and derivatives on the vacuum
......@@ -19,7 +19,6 @@ CONTAINS
USE m_types
IMPLICIT NONE
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
......@@ -28,31 +27,35 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp ,ivac,jsp1,jsp2
INTEGER, INTENT (IN) :: ivac,jspin1,jspin2
REAL, INTENT (OUT) :: wronk
! ..
! .. Array Arguments ..
INTEGER, INTENT (IN) :: nv2(input%jspins)
INTEGER, INTENT (IN) :: kvac1(dimension%nv2d,input%jspins),kvac2(dimension%nv2d,input%jspins)
COMPLEX, INTENT (IN) :: vxy(vacuum%nmzxyd,stars%ng2-1)
COMPLEX, INTENT (OUT):: tddv(dimension%nv2d,dimension%nv2d),tduv(dimension%nv2d,dimension%nv2d)
COMPLEX, INTENT (OUT):: tudv(dimension%nv2d,dimension%nv2d),tuuv(dimension%nv2d,dimension%nv2d)
REAL, INTENT (IN) :: vz(vacuum%nmzd,2,4) ,evac(2,input%jspins)
INTEGER, INTENT (IN) :: nv2(:)!(input%jspins)
INTEGER, INTENT (IN) :: kvac(:,:,:)!(2,dimension%nv2d,input%jspins)
COMPLEX, INTENT (IN) :: vxy(:,:,:,:) !(vacuum%nmzxyd,stars%ng2-1,nvac,:)
COMPLEX, INTENT (OUT):: tddv(:,:),tduv(:,:)!(dimension%nv2d,dimension%nv2d)
COMPLEX, INTENT (OUT):: tudv(:,:),tuuv(:,:)!(dimension%nv2d,dimension%nv2d)
REAL,ALLOCATABLE,INTENT (IN) :: vz(:,:,:) !(vacuum%nmzd,2,4) ,
REAL, INTENT (IN) :: evac(:,:)!(2,input%jspins)
REAL, INTENT (IN) :: bkpt(3)
REAL, INTENT (OUT):: udz(dimension%nv2d,input%jspins),uz(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: dudz(dimension%nv2d,input%jspins),duz(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: ddnv(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: udz(:,:),uz(:,:)!(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: dudz(:,:),duz(:,:)!(dimension%nv2d,input%jspins)
REAL, INTENT (OUT):: ddnv(:,:)!(dimension%nv2d,input%jspins)
! ..
! .. Local Scalars ..
REAL ev,scale,xv,yv,vzero
REAL ev,scale,xv,yv,vzero,fac
COMPLEX phase
INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin
INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,ipot
LOGICAL tail
! ..
! .. Local Arrays ..
REAL u(vacuum%nmzd,dimension%nv2d,input%jspins),ud(vacuum%nmzd,dimension%nv2d,input%jspins)
REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins)
REAL v(3),x(vacuum%nmzd), qssbti(2,2)
! ..
fac=MERGE(1.0,-1.0,jspin1<=jspin2)
ipot=MERGE(jspin1,3,jspin1==jspin2)
tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0
udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0.
tail = .true.
......@@ -60,15 +63,12 @@ CONTAINS
!---> wronksian for the schrodinger equation given by an identity
wronk = 2.0
!---> setup the spin-spiral q-vector
qssbti(1,1) = - noco%qss(1)/2
qssbti(2,1) = - noco%qss(2)/2
qssbti(1,2) = + noco%qss(1)/2
qssbti(2,2) = + noco%qss(2)/2
qssbti(1:2,1) = - noco%qss(1:2)/2
qssbti(1:2,2) = + noco%qss(1:2)/2
!---> generate basis functions for each 2-d k+g
DO jspin = 1,input%jspins
DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
DO ik = 1,nv2(jspin)
v(1) = bkpt(1) + kvac1(ik,jspin) + qssbti(1,jspin)
v(2) = bkpt(2) + kvac2(ik,jspin) + qssbti(2,jspin)
v(1:2) = bkpt(1:2) + kvac(:,ik,jspin) + qssbti(1:2,jspin)
v(3) = 0.0
ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat))
vzero = vz(vacuum%nmzd,ivac,jspin)
......@@ -87,12 +87,12 @@ CONTAINS
enddo
ENDDO
!---> set up the tuuv, etc. matrices
DO ik = 1,nv2(jsp1)
DO jk = 1,nv2(jsp2)
DO ik = 1,nv2(jspin1)
DO jk = 1,nv2(jspin2)
!---> determine the warping component of the potential
i1 = kvac1(ik,jsp1) - kvac1(jk,jsp2)
i2 = kvac2(ik,jsp1) - kvac2(jk,jsp2)
i1 = kvac(1,ik,jspin1) - kvac(1,jk,jspin2)
i2 = kvac(2,ik,jspin1) - kvac(2,jk,jspin2)
i3 = 0
ind3 = stars%ig(i1,i2,i3)
IF (ind3.EQ.0) CYCLE
......@@ -114,98 +114,96 @@ CONTAINS
!---> tuuv
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tuuv(ik,jk) = phase*cmplx(xv,yv)
!---> tddv
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) =ud(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) =ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tddv(ik,jk) = phase*cmplx(xv,yv)
!---> tudv
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tudv(ik,jk) = phase*cmplx(xv,yv)
!---> tduv
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
DO i = 1,vacuum%nmzxy
x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
enddo
CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
tduv(ik,jk) = phase*cmplx(xv,yv)
ELSE
!---> diagonal (film muffin-tin) terms
IF (jsp1==jsp2) THEN
tuuv(ik,ik) = cmplx(evac(ivac,jsp1),0.0)
tddv(ik,ik) = cmplx(evac(ivac,jsp1)*ddnv(ik,jsp1),0.0)
IF (jspin1==jspin2) THEN
tuuv(ik,ik) = cmplx(evac(ivac,jspin1),0.0)
tddv(ik,ik) = cmplx(evac(ivac,jspin1)*ddnv(ik,jspin1),0.0)
tudv(ik,ik) = cmplx(0.5,0.0)
tduv(ik,ik) = cmplx(0.5,0.0)
ELSE
!---> tuuv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tuuv(ik,jk) = cmplx(xv,yv)
!---> tddv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tddv(ik,jk) = cmplx(xv,yv)
!---> tudv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tudv(ik,jk) = cmplx(xv,yv)
!---> tduv
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
DO i = 1,vacuum%nmz
x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
ENDDO
CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
tduv(ik,jk) = cmplx(xv,yv)
......
force/CMakeLists.txt
View file @ addfc05e
set(fleur_F77 ${fleur_F77}
force/bfgs.f
force/bfgs0.f
force/relax.F
)
set(fleur_F90 ${fleur_F90}
force/bfgs0.f90
force/bfgs.f90
force/force_b8.f90
force/force_0.f90
force/force_a12.f90
......
force/bfgs.f deleted 100644 → 0
View file @ ebcccf8e
MODULE m_bfgs
use m_juDFT
CONTAINS
SUBROUTINE bfgs(
> ntype,istep,istep0,force,
> zat,xa,thetad,epsdisp,epsforce,tote,
X xold,y,h,tau0,
< lconv)
!*******************************************************************
! relaxes the forces using the BFGS quasi-Newton method.
! input:
! istep = atomic step in this run
! istep0 = number of atomic steps in previous runs used
! in updating hessian
!
! output:
! lconv = logical true if forces are converged to tolerance
! given in epsforce
!
! the positions and forces from this step are added
! to file force.dat
!*******************************************************************
IMPLICIT NONE
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ntype,istep,istep0
REAL, INTENT (IN) :: epsdisp,epsforce,tote
REAL, INTENT (INOUT) :: thetad,xa
LOGICAL, INTENT (OUT) :: lconv
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: force(3,ntype),zat(ntype)
REAL, INTENT (INOUT) :: tau0(3,ntype)
REAL, INTENT (INOUT) :: y(3*ntype),xold(3*ntype)
REAL, INTENT (INOUT) :: h(3*ntype,3*ntype)
! ..
! .. Local Scalars ..
INTEGER :: i,j,n,nn,ist,na
REAL :: py,yy,alpha,s,zatm,fmax,gamma,d2,dispmax
! ..
! .. Local Arrays ..
REAL :: f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype)
!
n=3*ntype
ist = istep+istep0
!--- > get positions and forces in correct form and output
DO na = 1,ntype
nn = 3*(na-1)
xnew(nn+1:nn+3) = tau0(:,na)
f(nn+1:nn+3) = force(:,na)
ENDDO
!Write new entry into forces.dat
OPEN(43,file ='forces.dat',status ='unknown',form='formatted'
$ ,position ='append')
WRITE (43,'(a,f20.10)') "energy =",tote
WRITE (43,'(3f16.9,3f14.9)')
> ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1)
CLOSE(43)
!--- > get maximum force
fmax = 0.0
DO na = 1,ntype
nn = 3*(na-1)
fmax = MAX( fmax, (f(nn+1)**2+f(nn+2)**2+f(nn+3)**2) )
ENDDO
fmax = SQRT(fmax)
WRITE (6,1000) istep,fmax
IF ( fmax<epsforce) THEN
lconv = .TRUE.
RETURN
ELSE
lconv = .FALSE.
ENDIF
1000 FORMAT (1x/,' atomic step',i4,': maximum force =',
+ 1p,e14.6,' hartrees/a.u.')
!------------------------------------------------------------
! ===> if first step, go along gradient
IF (ist==1) THEN
!--- > choose a reasonable first guess for scaling, but
!--- > limit displacement to a maximum of 0.25 a.u.
!--- > (may need to be changed for different systems)
!--- > this choice is based on a Debye temperature of 330K;
!--- > modify as needed (change thetad in param.8)
zatm = 0.0
DO i = 1,ntype
zatm = MAX(zatm,zat(i))
ENDDO
IF (ABS(xa)<1.0e-10) THEN
WRITE (6,*) 'WARNING, xa = 0.0 set to 2.0'
xa = 2.0
ENDIF
IF (ABS(thetad)<1.0e-10) THEN
WRITE (6,*) 'WARNING, thetad = 0.0 set to 330.0'
thetad = 330.0
ENDIF
alpha = (250.0/(zatm*xa))*((330./thetad)**2)
IF ( alpha*fmax*xa > 0.15 ) alpha = 0.25/(fmax*xa)
p(:) = alpha*f(:)
!--- > set h to identity matrix
h = 0.0
DO j = 1,n
h(j,j) = 1.0
ENDDO
ELSE
!------------------------------------------------------------
! ===> k-th step
!
!--- > determine p
p(:) = xnew(:)-xold(:)
!--- > update the change in gradients
y(:) = y(:)-f(:)
!--- > get necessary inner products and H|y>
py = dot_PRODUCT(p,y)
yy = 0.0
DO i = 1,n
s = 0.0
DO j = 1,n
s = s+y(j)*h(j,i)
ENDDO
v(i) = s
yy = yy+y(i)*s
ENDDO
!--- > check that update will leave h positive definite;
!--- > if not, then stop
IF (py<=0.0) THEN
WRITE (6,*) ' bfgs: <p|y> < 0'
WRITE (6,*) ' check convergence of forces'
CALL juDFT_error("bfgs: <p|y><0",calledby="bfgs")
ELSE
!--- > update h
IF (ist==2) THEN
gamma = py/yy
ELSE
gamma = 1.0
ENDIF
DO j = 1,n
DO i = 1,n
h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma
+ + (1.+gamma*yy/py)*p(i)*p(j)/py
ENDDO
ENDDO
!--- > generate p
DO i = 1,n
s = 0.0
DO j = 1,n
s = s+f(j)*h(j,i)
ENDDO
p(i) = s
ENDDO
ENDIF
ENDIF
!-------------------------------------------------------------
!--- > put xold and y in the correct form for the next step
DO i = 1,n
xold(i) = xnew(i)
y(i) = f(i)
ENDDO
!--- > if displacements are all less than epsdisp, then converged
dispmax = 0.0
DO na = 1,ntype
nn = 3*(na-1)
d2 = p(nn+1)**2 + p(nn+2)**2 + p(nn+3)**2
dispmax = MAX( dispmax, d2)