Commit dfec26e9 authored by Daniel Wortmann's avatar Daniel Wortmann

slowly reaching compilable version...

parent fd2b2b92
......@@ -15,33 +15,24 @@ MODULE m_mix
contains
SUBROUTINE mix( field, xcpot, DIMENSION, obsolete, sliceplot, mpi, &
SUBROUTINE mix( field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, hybrid, archiveType, inDen, outDen, results )
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_kerker
use m_types_mixvector
#ifdef CPP_MPI
use m_mpi_bc_potden
#endif
use m_distance
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
......@@ -50,10 +41,7 @@ contains
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
......@@ -63,8 +51,10 @@ contains
real :: fix
type(t_potden) :: resDen, vYukawa
TYPE(t_mixvector) :: sm, fsm, fmMet, smMet,fsm_mag
TYPE(t_mixvector),allocatable :: sm(:), fsm(:)
TYPE(t_mixvector) :: fmMet, smMet,fsm_mag
LOGICAL :: l_densitymatrix
integer :: it
MPI0_a: IF( mpi%irank == 0 ) THEN
......@@ -104,51 +94,24 @@ contains
ENDIF
ENDIF
CALL mixvector_init(mpi%mpi_comm,l_densitymatrix,oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms)
maxiter=merge(1,input%maxiter,input%imix==0)
CALL mixing_history(maxiter,inden,outden,sm,fsm,it)
CALL sm%alloc()
CALL fsm%alloc()
CALL fmMet%alloc()
CALL fsm_mag%alloc()
!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 sm%from_density(inDen)
CALL fsm%from_density(outDen)
!store the difference fsm - sm in fsm
fsm = fsm - sm
! calculate Magnetisation-difference
CALL fsm_mag%from_density(outden,swapspin=.true.)
fsm_mag=fsm_mag-sm
! Apply metric w to fsm and store in fmMet: w |fsm>
fmMet=fsm%apply_metric()
!CALL distance(fsm,fmMet)
CALL distance(mpi%irank,cell%vol,input%jspins,fsm,sm,inDen%iter,outDen,results)
! KERKER PRECONDITIONER
IF( input%preconditioning_param /= 0 ) THEN
call kerker(field, DIMENSION, mpi, &
IF( input%preconditioning_param /= 0 ) call kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm )
ENDIF
oneD, inDen, outDen, fsm(it) )
!mixing of the densities
SELECT CASE(input%imix)
CASE(0)
CALL stmix(atoms,input,noco,fsm,fsm_mag,sm)
CASE(9)
CALL pulay()
CASE default
PRINT *,"New Broyden"
!CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
! hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
END SELECT
if(input%imix==0.or.it==1) CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it))
!if(it>1.and.input%imix==9) CALL pulay(input%alpha,fsm,sm)
if(it>1.and.(input%imax==3.or.input%imix==5.or.input%imix==7)) Call broyden(input%alpha,fsm,sm)
!initiatlize mixed density and extract it
CALL sm%to_density(inDen)
CALL sm(it)%to_density(inDen)
!fix charge of the new density
CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.FALSE., fix)
......@@ -190,14 +153,7 @@ contains
inDen%iter = inDen%iter + 1
7900 FORMAT (/,'----> distance of charge densities for spin ',i2,' it=',i5,':',f13.6,' me/bohr**3')
7901 FORMAT (/,'----> HF 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')
8001 FORMAT (/,'----> HF 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')
8011 FORMAT (/,'----> HF distance of spin densities for it=',i5,':', f13.6,' me/bohr**3')
8020 FORMAT (4d25.14)
8030 FORMAT (10i10)
END SUBROUTINE mix
......
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 :: dfivi,fmvm,vmnorm
REAL,ALLOCATABLE :: am(:)
TYPE(t_mixvector) :: fm1,sm1,ui,um,vi,vm
Type(t_mixvector) :: u_store(:),v_store(:)
hlen=size(fm)
ALLOCATE(u_store(hlen-1),v_store(h_len-1))
do it=1,hlen-1
call u_store(it)%alloc()
call v_store(it)%alloc()
enddo
dfivi = 0.0
CALL fm1%alloc()
CALL sm1%alloc()
CALL ui%alloc()
CALL um%alloc()
CALL vi%alloc()
CALL vm%alloc()
ALLOCATE (am(hlen-1))
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(n)
vi=v_store(n)
am(it) = vi.dot.fm
! 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(n)
! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
vm=vm-am(it)*dfivi*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 = vmnorm
if (n<hlen) u_store(n)=um
if (n<hlen) v_store(n)=vm
enddo
! update rho(m+1)
! calculate <fm|w|vm>
fmvm = vm.dot.fm
! calculate sm(:) = (1.0-fmvm)*um(:) + sm
sm=sm+(1.0-fmvm)*um
END SUBROUTINE broyden
END MODULE m_broyden
......@@ -5,68 +5,82 @@
!--------------------------------------------------------------------------------
module m_distance
contains
SUBROUTINE distance(irank,fsm)
SUBROUTINE distance(irank,vol,jspins,fsm,sm,iter,outden,results)
use m_types
use m_types_mixvector
use m_xmlOutput
implicit none
integer,intent(in) :: irank,jspins,iter
real,intent(in) :: vol
type(t_mixvector),INTENT(IN) :: fsm,sm
TYPE(t_potden),INTENT(IN) :: outden
TYPE(t_results),INTENT(INOUT) :: results
REAL :: dist(6) !1:up,2:down,3:total,4:
integer ::js
REAL :: dist(6) !1:up,2:down,3:spinoff,4:total,5:magnet,6:noco
type(t_mixvector)::fmMet,fsm_mag
character(len=100)::attributes(2)
CALL fmMet%alloc()
if (input%jspins==2) CALL fsm_mag%alloc()
! calculate Magnetisation-difference
CALL fsm_mag%from_density(outden,swapspin=.true.)
fsm_mag=fsm_mag-sm(it)
CALL fsm_up%alloc()
IF (input%jspin
fsm_up=
! Apply metric w to fsm and store in fmMet: w |fsm>
fmMet=fsm%apply_metric()
dist(:) = 0.0
! Apply metric w to fsm and store in fmMet: w |fsm>
DO js = 1,input%jspins
dist(js) = fsm%multiply_dot_mask(fmMet,(/.true.,.true.,.true.,.false./),js)
END DO
dist(6) = fsm%multiply_dot_mask(fmMet,(/.true.,.true.,.true.,.false./),3)
IF (input%jspins.EQ.2) then
dist(3) = fsm_mag%multiply_dot_mask(fmMet,(/.true.,.true.,.true.,.false./),1)+&
fsm_mag%multiply_dot_mask(fmMet,(/.true.,.true.,.true.,.false./),2)
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>1) return
!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) = dot_product(fsm%vec(nmaph*(js-1)+1:nmaph*(js-1)+nmaph),fmMet%vec(nmaph*(js-1)+1:nmaph*(js-1)+nmaph))
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/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
WRITE ( 6,FMT=7900) js,iter,1000*SQRT(ABS(dist(js)/vol))
END DO
IF (noco%l_noco) dist(6) = dot_product(fsm%vec(nmaph*2+1:nmaph*2+nmap-2*nmaph),fmMet%vec(nmaph*2+1:nmaph*2+nmap-2*nmaph))
IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
IF (abs(dist(6))>1E-15) WRITE (6,FMT=7900) 3,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 (input%jspins.EQ.2) THEN
dist(3) = fsm.dot.fmMet
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/)))
(/1000*SQRT(ABS(dist(4)/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
(/1000*SQRT(ABS(dist(5)/vol))/),reshape((/19,20/),(/1,2/)))
WRITE ( 6,FMT=8000) iter,1000*SQRT(ABS(dist(4)/vol))
WRITE ( 6,FMT=8010) 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
results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
DEALLOCATE (smMet,fmMet)
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)
8030 FORMAT (10i10)
end SUBROUTINE distance
end module m_distance
!--------------------------------------------------------------------------------
! 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_mixing_history
use m_types_mixvector
implicit none
integer:: iter_stored=0
type(t_mixvector),allocatable::sm_store(:),fsm_store(:)
contains
subroutine mixing_history(maxiter,inden,outden,sm,fsm,it)
use m_types
implicit none
integer,intent(in)::maxiter
type(t_potden),intent(in)::inden,outden
type(t_mixvector),ALLOCATABLE::sm(:),fsm(:)
integer,intent(out)::it
if (.not.allocated(sm_store)) THEN
allocate(sm_store(maxiter),fsm_store(maxiter))
endif
it=min(iter_stored+1,maxiter)
allocate(sm(it),fsm(it))
CALL sm(it)%alloc()
CALL fsm(it)%alloc()
CALL sm(it)%from_density(inDen)
CALL fsm(it)%from_density(outDen)
!store the difference fsm - sm in fsm
fsm(it) = fsm(it) - sm(it)
do n=it-1,1,-1 !Copy from storage
sm(n)=sm_store(n)
fm(n)=fsm_store(n)
ENDDO
if(iter_stored<maxiter) THEN
iter_stored=iter_stored+1
sm_store(:iter_stored)=sm(:iter_stored)
fsm_store(:iter_stored)=fsm(:iter_stored)
else
sm_store(:maxiter-1)=sm(2:maxiter)
fsm_store(:maxiter-1)=sm(2:maxiter)
endif
end subroutine mixing_history
end MODULE m_mixing_history
......@@ -520,17 +520,27 @@ CONTAINS
#endif
END FUNCTION multiply_dot
FUNCTION multiply_dot_mask(vec1,vec2,mask)RESULT(dprod)
FUNCTION multiply_dot_mask(vec1,vec2,mask,spin)RESULT(dprod)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
LOGICAL,INTENT(IN) ::mask(4)
INTEGER,INTENT(IN) ::spin
REAL ::dprod
dprod=0.0
IF (mask(1)) dprod=dot_PRODUCT(vec1%vec_pw,vec2%vec_pw)
IF (mask(2)) dprod=dprod+dot_PRODUCT(vec1%vec_mt,vec2%vec_mt)
IF (mask(3)) dprod=dprod+dot_PRODUCT(vec1%vec_vac,vec2%vec_vac)
IF (mask(4)) dprod=dprod+dot_PRODUCT(vec1%vec_misc,vec2%vec_misc)
DO js=1,3
IF (mask(1).and.(spin==js.or.spin==0.and.start_pw(js)>0)) &
dprod=dprod+dot_PRODUCT(vec1%vec_pw(start_pw(js):stop_pw(js)),&
vec2%vec_pw(start_pw(js):stop_pw(js)))
IF (mask(2).and.(spin==js.or.spin==0.and.start_mt(js)>0)) &
dprod=dprod+dot_PRODUCT(vec1%vec_mt(start_mt(js):stop_mt(js)),&
vec2%vec_mt(start_mt(js):stop_mt(js)))
IF (mask(3).and.(spin==js.or.spin==0.and.start_vac(js)>0)) &
dprod=dprod+dot_PRODUCT(vec1%vec_vac(start_vac(js):stop_vac(js)),&
vec2%vec_vac(start_vac(js):stop_vac(js)))
IF (mask(4).and.(spin==js.or.spin==0.and.start_misc(js)>0)) &
dprod=dprod+dot_PRODUCT(vec1%vec_misc(start_misc(js):stop_misc(js)),&
vec2%vec_misc(start_misc(js):stop_misc(js)))
enddo
#ifdef CPP_MPI
CALL MPI_REDUCE_ALL()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment