Commit 130091cb authored by Victor Yu's avatar Victor Yu Committed by Victor Yu
Browse files

Add serial eigensolver with frozen core

parent 16076a04
......@@ -7,7 +7,7 @@ SET(elsi_URL "http://elsi-interchange.org")
SET(elsi_EMAIL "elsi-team@duke.edu")
SET(elsi_LICENSE "BSD 3")
SET(elsi_DESCRIPTION "Electronic Structure Infrastructure")
SET(elsi_DATESTAMP "20201117")
SET(elsi_DATESTAMP "20210124")
### CMake modules ###
LIST(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
......
......@@ -9,17 +9,21 @@
!!
module ELSI_LAPACK
use ELSI_CONSTANT, only: FC_BASIC,FC_PLUS_V
use ELSI_DATATYPE, only: elsi_param_t,elsi_basic_t
use ELSI_ELPA, only: elsi_elpa_tridiag
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_OUTPUT, only: elsi_say,elsi_get_time
use ELSI_PRECISION, only: r8,i4
use ELSI_SORT, only: elsi_heapsort,elsi_permute
implicit none
private
public :: elsi_solve_lapack
public :: elsi_do_fc_lapack
public :: elsi_undo_fc_lapack
interface elsi_solve_lapack
module procedure elsi_solve_lapack_real
......@@ -46,6 +50,16 @@ module ELSI_LAPACK
module procedure elsi_back_ev_sp_cmplx
end interface
interface elsi_do_fc_lapack
module procedure elsi_do_fc_lapack_real
module procedure elsi_do_fc_lapack_cmplx
end interface
interface elsi_undo_fc_lapack
module procedure elsi_undo_fc_lapack_real
module procedure elsi_undo_fc_lapack_cmplx
end interface
contains
!>
......@@ -283,6 +297,14 @@ subroutine elsi_solve_lapack_real(ph,bh,ham,ovlp,eval,evec)
call elsi_back_ev_sp(ph,bh,ovlp,evec)
end if
! Switch back to full dimension in case of frozen core
if(ph%n_basis_c > 0) then
ph%n_basis = ph%n_basis_v+ph%n_basis_c
ph%n_states = ph%n_states+ph%n_basis_c
ph%n_good = ph%n_good+ph%n_basis_c
ph%n_states_solve = ph%n_states_solve+ph%n_basis_c
end if
end subroutine
!>
......@@ -389,6 +411,172 @@ subroutine elsi_check_ovlp_sp_real(ph,bh,ovlp,eval,evec)
end subroutine
!>
!! Freeze core orbitals by transforming Hamiltonian and overlap.
!!
subroutine elsi_do_fc_lapack_real(ph,ham,ovlp,evec,perm)
implicit none
type(elsi_param_t), intent(inout) :: ph
real(kind=r8), intent(inout) :: ham(ph%n_basis,ph%n_basis)
real(kind=r8), intent(inout) :: ovlp(ph%n_basis,ph%n_basis)
real(kind=r8), intent(out) :: evec(ph%n_basis,ph%n_basis)
integer(kind=i4), intent(in) :: perm(ph%n_basis)
integer(kind=i4) :: i
character(len=*), parameter :: caller = "elsi_do_fc_lapack_real"
if(ph%fc_perm) then
evec(:,:) = ovlp
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(:,i) = ovlp(:,perm(i))
end if
end do
ovlp(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
ovlp(i,:) = evec(perm(i),:)
end if
end do
end if
evec(:,:) = ovlp
! S_vv = S_vv - S_vc * S_cv
call dsyrk("U","N",ph%n_basis_v,ph%n_basis_c,-1.0_r8,evec(ph%n_basis_c+1,1),&
ph%n_basis,1.0_r8,ovlp(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
if(ph%fc_perm) then
evec(:,:) = ham
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(:,i) = ham(:,perm(i))
end if
end do
ham(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
ham(i,:) = evec(perm(i),:)
end if
end do
end if
! Compute H_vv
call dgemm("N","N",ph%n_basis_v,ph%n_basis_c,ph%n_basis_c,1.0_r8,&
ovlp(ph%n_basis_c+1,1),ph%n_basis,ham,ph%n_basis,0.0_r8,evec,ph%n_basis)
if(ph%fc_method == FC_PLUS_V) then
! H_vv = H_vv + S_vc * H_cc * S_cv - H_vc * S_cv - S_vc * H_cv
! More accurate than H_vv = H_vv - S_vc * H_cc * S_cv
call dgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,1.0_r8,evec,&
ph%n_basis,ovlp(1,ph%n_basis_c+1),ph%n_basis,1.0_r8,&
ham(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
call dgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,1.0_r8,&
ham(ph%n_basis_c+1,1),ph%n_basis,ovlp(1,ph%n_basis_c+1),ph%n_basis,&
0.0_r8,evec(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
ham(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis) =&
ham(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis)&
-evec(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis)&
-transpose(evec(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis))
else
! H_vv = H_vv - S_vc * H_cc * S_cv
call dgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,-1.0_r8,evec,&
ph%n_basis,ovlp(1,ph%n_basis_c+1),ph%n_basis,1.0_r8,&
ham(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
end if
! Switch to valence dimension
ph%n_basis = ph%n_basis-ph%n_basis_c
ph%n_states = ph%n_states-ph%n_basis_c
ph%n_good = ph%n_good-ph%n_basis_c
ph%n_states_solve = ph%n_states_solve-ph%n_basis_c
end subroutine
!>
!! Transforming eigenvectors back to unfrozen eigenproblem.
!!
subroutine elsi_undo_fc_lapack_real(ph,bh,ham,ovlp,evec,perm,eval_c)
implicit none
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(inout) :: bh
real(kind=r8), intent(inout) :: ham(ph%n_basis,ph%n_basis)
real(kind=r8), intent(in) :: ovlp(ph%n_basis,ph%n_basis)
real(kind=r8), intent(out) :: evec(ph%n_basis,ph%n_basis)
integer(kind=i4), intent(in) :: perm(ph%n_basis)
real(kind=r8), intent(out) :: eval_c(ph%n_basis_c)
integer(kind=i4) :: i
integer(kind=i4) :: j
integer(kind=i4), allocatable :: idx(:)
integer(kind=i4), allocatable :: tmp(:)
character(len=*), parameter :: caller = "elsi_undo_fc_lapack_real"
call elsi_allocate(bh,idx,ph%n_basis_c,"idx",caller)
call elsi_allocate(bh,tmp,ph%n_basis_c,"tmp",caller)
eval_c(:) = 0.0_r8
do i = 1,ph%n_basis_c
idx(i) = i
if(ph%fc_method > FC_BASIC) then
eval_c(i) = ham(i,i)/ovlp(i,i)
else
eval_c(i) = ham(i,i)
end if
end do
call elsi_heapsort(ph%n_basis_c,eval_c,tmp)
call elsi_permute(ph%n_basis_c,tmp,idx)
evec(:,1:ph%n_basis_c) = 0.0_r8
do j = 1,ph%n_basis_c
i = idx(j)
if(ph%fc_method > FC_BASIC) then
evec(i,j) = 1.0_r8/sqrt(ovlp(i,i))
else
evec(i,j) = 1.0_r8
end if
end do
call elsi_deallocate(bh,idx,"idx")
call elsi_deallocate(bh,tmp,"tmp")
! C_cv = -S_cv * C_vv
call dgemm("N","N",ph%n_basis_c,ph%n_basis_v,ph%n_basis_v,-1.0_r8,&
ovlp(1,ph%n_basis_c+1),ph%n_basis,evec(ph%n_basis_c+1,ph%n_basis_c+1),&
ph%n_basis,0.0_r8,evec(1,ph%n_basis_c+1),ph%n_basis)
if(ph%fc_perm) then
ham(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(perm(i),:) = ham(i,:)
end if
end do
end if
end subroutine
!>
!! Cholesky factorize the overlap matrix in place.
!!
......@@ -632,6 +820,14 @@ subroutine elsi_solve_lapack_cmplx(ph,bh,ham,ovlp,eval,evec)
call elsi_back_ev_sp(ph,bh,ovlp,evec)
end if
! Switch back to full dimension in case of frozen core
if(ph%n_basis_c > 0) then
ph%n_basis = ph%n_basis_v+ph%n_basis_c
ph%n_states = ph%n_states+ph%n_basis_c
ph%n_good = ph%n_good+ph%n_basis_c
ph%n_states_solve = ph%n_states_solve+ph%n_basis_c
end if
end subroutine
!>
......@@ -741,4 +937,174 @@ subroutine elsi_check_ovlp_sp_cmplx(ph,bh,ovlp,eval,evec)
end subroutine
!>
!! Freeze core orbitals by transforming Hamiltonian and overlap.
!!
subroutine elsi_do_fc_lapack_cmplx(ph,ham,ovlp,evec,perm)
implicit none
type(elsi_param_t), intent(inout) :: ph
complex(kind=r8), intent(inout) :: ham(ph%n_basis,ph%n_basis)
complex(kind=r8), intent(inout) :: ovlp(ph%n_basis,ph%n_basis)
complex(kind=r8), intent(out) :: evec(ph%n_basis,ph%n_basis)
integer(kind=i4), intent(in) :: perm(ph%n_basis)
integer(kind=i4) :: i
character(len=*), parameter :: caller = "elsi_do_fc_lapack_cmplx"
if(ph%fc_perm) then
evec(:,:) = ovlp
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(:,i) = ovlp(:,perm(i))
end if
end do
ovlp(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
ovlp(i,:) = evec(perm(i),:)
end if
end do
end if
evec(:,:) = ovlp
! S_vv = S_vv - S_vc * S_cv
call zherk("U","N",ph%n_basis_v,ph%n_basis_c,(-1.0_r8,0.0_r8),&
evec(ph%n_basis_c+1,1),ph%n_basis,(1.0_r8,0.0_r8),&
ovlp(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
if(ph%fc_perm) then
evec(:,:) = ham
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(:,i) = ham(:,perm(i))
end if
end do
ham(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
ham(i,:) = evec(perm(i),:)
end if
end do
end if
! Compute H_vv
call zgemm("N","N",ph%n_basis_v,ph%n_basis_c,ph%n_basis_c,(1.0_r8,0.0_r8),&
ovlp(ph%n_basis_c+1,1),ph%n_basis,ham,ph%n_basis,(0.0_r8,0.0_r8),evec,&
ph%n_basis)
if(ph%fc_method == FC_PLUS_V) then
! H_vv = H_vv + S_vc * H_cc * S_cv - H_vc * S_cv - S_vc * H_cv
! More accurate than H_vv = H_vv - S_vc * H_cc * S_cv
call zgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,&
(1.0_r8,0.0_r8),evec,ph%n_basis,ovlp(1,ph%n_basis_c+1),ph%n_basis,&
(1.0_r8,0.0_r8),ham(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
call zgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,&
(1.0_r8,0.0_r8),ham(ph%n_basis_c+1,1),ph%n_basis,&
ovlp(1,ph%n_basis_c+1),ph%n_basis,(0.0_r8,0.0_r8),&
evec(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
ham(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis) =&
ham(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis)&
-evec(ph%n_basis_c+1:ph%n_basis,ph%n_basis_c+1:ph%n_basis)&
-conjg(transpose(evec(ph%n_basis_c+1:ph%n_basis,&
ph%n_basis_c+1:ph%n_basis)))
else
! H_vv = H_vv - S_vc * H_cc * S_cv
call zgemm("N","N",ph%n_basis_v,ph%n_basis_v,ph%n_basis_c,&
(-1.0_r8,0.0_r8),evec,ph%n_basis,ovlp(1,ph%n_basis_c+1),ph%n_basis,&
(1.0_r8,0.0_r8),ham(ph%n_basis_c+1,ph%n_basis_c+1),ph%n_basis)
end if
! Switch to valence dimension
ph%n_basis = ph%n_basis-ph%n_basis_c
ph%n_states = ph%n_states-ph%n_basis_c
ph%n_good = ph%n_good-ph%n_basis_c
ph%n_states_solve = ph%n_states_solve-ph%n_basis_c
end subroutine
!>
!! Transforming eigenvectors back to unfrozen eigenproblem.
!!
subroutine elsi_undo_fc_lapack_cmplx(ph,bh,ham,ovlp,evec,perm,eval_c)
implicit none
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(inout) :: bh
complex(kind=r8), intent(inout) :: ham(ph%n_basis,ph%n_basis)
complex(kind=r8), intent(in) :: ovlp(ph%n_basis,ph%n_basis)
complex(kind=r8), intent(out) :: evec(ph%n_basis,ph%n_basis)
integer(kind=i4), intent(in) :: perm(ph%n_basis)
real(kind=r8), intent(out) :: eval_c(ph%n_basis_c)
integer(kind=i4) :: i
integer(kind=i4) :: j
integer(kind=i4), allocatable :: idx(:)
integer(kind=i4), allocatable :: tmp(:)
character(len=*), parameter :: caller = "elsi_undo_fc_lapack_cmplx"
call elsi_allocate(bh,idx,ph%n_basis_c,"idx",caller)
call elsi_allocate(bh,tmp,ph%n_basis_c,"tmp",caller)
eval_c(:) = 0.0_r8
do i = 1,ph%n_basis_c
idx(i) = i
if(ph%fc_method > FC_BASIC) then
eval_c(i) = ham(i,i)/ovlp(i,i)
else
eval_c(i) = ham(i,i)
end if
end do
call elsi_heapsort(ph%n_basis_c,eval_c,tmp)
call elsi_permute(ph%n_basis_c,tmp,idx)
evec(:,1:ph%n_basis_c) = (0.0_r8,0.0_r8)
do j = 1,ph%n_basis_c
i = idx(j)
if(ph%fc_method > FC_BASIC) then
evec(i,j) = (1.0_r8,0.0_r8)/sqrt(real(ovlp(i,i),kind=r8))
else
evec(i,j) = (1.0_r8,0.0_r8)
end if
end do
call elsi_deallocate(bh,idx,"idx")
call elsi_deallocate(bh,tmp,"tmp")
! C_cv = -S_cv * C_vv
call zgemm("N","N",ph%n_basis_c,ph%n_basis_v,ph%n_basis_v,(-1.0_r8,0.0_r8),&
ovlp(1,ph%n_basis_c+1),ph%n_basis,evec(ph%n_basis_c+1,ph%n_basis_c+1),&
ph%n_basis,(0.0_r8,0.0_r8),evec(1,ph%n_basis_c+1),ph%n_basis)
if(ph%fc_perm) then
ham(:,:) = evec
do i = 1,ph%n_basis
if(i /= perm(i)) then
evec(perm(i),:) = ham(i,:)
end if
end do
end if
end subroutine
end module ELSI_LAPACK
......@@ -9,7 +9,7 @@
!!
module ELSI_OUTPUT
use ELSI_CONSTANT, only: MULTI_PROC,SINGLE_PROC,BLACS_DENSE,PEXSI_CSC,&
use ELSI_CONSTANT, only: UNSET,MULTI_PROC,SINGLE_PROC,BLACS_DENSE,PEXSI_CSC,&
SIESTA_CSC,GENERIC_COO,ELPA_SOLVER,OMM_SOLVER,PEXSI_SOLVER,&
EIGENEXA_SOLVER,SIPS_SOLVER,NTPOLY_SOLVER,MAGMA_SOLVER,BSEPACK_SOLVER,&
METHFESSEL_PAXTON
......@@ -276,6 +276,14 @@ subroutine elsi_print_handle_summary(ph,bh,jh,caller)
end select
end if
if(caller(6:6) == "e" .and. ph%solver == ELPA_SOLVER&
.and. ph%matrix_format == BLACS_DENSE .and. ph%n_basis_c > 0) then
call fjson_write_name_value(jh,"frozen_core_method",ph%fc_method)
call fjson_write_name_value(jh,"frozen_core_permute",ph%fc_perm)
call fjson_write_name_value(jh,"n_core",ph%n_basis_c)
call fjson_write_name_value(jh,"n_valence",ph%n_basis_v)
end if
end subroutine
!>
......@@ -319,10 +327,27 @@ subroutine elsi_print_elpa_settings(ph,jh)
character(len=*), parameter :: caller = "elsi_print_elpa_settings"
call fjson_start_name_object(jh,"solver_settings")
call fjson_write_name_value(jh,"elpa_solver",ph%elpa_solver)
if(ph%elpa_solver == UNSET) then
call fjson_write_name_value(jh,"elpa_solver",2)
else
call fjson_write_name_value(jh,"elpa_solver",ph%elpa_solver)
end if
call fjson_write_name_value(jh,"elpa_n_single_precision",ph%elpa_n_single)
call fjson_write_name_value(jh,"elpa_gpu",ph%elpa_gpu)
call fjson_write_name_value(jh,"elpa_autotune",ph%elpa_autotune)
if(ph%elpa_gpu == UNSET) then
call fjson_write_name_value(jh,"elpa_gpu",0)
else
call fjson_write_name_value(jh,"elpa_gpu",ph%elpa_gpu)
end if
if(ph%elpa_autotune == UNSET) then
call fjson_write_name_value(jh,"elpa_autotune",0)
else
call fjson_write_name_value(jh,"elpa_autotune",ph%elpa_autotune)
end if
call fjson_finish_object(jh)
end subroutine
......
......@@ -18,7 +18,8 @@ module ELSI_SOLVER
use ELSI_EIGENEXA, only: elsi_init_eigenexa,elsi_solve_eigenexa
use ELSI_ELPA, only: elsi_init_elpa,elsi_solve_elpa,elsi_do_fc_elpa,&
elsi_undo_fc_elpa
use ELSI_LAPACK, only: elsi_solve_lapack
use ELSI_LAPACK, only: elsi_solve_lapack,elsi_do_fc_lapack,&
elsi_undo_fc_lapack
use ELSI_MAGMA, only: elsi_init_magma,elsi_solve_magma
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_MPI
......@@ -204,7 +205,21 @@ subroutine elsi_ev_real(eh,ham,ovlp,eval,evec)
select case(solver)
case(ELPA_SOLVER)
if(eh%ph%parallel_mode == SINGLE_PROC) then
call elsi_solve_lapack(eh%ph,eh%bh,ham,ovlp,eval,evec)
if(eh%ph%n_basis_c > 0) then
call elsi_do_fc_lapack(eh%ph,ham,ovlp,evec,eh%perm_fc)
call elsi_solve_lapack(eh%ph,eh%bh,&
ham(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
ovlp(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
eval(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
evec(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v))
call elsi_undo_fc_lapack(eh%ph,eh%bh,ham,ovlp,evec,eh%perm_fc,&
eval(1:eh%ph%n_basis_c))
else
call elsi_solve_lapack(eh%ph,eh%bh,ham,ovlp,eval,evec)
end if
else
if(eh%ph%n_basis_c > 0) then
if(.not. allocated(eh%ham_real_v)) then
......@@ -313,7 +328,21 @@ subroutine elsi_ev_complex(eh,ham,ovlp,eval,evec)
select case(eh%ph%solver)
case(ELPA_SOLVER)
if(eh%ph%parallel_mode == SINGLE_PROC) then
call elsi_solve_lapack(eh%ph,eh%bh,ham,ovlp,eval,evec)
if(eh%ph%n_basis_c > 0) then
call elsi_do_fc_lapack(eh%ph,ham,ovlp,evec,eh%perm_fc)
call elsi_solve_lapack(eh%ph,eh%bh,&
ham(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
ovlp(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
eval(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v),&
evec(eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v,&
eh%ph%n_basis_c+1:eh%ph%n_basis_c+eh%ph%n_basis_v))
call elsi_undo_fc_lapack(eh%ph,eh%bh,ham,ovlp,evec,eh%perm_fc,&
eval(1:eh%ph%n_basis_c))
else
call elsi_solve_lapack(eh%ph,eh%bh,ham,ovlp,eval,evec)
end if
else
if(eh%ph%n_basis_c > 0) then
if(.not. allocated(eh%ham_cmplx_v)) then
......
......@@ -104,7 +104,7 @@ subroutine elsi_reset_param(ph)
ph%mu_max_steps = 100
ph%mu_mp_order = 1
ph%fc_method = 2
ph%n_basis_c = UNSET
ph%n_basis_c = 0
ph%n_basis_v = UNSET
ph%n_lrow_v = UNSET
ph%n_lcol_v = UNSET
......
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