Commit 18490aa8 authored by Henning Janssen's avatar Henning Janssen

Merge branch 'develop' into DFT+HubbardI

parents e544d13c 1e77c05a
......@@ -60,7 +60,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
......
......@@ -332,6 +332,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
......
......@@ -898,6 +898,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
......@@ -215,7 +215,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
......@@ -414,10 +414,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
......@@ -488,7 +487,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
......@@ -542,7 +542,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
......
This diff is collapsed.
......@@ -2,11 +2,17 @@ set(fleur_F77 ${fleur_F77}
mix/metr_z0.f
)
set(fleur_F90 ${fleur_F90}
mix/broyden.F90
mix/broyden2.F90
mix/brysh1.f90
mix/brysh2.f90
mix/metric.f90
mix/type_mixvector.F90
mix/kerker.F90
mix/broyden_history.F90
mix/mixing_history.F90
mix/distance.F90
mix/pulay.F90
mix/a_pulay.F90
#mix/broyden2.F90
#mix/brysh1.f90
#mix/brysh2.f90
#mix/metric.f90
mix/potdis.f90
mix/stdmix.f90
mix/u_mix.f90
......
MODULE m_a_pulay
USE m_juDFT
PRIVATE
REAL :: distance(3)
PUBLIC a_pulay
CONTAINS
SUBROUTINE a_pulay(alpha,fm,sm)
USE m_pulay
USE m_types_mat
USE m_types_mixvector
USE m_mixing_history
IMPLICIT NONE
REAL,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
REAL :: fac(2)=(/1.0,0.5/)
INTEGER,parameter :: parallel_images=2
INTEGER,PARAMETER :: maxiter=4,mode=2
INTEGER :: hlen,image,i,ii,local_length
INTEGER :: local_hist(maxiter*parallel_images)
hlen=SIZE(sm) !Number of results in history
image=MOD(hlen-1,parallel_images)+1
PRINT *,"HI:",hlen,image
IF (hlen<parallel_images) THEN
!create inital parallel images by simple mixing of first result
sm(hlen)=sm(1)+(alpha*fac(image))*fm(1)
RETURN
ENDIF
SELECT CASE(mode)
CASE(1)
local_length=0
DO i=image,hlen,parallel_images
local_length=local_length+1
local_hist(local_length)=i
ENDDO
IF (local_length>maxiter) THEN !recreate to enable cross-image mixing
local_length=hlen+1-parallel_images
local_hist(1)=image
local_hist(2:local_length)=(/(i,i=parallel_images+1,hlen)/)
END IF
sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)))
IF (hlen==parallel_images*(maxiter+1)) CALL mixing_history_limit(parallel_images)
case(2)
local_length=hlen-image
local_hist(1:local_length)=(/(i,i=1,local_length)/)
local_length=local_length+1
local_hist(local_length)=hlen
sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)))
IF (hlen==parallel_images*(maxiter)) CALL mixing_history_limit(parallel_images)
END SELECT
PRINT *,"P:",local_hist(:local_length)
END SUBROUTINE a_pulay
FUNCTION simple_pulay(alpha,fm,sm)RESULT(sm_out)
USE m_types_mixvector
IMPLICIT NONE
REAL,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(IN) :: sm(:)
TYPE(t_mixvector) :: sm_out
TYPE(t_mixvector) :: mdf
REAL,ALLOCATABLE :: a(:,:),b(:),work(:)
INTEGER,ALLOCATABLE :: ipiv(:)
INTEGER :: n,nn,info,lwork,hlen
hlen=SIZE(sm)
IF (hlen==1) THEN
sm_out=sm(1)+alpha*fm(1)
RETURN
ENDIF
ALLOCATE(a(hlen+1,hlen+1),b(hlen+1),ipiv(hlen+1),work((hlen+1)**2))
DO n=1,hlen
mdf=fm(n)%apply_metric()
DO nn=1,n
a(n,nn)=mdf.dot.fm(nn)
a(nn,n)=a(n,nn)
ENDDO
ENDDO
a(:,hlen+1)=1.0
a(hlen+1,:)=1.0
a(hlen+1,hlen+1)=0.0
b=0.0;b(hlen+1)=1.
CALL DSYSV( 'Upper', hlen+1, 1, a, SIZE(a,1), ipiv, b, SIZE(b,1), work, SIZE(work), INFO )
sm_out=0.0*mdf
DO n=1,hlen
sm_out=sm_out+b(n)*(sm(n)+alpha*fm(n))
ENDDO
END FUNCTION SIMPLE_PULAY
END MODULE m_a_pulay
......@@ -11,79 +11,46 @@ MODULE m_broyden
! fm1 : output minus inputcharge density of iteration m-1
!################################################################
CONTAINS
SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,lpot)
#include"cpp_double.h"
USE m_metric
SUBROUTINE broyden(input,fm,sm)
USE m_types
USE m_broyd_io
USE m_types_mixvector
IMPLICIT NONE
TYPE(t_oneD),INTENT(IN) :: oneD
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(IN) :: atoms
TYPE(t_hybrid),INTENT(IN) :: hybrid
! Scalar Arguments
INTEGER, INTENT (IN) :: mmap,nmap
INTEGER, INTENT (IN) :: mapmt,mapvac2
LOGICAL,OPTIONAL,INTENT(IN) :: lpot
! Array Arguments
REAL, INTENT (IN) :: fm(nmap)
REAL, INTENT (INOUT) :: sm(nmap)
TYPE(t_input),INTENT(IN) :: input
TYPE(t_mixvector),INTENT(IN) :: fm
TYPE(t_mixvector),INTENT(INOUT) :: sm
! Local Scalars
INTEGER :: i,it,k,nit,iread,nmaph, mit
REAL :: bm,dfivi,fmvm,smnorm,vmnorm,alphan
LOGICAL :: l_pot, l_exist
LOGICAL :: l_exist
REAL, PARAMETER :: one=1.0, zero=0.0
! Local Arrays
REAL, ALLOCATABLE :: am(:)
REAL, ALLOCATABLE :: fm1(:),sm1(:),ui(:),um(:),vi(:),vm(:)
! External Functions
REAL CPP_BLAS_sdot
EXTERNAL CPP_BLAS_sdot
! External Subroutines
EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal
REAL,ALLOCATABLE :: am(:)
TYPE(t_mixvector) :: fm1,sm1,ui,um,vi,vm
dfivi = zero
l_pot = .FALSE.
IF (PRESENT(lpot)) l_pot = lpot
ALLOCATE (fm1(mmap),sm1(mmap),ui(mmap),um(mmap),vi(mmap),vm(mmap))
CALL fm1%alloc()
CALL sm1%alloc()
CALL ui%alloc()
CALL um%alloc()
CALL vi%alloc()
CALL vm%alloc()
ALLOCATE (am(input%maxiter+1))
fm1 = 0.0
sm1 = 0.0
ui = 0.0
um = 0.0
vi = 0.0
vm = 0.0
am = 0.0
mit = 0
l_exist = initBroydenHistory(input,hybrid,nmap) ! returns true if there already exists a Broyden history
IF(.NOT.l_exist) mit = 1
IF (mit.NE.1) THEN
IF (l_exist) THEN
! load input charge density (sm1) and difference of
! in and out charge densities (fm1) from previous iteration (m-1)
CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))
CALL readLastIterInAndDiffDen(mit,alphan,sm1,fm1)
IF (ABS(input%alpha-alphan) > 0.0001) THEN
WRITE (6,*) 'mixing parameter has been changed; reset'
WRITE (6,*) 'broyden algorithm or set alpha to',alphan
......@@ -92,8 +59,10 @@ CONTAINS
! generate F_m - F_(m-1) ... sm1
! and rho_m - rho_(m-1) .. fm1
sm1(1:nmap) = sm(1:nmap) - sm1(1:nmap)
fm1(1:nmap) = fm(1:nmap) - fm1(1:nmap)
sm1 = sm - sm1
fm1 = fm - fm1
ELSE
mit=1
END IF
! save F_m and rho_m for next iteration
......@@ -101,95 +70,50 @@ CONTAINS
IF (nit > input%maxiter+1) nit = 1
CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm)
IF (mit.EQ.1) THEN
IF (.NOT.l_exist) THEN
! update for rho for mit=1 is straight mixing
! sm = sm + alpha*fm
CALL CPP_BLAS_saxpy(nmap,input%alpha,fm,1,sm,1)
sm=sm+input%alpha*fm
ELSE
! |vi> = w|vi>
! loop to generate um : um = sm1 + alpha*fm1 - \sum <fm1|w|vi> ui
um(:nmap) = input%alpha * fm1(:nmap) + sm1(:nmap)
um = input%alpha * fm1 + sm1
iread = MIN(mit-1,input%maxiter+1)
DO it = 2, iread
CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)
am(it) = vi.dot.fm
! calculate um(:) = -am(it)*ui(:) + um(:)
CALL CPP_BLAS_saxpy(nmap,-am(it),ui,1,um,1)
um=um-am(it)*ui
WRITE(6,FMT='(5x,"<vi|w|Fm> for it",i2,5x,f10.6)')it,am(it)
END DO
IF (input%imix.EQ.3) THEN
!****************************************
! broyden's first method
!****************************************
! convolute drho(m) with the metric: |fm1> = w|sm1>
CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
mmap,nmaph,mapmt,mapvac2,sm1,fm1,l_pot)
! calculate the norm of sm1 : <sm1|w|sm1>
smnorm = CPP_BLAS_sdot(nmap,sm1,1,fm1,1)
! generate vm = alpha*sm1 - \sum <ui|w|sm1> vi
vm(:) = input%alpha * fm1(:)
DO it = 2,iread
CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)
! calculate vm(:) = -bm*vi(:) + vm
CALL CPP_BLAS_saxpy(nmap,-bm,vi,1,vm,1)
!write(6,FMT='(5x,"<ui|w|Fm> for it",i2,5x,f10.6)') it, bm
END DO
! complete evaluation of vm
! vmnorm = <um|w|sm1>-<sm1|w|sm1>
vmnorm = CPP_BLAS_sdot(nmap,fm1,1,um,1) - smnorm
! if (vmnorm.lt.tol_10) stop
CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
ELSE IF (input%imix.EQ.5) THEN
!****************************************
! broyden's second method
!****************************************
! multiply fm1 with metric matrix and store in vm: w |fm1>
CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)
! calculate the norm of fm1 and normalize vm it: vm = wfm1 / <fm1|w|fm1>
vmnorm = one / CPP_BLAS_sdot(nmap,fm1,1,vm,1)
CALL CPP_BLAS_sscal(nmap,vmnorm,vm,1)
ELSE IF (input%imix.EQ.7) THEN
IF (input%imix.EQ.7) THEN
!****************************************
! generalized anderson method
!****************************************
! calculate vm = alpha*wfm1 -\sum <fm1|w|vi> <fi1|w|vi><vi|
! convolute fm1 with the metrik and store in vm
CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)
vm=fm1%apply_metric()
DO it = 2,iread
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
CALL CPP_BLAS_saxpy(nmap,-am(it)*dfivi,vi,1,vm,1)
vm=vm-am(it)*dfivi*vi
END DO
vmnorm = CPP_BLAS_sdot(nmap,fm1,1,vm,1)
vmnorm=fm1.dot.vm
! if (vmnorm.lt.tol_10) stop
! calculate vm(:) = (1.0/vmnorm)*vm(:)
CALL CPP_BLAS_sscal(nmap,one/vmnorm,vm,1)
vm=(1.0/vmnorm)*vm
! save dfivi(mit) for next iteration
dfivi = vmnorm
ELSE
CALL judft_error("Only generalized Anderson implemented")
END IF
! write um,vm and dfivi on file broyd.?
......@@ -199,12 +123,10 @@ CONTAINS
! update rho(m+1)
! calculate <fm|w|vm>
fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1)
fmvm = vm.dot.fm
! calculate sm(:) = (1.0-fmvm)*um(:) + sm
CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1)
sm=sm+(one-fmvm)*um
END IF
DEALLOCATE (fm1,sm1,ui,um,vi,vm,am)
END SUBROUTINE broyden
END MODULE m_broyden
MODULE m_broyden
USE m_juDFT
!################################################################
! IMIX = 3 : BROYDEN'S FIRST METHOD
! IMIX = 5 : BROYDEN'S SECOND METHOD
! IMIX = 7 : GENERALIZED ANDERSEN METHOD
! sm : input charge density of iteration m
! afterwards update rho(m+1)
! fm : output minus input charge density of iteration m
! sm1 : input charge density of iteration m-1
! fm1 : output minus inputcharge density of iteration m-1
!################################################################
CONTAINS
SUBROUTINE broyden(alpha,fm,sm)
USE m_types
USE m_types_mixvector
IMPLICIT NONE
real,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
! Locals
INTEGER :: n,it,hlen
REAL :: fmvm,vmnorm
REAL,ALLOCATABLE :: am(:),dfivi(:)
TYPE(t_mixvector) :: fm1,sm1,ui,um,vi,vm
TYPE(t_mixvector),allocatable :: u_store(:),v_store(:)
hlen=size(fm)
IF (hlen<2) THEN !Do a simple mixing step
sm(hlen)=sm(hlen)+alpha*fm(hlen)
RETURN
ENDIF
ALLOCATE(u_store(hlen-2),v_store(hlen-2))
do it=1,hlen-2
call u_store(it)%alloc()
call v_store(it)%alloc()
enddo
CALL fm1%alloc()
CALL sm1%alloc()
CALL ui%alloc()
CALL um%alloc()
CALL vi%alloc()
CALL vm%alloc()
ALLOCATE (am(hlen-1),dfivi(hlen-1))
dfivi = 0.0
am = 0.0
DO n=2,hlen
sm1 = sm(n) - sm(n-1)
fm1 = fm(n) - fm(n-1)
! |vi> = w|vi>
! loop to generate um : um = sm1 + alpha*fm1 - \sum <fm1|w|vi> ui
um = alpha * fm1 + sm1
DO it = n-2,1,-1
ui=u_store(it)
vi=v_store(it)
am(it) = vi.dot.fm1
! calculate um(:) = -am(it)*ui(:) + um(:)
um=um-am(it)*ui
WRITE(6,FMT='(5x,"<vi|w|Fm> for it",i2,5x,f10.6)')it,am(it)
END DO
! calculate vm = alpha*wfm1 -\sum <fm1|w|vi> <fi1|w|vi><vi|
! convolute fm1 with the metrik and store in vm
vm=fm1%apply_metric()
DO it = n-2,1,-1
vi=v_store(it)
! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
vm=vm-am(it)*dfivi(it)*vi
END DO
vmnorm=fm1.dot.vm
! if (vmnorm.lt.tol_10) stop
! calculate vm(:) = (1.0/vmnorm)*vm(:)
vm=(1.0/vmnorm)*vm
! save dfivi(mit) for next iteration
dfivi(n-1) = vmnorm
IF (n<hlen) u_store(n-1)=um
IF (n<hlen) v_store(n-1)=vm
enddo
! update rho(m+1)
! calculate <fm|w|vm>
fmvm = vm.dot.fm(hlen)
! calculate sm(:) = (1.0-fmvm)*um(:) + sm
sm(hlen)=sm(hlen)+(1.0-fmvm)*um
END SUBROUTINE broyden
END MODULE m_broyden
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! 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_distance
contains
SUBROUTINE distance(irank,vol,jspins,fsm,inden,outden,results,fsm_mag)
use m_types
use m_types_mixvector
use m_xmlOutput
implicit none
integer,intent(in) :: irank,jspins
real,intent(in) :: vol
type(t_mixvector),INTENT(IN) :: fsm
TYPE(t_potden),INTENT(INOUT) :: inden,outden
TYPE(t_results),INTENT(INOUT) :: results
type(t_mixvector),INTENT(OUT) :: fsm_mag
integer ::js
REAL :: dist(6) !1:up,2:down,3:spinoff,4:total,5:magnet,6:noco
TYPE(t_mixvector)::fmMet
character(len=100)::attributes(2)
CALL fmMet%alloc()
IF (jspins==2) THEN
CALL fsm_mag%alloc()
! calculate Magnetisation-difference
CALL fsm_mag%from_density(outden,swapspin=.TRUE.)
call fmMet%from_density(inden,swapspin=.true.)
fsm_mag=fsm_mag-fmMet
ENDIF
! Apply metric w to fsm and store in fmMet: w |fsm>
fmMet=fsm%apply_metric()
dist(:) = 0.0
DO js = 1,jspins
dist(js) = fsm%multiply_dot_mask(fmMet,(/.true.,.true.,.true.,.false./),js)
END DO
IF (SIZE(outden%pw,2)>2) dist(6) = fsm%multiply_dot_mask(fmMet,(/.TRUE.,.TRUE.,.TRUE.,.FALSE./),3)
IF (jspins.EQ.2) THEN
dist(3) = fmMet%multiply_dot_mask(fsm_mag,(/.true.,.true.,.true.,.false./),1)
dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
ENDIF
results%last_distance=maxval(1000*SQRT(ABS(dist/vol)))
if (irank>0) return
!calculate the distance of charge densities for each spin
CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
DO js = 1,jspins
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/vol))
CALL writeXMLElementForm('chargeDensity',(/'spin ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/vol))
END DO
IF (abs(dist(6))>1E-15) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/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 (jspins.EQ.2) THEN
CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(4)/vol))/),reshape((/10,20/),(/1,2/)))
CALL writeXMLElementFormPoly('spinDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(5)/vol))/),reshape((/19,20/),(/1,2/)))
WRITE ( 6,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/vol))
WRITE ( 6,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/vol))
!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
CALL closeXMLElement('densityConvergence')
7900 FORMAT (/,'----> distance of charge densities for spin ',i2,' it=',i5,':',f13.6,' me/bohr**3')
8000 FORMAT (/,'----> distance of charge densities for it=',i5,':', f13.6,' me/bohr**3')
8010 FORMAT (/,'----> distance of spin densities for it=',i5,':', f13.6,' me/bohr**3')
8020 FORMAT (4d25.14)