Commit 6645cae4 authored by Miriam Hinzen's avatar Miriam Hinzen

Merge branch 'kerker' into 'develop'

Merge branch kerker into develop

See merge request fleur/fleur!6
parents dbcba7f2 107ac463
......@@ -24,6 +24,7 @@ MODULE m_constants
INTEGER, PARAMETER :: POTDEN_TYPE_POTTOT = 1 ! 0 < POTDEN_TYPE <= 1000 ==> potential
INTEGER, PARAMETER :: POTDEN_TYPE_POTCOUL = 2
INTEGER, PARAMETER :: POTDEN_TYPE_POTX = 3
INTEGER, PARAMETER :: POTDEN_TYPE_POTYUK = 4
INTEGER, PARAMETER :: POTDEN_TYPE_DEN = 1001 ! 1000 < POTDEN_TYPE ==> density
CHARACTER(2),DIMENSION(0:103),PARAMETER :: namat_const=(/&
......
......@@ -146,7 +146,8 @@
input%pallst = .false. ; obsolete%lwb = .false. ; vacuum%starcoeff = .false.
input%strho = .false. ; input%l_f = .false. ; atoms%l_geo(:) = .true.
noco%l_noco = noco%l_ss ; input%jspins = 1
input%itmax = 9 ; input%maxiter = 99 ; input%imix = 7 ; input%alpha = 0.05 ; input%minDistance = 0.0
input%itmax = 9 ; input%maxiter = 99 ; input%imix = 7 ; input%alpha = 0.05
input%preconditioning_param = 0.0 ; input%minDistance = 0.0
input%spinf = 2.0 ; obsolete%lepr = 0 ; input%coretail_lmax = 0
sliceplot%kk = 0 ; sliceplot%nnne = 0 ; vacuum%nstars = 0 ; vacuum%nstm = 0
input%isec1 = 99 ; nu = 5 ; vacuum%layerd = 1 ; iofile = 6
......
......@@ -330,6 +330,7 @@ SUBROUTINE r_inpXML(&
END SELECT
input%alpha = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@alpha'))
input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@preconditioning_param'))
input%spinf = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@spinf'))
! Get parameters for core electrons
......
......@@ -563,17 +563,17 @@
IF( check .eq. ',' ) THEN
READ (UNIT=5,FMT=8060,END=99,ERR=99) &
& input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
WRITE (6,9180) input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
& input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
WRITE (6,9180) input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
8060 FORMAT (6x,i2,9x,i3,6x,i2,7x,f6.2,7x,f6.2)
ELSE
READ (UNIT=5,FMT=8061,END=99,ERR=99) &
& input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
WRITE (6,9180) input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
& input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
WRITE (6,9180) input%itmax,input%maxiter,input%imix,input%alpha,input%spinf
8061 FORMAT (6x,i3,9x,i3,6x,i2,7x,f6.2,7x,f6.2)
END IF
input%preconditioning_param = 0.0
chform = '(5x,l1,'//chntype//'f6.2)'
! chform = '(5x,l1,23f6.2)'
......
......@@ -168,8 +168,8 @@ SUBROUTINE w_inpXML(&
110 FORMAT(' <cutoffs Kmax="',f0.8,'" Gmax="',f0.8,'" GmaxXC="',f0.8,'" numbands="',i0,'"/>')
WRITE (fileNum,110) input%rkmax,stars%gmaxInit,xcpot%gmaxxc,input%gw_neigd
! <scfLoop itmax="9" maxIterBroyd="99" imix="Anderson" alpha="0.05" spinf="2.00"/>
120 FORMAT(' <scfLoop itmax="',i0,'" minDistance="',f0.8,'" maxIterBroyd="',i0,'" imix="',a,'" alpha="',f0.8,'" spinf="',f0.8,'"/>')
! <scfLoop itmax="9" maxIterBroyd="99" imix="Anderson" alpha="0.05" preconditioning_param="0.0" spinf="2.00"/>
120 FORMAT(' <scfLoop itmax="',i0,'" minDistance="',f0.8,'" maxIterBroyd="',i0,'" imix="',a,'" alpha="',f0.8,'" preconditioning_param="',f3.1,'" spinf="',f0.8,'"/>')
SELECT CASE (input%imix)
CASE (1)
mixingScheme='straight'
......@@ -182,7 +182,7 @@ SUBROUTINE w_inpXML(&
CASE DEFAULT
mixingScheme='errorUnknownMixing'
END SELECT
WRITE (fileNum,120) input%itmax,input%minDistance,input%maxiter,TRIM(mixingScheme),input%alpha,input%spinf
WRITE (fileNum,120) input%itmax,input%minDistance,input%maxiter,TRIM(mixingScheme),input%alpha,input%preconditioning_param,input%spinf
! <coreElectrons ctail="T" frcor="F" kcrel="0"/>
130 FORMAT(' <coreElectrons ctail="',l1,'" frcor="',l1,'" kcrel="',i0,'" coretail_lmax="',i0,'"/>')
......
......@@ -537,6 +537,7 @@
<xsd:attribute default="99" name="maxIterBroyd" type="xsd:nonNegativeInteger" use="optional"/>
<xsd:attribute name="imix" type="MixingEnum" use="required"/>
<xsd:attribute name="alpha" type="xsd:string" use="required"/>
<xsd:attribute default="0.0" name="preconditioning_param" type="xsd:string" use="optional"/>
<xsd:attribute default="2.0" name="spinf" type="xsd:string" use="optional"/>
<xsd:attribute default="0.0" name="minDistance" type="xsd:string" use="optional"/>
<xsd:attribute name="maxTimeToStartIter" type="xsd:string" use="optional"/>
......
......@@ -799,10 +799,11 @@
<xsd:attribute name="itmax" type="xsd:positiveInteger" use="required"/>
<xsd:attribute default="99" name="maxIterBroyd" type="xsd:nonNegativeInteger" use="optional"/>
<xsd:attribute name="imix" type="MixingEnum" use="required"/>
<xsd:attribute name="alpha" type="xsd:double" use="required"/>
<xsd:attribute default="2.0" name="spinf" type="xsd:double" use="optional"/>
<xsd:attribute default="0.0" name="minDistance" type="xsd:double" use="optional"/>
<xsd:attribute name="maxTimeToStartIter" type="xsd:double" use="optional"/>
<xsd:attribute name="alpha" type="xsd:string" use="required"/>
<xsd:attribute default="0.0" name="preconditioning_param" type="xsd:string" use="optional"/>
<xsd:attribute default="2.0" name="spinf" type="xsd:string" use="optional"/>
<xsd:attribute default="0.0" name="minDistance" type="xsd:string" use="optional"/>
<xsd:attribute name="maxTimeToStartIter" type="xsd:string" use="optional"/>
</xsd:complexType>
<xsd:complexType name="VacDOSType">
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -76,7 +76,7 @@ CONTAINS
! Types, these variables contain a lot of data!
TYPE(t_input) :: input
TYPE(t_field) :: field
TYPE(t_field) :: field, field2
TYPE(t_dimension):: DIMENSION
TYPE(t_atoms) :: atoms
TYPE(t_sphhar) :: sphhar
......@@ -112,13 +112,15 @@ CONTAINS
#endif
mpi%mpi_comm = mpi_comm
CALL timestart("Initialization")
CALL fleur_init(mpi,input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
sliceplot,banddos,obsolete,enpara,xcpot,results,kpts,hybrid,&
oneD,coreSpecInput,wann,l_opti)
CALL timestop("Initialization")
if( input%preconditioning_param /= 0 .and. input%film ) call juDFT_error('Currently no preconditioner for films', calledby = 'fleur' )
IF (l_opti) CALL optional(mpi,atoms,sphhar,vacuum,dimension,&
stars,input,sym,cell,sliceplot,obsolete,xcpot,noco,oneD)
......@@ -236,13 +238,13 @@ CONTAINS
!---< gwf
CALL timestart("generation of potential")
CALL vgen(hybrid,field,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,&
sym,obsolete,cell, oneD,sliceplot,mpi ,results,noco,inDen,vTot,vx,vCoul)
CALL vgen( hybrid, field, input, xcpot, DIMENSION, atoms, sphhar, stars, vacuum, &
sym, obsolete, cell, oneD, sliceplot, mpi, results, noco, inDen, vTot, vx, &
vCoul )
CALL timestop("generation of potential")
#ifdef CPP_MPI
CALL MPI_BARRIER(mpi%mpi_comm,ierr)
#endif
......@@ -251,7 +253,6 @@ CONTAINS
forcetheoloop:DO WHILE(forcetheo%next_job(it==input%itmax,noco))
CALL timestart("generation of hamiltonian and diagonalization (total)")
CALL timestart("eigen")
vTemp = vTot
......@@ -418,19 +419,21 @@ CONTAINS
CALL forcetheo%postprocess()
CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
IF (mpi%irank.EQ.0) THEN
! ----> mix input and output densities
CALL timestart("mixing")
CALL mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,hybrid,archiveType,inDen,outDen,results)
CALL timestop("mixing")
WRITE (6,FMT=8130) it
WRITE (16,FMT=8130) it
8130 FORMAT (/,5x,'******* it=',i3,' is completed********',/,/)
WRITE(*,*) "Iteration:",it," Distance:",results%last_distance
CALL timestop("Iteration")
!+t3e
ENDIF ! mpi%irank.EQ.0
field2 = field
! ----> mix input and output densities
CALL timestart("mixing")
CALL mix( field2, xcpot, dimension, obsolete, sliceplot, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, hybrid, archiveType, inDen, outDen, results )
CALL timestop("mixing")
if( mpi%irank == 0 ) then
WRITE (6,FMT=8130) it
WRITE (16,FMT=8130) it
8130 FORMAT (/,5x,'******* it=',i3,' is completed********',/,/)
WRITE(*,*) "Iteration:",it," Distance:",results%last_distance
CALL timestop("Iteration")
end if ! mpi%irank.EQ.0
#ifdef CPP_MPI
......
......@@ -250,6 +250,7 @@
CALL MPI_BCAST(input%jspins,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
CALL MPI_BCAST(atoms%n_u,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
CALL MPI_BCAST(atoms%lmaxd,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
call MPI_BCAST( input%preconditioning_param, 1, MPI_DOUBLE, 0, mpi%mpi_comm, ierr )
#endif
CALL ylmnorm_init(atoms%lmaxd)
!
......
......@@ -3,285 +3,331 @@
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mix
module m_mix
!*************************************************************************
!------------------------------------------------------------------------
! mixing of charge densities or potentials:
! IMIX= 0 : linear mixing
! IMIX = 3 : BROYDEN'S FIRST METHOD
! IMIX = 5 : BROYDEN'S SECOND METHOD
! IMIX = 7 : GENERALIZED ANDERSEN METHOD
!************************************************************************
! IMIX = 0 : linear mixing
! IMIX = 3 : Broyden's First method
! IMIX = 5 : Broyden's Second method
! IMIX = 7 : Generalized Anderson method
!------------------------------------------------------------------------
CONTAINS
contains
SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
hybrid,archiveType,inDen,outDen,results)
subroutine mix( field, xcpot, dimension, obsolete, sliceplot, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, hybrid, archiveType, inDen, outDen, results )
#include"cpp_double.h"
USE m_juDFT
USE m_constants
USE m_cdn_io
USE m_broyd_io
USE m_brysh1
USE m_stmix
USE m_broyden
USE m_broyden2
USE m_brysh2
USE m_metric
USE m_qfix
USE m_types
USE m_xmlOutput
USE m_umix
IMPLICIT NONE
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_hybrid),INTENT(IN) :: hybrid
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(INOUT) :: atoms !n_u is modified temporarily
TYPE(t_potden),INTENT(INOUT) :: outDen
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_potden),INTENT(INOUT) :: inDen
INTEGER, INTENT(IN) :: archiveType
!Local Scalars
REAL fix,intfac,vacfac
INTEGER i,imap,js
INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
INTEGER iofl,n_u_keep
LOGICAL l_exist,l_ldaU, l_densityMatrixPresent, l_pot
!Local Arrays
REAL dist(6)
REAL, ALLOCATABLE :: sm(:), fsm(:), fmMet(:), smMet(:)
CHARACTER(LEN=20) :: attributes(2)
COMPLEX :: n_mmpTemp(-3:3,-3:3,MAX(1,atoms%n_u),input%jspins)
!External functions
REAL CPP_BLAS_sdot
EXTERNAL CPP_BLAS_sdot
! YM: I have exported 'vol' from outside, be aware
! IF (film) THEN
! vol = 2.0 * z1 * area
! ELSE
! vol = omtil
! ENDIF
l_densityMatrixPresent = ANY(inDen%mmpMat(:,:,:,:).NE.0.0)
!In systems without inversions symmetry the interstitial star-
!coefficients are complex. Thus twice as many numbers have to be
!stored.
intfac = 2.0
IF (sym%invs) intfac = 1.0
!The corresponding is true for the coeff. of the warping vacuum
!density depending on the two dimensional inversion.
vacfac = 2.0
IF (sym%invs2) vacfac = 1.0
mmaph = intfac*stars%ng3 + atoms%ntype*(sphhar%nlhd+1)*atoms%jmtd +&
vacfac*vacuum%nmzxyd*(oneD%odi%n2d-1)*vacuum%nvac + vacuum%nmzd*vacuum%nvac
mmap = mmaph*input%jspins
!in a non-collinear calculations extra space is needed for the
!off-diag. part of the density matrix. these coeff. are generally
!complex independ of invs and invs2.
IF (noco%l_noco) THEN
mmap = mmap + 2*stars%ng3 + 2*vacuum%nmzxyd*(oneD%odi%n2d-1)*vacuum%nvac + &
2*vacuum%nmzd*vacuum%nvac
END IF
! LDA+U (start)
n_mmpTemp = inDen%mmpMat
n_u_keep=atoms%n_u
IF (atoms%n_u.GT.0) CALL u_mix(input,atoms,inDen%mmpMat,outDen%mmpMat)
IF (l_densityMatrixPresent) THEN
!In an LDA+U caclulation, also the density matrix is included in the
!supervectors (sm,fsm) if no linear mixing is performed on it.
IF (input%ldauLinMix) THEN
atoms%n_u = 0
use m_juDFT
use m_constants
use m_cdn_io
use m_broyd_io
use m_brysh1
use m_stmix
use m_broyden
use m_broyden2
use m_brysh2
use m_metric
use m_qfix
use m_types
use m_xmlOutput
use m_umix
use m_vgen_coulomb
#ifdef CPP_MPI
use m_mpi_bc_potden
#endif
implicit none
type(t_oneD), intent(in) :: oneD
type(t_hybrid), intent(in) :: hybrid
type(t_input), intent(in) :: input
type(t_vacuum), intent(in) :: vacuum
type(t_noco), intent(in) :: noco
type(t_sym), intent(in) :: sym
type(t_stars), intent(in) :: stars
type(t_cell), intent(in) :: cell
type(t_sphhar), intent(in) :: sphhar
type(t_field), intent(inout) :: field
class(t_xcpot), intent(in) :: xcpot
type(t_dimension), intent(in) :: dimension
type(t_obsolete), intent(in) :: obsolete
type(t_sliceplot), intent(in) :: sliceplot
type(t_mpi), intent(in) :: mpi
type(t_atoms), intent(inout) :: atoms !n_u is modified temporarily
type(t_potden), intent(inout) :: outDen
type(t_results), intent(inout) :: results
type(t_potden), intent(inout) :: inDen
integer, intent(in) :: archiveType
real :: fix, intfac, vacfac
integer :: i, imap, js, n, lh
integer :: mmap, mmaph, nmaph, nmap, mapmt, mapvac, mapvac2
integer :: iofl, n_u_keep
logical :: l_exist, l_ldaU, l_densityMatrixPresent, l_pot
real :: dist(6)
real, allocatable :: sm(:), fsm(:), fmMet(:), smMet(:)
character(len=20) :: attributes(2)
complex :: n_mmpTemp(-3:3,-3:3,max(1,atoms%n_u),input%jspins)
type(t_potden) :: resDen, vYukawa
integer :: ierr(2)
!External functions
real :: CPP_BLAS_sdot
external :: CPP_BLAS_sdot
! YM: I have exported 'vol' from outside, be aware
! IF (film) THEN
! vol = 2.0 * z1 * area
! ELSE
! vol = omtil
! ENDIF
MPI0_a: if( mpi%irank == 0 ) then
l_densityMatrixPresent = any( inDen%mmpMat(:,:,:,:) /= 0.0 )
!In systems without inversions symmetry the interstitial star-
!coefficients are complex. Thus twice as many numbers have to be
!stored.
intfac = 2.0
if ( sym%invs ) intfac = 1.0
!The corresponding is true for the coeff. of the warping vacuum
!density depending on the two dimensional inversion.
vacfac = 2.0
if ( sym%invs2 ) vacfac = 1.0
mmaph = intfac * stars%ng3 + atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd + &
vacfac * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + vacuum%nmzd * vacuum%nvac
mmap = mmaph * input%jspins
!in a non-collinear calculations extra space is needed for the
!off-diag. part of the density matrix. these coeff. are generally
!complex independ of invs and invs2.
if ( noco%l_noco ) then
mmap = mmap + 2 * stars%ng3 + 2 * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + &
2 * vacuum%nmzd * vacuum%nvac
end if
! LDA+U (start)
n_mmpTemp = inDen%mmpMat
n_u_keep = atoms%n_u
if ( atoms%n_u > 0 ) call u_mix( input, atoms, inDen%mmpMat, outDen%mmpMat )
if ( l_densityMatrixPresent ) then
!In an LDA+U caclulation, also the density matrix is included in the
!supervectors (sm,fsm) if no linear mixing is performed on it.
if ( input%ldauLinMix ) then
atoms%n_u = 0
else
mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
end if
else
atoms%n_u = 0
end if
! LDA+U (end)
allocate( sm(mmap), fsm(mmap) )
allocate( smMet(mmap), fmMet(mmap) )
dist(:) = 0.0
!determine type of mixing:
!imix=0:straight, imix=o broyden first, imix=5:broyden second
!imix=:generalozed anderson mixing
select case( input%imix )
case( 0 )
write( 16, fmt='(a,2f10.5)' ) 'STRAIGHT MIXING',input%alpha
case( 3 )
write( 16, fmt='(a,f10.5)' ) 'BROYDEN FIRST MIXING',input%alpha
case( 5 )
write( 16, fmt='(a,f10.5)' ) 'BROYDEN SECOND MIXING',input%alpha
case( 7 )
write( 16, fmt='(a,f10.5)' ) 'ANDERSON GENERALIZED',input%alpha
case default
call juDFT_error( "mix: input%imix =/= 0,3,5,7 ", calledby ="mix" )
end select
if ( input%jspins == 2 .and. input%imix /= 0 ) then
write( 6, '(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')' )
end if
!put input charge density into array sm
!(in the spin polarized case the arrays sm and fsm consist of spin up and spin down densities)
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, inDen, nmap, nmaph, mapmt, mapvac, mapvac2, sm )
!put output charge density into array fsm
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, outDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
!store the difference fsm - sm in fsm
fsm(:nmap) = fsm(:nmap) - sm(:nmap)
l_pot = .false.
! Apply metric w to fsm and store in fmMet: w |fsm>
call metric( cell, atoms, vacuum, sphhar, input, noco, stars, sym, oneD, &
mmap, nmaph, mapmt, mapvac2, fsm, fmMet, l_pot )
end if MPI0_a
! KERKER PRECONDITIONER
if( input%preconditioning_param /= 0 ) then
call resDen%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1001 )
call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 4 )
MPI0_b: if( mpi%irank == 0 ) then
call resDen%Residual( outDen, inDen )
if( input%jspins == 2 ) call resDen%SpinsToChargeAndMagnetisation()
end if MPI0_b
#ifdef CPP_MPI
call mpi_bc_potden( mpi, stars, sphhar, atoms, input, vacuum, oneD, noco, resDen )
#endif
call vgen_coulomb( 1, mpi, dimension, oneD, input, field, vacuum, sym, stars, cell, &
sphhar, atoms, resDen, vYukawa )
end if
MPI0_c: if( mpi%irank == 0 ) then
if( input%preconditioning_param /= 0 ) then
resDen%pw(1:stars%ng3,1) = resDen%pw(1:stars%ng3,1) - input%preconditioning_param ** 2 / fpi_const * vYukawa%pw(1:stars%ng3,1)
do n = 1, atoms%ntype
do lh = 0, sphhar%nlhd
resDen%mt(1:atoms%jri(n),lh,n,1) = resDen%mt(1:atoms%jri(n),lh,n,1) &
- input%preconditioning_param ** 2 / fpi_const &
* vYukawa%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2
end do
end do
if( input%jspins == 2 ) call resDen%ChargeAndMagnetisationToSpins()
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, resDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
end if
! end of preconditioner
!mixing of the densities
IF (input%imix.EQ.0) THEN
CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
ELSE
mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
! Replace the broyden call above by the commented metric and broyden2 calls
! below to switch on the continuous restart of the Broyden method.
! Apply metric w to sm and store in smMet: w |sm>
! CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
! mmap,nmaph,mapmt,mapvac2,sm,smMet,l_pot)
!
! CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
! hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm,fmMet,smMet)
END IF
ELSE
atoms%n_u = 0
END IF
! LDA+U (end)
ALLOCATE (sm(mmap),fsm(mmap))
ALLOCATE (smMet(mmap),fmMet(mmap))
dist(:) = 0.0
!determine type of mixing:
!imix=0:straight, imix=o broyden first, imix=5:broyden second
!imix=:generalozed anderson mixing
IF (input%imix.EQ.0) THEN
WRITE (16,FMT='(a,2f10.5)') 'STRAIGHT MIXING',input%alpha
ELSE IF (input%imix.EQ.3) THEN
WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
ELSE IF (input%imix.EQ.5) THEN
WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
ELSE IF (input%imix.EQ.7) THEN
WRITE (16,FMT='(a,f10.5)') 'ANDERSON GENERALIZED',input%alpha
ELSE
CALL juDFT_error("mix: input%imix =/= 0,3,5,7 ",calledby ="mix")
END IF
IF (input%jspins.EQ.2.AND.input%imix.NE.0) THEN
WRITE(6,'(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')')
END IF
!put input charge density into array sm
!(in the spin polarized case the arrays sm and fsm consist of spin up and spin down densities)
CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,inDen,nmap,nmaph,mapmt,mapvac,mapvac2,sm)
!put output charge density into array fsm
CALL brysh1(input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
intfac,vacfac,outDen,nmap,nmaph,mapmt,mapvac,mapvac2,fsm)
!store the difference fsm - sm in fsm
fsm(:nmap) = fsm(:nmap) - sm(:nmap)
l_pot = .FALSE.
! Apply metric w to fsm and store in fmMet: w |fsm>
CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
mmap,nmaph,mapmt,mapvac2,fsm,fmMet,l_pot)
!mixing of the densities
IF (input%imix.EQ.0) THEN
CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
ELSE
CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
! Replace the broyden call above by the commented metric and broyden2 calls
! below to switch on the continuous restart of the Broyden method.
! Apply metric w to sm and store in smMet: w |sm>
!initiatlize mixed density and extract it with brysh2 call
inDen%mmpMat = CMPLX(0.0,0.0)
CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen)
!calculate the distance of charge densities...
!induce metric in fsm use sm as an output array: |sm> = w |fsm>
! CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
! mmap,nmaph,mapmt,mapvac2,sm,smMet,l_pot)
!
! CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
! hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm,fmMet,smMet)
END IF
!initiatlize mixed density and extract it with brysh2 call
inDen%mmpMat = CMPLX(0.0,0.0)
CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen)
!calculate the distance of charge densities...
!induce metric in fsm use sm as an output array: |sm> = w |fsm>
! CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
! mmap,nmaph,mapmt,mapvac2,fsm, sm)
!calculate the charge density distance for each spin
IF(hybrid%l_calhf) THEN
CALL openXMLElement('densityConvergence',(/'units ','comment'/),(/'me/bohr^3','HF '/))
ELSE
CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
END IF
DO js = 1,input%jspins
dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
CALL writeXMLElementForm('chargeDensity',(/'spin ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
IF( hybrid%l_calhf ) THEN
WRITE (16,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
WRITE ( 6,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
ELSE
WRITE (16,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
END IF
END DO
IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
!calculate the distance of total charge and spin density
!|rho/m(o) - rho/m(i)| = |rh1(o) -rh1(i)|+ |rh2(o) -rh2(i)| +/_
! +/_2<rh2(o) -rh2(i)|rh1(o) -rh1(i)>
IF (input%jspins.EQ.2) THEN
dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(4)/cell%vol))/),reshape((/10,20/),(/1,2/)))
CALL writeXMLElementFormPoly('spinDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(5)/cell%vol))/),reshape((/19,20/),(/1,2/)))
IF( hybrid%l_calhf ) THEN
WRITE (16,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE (16,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
WRITE ( 6,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
! mmap,nmaph,mapmt,mapvac2,fsm, sm)
!calculate the charge density distance for each spin
IF(hybrid%l_calhf) THEN
CALL openXMLElement('densityConvergence',(/'units ','comment'/),(/'me/bohr^3','HF '/))
ELSE
WRITE (16,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE (16,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
WRITE ( 6,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
END IF
!dist/vol should always be >= 0 ,
!but for dist=0 numerically you might obtain dist/vol < 0
!(e.g. when calculating non-magnetic systems with jspins=2).
END IF
results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
DEALLOCATE (sm,fsm,smMet,fmMet)
CALL closeXMLElement('densityConvergence')
DO js = 1,input%jspins
dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
CALL writeXMLElementForm('chargeDensity',(/'spin ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
IF( hybrid%l_calhf ) THEN
WRITE (16,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol