Commit 1dc51216 authored by Victor Yu's avatar Victor Yu

Update DM extrapolation with NTPoly

parent 5c964d97
MODULE NTPoly
USE DensityMatrixSolversModule
USE GeometryOptimizationModule
USE LoggingModule
USE PermutationModule
USE PMatrixMemoryPoolModule
......
......@@ -9,19 +9,21 @@
!!
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,ProcessGrid_t,&
ConstructNewProcessGrid,DestructProcessGrid,ConstructRandomPermutation,&
DestructPermutation,InverseSquareRoot,SolverParameters_t,Triplet_r,&
Triplet_c,TripletList_r,TripletList_c,ConstructTripletList,&
AppendToTripletList,DestructTripletList,ActivateLogger,DeactivateLogger
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,&
ConstructTripletList,AppendToTripletList,DestructTripletList,&
ActivateLogger,DeactivateLogger
implicit none
......@@ -31,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
......@@ -208,6 +213,47 @@ subroutine elsi_compute_edm_ntpoly(ph,bh,ham,edm)
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 density matrix extrapolation"
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 cleans up NTPoly.
!!
......
......@@ -10,11 +10,12 @@
!!
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
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_solve_ntpoly
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_update_dm_ntpoly,Matrix_ps,&
CopyMatrix,DestructMatrix
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
......@@ -92,6 +93,8 @@ subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
real(kind=r8) :: dummy(1,1)
logical :: unit_ovlp_save
type(Matrix_ps) :: nt_ovlp0
character(len=*), parameter :: caller = "elsi_extrapolate_dm_real"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -130,14 +133,21 @@ subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
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)
if(eh%ph%solver == NTPOLY_SOLVER) then
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp0,dummy,nt_ovlp0,&
eh%ph%nt_dm)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,dummy,eh%ph%nt_ham,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
else
call elsi_init_ntpoly(eh%ph,eh%bh)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp0,eh%ph%nt_dm,nt_ovlp0)
end if
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
......@@ -145,8 +155,14 @@ subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
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)
call elsi_update_dm_ntpoly(eh%ph,eh%bh,nt_ovlp0,eh%ph%nt_ovlp,&
eh%ph%nt_dm,eh%ph%nt_ham)
if(eh%ph%solver == NTPOLY_SOLVER) then
call DestructMatrix(nt_ovlp0)
end if
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_ham,dm)
end select
eh%ph%solver = solver_save
......@@ -169,6 +185,8 @@ subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
complex(kind=r8) :: dummy(1,1)
logical :: unit_ovlp_save
type(Matrix_ps) :: nt_ovlp0
character(len=*), parameter :: caller = "elsi_extrapolate_dm_complex"
call elsi_check_init(eh%bh,eh%handle_init,caller)
......@@ -194,14 +212,21 @@ subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
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)
if(eh%ph%solver == NTPOLY_SOLVER) then
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp0,dummy,nt_ovlp0,&
eh%ph%nt_dm)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,dummy,eh%ph%nt_ham,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
else
call elsi_init_ntpoly(eh%ph,eh%bh)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp0,eh%ph%nt_dm,nt_ovlp0)
end if
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
......@@ -209,8 +234,14 @@ subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
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)
call elsi_update_dm_ntpoly(eh%ph,eh%bh,nt_ovlp0,eh%ph%nt_ovlp,&
eh%ph%nt_dm,eh%ph%nt_ham)
if(eh%ph%solver == NTPOLY_SOLVER) then
call DestructMatrix(nt_ovlp0)
end if
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_ham,dm)
end select
eh%ph%solver = solver_save
......
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