Commit f03baee7 authored by Victor Yu's avatar Victor Yu

Merge branch 'geometry_optimization' into 'master'

Geometry optimization helpers

See merge request elsi-devel/elsi-interface!84
parents 52e1318d ca6b706f
......@@ -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 "20181107")
SET(elsi_DATESTAMP "20181114")
### CMake modules ###
LIST(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
......
MODULE NTPoly
USE DensityMatrixSolversModule
USE GeometryOptimizationModule
USE LoggingModule
USE PermutationModule
USE PMatrixMemoryPoolModule
......
......@@ -26,7 +26,7 @@ CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> Deactivate the logger.
SUBROUTINE DeactivateLogger
IsActive = .TRUE.
IsActive = .FALSE.
END SUBROUTINE DeactivateLogger
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> Call this subroutine when you enter into a section with verbose output
......
......@@ -302,13 +302,11 @@ void c_elsi_orthonormalize_ev_complex(elsi_handle handle_c,
double _Complex *evec);
void c_elsi_extrapolate_dm_real(elsi_handle handle_c,
double *ovlp0,
double *ovlp1,
double *ovlp,
double *dm);
void c_elsi_extrapolate_dm_complex(elsi_handle handle_c,
double _Complex *ovlp0,
double _Complex *ovlp1,
double _Complex *ovlp,
double _Complex *dm);
void c_elsi_construct_dm_real(elsi_handle handle_c,
......
......@@ -398,13 +398,11 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -750,13 +748,11 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......
......@@ -1545,19 +1545,17 @@ subroutine elsi_orthonormalize_ev_complex_c_wrapper(h_c,ovlp_c,evec_c)&
end subroutine
subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp_c,dm_c)&
bind(C,name="c_elsi_extrapolate_dm_real")
implicit none
type(c_ptr), value, intent(in) :: h_c
type(c_ptr), value, intent(in) :: ovlp0_c
type(c_ptr), value, intent(in) :: ovlp1_c
type(c_ptr), value, intent(in) :: ovlp_c
type(c_ptr), value, intent(in) :: dm_c
type(elsi_handle), pointer :: h_f
real(kind=c_double), pointer :: ovlp0_f(:,:)
real(kind=c_double), pointer :: ovlp1_f(:,:)
real(kind=c_double), pointer :: ovlp_f(:,:)
real(kind=c_double), pointer :: dm_f(:,:)
integer(kind=c_int) :: lrow
......@@ -1568,27 +1566,24 @@ subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
lrow = h_f%bh%n_lrow
lcol = h_f%bh%n_lcol
call c_f_pointer(ovlp0_c,ovlp0_f,shape=[lrow,lcol])
call c_f_pointer(ovlp1_c,ovlp1_f,shape=[lrow,lcol])
call c_f_pointer(ovlp_c,ovlp_f,shape=[lrow,lcol])
call c_f_pointer(dm_c,dm_f,shape=[lrow,lcol])
call elsi_extrapolate_dm_real(h_f,ovlp0_f,ovlp1_f,dm_f)
call elsi_extrapolate_dm_real(h_f,ovlp_f,dm_f)
end subroutine
subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp_c,dm_c)&
bind(C,name="c_elsi_extrapolate_dm_complex")
implicit none
type(c_ptr), value, intent(in) :: h_c
type(c_ptr), value, intent(in) :: ovlp0_c
type(c_ptr), value, intent(in) :: ovlp1_c
type(c_ptr), value, intent(in) :: ovlp_c
type(c_ptr), value, intent(in) :: dm_c
type(elsi_handle), pointer :: h_f
complex(kind=c_double), pointer :: ovlp0_f(:,:)
complex(kind=c_double), pointer :: ovlp1_f(:,:)
complex(kind=c_double), pointer :: ovlp_f(:,:)
complex(kind=c_double), pointer :: dm_f(:,:)
integer(kind=c_int) :: lrow
......@@ -1599,11 +1594,10 @@ subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
lrow = h_f%bh%n_lrow
lcol = h_f%bh%n_lcol
call c_f_pointer(ovlp0_c,ovlp0_f,shape=[lrow,lcol])
call c_f_pointer(ovlp1_c,ovlp1_f,shape=[lrow,lcol])
call c_f_pointer(ovlp_c,ovlp_f,shape=[lrow,lcol])
call c_f_pointer(dm_c,dm_f,shape=[lrow,lcol])
call elsi_extrapolate_dm_complex(h_f,ovlp0_f,ovlp1_f,dm_f)
call elsi_extrapolate_dm_complex(h_f,ovlp_f,dm_f)
end subroutine
......
......@@ -84,6 +84,7 @@ module ELSI_DATATYPE
integer(kind=i4) :: n_calls_all ! Total number of calls
! Overlap
logical :: save_ovlp
logical :: unit_ovlp
logical :: ill_ovlp
logical :: ill_check
......@@ -193,6 +194,7 @@ module ELSI_DATATYPE
logical :: nt_started = .false.
type(Matrix_ps) :: nt_ham
type(Matrix_ps) :: nt_ovlp
type(Matrix_ps) :: nt_ovlp_copy
type(Matrix_ps) :: nt_dm
type(SolverParameters_t) :: nt_options
type(Permutation_t) :: nt_perm
......
......@@ -413,13 +413,11 @@ subroutine elsi_update_dm_elpa_real(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......@@ -767,13 +765,11 @@ subroutine elsi_update_dm_elpa_cmplx(ph,bh,row_map,col_map,ovlp0,ovlp1,dm)
call elsi_get_time(t0)
if(ph%elpa_first) then
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1
call elsi_elpa_cholesky(ph,bh,ovlp1)
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
end if
! ovlp1 = U_1^(-1)
call elsi_elpa_invert(ph,bh,ovlp1)
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
......
......@@ -43,6 +43,7 @@ module ELSI
public :: elsi_set_output_log
public :: elsi_set_output_tag
public :: elsi_set_uuid
public :: elsi_set_save_ovlp
public :: elsi_set_unit_ovlp
public :: elsi_set_zero_def
public :: elsi_set_illcond_check
......
......@@ -154,6 +154,7 @@ subroutine elsi_to_original_ev_sp_real(ph,bh,ovlp,evec)
if(ph%ill_ovlp) then
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
tmp = evec
call dgemm("N","N",ph%n_basis,ph%n_states_solve,ph%n_good,1.0_r8,ovlp,&
......@@ -485,6 +486,7 @@ subroutine elsi_to_original_ev_sp_cmplx(ph,bh,ovlp,evec)
if(ph%ill_ovlp) then
call elsi_allocate(bh,tmp,bh%n_lrow,bh%n_lcol,"tmp",caller)
tmp = evec
call zgemm("N","N",ph%n_basis,ph%n_states_solve,ph%n_good,&
......
......@@ -33,6 +33,7 @@ module ELSI_MUTATOR
public :: elsi_set_output_log
public :: elsi_set_output_tag
public :: elsi_set_uuid
public :: elsi_set_save_ovlp
public :: elsi_set_unit_ovlp
public :: elsi_set_zero_def
public :: elsi_set_illcond_check
......@@ -234,6 +235,29 @@ subroutine elsi_set_uuid(eh,uuid)
end subroutine
!>
!! This routine sets whether to save the overlap matrix or its Cholesky factor
!! for density matrix extrapolation.
!!
subroutine elsi_set_save_ovlp(eh,save_ovlp)
implicit none
type(elsi_handle), intent(inout) :: eh !< Handle
integer(kind=i4), intent(in) :: save_ovlp !< Save overlap?
character(len=*), parameter :: caller = "elsi_set_save_ovlp"
call elsi_check_init(eh%bh,eh%handle_init,caller)
if(save_ovlp == 0) then
eh%ph%save_ovlp = .false.
else
eh%ph%save_ovlp = .true.
end if
end subroutine
!>
!! This routine sets the overlap matrix to be identity.
!!
......
......@@ -9,15 +9,16 @@
!!
module ELSI_NTPOLY
use ELSI_CONSTANTS, only: NTPOLY_PM,NTPOLY_TC2,NTPOLY_TRS4,NTPOLY_HPCP
use ELSI_CONSTANTS, only: NTPOLY_SOLVER,NTPOLY_PM,NTPOLY_TC2,NTPOLY_TRS4,&
NTPOLY_HPCP
use ELSI_DATATYPE, only: elsi_param_t,elsi_basic_t
use ELSI_IO, only: elsi_say,elsi_get_time
use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
use ELSI_MPI, only: elsi_check_mpi,mpi_logical
use ELSI_PRECISION, only: r8,i4
use NTPOLY, only: PM,TRS2,TRS4,HPCP,EnergyDensityMatrix,Matrix_ps,&
ConstructEmptyMatrix,DestructMatrix,CopyMatrix,ScaleMatrix,&
FillMatrixFromTripletList,GetMatrixTripletList,FilterMatrix,&
use NTPOLY, only: PM,TRS2,TRS4,HPCP,EnergyDensityMatrix,LowdinExtrapolate,&
PurificationExtrapolate,Matrix_ps,ConstructEmptyMatrix,DestructMatrix,&
CopyMatrix,ScaleMatrix,FillMatrixFromTripletList,GetMatrixTripletList,&
ProcessGrid_t,ConstructNewProcessGrid,DestructProcessGrid,&
ConstructRandomPermutation,DestructPermutation,InverseSquareRoot,&
SolverParameters_t,Triplet_r,Triplet_c,TripletList_r,TripletList_c,&
......@@ -32,8 +33,11 @@ module ELSI_NTPOLY
public :: elsi_cleanup_ntpoly
public :: elsi_solve_ntpoly
public :: elsi_compute_edm_ntpoly
public :: elsi_update_dm_ntpoly
public :: Matrix_ps
public :: ConstructEmptyMatrix
public :: CopyMatrix
public :: DestructMatrix
public :: FillMatrixFromTripletList
public :: GetMatrixTripletList
public :: TripletList_r
......@@ -83,6 +87,8 @@ subroutine elsi_init_ntpoly(ph,bh)
call elsi_check_mpi(bh,"MPI_Bcast",ierr,caller)
ph%nt_filter = max(ph%nt_filter,bh%def0)
if(ph%nt_output .and. bh%myid_all == 0) then
call ActivateLogger()
end if
......@@ -117,6 +123,7 @@ subroutine elsi_solve_ntpoly(ph,bh,ham,ovlp,dm)
if(ph%nt_first) then
call elsi_get_time(t0)
call DestructPermutation(ph%nt_perm)
call ConstructRandomPermutation(ph%nt_perm,ovlp%logical_matrix_dimension)
ph%nt_options = SolverParameters_t(ph%nt_tol,ph%nt_filter,ph%nt_max_iter,&
......@@ -141,8 +148,6 @@ subroutine elsi_solve_ntpoly(ph,bh,ham,ovlp,dm)
call elsi_get_time(t0)
call ConstructEmptyMatrix(dm,ham)
ne = int(ph%n_electrons,kind=i4)
select case(ph%nt_method)
......@@ -158,10 +163,6 @@ subroutine elsi_solve_ntpoly(ph,bh,ham,ovlp,dm)
call ScaleMatrix(dm,ph%spin_degen)
if(ph%nt_filter < bh%def0) then
call FilterMatrix(dm,bh%def0)
end if
call elsi_get_time(t1)
write(msg,"(2X,A)") "Finished density matrix purification"
......@@ -198,19 +199,55 @@ subroutine elsi_compute_edm_ntpoly(ph,bh,ham,edm)
factor = 1.0_r8/ph%spin_degen
call ConstructEmptyMatrix(tmp,ham)
call CopyMatrix(edm,tmp)
call EnergyDensityMatrix(ham,tmp,edm,ph%nt_filter)
call DestructMatrix(tmp)
call ScaleMatrix(edm,factor)
if(ph%nt_filter < bh%def0) then
call FilterMatrix(edm,bh%def0)
call elsi_get_time(t1)
write(msg,"(2X,A)") "Finished energy density matrix calculation"
call elsi_say(bh,msg)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
end subroutine
!>
!! This routine extrapolates density matrix using Lowdin decomposition of the
!! old and new overlap matrices or 2nd order trace resetting purification.
!!
subroutine elsi_update_dm_ntpoly(ph,bh,ovlp0,ovlp1,dm0,dm1)
implicit none
type(elsi_param_t), intent(inout) :: ph
type(elsi_basic_t), intent(in) :: bh
type(Matrix_ps), intent(inout) :: ovlp0
type(Matrix_ps), intent(inout) :: ovlp1
type(Matrix_ps), intent(inout) :: dm0
type(Matrix_ps), intent(inout) :: dm1
integer(kind=i4) :: ne
real(kind=r8) :: t0
real(kind=r8) :: t1
character(len=200) :: msg
character(len=*), parameter :: caller = "elsi_update_dm_ntpoly"
call elsi_get_time(t0)
ne = int(ph%n_electrons,kind=i4)
if(ph%solver /= NTPOLY_SOLVER) then
call PurificationExtrapolate(dm0,ovlp1,ne,dm1,ph%nt_options)
else
call LowdinExtrapolate(dm0,ovlp0,ovlp1,dm1,ph%nt_options)
end if
call elsi_get_time(t1)
write(msg,"(2X,A)") "Finished energy density matrix calculation"
write(msg,"(2X,A)") "Finished density matrix extrapolation"
call elsi_say(bh,msg)
write(msg,"(2X,A,F10.3,A)") "| Time :",t1-t0," s"
call elsi_say(bh,msg)
......@@ -232,6 +269,7 @@ subroutine elsi_cleanup_ntpoly(ph)
call DeactivateLogger()
call DestructMatrix(ph%nt_ham)
call DestructMatrix(ph%nt_ovlp)
call DestructMatrix(ph%nt_ovlp_copy)
call DestructMatrix(ph%nt_dm)
call DestructPermutation(ph%nt_perm)
call DestructProcessGrid(ph%nt_pgrid)
......
......@@ -977,9 +977,9 @@ subroutine elsi_write_mat_real_mp(rwh,f_name,mat)
integer(kind=i4) :: n_lcol0
integer(kind=i4) :: prev_nnz
integer(kind=i8) :: offset
real(kind=r8) :: dummy1(1,1)
real(kind=r8) :: dummy2(1)
real(kind=r8), allocatable :: dummy1(:,:)
real(kind=r8), allocatable :: dummy2(:)
real(kind=r8), allocatable :: mat_csc(:)
integer(kind=i4), allocatable :: row_ind(:)
integer(kind=i4), allocatable :: col_ptr(:)
......@@ -1001,21 +1001,15 @@ subroutine elsi_write_mat_real_mp(rwh,f_name,mat)
eh%bh%n_lcol_sp = rwh%n_basis-(rwh%bh%n_procs-1)*eh%bh%n_lcol_sp
end if
call elsi_allocate(rwh%bh,dummy1,1,1,"dummy1",caller)
call elsi_blacs_to_sips_hs_dim(eh%ph,eh%bh,mat,dummy1)
call elsi_allocate(rwh%bh,mat_csc,eh%bh%nnz_l_sp,"mat_csc",caller)
call elsi_allocate(rwh%bh,dummy2,1,"dummy2",caller)
call elsi_allocate(rwh%bh,row_ind,eh%bh%nnz_l_sp,"row_ind",caller)
call elsi_allocate(rwh%bh,col_ptr,eh%bh%n_lcol_sp+1,"col_ptr",caller)
call elsi_blacs_to_sips_hs(eh%ph,eh%bh,mat,dummy1,mat_csc,dummy2,row_ind,&
col_ptr)
call elsi_deallocate(rwh%bh,dummy1,"dummy1")
call elsi_deallocate(rwh%bh,dummy2,"dummy2")
! Open file
f_mode = mpi_mode_wronly+mpi_mode_create
......@@ -1109,9 +1103,9 @@ subroutine elsi_write_mat_complex_mp(rwh,f_name,mat)
integer(kind=i4) :: n_lcol0
integer(kind=i4) :: prev_nnz
integer(kind=i8) :: offset
complex(kind=r8) :: dummy1(1,1)
complex(kind=r8) :: dummy2(1)
complex(kind=r8), allocatable :: dummy1(:,:)
complex(kind=r8), allocatable :: dummy2(:)
complex(kind=r8), allocatable :: mat_csc(:)
integer(kind=i4), allocatable :: row_ind(:)
integer(kind=i4), allocatable :: col_ptr(:)
......@@ -1133,21 +1127,15 @@ subroutine elsi_write_mat_complex_mp(rwh,f_name,mat)
eh%bh%n_lcol_sp = rwh%n_basis-(rwh%bh%n_procs-1)*eh%bh%n_lcol_sp
end if
call elsi_allocate(rwh%bh,dummy1,1,1,"dummy1",caller)
call elsi_blacs_to_sips_hs_dim(eh%ph,eh%bh,mat,dummy1)
call elsi_allocate(rwh%bh,mat_csc,eh%bh%nnz_l_sp,"mat_csc",caller)
call elsi_allocate(rwh%bh,dummy2,1,"dummy2",caller)
call elsi_allocate(rwh%bh,row_ind,eh%bh%nnz_l_sp,"row_ind",caller)
call elsi_allocate(rwh%bh,col_ptr,eh%bh%n_lcol_sp+1,"col_ptr",caller)
call elsi_blacs_to_sips_hs(eh%ph,eh%bh,mat,dummy1,mat_csc,dummy2,row_ind,&
col_ptr)
call elsi_deallocate(rwh%bh,dummy1,"dummy1")
call elsi_deallocate(rwh%bh,dummy2,"dummy2")
! Open file
f_mode = mpi_mode_wronly+mpi_mode_create
......
......@@ -9,8 +9,8 @@
!!
module ELSI_SETUP
use ELSI_CONSTANTS, only: PEXSI_SOLVER,SINGLE_PROC,MULTI_PROC,PEXSI_CSC,&
SIESTA_CSC,UNSET
use ELSI_CONSTANTS, only: ELPA_SOLVER,PEXSI_SOLVER,SINGLE_PROC,MULTI_PROC,&
PEXSI_CSC,SIESTA_CSC,UNSET
use ELSI_DATATYPE, only: elsi_handle,elsi_param_t,elsi_basic_t
use ELSI_ELPA, only: elsi_cleanup_elpa
use ELSI_IO, only: elsi_say,elsi_final_print,fjson_close_file,&
......@@ -349,7 +349,6 @@ subroutine elsi_reinit(eh)
eh%ph%sips_first = .true.
eh%ph%nt_first = .true.
call elsi_cleanup_ntpoly(eh%ph)
call elsi_cleanup_pexsi(eh%ph)
call elsi_cleanup_sips(eh%ph)
......@@ -400,6 +399,16 @@ subroutine elsi_reinit(eh)
if(allocated(eh%col_ptr_sp2)) then
call elsi_deallocate(eh%bh,eh%col_ptr_sp2,"col_ptr_sp2")
end if
if(.not. eh%ph%solver == ELPA_SOLVER) then
if(allocated(eh%ovlp_real_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_real_copy,"ovlp_real_copy")
end if
if(allocated(eh%ovlp_cmplx_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_cmplx_copy,"ovlp_cmplx_copy")
end if
end if
end if
end subroutine
......
This diff is collapsed.
......@@ -10,11 +10,11 @@
!!
module ELSI_TOOLS
use ELSI_CONSTANTS, only: ELPA_SOLVER,OMM_SOLVER,SIPS_SOLVER
use ELSI_CONSTANTS, only: ELPA_SOLVER,OMM_SOLVER,SIPS_SOLVER,NTPOLY_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_MALLOC, only: elsi_allocate
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_update_dm_ntpoly,CopyMatrix
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
......@@ -79,120 +79,89 @@ end subroutine
!>
!! This routine extrapolates density matrix for a new overlap.
!!
subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
subroutine elsi_extrapolate_dm_real(eh,ovlp,dm)
implicit none
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) :: ovlp(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)
solver_save = eh%ph%solver
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)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,&
eh%ovlp_real_copy,ovlp,dm)
eh%ovlp_real_copy = ovlp
case default
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)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp,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,&
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp,dm,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_update_dm_ntpoly(eh%ph,eh%bh,eh%ph%nt_ovlp_copy,eh%ph%nt_ovlp,&
eh%ph%nt_ham,eh%ph%nt_dm)
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_dm,dm)
call CopyMatrix(eh%ph%nt_ovlp,eh%ph%nt_ovlp_copy)
end select
eh%ph%solver = solver_save
end subroutine
!>
!! This routine extrapolates density matrix for a new overlap.
!!
subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
subroutine elsi_extrapolate_dm_complex(eh,ovlp,dm)
implicit none
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) :: ovlp(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)
solver_save = eh%ph%solver
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)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,&
eh%ovlp_cmplx_copy,ovlp,dm)
eh%ovlp_cmplx_copy = ovlp
case default
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)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp,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,&
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp,dm,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_update_dm_ntpoly(eh%ph,eh%bh,eh%ph%nt_ovlp_copy,eh%ph%nt_ovlp,&