Commit 32055251 authored by Victor Yu's avatar Victor Yu
Browse files

Simplify density matrix build

Next will move part of density matrix build to GPUs.
parent 0b4901b5
......@@ -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 "20200524")
SET(elsi_DATESTAMP "20200525")
### CMake modules ###
LIST(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
......
......@@ -731,7 +731,7 @@ CONTAINS
END SUBROUTINE
SUBROUTINE sips_get_dm(ncol_l,nnz_l,col_idx,row_ptr,nev,occ,dm_val)
SUBROUTINE sips_get_dm_edm(ncol_l,nnz_l,col_idx,row_ptr,nev,occ,dm_val,i_dm)
IMPLICIT NONE
......@@ -742,6 +742,7 @@ CONTAINS
PetscInt, INTENT(IN) :: nev ! Number of eigenvalues
PetscReal, INTENT(IN) :: occ(nev) ! Occupation numbers
PetscReal, INTENT(OUT) :: dm_val(nnz_l) ! Non-zero values
PetscInt, INTENT(IN) :: i_dm ! 0: dm, 1: edm
Vec :: xr
Vec :: xrseq
......@@ -752,6 +753,7 @@ CONTAINS
PetscInt :: i_val
PetscInt :: nnz_this_row
PetscReal :: tmp
PetscReal :: eval
PetscReal, POINTER :: vec_tmp(:)
! Index conversion
......@@ -767,8 +769,16 @@ CONTAINS
dm_val = 0.0_dp
DO i = 1,nev
CALL EPSGetEigenvector(eps,i-1,xr,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
IF (i_dm == 0) THEN
CALL EPSGetEigenvector(eps,i-1,xr,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
eval = 1.0_dp
ELSE
CALL EPSGetEigenpair(eps,i-1,eval,PETSC_NULL_SCALAR,xr,&
PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
END IF
CALL VecScatterBegin(ctx,xr,xrseq,INSERT_VALUES,SCATTER_FORWARD,&
ierr)
......@@ -787,7 +797,8 @@ CONTAINS
DO k = 1,nnz_this_row
i_val = i_val+1
tmp = occ(i)*vec_tmp(j)*vec_tmp(col_idx(i_val)+1)
tmp = vec_tmp(j)*vec_tmp(col_idx(i_val)+1)
tmp = eval*occ(i)*tmp
dm_val(i_val) = dm_val(i_val)+tmp
END DO
END DO
......@@ -811,87 +822,4 @@ CONTAINS
END SUBROUTINE
SUBROUTINE sips_get_edm(ncol_l,nnz_l,col_idx,row_ptr,nev,occ,edm_val)
IMPLICIT NONE
PetscInt, INTENT(IN) :: ncol_l ! Local size
PetscInt, INTENT(IN) :: nnz_l ! Local non-zeros
PetscInt, INTENT(INOUT) :: col_idx(nnz_l) ! Column index
PetscInt, INTENT(INOUT) :: row_ptr(ncol_l+1) ! Row pointer
PetscInt, INTENT(IN) :: nev ! Number of eigenvalues
PetscReal, INTENT(IN) :: occ(nev) ! Occupation numbers
PetscReal, INTENT(OUT) :: edm_val(nnz_l) ! Non-zero values
Vec :: xr
Vec :: xrseq
VecScatter :: ctx
PetscInt :: i
PetscInt :: j
PetscInt :: k
PetscInt :: i_val
PetscInt :: nnz_this_row
PetscReal :: tmp
PetscReal :: eval
PetscReal, POINTER :: vec_tmp(:)
! Index conversion
col_idx = col_idx-1
row_ptr = row_ptr-1
CALL MatCreateVecs(math,xr,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
CALL VecScatterCreateToAll(xr,ctx,xrseq,ierr)
CHKERRQ(ierr)
edm_val = 0.0_dp
DO i = 1,nev
CALL EPSGetEigenpair(eps,i-1,eval,PETSC_NULL_SCALAR,xr,&
PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
CALL VecScatterBegin(ctx,xr,xrseq,INSERT_VALUES,SCATTER_FORWARD,&
ierr)
CHKERRQ(ierr)
CALL VecScatterEnd(ctx,xr,xrseq,INSERT_VALUES,SCATTER_FORWARD,ierr)
CHKERRQ(ierr)
CALL VecGetArrayReadF90(xrseq,vec_tmp,ierr)
CHKERRQ(ierr)
i_val = 0
DO j = istart+1,iend
nnz_this_row = row_ptr(j-istart+1)-row_ptr(j-istart)
DO k = 1,nnz_this_row
i_val = i_val+1
tmp = eval*occ(i)*vec_tmp(j)*&
vec_tmp(col_idx(i_val)+1)
edm_val(i_val) = edm_val(i_val)+tmp
END DO
END DO
CALL VecRestoreArrayReadF90(xrseq,vec_tmp,ierr)
CHKERRQ(ierr)
END DO
! Index conversion
col_idx = col_idx+1
row_ptr = row_ptr+1
CALL VecDestroy(xr,ierr)
CHKERRQ(ierr)
CALL VecDestroy(xrseq,ierr)
CHKERRQ(ierr)
CALL VecScatterDestroy(ctx,ierr)
CHKERRQ(ierr)
END SUBROUTINE
END MODULE M_SIPS
......@@ -63,10 +63,10 @@ module ELSI_CONSTANT
integer(kind=i4), parameter :: CUBIC = 3
integer(kind=i4), parameter :: COLD = 4
! PEXSI result
integer(kind=i4), parameter :: PEXSI_DM = 0
integer(kind=i4), parameter :: PEXSI_EDM = 1
integer(kind=i4), parameter :: PEXSI_FDM = 2
! Density matrix
integer(kind=i4), parameter :: GET_DM = 0
integer(kind=i4), parameter :: GET_EDM = 1
integer(kind=i4), parameter :: GET_FDM = 2
! Density matrix purification
integer(kind=i4), parameter :: NTPOLY_PM = 0
......
......@@ -10,8 +10,9 @@
module ELSI_GET
use ELSI_CONSTANT, only: PEXSI_CSC,SIESTA_CSC,GENERIC_COO,ELPA_SOLVER,&
OMM_SOLVER,PEXSI_SOLVER,EIGENEXA_SOLVER,SIPS_SOLVER,NTPOLY_SOLVER
OMM_SOLVER,PEXSI_SOLVER,EIGENEXA_SOLVER,SIPS_SOLVER,NTPOLY_SOLVER,GET_EDM
use ELSI_DATATYPE, only: elsi_handle
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_MPI, only: elsi_stop
use ELSI_NTPOLY, only: elsi_compute_edm_ntpoly
use ELSI_OMM, only: elsi_compute_edm_omm
......@@ -22,8 +23,8 @@ module ELSI_GET
elsi_ntpoly_to_siesta_dm,elsi_ntpoly_to_sips_dm,elsi_pexsi_to_blacs_dm,&
elsi_pexsi_to_generic_dm,elsi_pexsi_to_siesta_dm,elsi_sips_to_blacs_dm,&
elsi_sips_to_generic_dm,elsi_sips_to_siesta_dm
use ELSI_SIPS, only: elsi_build_edm_sips
use ELSI_UTIL, only: elsi_check_init,elsi_reduce_energy,elsi_build_edm
use ELSI_SIPS, only: elsi_build_dm_edm_sips
use ELSI_UTIL, only: elsi_check_init,elsi_reduce_energy,elsi_build_dm_edm
implicit none
......@@ -277,6 +278,8 @@ subroutine elsi_get_edm_real(eh,edm)
integer(kind=i4) :: solver_save
character(len=200) :: msg
real(kind=r8), allocatable :: factor(:)
character(len=*), parameter :: caller = "elsi_get_edm_real"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -301,8 +304,13 @@ subroutine elsi_get_edm_real(eh,edm)
select case(eh%ph%solver)
case(ELPA_SOLVER,EIGENEXA_SOLVER)
call elsi_build_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%eval(1:eh%ph%n_states),eh%evec_real,edm)
call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
factor(:) = -eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt)*eh%eval(1:eh%ph%n_states)
call elsi_build_dm_edm(eh%ph,eh%bh,factor,eh%evec_real,edm,GET_EDM)
call elsi_deallocate(eh%bh,factor,"factor")
case(OMM_SOLVER)
call elsi_compute_edm_omm(eh%ph,eh%bh,eh%omm_c_real,edm)
case(PEXSI_SOLVER)
......@@ -310,8 +318,8 @@ subroutine elsi_get_edm_real(eh,edm)
call elsi_pexsi_to_blacs_dm(eh%ph,eh%bh,eh%dm_real_sp,eh%row_ind_sp1,&
eh%col_ptr_sp1,edm)
case(SIPS_SOLVER)
call elsi_build_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp)
call elsi_build_dm_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp,GET_EDM)
call elsi_sips_to_blacs_dm(eh%ph,eh%bh,eh%dm_real_sp,eh%row_ind_sp1,&
eh%col_ptr_sp1,edm)
case(NTPOLY_SOLVER)
......@@ -340,6 +348,8 @@ subroutine elsi_get_edm_real_sparse(eh,edm)
integer(kind=i4) :: solver_save
character(len=200) :: msg
real(kind=r8), allocatable :: factor(:)
character(len=*), parameter :: caller = "elsi_get_edm_real_sparse"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -364,8 +374,14 @@ subroutine elsi_get_edm_real_sparse(eh,edm)
select case(eh%ph%solver)
case(ELPA_SOLVER,EIGENEXA_SOLVER)
call elsi_build_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%eval(1:eh%ph%n_states),eh%evec_real,eh%dm_real_den)
call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
factor(:) = -eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt)*eh%eval(1:eh%ph%n_states)
call elsi_build_dm_edm(eh%ph,eh%bh,factor,eh%evec_real,eh%dm_real_den,&
GET_EDM)
call elsi_deallocate(eh%bh,factor,"factor")
select case(eh%ph%matrix_format)
case(PEXSI_CSC)
......@@ -417,16 +433,16 @@ subroutine elsi_get_edm_real_sparse(eh,edm)
case(SIPS_SOLVER)
select case(eh%ph%matrix_format)
case(PEXSI_CSC)
call elsi_build_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),edm)
call elsi_build_dm_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),edm,GET_EDM)
case(SIESTA_CSC)
call elsi_build_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp)
call elsi_build_dm_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp,GET_EDM)
call elsi_sips_to_siesta_dm(eh%ph,eh%bh,eh%dm_real_sp,eh%row_ind_sp1,&
eh%col_ptr_sp1,edm,eh%row_ind_sp2,eh%col_ptr_sp2)
case(GENERIC_COO)
call elsi_build_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp)
call elsi_build_dm_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp,GET_EDM)
call elsi_sips_to_generic_dm(eh%ph,eh%bh,eh%dm_real_sp,eh%row_ind_sp1,&
eh%col_ptr_sp1,eh%map_sp1,edm,eh%perm_sp3)
case default
......@@ -470,6 +486,8 @@ subroutine elsi_get_edm_complex(eh,edm)
integer(kind=i4) :: solver_save
character(len=200) :: msg
real(kind=r8), allocatable :: factor(:)
character(len=*), parameter :: caller = "elsi_get_edm_complex"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -489,8 +507,13 @@ subroutine elsi_get_edm_complex(eh,edm)
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_build_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%eval(1:eh%ph%n_states),eh%evec_cmplx,edm)
call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
factor(:) = -eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt)*eh%eval(1:eh%ph%n_states)
call elsi_build_dm_edm(eh%ph,eh%bh,factor,eh%evec_cmplx,edm,GET_EDM)
call elsi_deallocate(eh%bh,factor,"factor")
case(OMM_SOLVER)
call elsi_compute_edm_omm(eh%ph,eh%bh,eh%omm_c_cmplx,edm)
case(PEXSI_SOLVER)
......@@ -523,6 +546,8 @@ subroutine elsi_get_edm_complex_sparse(eh,edm)
integer(kind=i4) :: solver_save
character(len=200) :: msg
real(kind=r8), allocatable :: factor(:)
character(len=*), parameter :: caller = "elsi_get_edm_complex_sparse"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -542,8 +567,14 @@ subroutine elsi_get_edm_complex_sparse(eh,edm)
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_build_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%eval(1:eh%ph%n_states),eh%evec_cmplx,eh%dm_cmplx_den)
call elsi_allocate(eh%bh,factor,eh%ph%n_states,"factor",caller)
factor(:) = -eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt)*eh%eval(1:eh%ph%n_states)
call elsi_build_dm_edm(eh%ph,eh%bh,factor,eh%evec_cmplx,eh%dm_cmplx_den,&
GET_EDM)
call elsi_deallocate(eh%bh,factor,"factor")
select case(eh%ph%matrix_format)
case(PEXSI_CSC)
......
......@@ -9,7 +9,7 @@
!!
module ELSI_PEXSI
use ELSI_CONSTANT, only: UNSET,PEXSI_CSC,PEXSI_DM,PEXSI_EDM,PEXSI_FDM
use ELSI_CONSTANT, only: UNSET,PEXSI_CSC,GET_DM,GET_EDM,GET_FDM
use ELSI_DATATYPE, only: elsi_param_t,elsi_basic_t
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_MPI, only: MPI_SUM,MPI_REAL8,MPI_COMPLEX16
......@@ -408,7 +408,7 @@ subroutine elsi_solve_pexsi_real(ph,bh,row_ind,col_ptr,ne_vec,ham,ovlp,dm)
! Get free energy density matrix
call elsi_get_time(t0)
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_FDM,ne_vec,dm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_FDM,ne_vec,dm)
! Compute energy = Tr(S*FDM)
if(ph%pexsi_my_prow == 0) then
......@@ -428,7 +428,7 @@ subroutine elsi_solve_pexsi_real(ph,bh,row_ind,col_ptr,ne_vec,ham,ovlp,dm)
! Get density matrix
call elsi_get_time(t0)
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_DM,ne_vec,dm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_DM,ne_vec,dm)
! Compute energy = Tr(H*DM)
if(ph%pexsi_my_prow == 0) then
......@@ -481,7 +481,7 @@ subroutine elsi_compute_edm_pexsi_real(ph,bh,ne_vec,edm)
call elsi_get_time(t0)
! Get energy density matrix
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_EDM,ne_vec,edm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_EDM,ne_vec,edm)
call elsi_get_time(t1)
......@@ -527,12 +527,12 @@ subroutine elsi_retrieve_dm_pexsi_real(ph,bh,which,ne_vec,dm)
call elsi_allocate(bh,tmp,bh%nnz_l_sp1,"tmp",caller)
select case(which)
case(PEXSI_DM)
case(GET_DM)
call f_ppexsi_retrieve_real_dm(ph%pexsi_plan,tmp,local_energy,ierr)
case(PEXSI_EDM)
case(GET_EDM)
call f_ppexsi_retrieve_real_edm(ph%pexsi_plan,ph%pexsi_options,tmp,&
local_energy,ierr)
case(PEXSI_FDM)
case(GET_FDM)
call f_ppexsi_retrieve_real_fdm(ph%pexsi_plan,ph%pexsi_options,tmp,&
local_energy,ierr)
end select
......@@ -884,7 +884,7 @@ subroutine elsi_solve_pexsi_cmplx(ph,bh,row_ind,col_ptr,ne_vec,ham,ovlp,dm)
! Get free energy density matrix
call elsi_get_time(t0)
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_FDM,ne_vec,dm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_FDM,ne_vec,dm)
! Compute energy = Tr(S*FDM)
if(ph%pexsi_my_prow == 0) then
......@@ -905,7 +905,7 @@ subroutine elsi_solve_pexsi_cmplx(ph,bh,row_ind,col_ptr,ne_vec,ham,ovlp,dm)
! Get density matrix
call elsi_get_time(t0)
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_DM,ne_vec,dm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_DM,ne_vec,dm)
! Compute energy = Tr(H*DM)
if(ph%pexsi_my_prow == 0) then
......@@ -959,7 +959,7 @@ subroutine elsi_compute_edm_pexsi_cmplx(ph,bh,ne_vec,edm)
call elsi_get_time(t0)
! Get energy density matrix
call elsi_retrieve_dm_pexsi(ph,bh,PEXSI_EDM,ne_vec,edm)
call elsi_retrieve_dm_pexsi(ph,bh,GET_EDM,ne_vec,edm)
call elsi_get_time(t1)
......@@ -1005,12 +1005,12 @@ subroutine elsi_retrieve_dm_pexsi_cmplx(ph,bh,which,ne_vec,dm)
call elsi_allocate(bh,tmp,bh%nnz_l_sp1,"tmp",caller)
select case(which)
case(PEXSI_DM)
case(GET_DM)
call f_ppexsi_retrieve_complex_dm(ph%pexsi_plan,tmp,local_energy,ierr)
case(PEXSI_EDM)
case(GET_EDM)
call f_ppexsi_retrieve_complex_edm(ph%pexsi_plan,ph%pexsi_options,tmp,&
local_energy,ierr)
case(PEXSI_FDM)
case(GET_FDM)
call f_ppexsi_retrieve_complex_fdm(ph%pexsi_plan,ph%pexsi_options,tmp,&
local_energy,ierr)
end select
......
......@@ -9,7 +9,7 @@
!!
module ELSI_SIPS
use ELSI_CONSTANT, only: UNSET
use ELSI_CONSTANT, only: UNSET,GET_DM
use ELSI_DATATYPE, only: elsi_param_t,elsi_basic_t
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_MPI, only: elsi_stop
......@@ -18,8 +18,8 @@ module ELSI_SIPS
use M_SIPS, only: sips_initialize,sips_load_ham_ovlp,sips_load_ham,&
sips_update_ham,sips_set_eps,sips_update_eps,sips_set_slices,&
sips_get_slices,sips_get_inertias,sips_get_slices_from_inertias,&
sips_solve_eps,sips_get_eigenvalues,sips_get_eigenvectors,sips_get_dm,&
sips_get_edm,sips_finalize
sips_solve_eps,sips_get_eigenvalues,sips_get_eigenvectors,&
sips_get_dm_edm,sips_finalize
implicit none
......@@ -28,19 +28,14 @@ module ELSI_SIPS
public :: elsi_init_sips
public :: elsi_cleanup_sips
public :: elsi_solve_sips
public :: elsi_build_dm_sips
public :: elsi_build_edm_sips
public :: elsi_build_dm_edm_sips
interface elsi_solve_sips
module procedure elsi_solve_sips_real
end interface
interface elsi_build_dm_sips
module procedure elsi_build_dm_sips_real
end interface
interface elsi_build_edm_sips
module procedure elsi_build_edm_sips_real
interface elsi_build_dm_edm_sips
module procedure elsi_build_dm_edm_sips_real
end interface
contains
......@@ -262,9 +257,9 @@ subroutine elsi_solve_sips_real(ph,bh,row_ind,col_ptr,ham,ovlp,eval,evec)
end subroutine
!>
!! Construct the density matrix.
!! Construct density matrix or energy-weighted density matrix from eigenvectors.
!!
subroutine elsi_build_dm_sips_real(ph,bh,row_ind,col_ptr,occ,dm)
subroutine elsi_build_dm_edm_sips_real(ph,bh,row_ind,col_ptr,occ,dm,which)
implicit none
......@@ -274,55 +269,27 @@ subroutine elsi_build_dm_sips_real(ph,bh,row_ind,col_ptr,occ,dm)
integer(kind=i4), intent(inout) :: col_ptr(bh%n_lcol_sp1+1)
real(kind=r8), intent(in) :: occ(ph%n_states)
real(kind=r8), intent(out) :: dm(bh%nnz_l_sp1)
integer(kind=i4), intent(in) :: which
real(kind=r8) :: t0
real(kind=r8) :: t1
character(len=200) :: msg
character(len=*), parameter :: caller = "elsi_build_dm_sips_real"
character(len=*), parameter :: caller = "elsi_build_dm_edm_sips_real"
call elsi_get_time(t0)
call sips_get_dm(bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,col_ptr,ph%n_states,occ,&
dm)
call sips_get_dm_edm(bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,col_ptr,ph%n_states,&
occ,dm,which)
call elsi_get_time(t1)
write(msg,"(A)") "Finished density matrix calculation"
call elsi_say(bh,msg)
write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
end subroutine
!>
!! Construct the energy-weighted density matrix.
!!
subroutine elsi_build_edm_sips_real(ph,bh,row_ind,col_ptr,occ,edm)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_basic_t), intent(in) :: bh
integer(kind=i4), intent(inout) :: row_ind(bh%nnz_l_sp1)
integer(kind=i4), intent(inout) :: col_ptr(bh%n_lcol_sp1+1)
real(kind=r8), intent(in) :: occ(ph%n_states)
real(kind=r8), intent(out) :: edm(bh%nnz_l_sp1)
real(kind=r8) :: t0
real(kind=r8) :: t1
character(len=200) :: msg
character(len=*), parameter :: caller = "elsi_build_edm_sips_real"
call elsi_get_time(t0)
call sips_get_edm(bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,col_ptr,ph%n_states,&
occ,edm)
call elsi_get_time(t1)
if(which == GET_DM) then
write(msg,"(A)") "Finished density matrix calculation"
else
write(msg,"(A)") "Finished energy density matrix calculation"
end if
write(msg,"(A)") "Finished energy density matrix calculation"
call elsi_say(bh,msg)
write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
......
......@@ -12,7 +12,7 @@ module ELSI_SOLVER
use ELSI_BSEPACK, only: elsi_solve_bsepack
use ELSI_CONSTANT, only: ELPA_SOLVER,OMM_SOLVER,PEXSI_SOLVER,SIPS_SOLVER,&
NTPOLY_SOLVER,EIGENEXA_SOLVER,MAGMA_SOLVER,BSEPACK_SOLVER,MULTI_PROC,&
SINGLE_PROC,PEXSI_CSC,SIESTA_CSC,GENERIC_COO
SINGLE_PROC,PEXSI_CSC,SIESTA_CSC,GENERIC_COO,GET_DM
use ELSI_DATATYPE, only: elsi_handle,elsi_param_t,elsi_basic_t
use ELSI_DECISION, only: elsi_decide_ev,elsi_decide_dm
use ELSI_EIGENEXA, only: elsi_init_eigenexa,elsi_solve_eigenexa
......@@ -43,9 +43,9 @@ module ELSI_SOLVER
elsi_sips_to_blacs_ev,elsi_sips_to_blacs_hs,elsi_sips_to_generic_dm,&
elsi_sips_to_ntpoly_hs,elsi_sips_to_siesta_dm
use ELSI_SETUP, only: elsi_set_blacs
use ELSI_SIPS, only: elsi_init_sips,elsi_solve_sips,elsi_build_dm_sips
use ELSI_SIPS, only: elsi_init_sips,elsi_solve_sips,elsi_build_dm_edm_sips
use ELSI_UTIL, only: elsi_check,elsi_check_init,elsi_reduce_energy,&
elsi_build_dm
elsi_build_dm_edm
implicit none
......@@ -652,8 +652,8 @@ subroutine elsi_dm_real(eh,ham,ovlp,dm,ebs)
call elsi_solve_elpa(eh%ph,eh%bh,ham,ovlp,eh%eval,eh%evec_real)
call elsi_get_occ_for_dm(eh%ph,eh%bh,eh%eval,eh%occ)
call elsi_build_dm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_real,dm)
call elsi_build_dm_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_real,dm,GET_DM)
call elsi_get_band_energy(eh%ph,eh%bh,ebs,solver)
eh%ph%eval_ready = .true.
......@@ -760,8 +760,8 @@ subroutine elsi_dm_real(eh,ham,ovlp,dm,ebs)
call elsi_solve_eigenexa(eh%ph,eh%bh,ham,ovlp,eh%eval,eh%evec_real)
call elsi_get_occ_for_dm(eh%ph,eh%bh,eh%eval,eh%occ)
call elsi_build_dm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_real,dm)
call elsi_build_dm_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_real,dm,GET_DM)
call elsi_get_band_energy(eh%ph,eh%bh,ebs,solver)
eh%ph%eval_ready = .true.
......@@ -826,8 +826,8 @@ subroutine elsi_dm_real(eh,ham,ovlp,dm,ebs)
call elsi_solve_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%ham_real_sp,eh%ovlp_real_sp,eh%eval,eh%evec_real)
call elsi_get_occ_for_dm(eh%ph,eh%bh,eh%eval,eh%occ)
call elsi_build_dm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp)
call elsi_build_dm_edm_sips(eh%ph,eh%bh,eh%row_ind_sp1,eh%col_ptr_sp1,&
eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),eh%dm_real_sp,GET_DM)
call elsi_sips_to_blacs_dm(eh%ph,eh%bh,eh%dm_real_sp,eh%row_ind_sp1,&
eh%col_ptr_sp1,dm)
call elsi_get_band_energy(eh%ph,eh%bh,ebs,solver)
......@@ -935,8 +935,8 @@ subroutine elsi_dm_complex(eh,ham,ovlp,dm,ebs)
call elsi_solve_elpa(eh%ph,eh%bh,ham,ovlp,eh%eval,eh%evec_cmplx)
call elsi_get_occ_for_dm(eh%ph,eh%bh,eh%eval,eh%occ)
call elsi_build_dm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_cmplx,dm)
call elsi_build_dm_edm(eh%ph,eh%bh,eh%occ(:,eh%ph%i_spin,eh%ph%i_kpt),&
eh%evec_cmplx,dm,GET_DM)
call elsi_get_band_energy(eh%ph,eh%bh,ebs,