Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[feature] restore function to set density matrix #884

Merged
merged 1 commit into from
Aug 21, 2023
Merged

[feature] restore function to set density matrix #884

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 58 additions & 9 deletions src/api/sirius.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5938,33 +5938,33 @@ end subroutine sirius_save_state

!
!> @brief Save DFT ground state (density and potential)
!> @param [in] gs_handler Ground-state handler.
!> @param [in] handler Ground-state handler.
!> @param [in] file_name Name of the file that stores the saved data.
!> @param [out] error_code Error code
subroutine sirius_load_state(gs_handler,file_name,error_code)
subroutine sirius_load_state(handler,file_name,error_code)
implicit none
!
type(sirius_ground_state_handler), target, intent(in) :: gs_handler
type(sirius_ground_state_handler), target, intent(in) :: handler
character(*), target, intent(in) :: file_name
integer, optional, target, intent(out) :: error_code
!
type(C_PTR) :: gs_handler_ptr
type(C_PTR) :: handler_ptr
type(C_PTR) :: file_name_ptr
character(C_CHAR), target, allocatable :: file_name_c_type(:)
type(C_PTR) :: error_code_ptr
!
interface
subroutine sirius_load_state_aux(gs_handler,file_name,error_code)&
subroutine sirius_load_state_aux(handler,file_name,error_code)&
&bind(C, name="sirius_load_state")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: gs_handler
type(C_PTR), value :: handler
type(C_PTR), value :: file_name
type(C_PTR), value :: error_code
end subroutine
end interface
!
gs_handler_ptr = C_NULL_PTR
gs_handler_ptr = C_LOC(gs_handler%handler_ptr_)
handler_ptr = C_NULL_PTR
handler_ptr = C_LOC(handler%handler_ptr_)
file_name_ptr = C_NULL_PTR
allocate(file_name_c_type(len(file_name)+1))
file_name_c_type = string_f2c(file_name)
Expand All @@ -5973,10 +5973,59 @@ subroutine sirius_load_state_aux(gs_handler,file_name,error_code)&
if (present(error_code)) then
error_code_ptr = C_LOC(error_code)
endif
call sirius_load_state_aux(gs_handler_ptr,file_name_ptr,error_code_ptr)
call sirius_load_state_aux(handler_ptr,file_name_ptr,error_code_ptr)
deallocate(file_name_c_type)
end subroutine sirius_load_state

!
!> @brief Set density matrix.
!> @param [in] handler Ground-state handler.
!> @param [in] ia Index of atom.
!> @param [in] dm Input density matrix.
!> @param [in] ld Leading dimension of the density matrix.
!> @param [out] error_code Error code.
subroutine sirius_set_density_matrix(handler,ia,dm,ld,error_code)
implicit none
!
type(sirius_ground_state_handler), target, intent(in) :: handler
integer, target, intent(in) :: ia
complex(8), target, intent(in) :: dm
integer, target, intent(in) :: ld
integer, optional, target, intent(out) :: error_code
!
type(C_PTR) :: handler_ptr
type(C_PTR) :: ia_ptr
type(C_PTR) :: dm_ptr
type(C_PTR) :: ld_ptr
type(C_PTR) :: error_code_ptr
!
interface
subroutine sirius_set_density_matrix_aux(handler,ia,dm,ld,error_code)&
&bind(C, name="sirius_set_density_matrix")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: handler
type(C_PTR), value :: ia
type(C_PTR), value :: dm
type(C_PTR), value :: ld
type(C_PTR), value :: error_code
end subroutine
end interface
!
handler_ptr = C_NULL_PTR
handler_ptr = C_LOC(handler%handler_ptr_)
ia_ptr = C_NULL_PTR
ia_ptr = C_LOC(ia)
dm_ptr = C_NULL_PTR
dm_ptr = C_LOC(dm)
ld_ptr = C_NULL_PTR
ld_ptr = C_LOC(ld)
error_code_ptr = C_NULL_PTR
if (present(error_code)) then
error_code_ptr = C_LOC(error_code)
endif
call sirius_set_density_matrix_aux(handler_ptr,ia_ptr,dm_ptr,ld_ptr,error_code_ptr)
end subroutine sirius_set_density_matrix


subroutine sirius_free_handler_ctx(handler, error_code)
implicit none
Expand Down
56 changes: 55 additions & 1 deletion src/api/sirius_api.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6071,7 +6071,7 @@ sirius_save_state(void** handler__, const char* file_name__, int* error_code__)
sirius_load_state:
doc: Save DFT ground state (density and potential)
arguments:
gs_handler:
handler:
type: gs_handler
attr: in, required
doc: Ground-state handler.
Expand All @@ -6098,4 +6098,58 @@ sirius_load_state(void** handler__, const char* file_name__, int* error_code__)
error_code__);
}

/*
@api begin
sirius_set_density_matrix:
doc: Set density matrix.
arguments:
handler:
type: gs_handler
attr: in, required
doc: Ground-state handler.
ia:
type: int
attr: in, required
doc: Index of atom.
dm:
type: complex
attr: in, required
doc: Input density matrix.
ld:
type: int
attr: in, required
doc: Leading dimension of the density matrix.
error_code:
type: int
attr: out, optional
doc: Error code.
@api end
*/
void
sirius_set_density_matrix(void** handler__, int const* ia__, std::complex<double>* dm__, int const* ld__, int* error_code__)
{
call_sirius(
[&]() {
auto& gs = get_gs(handler__);
sddk::mdarray<std::complex<double>, 3> dm(dm__, *ld__, *ld__, 3);
int ia = *ia__ - 1;
auto& atom = gs.ctx().unit_cell().atom(ia);
auto idx_map = atomic_orbital_index_map_QE(atom.type());
int nbf = atom.mt_basis_size();
RTE_ASSERT(nbf <= *ld__);

for (int icomp = 0; icomp < gs.ctx().num_mag_comp(); icomp++) {
for (int xi1 = 0; xi1 < nbf; xi1++) {
int p1 = phase_Rlm_QE(atom.type(), xi1);
for (int xi2 = 0; xi2 < nbf; xi2++) {
int p2 = phase_Rlm_QE(atom.type(), xi2);
gs.density().density_matrix(ia)(xi1, xi2, icomp) = dm(idx_map[xi1], idx_map[xi2], icomp) *
static_cast<double>(p1 * p2);
}
}
}
},
error_code__);
}

} // extern "C"
Loading