Commit 330d2730 authored by Victor Yu's avatar Victor Yu

Add density matrix extrapolation by purification

parent 3f26170a
......@@ -380,7 +380,7 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(in) :: bh
integer(kind=i4), intent(in) :: row_map(ph%n_basis)
integer(kind=i4), intent(in) :: col_map(ph%n_basis)
......@@ -398,11 +398,13 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -448,6 +450,8 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
ph%elpa_first = .false.
end subroutine
!>
......@@ -728,7 +732,7 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(in) :: bh
integer(kind=i4), intent(in) :: row_map(ph%n_basis)
integer(kind=i4), intent(in) :: col_map(ph%n_basis)
......@@ -746,11 +750,13 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -796,6 +802,8 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
ph%elpa_first = .false.
end subroutine
!>
......
......@@ -395,7 +395,7 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(in) :: bh
integer(kind=i4), intent(in) :: row_map(ph%n_basis)
integer(kind=i4), intent(in) :: col_map(ph%n_basis)
......@@ -413,11 +413,13 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -463,6 +465,8 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
ph%elpa_first = .false.
end subroutine
!>
......@@ -745,7 +749,7 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
implicit none
type(elsi_param_t), intent(in) :: ph
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(in) :: bh
integer(kind=i4), intent(in) :: row_map(ph%n_basis)
integer(kind=i4), intent(in) :: col_map(ph%n_basis)
......@@ -763,11 +767,13 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -813,6 +819,8 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
ph%elpa_first = .false.
end subroutine
!>
......
......@@ -10,10 +10,14 @@
!!
module ELSI_TOOLS
use ELSI_CONSTANTS, only: ELPA_SOLVER,OMM_SOLVER,SIPS_SOLVER
use ELSI_DATATYPE, only: elsi_handle
use ELSI_ELPA, only: elsi_update_dm_elpa
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_solve_ntpoly
use ELSI_OCC, only: elsi_mu_and_occ,elsi_entropy
use ELSI_PRECISION, only: i4,r8
use ELSI_REDIST, only: elsi_blacs_to_ntpoly_hs,elsi_ntpoly_to_blacs_dm
use ELSI_UTILS, only: elsi_check_init,elsi_set_full_mat,elsi_build_dm,&
elsi_build_edm,elsi_gram_schmidt
......@@ -79,16 +83,56 @@ subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
implicit none
type(elsi_handle), intent(in) :: eh !< Handle
type(elsi_handle), intent(inout) :: eh !< Handle
real(kind=r8), intent(inout) :: ovlp0(eh%bh%n_lrow,eh%bh%n_lcol) !< Old overlap
real(kind=r8), intent(inout) :: ovlp1(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
real(kind=r8), intent(inout) :: dm(eh%bh%n_lrow,eh%bh%n_lcol) !< Density matrix
integer(kind=i4) :: solver_save
logical :: unit_ovlp_save
real(kind=r8), allocatable :: dummy(:,:)
character(len=*), parameter :: caller = "elsi_extrapolate_dm_real"
call elsi_check_init(eh%bh,eh%handle_init,caller)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
if(eh%ph%solver == SIPS_SOLVER .and. eh%ph%sips_n_elpa > 0) then
solver_save = SIPS_SOLVER
eh%ph%solver = ELPA_SOLVER
end if
if(eh%ph%solver == OMM_SOLVER .and. eh%ph%omm_n_elpa > 0) then
solver_save = OMM_SOLVER
eh%ph%solver = ELPA_SOLVER
end if
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
case default
call elsi_init_ntpoly(eh%ph,eh%bh)
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_allocate(eh%bh,dummy,1,1,"dummy",caller)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,dummy,eh%ph%nt_ham,&
eh%ph%nt_dm)
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
eh%ph%nt_dm)
call elsi_deallocate(eh%bh,dummy,"dummy")
eh%ph%unit_ovlp = unit_ovlp_save
call elsi_solve_ntpoly(eh%ph,eh%bh,eh%ph%nt_ham,eh%ph%nt_ovlp,eh%ph%nt_dm)
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_dm,dm)
end select
eh%ph%solver = solver_save
end subroutine
......@@ -99,16 +143,51 @@ subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
implicit none
type(elsi_handle), intent(in) :: eh !< Handle
type(elsi_handle), intent(inout) :: eh !< Handle
complex(kind=r8), intent(inout) :: ovlp0(eh%bh%n_lrow,eh%bh%n_lcol) !< Old overlap
complex(kind=r8), intent(inout) :: ovlp1(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
complex(kind=r8), intent(inout) :: dm(eh%bh%n_lrow,eh%bh%n_lcol) !< Density matrix
integer(kind=i4) :: solver_save
logical :: unit_ovlp_save
complex(kind=r8), allocatable :: dummy(:,:)
character(len=*), parameter :: caller = "elsi_extrapolate_dm_complex"
call elsi_check_init(eh%bh,eh%handle_init,caller)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
if(eh%ph%solver == OMM_SOLVER .and. eh%ph%omm_n_elpa > 0) then
solver_save = OMM_SOLVER
eh%ph%solver = ELPA_SOLVER
end if
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
case default
call elsi_init_ntpoly(eh%ph,eh%bh)
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_allocate(eh%bh,dummy,1,1,"dummy",caller)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,dummy,eh%ph%nt_ham,&
eh%ph%nt_dm)
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
eh%ph%nt_dm)
call elsi_deallocate(eh%bh,dummy,"dummy")
eh%ph%unit_ovlp = unit_ovlp_save
call elsi_solve_ntpoly(eh%ph,eh%bh,eh%ph%nt_ham,eh%ph%nt_ovlp,eh%ph%nt_dm)
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_dm,dm)
end select
eh%ph%solver = solver_save
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