Commit 8092b66f authored by Nick R. Papior's avatar Nick R. Papior
Browse files

Enabled EDM convergence criteria, controls the forces

- The convergence criteria now enables converging the
  energy density matrix.

  This is a simple addition which does not introduce
  any complexity to the code.
  I expect it may be used to control the fineness of
  the force calculations.
  It does, however, not seem too much different than
  the Hamiltonian convergence.
parent 0e98f1ad
......@@ -5032,6 +5032,25 @@ simultaneously for complete control.
\end{fdfentry}
\begin{fdflogicalT}{SCF.Converge!EDM}
\index{SCF!mixing!energy density matrix convergence}
Logical variable to use the energy density matrix elements as monitor
of self-consistency: this is considered achieved when the maximum
absolute change (dEmax) in the energy density matrix elements is below
\fdf{SCF.Tolerance!EDM}. The meaning of dEmax is equivalent to that
of \fdf{SCF.Tolerance!DM}.
\end{fdflogicalT}
\begin{fdfentry}{SCF.Tolerance!EDM}[energy]<$10^{-3}\,\mathrm{eV}$>
If \fdf{SCF.Converge!EDM} is \fdftrue, then self-consistency is
achieved when the maximum absolute change in the energy density
matrix elements is below this value.
\end{fdfentry}
\begin{fdflogicalF}{SCF.Converge!FreeE}
\index{SCF!mixing!energy convergence}
......
......@@ -81,12 +81,20 @@
#endif
if (SIESTA_worker) then
! Save present density matrix
!$OMP parallel workshare default(shared)
!! Eold(:,:) = Escf(:,:)
Dold(:,:) = Dscf(:,:)
!$OMP end parallel workshare
endif
! Save present density matrix
!$OMP parallel default(shared)
if ( converge_EDM ) then
!$OMP workshare
Eold(:,:) = Escf(:,:)
Dold(:,:) = Dscf(:,:)
!$OMP end workshare
else
!$OMP workshare
Dold(:,:) = Dscf(:,:)
!$OMP end workshare
end if
!$OMP end parallel
end if
! Compute shift in Tr(H*DM) for fermi-level bracketting
! Use the current H, the previous iteration H, and the
......
......@@ -549,18 +549,35 @@ subroutine read_options( na, ns, nspin )
dDtol = fdf_get('SCF.Tolerance.DM',dDtol)
if ( IONode ) then
write(6,1) 'redata: Require DM convergence for SCF', converge_DM
write(6,7) 'redata: DM tolerance for SCF',dDtol
write(6,9) 'redata: DM tolerance for SCF',dDtol
end if
if (cml_p) then
call cmlAddParameter( xf=mainXML, name='SCF.Converge.DM', &
value=converge_DM, &
dictRef='siesta:ReqDMConv' )
call cmlAddParameter( xf=mainXML, name='SCF.Tolerance.DM', &
call cmlAddParameter( xf=mainXML, name='SCF.Tolerance.DM', &
value=dDtol, dictRef='siesta:dDtol', &
units='siestaUnits:eAng_3' )
end if
! Energy-density matrix convergence
converge_EDM = fdf_get('SCF.Converge.EDM', .false.)
tolerance_EDM = fdf_get('SCF.Tolerance.EDM',1.e-3_dp*eV, 'Ry')
if ( IONode ) then
write(6,1) 'redata: Require EDM convergence for SCF', converge_EDM
write(6,7) 'redata: EDM tolerance for SCF',tolerance_EDM/eV, ' eV'
end if
if (cml_p) then
call cmlAddParameter( xf=mainXML, name='SCF.Converge.EDM', &
value=converge_EDM, &
dictRef='siesta:ReqEDMConv' )
call cmlAddParameter( xf=mainXML, name='SCF.Tolerance.EDM', &
value=tolerance_EDM/eV, dictRef='siesta:EDM_tolerance', &
units='siestaUnits:eV' )
end if
! Hamiltonian convergence
if ( compat_pre_v4_DM_H ) then
tBool = .false.
......@@ -607,6 +624,7 @@ subroutine read_options( na, ns, nspin )
tBool = .false.
tBool = tBool .or. converge_Eharr
tBool = tBool .or. converge_FreeE
tBool = tBool .or. converge_EDM
tBool = tBool .or. converge_DM
tBool = tBool .or. converge_H
if ( .not. tBool ) then
......
......@@ -10,7 +10,8 @@
public :: scfconvergence_test
CONTAINS
subroutine scfconvergence_test( first, iscf, dDmax, dHmax,
subroutine scfconvergence_test( first, iscf,
& dDmax, dHmax, dEmax,
$ conv_harris, conv_freeE,
$ converged)
use precision, only: dp
......@@ -31,8 +32,10 @@
integer :: iscf
logical :: first
real(dp), intent(in) :: dDmax ! Max. change in dens. matrix
real(dp), intent(in) :: dHmax ! Max. change in H
real(dp), intent(in) :: dDmax ! Max. change in dens. matrix
real(dp), intent(in) :: dHmax ! Max. change in H
real(dp), intent(in) :: dEmax ! Max. change in EDM
type(converger_t), intent(inout) :: conv_harris, conv_freeE
logical, intent(out) :: converged
......@@ -85,6 +88,10 @@
conv_text = trim(conv_text)//'+FreeE'
converged = converged .and. is_converged(conv_FreeE)
end if
if ( converge_EDM ) then
conv_text = trim(conv_text)//'+EDM'
converged = converged .and. dEmax < tolerance_EDM
end if
if ( converge_DM ) then
conv_text = trim(conv_text)//'+DM'
converged = converged .and. dDmax < dDtol
......@@ -114,12 +121,21 @@
! No matter what we print out the differences of DM and H
if ( mix_charge ) then
write(6,"(a,f16.10)") "max |DM_i - DM_(i-1)| : ", dDmax
write(6,"(a,f16.10)") "max |DM_i - DM_(i-1)| : ",
& dDmax
if ( dEmax >= 0._dp )
& write(6,"(a,f16.10)") "max |EDM_i - EDM_(i-1)| (eV) : ",
& dEmax/eV
else
write(6,"(a,f16.10)") "max |DM_out - DM_in| : ", dDmax
write(6,"(a,f16.10)") "max |DM_out - DM_in| : ",
& dDmax
if ( dEmax >= 0._dp )
& write(6,"(a,f16.10)") "max |EDM_out - EDM_in| (eV) : ",
& dEmax/eV
endif
if ( dHmax >= 0._dp ) then
write(6,"(a,f16.10)") "max |H_out - H_in| (eV): ", dHmax/eV
write(6,"(a,f16.10)") "max |H_out - H_in| (eV) : ",
& dHmax/eV
end if
if ( is_LDAU_conv ) then
......
......@@ -40,7 +40,7 @@
use m_state_analysis
use m_steps
use sys, only : die, bye
use sparse_matrices, only: H, Hold
use sparse_matrices, only: H, Hold, Dold, Dscf, Eold, Escf
use m_convergence, only: converger_t
use m_convergence, only: reset, set_tolerance
use siesta_geom, only: na_u ! Number of atoms in unit cell
......@@ -53,7 +53,6 @@
use m_iodm_old, only: write_spmatrix
!
use units, only: eV, Ang
use sparse_matrices, only: H, Hold, Dold, Dscf
use m_pexsi_solver, only: prevDmax
use m_forces, only: fa
use units, only: eV, Ang
......@@ -97,6 +96,7 @@
logical :: SCFconverged
real(dp) :: dDmax ! Max. change in DM elements
real(dp) :: dHmax ! Max. change in H elements
real(dp) :: dEmax ! Max. change in EDM elements
real(dp) :: drhog ! Max. change in rho(G) (experimental)
type(converger_t) :: conv_harris, conv_freeE
!---------------------------------------------------------------------- BEGIN
......@@ -155,6 +155,7 @@
call dict_variable_add('SCF.iteration',iscf)
call dict_variable_add('SCF.dD',dDmax)
call dict_variable_add('SCF.dH',dHmax)
call dict_variable_add('SCF.dE',dEmax)
call dict_variable_add('SCF.drhoG',drhog)
#endif
......@@ -204,6 +205,7 @@
! The dHmax variable only has meaning for Hamiltonian
! mixing, or when requiring the Hamiltonian to be converged.
dHmax = -1._dp
dEmax = -1._dp
DO
! Conditions of exit:
......@@ -256,6 +258,8 @@
if (SIESTA_worker) then
! Maybe set Dold to zero if reading charge or H...
call compute_max_diff(Dold,Dscf,dDmax)
if ( converge_EDM )
& call compute_max_diff(Eold,Escf,dEmax)
call setup_hamiltonian( iscf )
call compute_max_diff(Hold,H,dHmax)
end if
......@@ -265,7 +269,11 @@
call compute_max_diff(Hold,H,dHmax)
end if
call compute_dm( iscf )
if (SIESTA_worker) call compute_max_diff(Dold,Dscf,dDmax)
if (SIESTA_worker) then
call compute_max_diff(Dold,Dscf,dDmax)
if ( converge_EDM )
& call compute_max_diff(Eold,Escf,dEmax)
end if
endif
......@@ -284,7 +292,8 @@
! Dscf=DM_out, Dold=DM_in, H=H_(DM_out), Hold=H_in(mixed)
! dDmax=maxdiff(DM_out,DM_in)
! dHmax=maxdiff(H(DM_out),H_in)
call scfconvergence_test( first, iscf, dDmax, dHmax,
call scfconvergence_test( first, iscf,
& dDmax, dHmax, dEmax,
& conv_harris, conv_freeE,
& SCFconverged )
!
......
......@@ -128,14 +128,16 @@ MODULE siesta_options
logical :: atmonly ! Set up pseudoatom information only?
logical :: harrisfun ! Use Harris functional?
logical :: muldeb ! Write Mulliken polpulations at every SCF step?
logical :: converge_FreeE ! free Energy conv. to finish SCF iteration?
real(dp) :: tolerance_FreeE ! Free-energy tolerance
logical :: converge_Eharr ! to finish SCF iteration?
real(dp) :: tolerance_Eharr ! Harris tolerance
logical :: converge_DM ! to finish SCF iteration?
real(dp) :: dDtol ! Tolerance in change of DM elements to finish SCF iteration
logical :: converge_H ! to finish SCF iteration?
real(dp) :: dHtol ! Tolerance in change of H elements to finish SCF iteration
logical :: converge_FreeE ! free Energy conv. to finish SCF iteration?
real(dp):: tolerance_FreeE ! Free-energy tolerance
logical :: converge_Eharr ! to finish SCF iteration?
real(dp):: tolerance_Eharr ! Harris tolerance
logical :: converge_EDM ! to finish SCF iteration?
real(dp):: tolerance_EDM ! Tolerance in change of EDM elements to finish SCF iteration
logical :: converge_DM ! to finish SCF iteration?
real(dp):: dDtol ! Tolerance in change of DM elements to finish SCF iteration
logical :: converge_H ! to finish SCF iteration?
real(dp):: dHtol ! Tolerance in change of H elements to finish SCF iteration
logical :: broyden_optim ! Use Broyden method to optimize geometry?
logical :: fire_optim ! Use FIRE method to optimize geometry?
logical :: struct_only ! Output initial structure only?
......
......@@ -23,12 +23,14 @@
integer, public, pointer :: listh(:)=>null(), listhptr(:)=>null(),
$ numh(:)=>null()
real(dp), public, pointer :: Dold(:,:)=>null(), Dscf(:,:)=>null(),
& Escf(:,:)=>null(), H(:,:)=>null(),
& xijo(:,:)=>null(), S(:)=>null()
real(dp), public, pointer :: Dold(:,:)=>null(), Dscf(:,:)=>null()
real(dp), public, pointer :: Eold(:,:)=>null(), Escf(:,:)=>null()
real(dp), public, pointer :: Hold(:,:)=>null(), H(:,:)=>null()
real(dp), public, pointer :: S(:)=>null()
real(dp), public, pointer :: xijo(:,:)=>null()
! Used to hold the "input" H when mixing the Hamiltonian
real(dp), dimension(:,:), public, pointer :: Hold => null()
! Pieces of H that do not depend on the SCF density matrix
! Formerly there was a single array H0 for this
......@@ -80,6 +82,7 @@
call de_alloc( Hold, 'Hold', 'sparseMat' )
call de_alloc( Dold, 'Dold', 'sparseMat' )
call de_alloc( Eold, 'Eold', 'sparseMat' )
end subroutine resetSparseMatrices
......
......@@ -20,15 +20,18 @@
use siesta_options
use units, only: Ang
use sparse_matrices, only: maxnh, numh, listh, listhptr
use sparse_matrices, only: Dscf, Dold, Escf, Hold
use sparse_matrices, only: xijo, H, S_1D, S
use sparse_matrices, only: Dold, Dscf, DM_2D
use sparse_matrices, only: Eold, Escf, EDM_2D
use sparse_matrices, only: Hold, H , H_2D
use sparse_matrices, only: xijo, xij_2D
use sparse_matrices, only: S , S_1D
use sparse_matrices, only: H_kin_1D, H_vkb_1D
use sparse_matrices, only: H_ldau_2D, H_so_2D
use sparse_matrices, only: H_2D, xij_2D
use sparse_matrices, only: sparse_pattern
use sparse_matrices, only: block_dist, single_dist
use sparse_matrices, only: DM_2D, EDM_2D, DM_history
use sparse_matrices, only: DM_history
use m_sparsity_handling, only: SpOrb_to_SpAtom
......@@ -438,10 +441,14 @@
!if (ionode) call print_type(EDM_2D)
Escf => val(EDM_2D)
call re_alloc(Dold,1,maxnh,1,spin%H,name='Dold',
call re_alloc(Dold,1,maxnh,1,spin%DM,name='Dold',
. routine='sparseMat',copy=.false.,shrink=.true.)
call re_alloc(Hold,1,maxnh,1,spin%H,name='Hold',
. routine='sparseMat',copy=.false.,shrink=.true.)
if ( converge_EDM ) then
call re_alloc(Eold,1,maxnh,1,spin%EDM,name='Eold',
. routine='sparseMat',copy=.false.,shrink=.true.)
end if
! Allocate/reallocate storage associated with Hamiltonian/Overlap matrix
write(oname,"(a,i0)") "H at geom step ", istep
......
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