Commit c2546518 authored by Yi Yao's avatar Yi Yao
Browse files

Merge branch 'occ_number_non_aufbau' into 'master'

Occ number non aufbau

See merge request elsi_project/elsi_interface!270
parents fc389db0 bd7bfc89
......@@ -116,8 +116,9 @@ module ELSI_DATATYPE
logical :: eval_ready = .false.
logical :: evec_ready = .false.
logical :: occ_ready = .false.
logical :: occ_non_aufbau = .false. ! occ for lowest excited state ( ...222000... => ...221100... )
! Chemical potential
! Chemical potential
real(kind=r8) :: mu ! Fermi level
real(kind=r8) :: ts ! Entropy
integer(kind=i4) :: mu_scheme
......
......@@ -58,6 +58,7 @@ module ELSI
public :: elsi_set_elpa_n_single
public :: elsi_set_elpa_gpu
public :: elsi_set_elpa_autotune
public :: elsi_set_occ_non_aufbau
public :: elsi_set_omm_flavor
public :: elsi_set_omm_n_elpa
public :: elsi_set_omm_tol
......
......@@ -48,6 +48,48 @@ subroutine elsi_mu_and_occ(ph,bh,n_electron,n_state,n_spin,n_kpt,k_wt,eval,occ,&
real(kind=r8), intent(out) :: occ(n_state,n_spin,n_kpt)
real(kind=r8), intent(out) :: mu
real(kind=r8) :: occ1(n_state,n_spin,n_kpt)
real(kind=r8) :: occ2(n_state,n_spin,n_kpt)
real(kind=r8) :: mu1
real(kind=r8) :: mu2
if (ph%occ_non_aufbau) then
! YY: This is for non Aufbau distribution occ number.
! YY: Used for (lowest excited state) delta SCF
! YY: normal occ_number : 2 2 2 2 0 0 0 0
! YY: non-Aufbau occ_number: 2 2 2 1 1 0 0 0
call elsi_mu_and_occ_normal(ph,bh,n_electron-2,n_state,n_spin,n_kpt,k_wt,&
eval,occ1,mu1)
call elsi_mu_and_occ_normal(ph,bh,n_electron+2,n_state,n_spin,n_kpt,k_wt,&
eval,occ2,mu2)
occ = (occ1 + occ2) / 2
mu = (mu1 + mu2) / 2
else
call elsi_mu_and_occ_normal(ph,bh,n_electron,n_state,n_spin,n_kpt,k_wt,&
eval,occ,mu)
end if
end subroutine
!>
!! Compute the chemical potential and occupation numbers normal distribution.
!!
subroutine elsi_mu_and_occ_normal(ph,bh,n_electron,n_state,n_spin,n_kpt,k_wt,eval,occ,&
mu)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_basic_t), intent(in) :: bh
real(kind=r8), intent(in) :: n_electron
integer(kind=i4), intent(in) :: n_state
integer(kind=i4), intent(in) :: n_spin
integer(kind=i4), intent(in) :: n_kpt
real(kind=r8), intent(in) :: k_wt(n_kpt)
real(kind=r8), intent(in) :: eval(n_state,n_spin,n_kpt)
real(kind=r8), intent(out) :: occ(n_state,n_spin,n_kpt)
real(kind=r8), intent(out) :: mu
real(kind=r8) :: mu_min
real(kind=r8) :: mu_max
real(kind=r8) :: buf
......@@ -58,7 +100,7 @@ subroutine elsi_mu_and_occ(ph,bh,n_electron,n_state,n_spin,n_kpt,k_wt,eval,occ,&
logical :: found_interval
character(len=200) :: msg
character(len=*), parameter :: caller = "elsi_mu_and_occ"
character(len=*), parameter :: caller = "elsi_mu_and_occ_normal"
! Determine upper and lower bounds of mu
mu_min = minval(eval)
......
......@@ -38,6 +38,7 @@ module ELSI_SET
public :: elsi_set_elpa_n_single
public :: elsi_set_elpa_gpu
public :: elsi_set_elpa_autotune
public :: elsi_set_occ_non_aufbau
public :: elsi_set_omm_flavor
public :: elsi_set_omm_n_elpa
public :: elsi_set_omm_tol
......@@ -1174,6 +1175,25 @@ subroutine elsi_set_mu_tol(eh,tol)
eh%ph%mu_tol = tol
end subroutine
!>
!! Set the non Aufbau occupation
!!
subroutine elsi_set_occ_non_aufbau(eh,occ_non_aufbau)
implicit none
type(elsi_handle), intent(inout) :: eh !< Handle
logical, intent(in) :: occ_non_aufbau !< whether to use non Aufbau occupation
character(len=200) :: msg
character(len=*), parameter :: caller = "elsi_set_occ_non_aufbau"
call elsi_check_init(eh%bh,eh%handle_init,caller)
eh%ph%occ_non_aufbau = occ_non_aufbau
end subroutine
!>
......
......@@ -18,7 +18,9 @@ LIST(APPEND ftest_src
Fortran/test_ev_real_den.f90
Fortran/test_ev_real_coo.f90
Fortran/test_ev_real_csc1.f90
Fortran/test_ev_real_csc2.f90)
Fortran/test_ev_real_csc2.f90
Fortran/test_occ_normal.f90
Fortran/test_occ_non_aufbau.f90)
MACRO(add_test_target arg1)
ADD_DEPENDENCIES(${arg1} elsi)
......
......@@ -10,7 +10,7 @@
program elsi_test
use ELSI_MPI
use ELSI_PRECISION, only: i4
use ELSI_PRECISION, only: i4, r8
implicit none
......@@ -25,6 +25,8 @@ program elsi_test
integer(kind=i4) :: myid
integer(kind=i4) :: ierr
real(kind=r8) :: mu_width
call MPI_Init(ierr)
call MPI_Comm_rank(MPI_COMM_WORLD,myid,ierr)
......@@ -133,6 +135,16 @@ program elsi_test
case default
call test_die()
end select
case("o") ! occupation number
read(arg3,*) mu_width
select case(arg2(1:1))
case("1") ! normal occupation distribution
call test_occ_normal(MPI_COMM_WORLD,mu_width)
case("2") ! Non-Aufbau occupation distribution
call test_occ_non_aufbau(MPI_COMM_WORLD,mu_width)
case default
call test_die()
end select
case default
call test_die()
end select
......
! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.
!>
!! This subroutine tests occ number normal occupation.
!!
subroutine test_occ_non_aufbau(comm,mu_width)
use ELSI
use ELSI_MPI
use ELSI_PRECISION, only: r8,i4
implicit none
integer(kind=i4), intent(in) :: comm
real(kind=r8), intent(in) :: mu_width
integer(kind=i4) :: n_proc
integer(kind=i4) :: nprow
integer(kind=i4) :: npcol
integer(kind=i4) :: myid
integer(kind=i4) :: myprow
integer(kind=i4) :: mypcol
integer(kind=i4) :: ierr
integer(kind=i4) :: blk
integer(kind=i4) :: blacs_ctxt
integer(kind=i4) :: n_basis
integer(kind=i4) :: l_rows
integer(kind=i4) :: l_cols
integer(kind=i4) :: l_rows2
integer(kind=i4) :: l_cols2
integer(kind=i4) :: header(8)
integer(kind=i4) :: i_state
integer(kind=i4) :: n_kpt
real(kind=r8) :: k_wt(1)
real(kind=r8) :: occ(100,1,1) !< Occupation members
real(kind=r8) :: mu !< Chemical potential
real(kind=r8) :: n_electrons
real(kind=r8) :: tol
real(kind=r8) :: t1
real(kind=r8) :: t2
logical :: file_exist
real(kind=r8), allocatable :: eval(:,:,:)
type(elsi_handle) :: eh
integer(kind=i4), external :: numroc
! Reference values
real(kind=r8), parameter :: e_ref = 0.31442451613_r8
character(len=*), parameter :: file_name = "elsi.in"
tol = 1.0e-8_r8
header(:) = 0
if(myid == 0) then
write(*,"(2X,A)") "################################"
write(*,"(2X,A)") "## ELSI TEST PROGRAMS ##"
write(*,"(2X,A)") "################################"
write(*,*)
write(*,*)
end if
n_kpt = 1
n_basis = 100
n_electrons = 100
allocate(eval(n_basis,1,1))
do i_state = 1, 100
eval(i_state,1,1) = dble(i_state)
end do
k_wt = 1.0
! Initialize ELSI
call elsi_init(eh,1,1,0,n_basis,n_electrons,n_basis)
call elsi_set_mpi(eh,comm)
eh%ph%occ_non_aufbau = .true.
eh%ph%mu_width = mu_width
! Customize ELSI
call elsi_set_output(eh,2)
call elsi_set_output_log(eh,1)
call elsi_compute_mu_and_occ(eh,n_electrons,n_basis,1,n_kpt,k_wt,eval,occ,mu)
if(myid == 0) then
write(*,"(2X,A)") "Finished occ normal"
write(*,*) "#chemical potential mu", mu
write(*,*) "# state eigenvalue occ"
do i_state = 1, 100
write(*,*) i_state, eval(i_state,1,1), occ(i_state,1,1)
end do
write(*,*)
end if
! Finalize ELSI
call elsi_finalize(eh)
deallocate(eval)
end subroutine
! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.
!>
!! This subroutine tests occ number normal occupation.
!!
subroutine test_occ_normal(comm,mu_width)
use ELSI
use ELSI_MPI
use ELSI_PRECISION, only: r8,i4
implicit none
integer(kind=i4), intent(in) :: comm
real(kind=r8), intent(in) :: mu_width
integer(kind=i4) :: n_proc
integer(kind=i4) :: nprow
integer(kind=i4) :: npcol
integer(kind=i4) :: myid
integer(kind=i4) :: myprow
integer(kind=i4) :: mypcol
integer(kind=i4) :: ierr
integer(kind=i4) :: blk
integer(kind=i4) :: blacs_ctxt
integer(kind=i4) :: n_basis
integer(kind=i4) :: l_rows
integer(kind=i4) :: l_cols
integer(kind=i4) :: l_rows2
integer(kind=i4) :: l_cols2
integer(kind=i4) :: header(8)
integer(kind=i4) :: i_state
integer(kind=i4) :: n_kpt
real(kind=r8) :: k_wt(1)
real(kind=r8) :: occ(100,1,1) !< Occupation members
real(kind=r8) :: mu !< Chemical potential
real(kind=r8) :: n_electrons
real(kind=r8) :: tol
real(kind=r8) :: t1
real(kind=r8) :: t2
logical :: file_exist
real(kind=r8), allocatable :: eval(:,:,:)
type(elsi_handle) :: eh
integer(kind=i4), external :: numroc
! Reference values
real(kind=r8), parameter :: e_ref = 0.31442451613_r8
character(len=*), parameter :: file_name = "elsi.in"
tol = 1.0e-8_r8
header(:) = 0
if(myid == 0) then
write(*,"(2X,A)") "################################"
write(*,"(2X,A)") "## ELSI TEST PROGRAMS ##"
write(*,"(2X,A)") "################################"
write(*,*)
write(*,*)
end if
n_kpt = 1
n_basis = 100
n_electrons = 100
allocate(eval(n_basis,1,1))
do i_state = 1, 100
eval(i_state,1,1) = dble(i_state)
end do
k_wt = 1.0
eh%ph%mu_width = mu_width
! Initialize ELSI
call elsi_init(eh,1,1,0,n_basis,n_electrons,n_basis)
call elsi_set_mpi(eh,comm)
! Customize ELSI
call elsi_set_output(eh,2)
call elsi_set_output_log(eh,1)
call elsi_compute_mu_and_occ(eh,n_electrons,n_basis,1,n_kpt,k_wt,eval,occ,mu)
if(myid == 0) then
write(*,"(2X,A)") "Finished occ normal"
write(*,*) "#chemical potential mu", mu
write(*,*) "# state eigenvalue occ"
do i_state = 1, 100
write(*,*) i_state, eval(i_state,1,1), occ(i_state,1,1)
end do
write(*,*)
end if
! Finalize ELSI
call elsi_finalize(eh)
deallocate(eval)
end subroutine
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