Commit a023d9f4 authored by Irina Lebedeva's avatar Irina Lebedeva
Browse files

Update omm_callback

parent 33f6d510
......@@ -67,6 +67,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
!**** INOUT ***********************************!
type(matrix), intent(inout) :: S ! overlap matrix (or its Cholesky-factorized upper
! triangular matrix) in orbital basis (m x m)
type(matrix), intent(inout) :: D_min ! density (or energy-weighted density) matrix
......@@ -124,13 +125,11 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
type(matrix), save :: HC ! H*C matrix (n x m)
type(matrix), save :: HG ! H*G matrix (n x m)
type(matrix), save :: work1 ! work matrix (n x n)
type(matrix), allocatable, save :: HW(:) ! Hamiltonian matrix in WF basis (n x n)
! for each value of ip
type(matrix), save :: HW ! Hamiltonian matrix in WF basis (n x n)
type(matrix), save :: HWd ! G^dag*H*C matrix (n x n)
type(matrix), save :: HWdd ! G^dag*H*G matrix (n x n)
type(matrix), allocatable, save :: SW(:) ! overlap matrix in WF basis (n x n)
! for each value of ip
type(matrix), allocatable, save :: SC(:) ! S*C matrix (n x m) for each value of ip
type(matrix), save :: SW ! overlap matrix in WF basis (n x n)
type(matrix), save :: SC ! S*C matrix (n x m)
type(matrix), save :: SG ! S*G matrix (n x m)
type(matrix), allocatable, save :: C_Chl(:) ! Cholesky-transformed WF coeffs.
! in orbital basis (n x m) for each value of ip
......@@ -174,6 +173,10 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
use_Cholesky = .false.
S_is_U = .false.
use_precon = .true.
! else if (flavour == 4) then
! use_Cholesky = .false.
! S_is_U = .false.
......@@ -222,7 +225,6 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (dealloc) then
if (work1%is_initialized) call m_deallocate(work1)
if (HG%is_initialized) call m_deallocate(HG)
if (HC%is_initialized) call m_deallocate(HC)
if (D%is_initialized) call m_deallocate(D)
if (PG_p%is_initialized) call m_deallocate(PG_p)
if (PG%is_initialized) call m_deallocate(PG)
......@@ -245,22 +247,59 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end do
deallocate(C_Chl)
end if
if (allocated(SC)) then
do i = 1, np
if (SC(i)%is_initialized) call m_deallocate(SC(i))
end do
deallocate(SC)
end if
if(present_eta) then
if (.not. SW%is_initialized) call m_allocate(SW, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
call m_set(SW, 'a', 0.0_dp, 0.0_dp, m_operation)
if (.not. (use_Cholesky .or. no_S)) then
if (.not. SC%is_initialized) call m_allocate(SC, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
call m_set(SC, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
! -calculate the overlap matrix in WF basis: SW=C^dag*S*C
if (use_Cholesky) then
call mm_multiply(C_Chl(ip), 'n', C_Chl(ip), 'c', SW, 1.0_dp, 0.0_dp, m_operation)
else if (no_S) then
call mm_multiply(C_min, 'n', C_min, 'c', SW, 1.0_dp, 0.0_dp, m_operation)
else
call calc_AW(S, C_min, SW, SC, m_operation)
end if
if (.not. QW%is_initialized) call m_allocate(QW, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
call m_set(QW, 'a', 0.0_dp, 1.0_dp, m_operation)
call m_add(SW, 'n', QW, -1.0_dp, 1.0_dp, m_operation)
end if
if (.not. work2%is_initialized) call m_allocate(work2, n, m, label=m_storage,&
if(dealloc) then
if (SC%is_initialized) call m_deallocate(SC)
if (SW%is_initialized) call m_deallocate(SW)
end if
if (.not. HW%is_initialized) call m_allocate(HW, n, n, label=m_storage,&
blocksize1 = bl_n,blocksize2 = bl_n)
if (.not. HC%is_initialized) call m_allocate(HC, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
if (.not. QW%is_initialized) call m_allocate(QW, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
call m_set(HW, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HC, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(QW, 'a', 0.0_dp, 2.0_dp, m_operation)
call m_add(SW(ip), 'n', QW, -1.0_dp, 1.0_dp, m_operation)
if(present_eta) call m_add(QW, 'n', HW(ip), eta, 1.0_dp, m_operation)
! -calculate the hamiltonian in WF basis: HW=C^dag*H*C for eta=0
if (use_Cholesky) then
call calc_AW(H, C_Chl(ip), HW, HC, m_operation)
else
call calc_AW(H, C_min, HW, HC, m_operation)
end if
if(dealloc) then
if (HC%is_initialized) call m_deallocate(HC)
end if
if(present_eta) call m_add(QW, 'n', HW, 2.0_dp * eta, 1.0_dp, m_operation)
if(dealloc) then
if (QW%is_initialized) call m_deallocate(QW)
end if
! impose sparsity on D_min
if(sparse) then
......@@ -268,27 +307,14 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call m_add(S, 'n', D_min, 1.0_dp, 0.0_dp)
end if
! calculate the energy-density matrix: ED=C*[(HW+eta*(2IW-SW))]*C^dag
call calc_A(HW(ip), C_min, D_min, work2, m_operation, keep_sparsity=sparse)
! shift the Hamiltonian back
if((.not. dealloc) .and. present_eta) call m_add(QW, 'n', HW(ip), -eta, 1.0_dp, m_operation)
if (.not. work2%is_initialized) call m_allocate(work2, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
! calculate the energy-density matrix: ED=C*[(HW(eta=0)+2*eta*(IW-SW))]*C^dag
call calc_A(HW, C_min, D_min, work2, m_operation, keep_sparsity=sparse)
if (dealloc .and. (ip == np)) then
if (dealloc) then
if (work2%is_initialized) call m_deallocate(work2)
if (allocated(HW)) then
do i = 1, np
if (HW(i)%is_initialized) call m_deallocate(HW(i))
end do
deallocate(HW)
end if
if (allocated(SW)) then
do i = 1, np
if (SW(i)%is_initialized) call m_deallocate(SW(i))
end do
deallocate(SW)
end if
if (QW%is_initialized) call m_deallocate(QW)
if (HW%is_initialized) call m_deallocate(HW)
end if
if(mpi_rank == 0) close(log_unit)
......@@ -303,11 +329,6 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end if
! allocate matrices involved
if (.not. allocated(HW)) allocate(HW(np))
if (.not. allocated(SW)) allocate(SW(np))
if (.not. (use_Cholesky .or. no_S)) then
if (.not. allocated(SC)) allocate(SC(np))
end if
if (use_Cholesky) then
if (.not. allocated(C_Chl)) allocate(C_Chl(np))
end if
......@@ -315,16 +336,16 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (.not. allocated(P)) allocate(P(np))
end if
if (.not. HW(ip)%is_initialized) call m_allocate(HW(ip), n, n, label=m_storage,&
if (.not. HW%is_initialized) call m_allocate(HW, n, n, label=m_storage,&
blocksize1 = bl_n,blocksize2 = bl_n)
if (.not. HWd%is_initialized) call m_allocate(HWd, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
if (.not. HWdd%is_initialized) call m_allocate(HWdd, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
if (.not. SW(ip)%is_initialized) call m_allocate(SW(ip), n, n, label=m_storage,&
if (.not. SW%is_initialized) call m_allocate(SW, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
if (.not. (use_Cholesky .or. no_S)) then
if (.not. SC(ip)%is_initialized) call m_allocate(SC(ip), n, m, label=m_storage,&
if (.not. SC%is_initialized) call m_allocate(SC, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
if (.not. SG%is_initialized) call m_allocate(SG, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
......@@ -360,7 +381,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (.not. work1%is_initialized) call m_allocate(work1, n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
call m_set(HW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HW, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HWd, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HWdd, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(work2, 'a', 0.0_dp, 0.0_dp, m_operation)
......@@ -371,18 +392,15 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call m_set(HC, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HG, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(work1, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(SW, 'a', 0.0_dp, 0.0_dp, m_operation)
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SG, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
if(new_S .or. first_call(ip)) then
call m_set(SW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SC(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
end if
call m_set(SC, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
if(sparse) then
if((.not. G%is_initialized) .or. (new_S .and. ip == 1)) then
! Creating a gradient matrix of proper sparsity
if (.not. G%is_initialized) call m_allocate(G, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
call m_set(G, 'a', 0.0_dp, 0.0_dp, m_operation)
......@@ -394,6 +412,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call m_set(G, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
! shift H for Kim's functional
! if(present_eta) then
! call m_add(S, 'n', H, -eta, 1.0_dp, m_operation)
......@@ -452,30 +471,28 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! first we calculate the energy and gradient for our initial guess, with the following steps:
! -calculate the hamiltonian in WF basis: HW=C^dag*H*C
if (use_Cholesky) then
call calc_HW_callback(H, C_Chl(ip), HW(ip), HC, m_operation)
call calc_HW_callback(H, C_Chl(ip), HW, HC, m_operation)
else
call calc_HW_callback(H, C_min, HW(ip), HC, m_operation)
call calc_HW_callback(H, C_min, HW, HC, m_operation)
end if
! -calculate the overlap matrix in WF basis: SW=C^dag*S*C
if (use_Cholesky) then
call mm_multiply(C_Chl(ip), 'n',&
C_Chl(ip), 'c', SW(ip), 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(C_Chl(ip), 'n', C_Chl(ip), 'c', SW, 1.0_dp, 0.0_dp, m_operation)
else if (no_S) then
call mm_multiply(C_min, 'n', C_min, 'c',&
SW(ip), 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(C_min, 'n', C_min, 'c', SW, 1.0_dp, 0.0_dp, m_operation)
else
call calc_AW(S, C_min, SW(ip), SC(ip), m_operation)
call calc_AW(S, C_min, SW, SC, m_operation)
end if
! calculate the gradient: G=2*(2*H*C-S*C*HW-H*C*SW)
! (note that we *reuse* H*C and S*C contained in HC and SC from the previous calls to calc_AW)
if (use_Cholesky) then
call calc_G(HW(ip), SW(ip), G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
call calc_G(HW, SW, G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
call calc_G(HW, SW, G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call calc_G(HW(ip), SW(ip), G, HC, SC(ip), m_operation, keep_sparsity=sparse)
call calc_G(HW, SW, G, HC, SC, m_operation, keep_sparsity=sparse)
end if
! calculate the preconditioned gradient by premultiplying G by P
......@@ -490,7 +507,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! (the energy at C is given by the zeroth-order coeff. alpha(0))
if (use_precon) then
call mm_multiply(HC, 'n', PG, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC, 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(HC, 'n', G, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
if (use_Cholesky) then
......@@ -498,7 +515,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else if (no_S) then
call mm_multiply(C_min, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC(ip), 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
end if
......@@ -515,7 +532,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end if
! calculate coefficients of the quartic energy function
call calc_coeff(HW(ip), SW(ip), HWd, SWd, HWdd, SWdd, work1, coeff, m_operation)
call calc_coeff(HW, SW, HWd, SWd, HWdd, SWdd, work1, coeff, m_operation)
e_min = coeff(0)
! if(present_eta) e_min = e_min + n_occ_i * eta
......@@ -559,7 +576,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else if (no_S) then
call mm_multiply(C_min, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC(ip), 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
call calc_HW_callback(H, D, HWdd, HG, m_operation)
if (use_Cholesky .or. no_S) then
......@@ -567,7 +584,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else
call calc_AW(S, D, SWdd, SG, m_operation)
end if
call calc_coeff(HW(ip), SW(ip), HWd, SWd, HWdd, SWdd, work1, coeff, m_operation)
call calc_coeff(HW, SW, HWd, SWd, HWdd, SWdd, work1, coeff, m_operation)
end if
! using the coeffs. calculated anlytically, we can find the minimum of the functional in the
! search direction, and calculate the energy at that minimum
......@@ -596,19 +613,19 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end if
! recalculate SW at the minimum (or for the rescaled coeffs.)
call m_add(SWd, 'c', SW(ip), x_min, 1.0_dp, m_operation)
call m_add(SWd, 'n', SW(ip), x_min, 1.0_dp, m_operation)
call m_add(SWdd, 'n', SW(ip), x_min**2, 1.0_dp, m_operation)
call m_add(SWd, 'c', SW, x_min, 1.0_dp, m_operation)
call m_add(SWd, 'n', SW, x_min, 1.0_dp, m_operation)
call m_add(SWdd, 'n', SW, x_min**2, 1.0_dp, m_operation)
! recalculate HW at the minimum (or for the rescaled coeffs.)
call m_add(HWd, 'c', HW(ip), x_min, 1.0_dp, m_operation)
call m_add(HWd, 'n', HW(ip), x_min, 1.0_dp, m_operation)
call m_add(HWdd, 'n', HW(ip), x_min**2, 1.0_dp, m_operation)
call m_add(HWd, 'c', HW, x_min, 1.0_dp, m_operation)
call m_add(HWd, 'n', HW, x_min, 1.0_dp, m_operation)
call m_add(HWdd, 'n', HW, x_min**2, 1.0_dp, m_operation)
! recalculate HC ans SC
call m_add(HG, 'n', HC, x_min, 1.0_dp, m_operation)
if(.not. (use_Cholesky .or. no_S)) then
call m_add(SG, 'n', SC(ip), x_min, 1.0_dp, m_operation)
call m_add(SG, 'n', SC, x_min, 1.0_dp, m_operation)
end if
e_diff = 2.0_dp * abs((e_min - e_min_old)/(e_min + e_min_old))
......@@ -623,14 +640,15 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! recalculate G at the minimum (or for the rescaled coeffs.)
if (use_Cholesky) then
call calc_G(HW(ip), SW(ip), G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
call calc_G(HW, SW, G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
call calc_G(HW, SW, G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call m_add(SG, 'n', SC(ip), x_min, 1.0_dp, m_operation)
call calc_G(HW, SW, G, HC, SC, m_operation, keep_sparsity=sparse)
end if
if (use_precon) call mm_multiply(G, 'n', P(ip), 'n', PG, 1.0_dp, 0.0_dp, m_operation)
if (ls_conv) then
call m_add(G, 'n', work2, 1.0_dp, 0.0_dp, m_operation)
call m_add(G_p, 'n', work2, -1.0_dp, 1.0_dp, m_operation)
......@@ -664,10 +682,12 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (G_p%is_initialized) call m_deallocate(G_p)
if (G%is_initialized) call m_deallocate(G)
if (SG%is_initialized) call m_deallocate(SG)
if (SC%is_initialized) call m_deallocate(SC)
if (SWdd%is_initialized) call m_deallocate(SWdd)
if (SWd%is_initialized) call m_deallocate(SWd)
if (HWdd%is_initialized) call m_deallocate(HWdd)
if (HWd%is_initialized) call m_deallocate(HWd)
if (HW%is_initialized) call m_deallocate(HW)
if (work1%is_initialized) call m_deallocate(work1)
end if
......@@ -677,7 +697,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! calculate the density matrix: D=C*(2*IW-SW)*C^dag
call m_set(QW, 'a', 0.0_dp, 2.0_dp, m_operation)
call m_add(SW(ip), 'n', QW, -1.0_dp, 1.0_dp, m_operation)
call m_add(SW, 'n', QW, -1.0_dp, 1.0_dp, m_operation)
if (use_Cholesky) then
call m_add(C_Chl(ip), 'n', C_min, 1.0_dp, 0.0_dp, m_operation)
call m_back_transform(S, C_min, m_operation)
......@@ -692,7 +712,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! calculate the trace of (2*IW-SW) to make sure we are occupying the right number
! of eigenstates in our solution
call mm_trace(SW(ip), QW, TrQS, m_operation)
call mm_trace(SW, QW, TrQS, m_operation)
if (mpi_rank == 0) then
write(log_unit,'(a,i5,a)') '| minim: icg = ', icg, ' |'
......@@ -702,11 +722,12 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
write(log_unit,'(a)') '+-----------------------------------------------+'
end if
if (dealloc) then
first_call(ip) = .false.
if(dealloc) then
if (QW%is_initialized) call m_deallocate(QW)
if (work2%is_initialized) call m_deallocate(work2)
if (SW%is_initialized) call m_deallocate(SW)
end if
first_call(ip) = .false.
if(mpi_rank == 0) close(log_unit)
......
Supports Markdown
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