Commit 1e77c05a authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'mixer' into 'develop'

Refactoring of Mixer finished

See merge request fleur/fleur!25
parents 13adb7a0 addfc05e
......@@ -59,7 +59,7 @@ kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F glo
)
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90
eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90 mpi/mpi_bc_tool.F90
global/sort.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
inpgen/closure.f90 inpgen/inpgen_arguments.F90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 juDFT/string.f90 io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90
......
......@@ -4,6 +4,8 @@ Contributors guide
Everyone is very welcome to contribute to the enhancement of FLEUR.
Please use the [gitlab service] (https://iffgit.fz-juelich.de/fleur/fleur) to obtain the
latest development version of FLEUR.
- The
##Coding rules for FLEUR:
In no particular order we try to collect items to consider when writing code for FLEUR
......
MODULE m_tlmplm_cholesky
use m_judft
IMPLICIT NONE
!*********************************************************************
! sets up the local Hamiltonian, i.e. the Hamiltonian in the
......
......@@ -21,7 +21,6 @@ io/cdnpot_io_common.F90
io/cdn_io.F90
io/pot_io.F90
io/banddos_io.F90
io/broyd_io.F90
io/rw_inp.f90
io/rw_noco.f90
io/r_inpXML.F90
......
......@@ -6,6 +6,7 @@
MODULE m_io_matrix
USE m_types
USE m_judft
IMPLICIT NONE
private
INTEGER:: mode=1
......
......@@ -330,6 +330,14 @@ CONTAINS
input%imix = 5
CASE ('Anderson')
input%imix = 7
CASE ("Pulay")
input%imix = 9
CASE ("pPulay")
input%imix = 11
CASE ("rPulay")
input%imix = 13
CASE ("aPulay")
input%imix = 15
CASE DEFAULT
STOP 'Error: unknown mixing scheme selected!'
END SELECT
......
......@@ -885,6 +885,10 @@
<xsd:enumeration value="Broyden1"/>
<xsd:enumeration value="Broyden2"/>
<xsd:enumeration value="Anderson"/>
<xsd:enumeration value="Pulay"/>
<xsd:enumeration value="rPulay"/>
<xsd:enumeration value="pPulay"/>
<xsd:enumeration value="aPulay"/>
</xsd:restriction>
</xsd:simpleType>
......
......@@ -41,7 +41,7 @@ CONTAINS
USE m_fleur_init
USE m_optional
USE m_cdn_io
USE m_broyd_io
USE m_mixing_history
USE m_qfix
USE m_vgen
USE m_writexcstuff
......@@ -212,7 +212,7 @@ CONTAINS
cell,oneD,enpara,results,sym,xcpot,vTot,iter,iterHF)
END SELECT
IF(hybrid%l_calhf) THEN
CALL system("rm broyd*")
call mixing_history_reset(mpi)
iter = 0
END IF
ENDIF
......@@ -411,10 +411,9 @@ CONTAINS
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")
CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
stars,atoms,sphhar,vacuum,input,&
sym,cell,noco,oneD,archiveType,inDen,outDen,results)
IF(mpi%irank == 0) THEN
WRITE (6,FMT=8130) iter
......@@ -485,7 +484,7 @@ CONTAINS
CLOSE(2)
PRINT *,"qfix set to F"
ENDIF
CALL resetBroydenHistory()
call mixing_history_reset(mpi)
ENDIF
CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank)
END SUBROUTINE priv_geo_end
......
......@@ -13,7 +13,7 @@ MODULE m_fleur_arguments
CHARACTER(len=200) :: values
END TYPE t_fleur_param
INTEGER,PARAMETER:: no_params=24
INTEGER,PARAMETER:: no_params=25
TYPE(t_fleur_param) :: fleur_param(no_params)=(/&
!Input options
t_fleur_param(0,"-toXML","Convert an old 'inp' file into the new XML format",""),&
......@@ -58,6 +58,7 @@ MODULE m_fleur_arguments
t_fleur_param(0,"-trace","Try to generate a stacktrace in case of an error",""),&
t_fleur_param(0,"-debugtime","Write the start/stop of all timers to the console",""),&
!Output
t_fleur_param(0,"-mix_io","Do not store mixing history in memory but do IO in each iteration",""),&
t_fleur_param(0,"-no_out","Do not open the 'out' file but write to stdout",""),&
t_fleur_param(0,"-genEnpara","Generate an 'enpara' file for the energy parameters",""),&
t_fleur_param(0,"-kpts_gw","add alternative k point set for GW in all outputs for the XML input file",""),&
......
......@@ -28,7 +28,7 @@
USE m_setupMPI
USE m_cdn_io
USE m_fleur_info
USE m_broyd_io
USE m_mixing_history
USE m_checks
USE m_prpqfftmap
USE m_writeOutHeader
......@@ -531,7 +531,7 @@
CALL results%init(dimension,input,atoms,kpts,noco)
IF (mpi%irank.EQ.0) THEN
IF(input%gw.NE.0) CALL resetBroydenHistory()
IF(input%gw.NE.0) CALL mixing_history_reset(mpi)
CALL setStartingDensity(noco%l_noco)
END IF
......
......@@ -3,7 +3,7 @@
! 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:
......@@ -15,332 +15,172 @@ module m_mix
contains
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"
SUBROUTINE mix_charge( field, DIMENSION, mpi, l_writehistory,&
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, archiveType, inDen, outDen, results )
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
use m_VYukawaFilm
#ifdef CPP_MPI
use m_mpi_bc_potden
#endif
USE m_kerker
use m_pulay
use m_a_pulay
use m_types_mixvector
USE m_distance
use m_mixing_history
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_stars),TARGET,INTENT(in) :: stars
TYPE(t_cell),TARGET,INTENT(in) :: cell
TYPE(t_sphhar),TARGET,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_atoms),TARGET,INTENT(in) :: atoms
type(t_potden), intent(inout) :: outDen
type(t_results), intent(inout) :: results
type(t_potden), intent(inout) :: inDen
integer, intent(in) :: archiveType
LOGICAL, INTENT(IN) :: l_writehistory
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)
real :: fix
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
IF (noco%l_mtnocopot) mmap= mmap+ 2*atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd
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( 6, fmt='(a,2f10.5)' ) 'STRAIGHT MIXING',input%alpha
case( 3 )
write( 6, fmt='(a,f10.5)' ) 'BROYDEN FIRST MIXING',input%alpha
case( 5 )
write( 6, fmt='(a,f10.5)' ) 'BROYDEN SECOND MIXING',input%alpha
case( 7 )
write( 6, 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 )
!calculate the distance of charge densities 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 ( 6,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
ELSE
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 ( 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))
ELSE
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))
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 (smMet,fmMet)
CALL closeXMLElement('densityConvergence')
end if MPI0_a
TYPE(t_mixvector),ALLOCATABLE :: sm(:), fsm(:)
TYPE(t_mixvector) :: fsm_mag
LOGICAL :: l_densitymatrix
INTEGER :: it,maxiter
CALL timestart("Charge Density Mixing")
l_densitymatrix=.FALSE.
IF (atoms%n_u>0) THEN
l_densitymatrix=.NOT.input%ldaulinmix
IF (mpi%irank==0) CALL u_mix(input,atoms,inDen%mmpMat,outDen%mmpMat)
IF (ALL(inDen%mmpMat==0.0)) THEN
l_densitymatrix=.FALSE.
inDen%mmpMat=outDen%mmpMat
if (mpi%irank.ne.0) inden%mmpmat=0.0
ENDIF
ENDIF
CALL timestart("Reading of distances")
CALL mixvector_init(mpi%mpi_comm,l_densitymatrix,oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms)
CALL mixing_history_open(mpi,input%maxiter)
maxiter=MERGE(1,input%maxiter,input%imix==0)
CALL mixing_history(input%imix,maxiter,inden,outden,sm,fsm,it)
CALL distance(mpi%irank,cell%vol,input%jspins,fsm(it),inDen,outDen,results,fsm_Mag)
CALL timestop("Reading of distances")
! KERKER PRECONDITIONER
if( input%preconditioning_param /= 0 ) then
CALL resDen%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
CALL vYukawa%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 )
MPI0_b: if( mpi%irank == 0 ) then
call resDen%subPotDen( 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
if ( .not. input%film ) then
call vgen_coulomb( 1, mpi, dimension, oneD, input, field, vacuum, sym, stars, cell, &
sphhar, atoms, resDen, vYukawa )
else
vYukawa%iter = resDen%iter
call VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, resDen, &
vYukawa )
end if
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
resDen%vacz = resDen%vacz - input%preconditioning_param ** 2 / fpi_const * vYukawa%vacz
resDen%vacxy = resDen%vacxy - input%preconditioning_param ** 2 / fpi_const * vYukawa%vacxy
if( input%jspins == 2 ) call resDen%ChargeAndMagnetisationToSpins()
! fix the preconditioned density
call outDen%addPotDen( resDen, inDen )
call qfix(mpi,stars, atoms, sym, vacuum, sphhar, input, cell, oneD, outDen, noco%l_noco, .false., .true., fix )
call resDen%subPotDen( outDen, inDen )
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, resDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
end if
! end of preconditioner
IF( input%preconditioning_param /= 0 ) THEN
CALL timestart("Preconditioner")
CALL kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm(it) )
!Store modified density in history
CALL mixing_history_store(fsm(it))
CALL timestop("Preconditioner")
END IF
CALL timestart("Mixing")
!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>
! 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)
DEALLOCATE (sm,fsm)
!fix charge of the new density
CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false., fix)
IF(atoms%n_u.NE.n_u_keep) THEN
inDen%mmpMat = n_mmpTemp
END IF
atoms%n_u=n_u_keep
IF(vacuum%nvac.EQ.1) THEN
inDen%vacz(:,2,:) = inDen%vacz(:,1,:)
IF (sym%invs) THEN
inDen%vacxy(:,:,2,:) = CONJG(inDen%vacxy(:,:,1,:))
ELSE
inDen%vacxy(:,:,2,:) = inDen%vacxy(:,:,1,:)
END IF
END IF
IF (atoms%n_u > 0) THEN
IF (.NOT.l_densityMatrixPresent) THEN
inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
CALL resetBroydenHistory()
END IF
ENDIF
!write out mixed density
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen)
SELECT CASE(input%imix)
CASE(0)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,f10.5)' ) &
'STRAIGHT MIXING: alpha=',input%alpha," spin-mixing=",MERGE(input%alpha*input%spinf,0.,input%jspins>1)
CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it))
CASE(3,5)
CALL judft_error("Broyden 1/2 method not implemented")
CASE(7)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0)' ) &
'GENERALIZED ANDERSON MIXING: alpha=',input%alpha," History-length=",it-1
Call broyden(input%alpha,fsm,sm)
CASE(9)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,0)
CASE(11)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'PERIODIC PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,input%maxiter)
CASE(13)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'RESTARTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,0)
IF (it==input%maxiter) CALL mixing_history_limit(0) !Restarting Pulay
CASE(15)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'ADAPTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL a_pulay(input%alpha,fsm,sm)
CASE DEFAULT
CALL judft_error("Unknown Mixing schema")
END SELECT
CALL timestop("Mixing")
CALL timestart("Postprocessing")
!extracte mixed density
inDen%pw=0.0;inDen%mt=0.0
IF (ALLOCATED(inDen%vacz)) inden%vacz=0.0
IF (ALLOCATED(inDen%vacxy)) inden%vacxy=0.0
IF (ALLOCATED(inDen%mmpMat).AND.l_densitymatrix) inden%mmpMat=0.0
CALL sm(it)%to_density(inDen)
IF (atoms%n_u>0.AND..NOT.l_densitymatrix.AND..NOT.input%ldaulinmix) THEN
!No density matrix was present
!but is now created...
inden%mmpMAT=outden%mmpMat
CALL mixing_history_reset(mpi)
CALL mixvector_reset()
ENDIF
!fix charge of the new density
IF (mpi%irank==0) CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.FALSE., fix)
IF(vacuum%nvac.EQ.1) THEN
inDen%vacz(:,2,:) = inDen%vacz(:,1,:)
IF (sym%invs) THEN
inDen%vacxy(:,:,2,:) = CONJG(inDen%vacxy(:,:,1,:))
ELSE
inDen%vacxy(:,:,2,:) = inDen%vacxy(:,:,1,:)
END IF
END IF
!write out mixed density
IF (mpi%irank==0) CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen)
#ifdef CPP_HDF
IF (judft_was_argument("-last_extra")) THEN
CALL system("rm cdn_last.hdf")
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen,'cdn_last')
IF (mpi%irank==0.and.judft_was_argument("-last_extra")) THEN
CALL system("rm cdn_last.hdf")
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen,'cdn_last')
END IF
END IF