diff --git a/.gitignore b/.gitignore index da51fad..d8acdc7 100644 --- a/.gitignore +++ b/.gitignore @@ -110,6 +110,9 @@ auto-save-list tramp .\#* +### VS Code ### +.vscode + # Org-mode .org-id-locations *_archive diff --git a/src/admm_dm_methods.cpp b/src/admm_dm_methods.cpp new file mode 100644 index 0000000..c90f367 --- /dev/null +++ b/src/admm_dm_methods.cpp @@ -0,0 +1,533 @@ +/*-------------------------------------------------------------------------------------------------- + CP2K: A general program to perform molecular dynamics simulations + Copyright 2000-2023 CP2K developers group + + SPDX-License-Identifier: GPL-2.0-or-later +--------------------------------------------------------------------------------------------------*/ + +/*************************************************************************************************** + \brief Contains ADMM methods which only require the density matrix + \par History + 11.2014 created [Ole Schuett] + \author Ole Schuett +***************************************************************************************************/ + +#include "./base/base_uses.h" + +class admm_dm_methods { + +public: + void admm_dm_calc_rho_aux(struct qs_environment_type); + void admm_dm_merge_ks_matrix(struct qs_environment_type); + +private: + char *moduleN = "admm_dm_methods"; + +}; + +/*************************************************************************************************** + \brief Entry methods: Calculates auxiliary density matrix from primary one. + \param qs_env ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::admm_dm_calc_rho_aux(struct qs_environment_type *qs_env) +{ + + char *routineN = "admm_dm_calc_rho_aux"; + + int handle; + struct admm_dm_type *admm_dm; + + NULLIFY (admm_dm); + timeset(routineN, handle); + get_admm_env(qs_env.admm_env, admm_dm=admm_dm); + + SELECT CASE (admm_dm.method) + CASE (do_admm_basis_projection) + CALL map_dm_projection(qs_env) + + CASE (do_admm_blocked_projection) + CALL map_dm_blocked(qs_env) + + CASE DEFAULT + CPABORT("admm_dm_calc_rho_aux: unknown method") + END SELECT + + if (admm_dm.purify) purify_mcweeny(qs_env); + + update_rho_aux(qs_env); + + timestop(handle); +} + +/*************************************************************************************************** + \brief Entry methods: Merges auxiliary Kohn-Sham matrix into primary one. + \param qs_env ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::admm_dm_merge_ks_matrix(struct qs_environment_type *qs_env) +{ + char *routineN = "admm_dm_merge_ks_matrix"; + + int handle; + struct admm_dm_type *admm_dm; + struct dbcsr_p_type *matrix_ks_merge; + + timeset(routineN, handle); + NULLIFY (admm_dm, matrix_ks_merge); + + get_admm_env(qs_env.dmm_env, admm_dm=admm_dm); + + IF (admm_dm.purify) THEN + CALL revert_purify_mcweeny(qs_env, matrix_ks_merge) + ELSE + CALL get_admm_env(qs_env.admm_env, matrix_ks_aux_fit=matrix_ks_merge) + END IF + + SELECT CASE (admm_dm.method) + CASE (do_admm_basis_projection) + CALL merge_dm_projection(qs_env, matrix_ks_merge) + + CASE (do_admm_blocked_projection) + CALL merge_dm_blocked(qs_env, matrix_ks_merge) + + CASE DEFAULT + CPABORT("admm_dm_merge_ks_matrix: unknown method") + END SELECT + + IF (admm_dm.purify) & + CALL dbcsr_deallocate_matrix_set(matrix_ks_merge) + + CALL timestop(handle) + +} + +/*************************************************************************************************** + \brief Calculates auxiliary density matrix via basis projection. + \param qs_env ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::map_dm_projection(qs_env) +{ + TYPE(qs_environment_type), POINTER :: qs_env + + INTEGER :: ispin + LOGICAL :: s_mstruct_changed + REAL(KIND=dp) :: threshold + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux, matrix_s_mixed, rho_ao, & + rho_ao_aux + TYPE(dbcsr_type) :: matrix_s_aux_inv, matrix_tmp + TYPE(dft_control_type), POINTER :: dft_control + TYPE(qs_rho_type), POINTER :: rho, rho_aux + + NULLIFY (dft_control, admm_dm, matrix_s_aux, matrix_s_mixed, rho, rho_aux) + NULLIFY (rho_ao, rho_ao_aux) + + CALL get_qs_env(qs_env, dft_control=dft_control, s_mstruct_changed=s_mstruct_changed, +rho=rho) + CALL get_admm_env(qs_env.admm_env, matrix_s_aux_fit=matrix_s_aux, rho_aux_fit=rho_aux, & + matrix_s_aux_fit_vs_orb=matrix_s_mixed, admm_dm=admm_dm) + + CALL qs_rho_get(rho, rho_ao=rho_ao) + CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux) + + IF (s_mstruct_changed) THEN + ! Calculate A = S_aux^(-1) * S_mixed + CALL dbcsr_create(matrix_s_aux_inv, template=matrix_s_aux[0].matrix, matrix_type="N") + threshold = MAX(admm_dm.eps_filter, 1.0e-12_dp) + CALL invert_Hotelling(matrix_s_aux_inv, matrix_s_aux[0].matrix, threshold) + + IF (.NOT. ASSOCIATED(admm_dm.matrix_A)) THEN + ALLOCATE (admm_dm.matrix_A) + CALL dbcsr_create(admm_dm.matrix_A, template=matrix_s_mixed[0].matrix, matrix_type="N") + END IF + CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_aux_inv, matrix_s_mixed[0].matrix, & + 0.0_dp, admm_dm.matrix_A) + CALL dbcsr_release(matrix_s_aux_inv) + END IF + + ! Calculate P_aux = A * P * A^T + CALL dbcsr_create(matrix_tmp, template=admm_dm.matrix_A) + DO ispin = 1, dft_control.nspins + CALL dbcsr_multiply("N", "N", 1.0_dp, admm_dm.matrix_A, rho_ao(ispin).matrix, & + 0.0_dp, matrix_tmp) + CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, admm_dm.matrix_A, & + 0.0_dp, rho_ao_aux(ispin).matrix) + END DO + CALL dbcsr_release(matrix_tmp); + +} + +/*************************************************************************************************** + \brief Calculates auxiliary density matrix via blocking. + \param qs_env ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::map_dm_blocked(qs_env) +{ + TYPE(qs_environment_type), POINTER :: qs_env + + INTEGER :: blk, iatom, ispin, jatom + LOGICAL :: found + REAL(dp), DIMENSION(:, :), POINTER :: sparse_block, sparse_block_aux + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_iterator_type) :: iter + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_aux + TYPE(dft_control_type), POINTER :: dft_control + TYPE(qs_rho_type), POINTER :: rho, rho_aux + + NULLIFY (dft_control, admm_dm, rho, rho_aux, rho_ao, rho_ao_aux) + + CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho) + CALL get_admm_env(qs_env.admm_env, rho_aux_fit=rho_aux, admm_dm=admm_dm) + + CALL qs_rho_get(rho, rho_ao=rho_ao) + CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux) + + ! ** set blocked density matrix to 0 + DO ispin = 1, dft_control.nspins + CALL dbcsr_set(rho_ao_aux(ispin).matrix, 0.0_dp) + ! ** now loop through the list and copy corresponding blocks + CALL dbcsr_iterator_start(iter, rho_ao(ispin).matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) + IF (admm_dm.block_map(iatom, jatom) == 1) THEN + CALL dbcsr_get_block_p(rho_ao_aux(ispin).matrix, & + row=iatom, col=jatom, BLOCK=sparse_block_aux, found=found) + IF (found) & + sparse_block_aux = sparse_block + END IF + END DO + CALL dbcsr_iterator_stop(iter) + END DO + +} + +/*************************************************************************************************** + \brief Call calculate_rho_elec() for auxiliary density + \param qs_env ... +***************************************************************************************************/ +void admm_dm_methods::update_rho_aux(qs_env) +{ + TYPE(qs_environment_type), POINTER :: qs_env + + INTEGER :: ispin + REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_aux + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux + TYPE(dft_control_type), POINTER :: dft_control + TYPE(pw_type), DIMENSION(:), POINTER :: rho_g_aux, rho_r_aux + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(qs_rho_type), POINTER :: rho_aux + TYPE(task_list_type), POINTER :: task_list_aux_fit + + NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, & + task_list_aux_fit, ks_env) + + CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control) + CALL get_admm_env(qs_env.admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux, +& + admm_dm=admm_dm) + + CALL qs_rho_get(rho_aux, & + rho_ao=rho_ao_aux, & + rho_r=rho_r_aux, & + rho_g=rho_g_aux, & + tot_rho_r=tot_rho_r_aux) + + DO ispin = 1, dft_control.nspins + CALL calculate_rho_elec(ks_env=ks_env, & + matrix_p=rho_ao_aux(ispin).matrix, & + rho=rho_r_aux(ispin), & + rho_gspace=rho_g_aux(ispin), & + total_rho=tot_rho_r_aux(ispin), & + soft_valid=.FALSE., & + basis_type="AUX_FIT", & + task_list_external=task_list_aux_fit) + END DO + + CALL qs_rho_set(rho_aux, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) + +} + +/*************************************************************************************************** + \brief Merges auxiliary Kohn-Sham matrix via basis projection. + \param qs_env ... + \param matrix_ks_merge Input: The KS matrix to be merged + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::merge_dm_projection(qs_env, matrix_ks_merge) +{ + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge + + INTEGER :: ispin + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks + TYPE(dbcsr_type) :: matrix_tmp + TYPE(dft_control_type), POINTER :: dft_control + + NULLIFY (admm_dm, dft_control, matrix_ks) + + CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks) + CALL get_admm_env(qs_env.admm_env, admm_dm=admm_dm) + + ! Calculate K += A^T * K_aux * A + CALL dbcsr_create(matrix_tmp, template=admm_dm.matrix_A, matrix_type="N") + + DO ispin = 1, dft_control.nspins + CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks_merge(ispin).matrix, admm_dm.matrix_A, & + 0.0_dp, matrix_tmp) + CALL dbcsr_multiply("T", "N", 1.0_dp, admm_dm.matrix_A, matrix_tmp, & + 1.0_dp, matrix_ks(ispin).matrix) + END DO + + CALL dbcsr_release(matrix_tmp) + +} + +/*************************************************************************************************** + \brief Merges auxiliary Kohn-Sham matrix via blocking. + \param qs_env ... + \param matrix_ks_merge Input: The KS matrix to be merged + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::merge_dm_blocked(qs_env, matrix_ks_merge) +{ + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge + + INTEGER :: blk, iatom, ispin, jatom + REAL(dp), DIMENSION(:, :), POINTER :: sparse_block + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_iterator_type) :: iter + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks + TYPE(dft_control_type), POINTER :: dft_control + + NULLIFY (admm_dm, dft_control, matrix_ks) + + CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks) + CALL get_admm_env(qs_env.admm_env, admm_dm=admm_dm) + + DO ispin = 1, dft_control.nspins + CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin).matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block, blk) + IF (admm_dm.block_map(iatom, jatom) == 0) & + sparse_block = 0.0_dp + END DO + CALL dbcsr_iterator_stop(iter) + CALL dbcsr_add(matrix_ks(ispin).matrix, matrix_ks_merge(ispin).matrix, 1.0_dp, 1.0_dp) + END DO + +} + +/*************************************************************************************************** + \brief Apply McWeeny purification to auxiliary density matrix + \param qs_env ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::purify_mcweeny(qs_env) +{ + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'purify_mcweeny' + + INTEGER :: handle, ispin, istep, nspins, unit_nr + REAL(KIND=dp) :: frob_norm + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, rho_ao_aux + TYPE(dbcsr_type) :: matrix_ps, matrix_psp, matrix_test + TYPE(dbcsr_type), POINTER :: matrix_p, matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mcweeny_history_type), POINTER :: history, new_hist_entry + TYPE(qs_rho_type), POINTER :: rho_aux_fit + + CALL timeset(routineN, handle) + NULLIFY (dft_control, admm_dm, matrix_s_aux_fit, rho_aux_fit, new_hist_entry, & + matrix_p, matrix_s, rho_ao_aux) + + unit_nr = cp_logger_get_default_unit_nr() + CALL get_qs_env(qs_env, dft_control=dft_control) + CALL get_admm_env(qs_env.admm_env, matrix_s_aux_fit=matrix_s_aux_fit, & + rho_aux_fit=rho_aux_fit, admm_dm=admm_dm) + + CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux) + + matrix_p => rho_ao_aux[0].matrix + CALL dbcsr_create(matrix_PS, template=matrix_p, matrix_type="N") + CALL dbcsr_create(matrix_PSP, template=matrix_p, matrix_type="S") + CALL dbcsr_create(matrix_test, template=matrix_p, matrix_type="S") + + nspins = dft_control.nspins + DO ispin = 1, nspins + matrix_p => rho_ao_aux(ispin).matrix + matrix_s => matrix_s_aux_fit[0].matrix + history => admm_dm.mcweeny_history(ispin).p + IF (ASSOCIATED(history)) CPABORT("purify_dm_mcweeny: history already associated") + IF (nspins == 1) CALL dbcsr_scale(matrix_p, 0.5_dp) + + DO istep = 1, admm_dm.mcweeny_max_steps + ! allocate new element in linked list + ALLOCATE (new_hist_entry) + new_hist_entry.next => history + history => new_hist_entry + history.count = istep + NULLIFY (new_hist_entry) + CALL dbcsr_create(history.m, template=matrix_p, matrix_type="N") + CALL dbcsr_copy(history.m, matrix_p, name="P from McWeeny") + + ! calc PS and PSP + CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, & + 0.0_dp, matrix_ps) + + CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ps, matrix_p, & + 0.0_dp, matrix_psp) + + !test convergence + CALL dbcsr_copy(matrix_test, matrix_psp) + CALL dbcsr_add(matrix_test, matrix_p, 1.0_dp, -1.0_dp) + frob_norm = dbcsr_frobenius_norm(matrix_test) + IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5,a,f16.8)') "McWeeny-Step", istep, & + ": Deviation of idempotency", frob_norm + IF (frob_norm < 1000_dp*admm_dm.eps_filter .AND. istep > 1) EXIT + + ! build next P matrix + CALL dbcsr_copy(matrix_p, matrix_PSP, name="P from McWeeny") + CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PS, matrix_PSP, & + 3.0_dp, matrix_p) + END DO + admm_dm.mcweeny_history(ispin).p => history + IF (nspins == 1) CALL dbcsr_scale(matrix_p, 2.0_dp) + END DO + + ! clean up + CALL dbcsr_release(matrix_PS) + CALL dbcsr_release(matrix_PSP) + CALL dbcsr_release(matrix_test) + CALL timestop(handle) +} + +/*************************************************************************************************** + \brief Prepare auxiliary KS-matrix for merge using reverse McWeeny + \param qs_env ... + \param matrix_ks_merge Output: The KS matrix for the merge + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::revert_purify_mcweeny(qs_env, matrix_ks_merge) +{ + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_merge + + CHARACTER(LEN=*), PARAMETER :: routineN = 'revert_purify_mcweeny' + + INTEGER :: handle, ispin, nspins, unit_nr + TYPE(admm_dm_type), POINTER :: admm_dm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, & + matrix_s_aux_fit, & + matrix_s_aux_fit_vs_orb + TYPE(dbcsr_type), POINTER :: matrix_k + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mcweeny_history_type), POINTER :: history_curr, history_next + + CALL timeset(routineN, handle) + unit_nr = cp_logger_get_default_unit_nr() + NULLIFY (admm_dm, dft_control, matrix_ks, matrix_ks_aux_fit, & + matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, & + history_next, history_curr, matrix_k) + + CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks) + CALL get_admm_env(qs_env.admm_env, matrix_s_aux_fit=matrix_s_aux_fit, admm_dm=admm_dm, & + matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, +matrix_ks_aux_fit=matrix_ks_aux_fit) + + nspins = dft_control.nspins + ALLOCATE (matrix_ks_merge(nspins)) + + DO ispin = 1, nspins + ALLOCATE (matrix_ks_merge(ispin).matrix) + matrix_k => matrix_ks_merge(ispin).matrix + CALL dbcsr_copy(matrix_k, matrix_ks_aux_fit(ispin).matrix, name="K") + history_curr => admm_dm.mcweeny_history(ispin).p + NULLIFY (admm_dm.mcweeny_history(ispin).p) + + ! reverse McWeeny iteration + DO WHILE (ASSOCIATED(history_curr)) + IF (unit_nr > 0) WRITE (unit_nr, '(t3,a,i5)') "Reverse McWeeny-Step ", history_curr.count + CALL reverse_mcweeny_step(matrix_k=matrix_k, & + matrix_s=matrix_s_aux_fit[0].matrix, & + matrix_p=history_curr.m) + CALL dbcsr_release(history_curr.m) + history_next => history_curr.next + DEALLOCATE (history_curr) + history_curr => history_next + NULLIFY (history_next) + END DO + + END DO + + ! clean up + timestop(handle); + +} + +/*************************************************************************************************** + \brief Multiply matrix_k with partial derivative of McWeeny by reversing it. + \param matrix_k ... + \param matrix_s ... + \param matrix_p ... + \author Ole Schuett +***************************************************************************************************/ +void admm_dm_methods::reverse_mcweeny_step(matrix_k, matrix_s, matrix_p) +{ + TYPE(dbcsr_type) :: matrix_k, matrix_s, matrix_p + + char *routineN = "reverse_mcweeny_step"; + + int handle; + TYPE(dbcsr_type) :: matrix_ps, matrix_sp, matrix_sum, & + matrix_tmp + + timeset(routineN, handle); + dbcsr_create(matrix_ps, template=matrix_p, matrix_type="N"); + dbcsr_create(matrix_sp, template=matrix_p, matrix_type="N"); + dbcsr_create(matrix_tmp, template=matrix_p, matrix_type="N"); + dbcsr_create(matrix_sum, template=matrix_p, matrix_type="N"); + + dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s, + 0.0_dp, matrix_ps); + dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_p, + 0.0_dp, matrix_sp); + + // TODO: can we exploid more symmetry? + dbcsr_multiply("N", "N", 3.0_dp, matrix_k, matrix_ps, + 0.0_dp, matrix_sum); + dbcsr_multiply("N", "N", 3.0_dp, matrix_sp, matrix_k, + 1.0_dp, matrix_sum); + + // matrix_tmp = KPS + dbcsr_multiply("N", "N", 1.0_dp, matrix_k, matrix_ps, + 0.0_dp, matrix_tmp); + dbcsr_multiply("N", "N", -2.0_dp, matrix_tmp, matrix_ps, + 1.0_dp, matrix_sum); + dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, + 1.0_dp, matrix_sum); + + // matrix_tmp = SPK + dbcsr_multiply("N", "N", 1.0_dp, matrix_sp, matrix_k, + 0.0_dp, matrix_tmp); + dbcsr_multiply("N", "N", -2.0_dp, matrix_sp, matrix_tmp, + 1.0_dp, matrix_sum); + + // overwrite matrix_k + dbcsr_copy(matrix_k, matrix_sum, name="K from reverse McWeeny"); + + // clean up + dbcsr_release(matrix_sum); + dbcsr_release(matrix_tmp); + dbcsr_release(matrix_ps); + dbcsr_release(matrix_sp); + timestop(handle); + +} diff --git a/src/admm_utils.c b/src/admm_utils.c new file mode 100644 index 0000000..56af759 --- /dev/null +++ b/src/admm_utils.c @@ -0,0 +1,179 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Contains methods used in the context of density fitting +//> \par History +//> 04.2008 created [Manuel Guidon] +//> 02.2013 moved from admm_methods +//> \author Manuel Guidon +// ************************************************************************************************** +#include "./base/base_uses.h" + +void admm_utils() { + USE admm_types, ONLY: admm_type + USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr + USE cp_gemm_interface, ONLY: cp_gemm + USE dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_create,& + dbcsr_deallocate_matrix,& + dbcsr_set,& + dbcsr_type,& + dbcsr_type_symmetric + USE input_constants, ONLY: do_admm_purify_cauchy,& + do_admm_purify_cauchy_subspace,& + do_admm_purify_mo_diag,& + do_admm_purify_mo_no_diag,& + do_admm_purify_none + USE kinds, ONLY: dp + + IMPLICIT NONE + PRIVATE + + PUBLIC :: admm_correct_for_eigenvalues, & + admm_uncorrect_for_eigenvalues + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_utils' + +//*** + +CONTAINS + +// ************************************************************************************************** +//> \brief ... +//> \param ispin ... +//> \param admm_env ... +//> \param ks_matrix ... +// ************************************************************************************************** + SUBROUTINE admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix) + INTEGER, INTENT(IN) :: ispin + TYPE(admm_type), POINTER :: admm_env + TYPE(dbcsr_type), POINTER :: ks_matrix + + INTEGER :: nao_aux_fit, nao_orb + TYPE(dbcsr_type), POINTER :: work + + nao_aux_fit = admm_env%nao_aux_fit + nao_orb = admm_env%nao_orb + + IF (.NOT. admm_env%block_dm) THEN + SELECT CASE (admm_env%purification_method) + CASE (do_admm_purify_cauchy_subspace) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + + // ** calculate A^T*H_tilde*A + CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%K(ispin)%matrix, admm_env%A, 0.0_dp, & + admm_env%work_aux_orb) + CALL cp_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, & + admm_env%H_corr(ispin)%matrix) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_diag) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + // ** calculate A^T*H_tilde*A + CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%K(ispin)%matrix, admm_env%A, 0.0_dp, & + admm_env%work_aux_orb) + CALL cp_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, & + admm_env%H_corr(ispin)%matrix) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_cauchy) + // do nothing + END SELECT + END IF + + END SUBROUTINE admm_correct_for_eigenvalues + +// ************************************************************************************************** +//> \brief ... +//> \param ispin ... +//> \param admm_env ... +//> \param ks_matrix ... +// ************************************************************************************************** + SUBROUTINE admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix) + INTEGER, INTENT(IN) :: ispin + TYPE(admm_type), POINTER :: admm_env + TYPE(dbcsr_type), POINTER :: ks_matrix + + INTEGER :: nao_aux_fit, nao_orb + TYPE(dbcsr_type), POINTER :: work + + nao_aux_fit = admm_env%nao_aux_fit + nao_orb = admm_env%nao_orb + + IF (.NOT. admm_env%block_dm) THEN + SELECT CASE (admm_env%purification_method) + CASE (do_admm_purify_cauchy_subspace) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_diag) + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_cauchy) + // do nothing + END SELECT + END IF + END SUBROUTINE admm_uncorrect_for_eigenvalues + +} diff --git a/src/base/PACKAGE b/src/base/PACKAGE new file mode 100644 index 0000000..397ef51 --- /dev/null +++ b/src/base/PACKAGE @@ -0,0 +1,4 @@ +{ + "description": "base routines needed to abstract away some machine/compiler dependent functionality", + "requires": [], +} diff --git a/src/base/base.h b/src/base/base.h new file mode 100644 index 0000000..dedff48 --- /dev/null +++ b/src/base/base.h @@ -0,0 +1,8 @@ +#ifndef _BASE_H +#define _BASE_H + + #define default_string_length 80 + #define default_path_length 1024 + const int max_line_length = 2*default_path_length; + +#endif diff --git a/src/base/base_hooks.cpp b/src/base/base_hooks.cpp new file mode 100644 index 0000000..8b099f0 --- /dev/null +++ b/src/base/base_hooks.cpp @@ -0,0 +1,238 @@ +/*-------------------------------------------------------------------------------------------------- + CP2K: A general program to perform molecular dynamics simulations + Copyright 2000-2023 CP2K developers group + + SPDX-License-Identifier: GPL-2.0-or-later +--------------------------------------------------------------------------------------------------*/ + +/*************************************************************************************************** + \brief Central dispatch for basic hooks + \author Ole Schuett +***************************************************************************************************/ +class base_hooks { + +private: + +public: + // API + void cp_abort(); + void cp_warn(); + void cp_hint(); + void timeset(); + void timestop(); + void cp_abort_hook(); + void cp_warn_hook(); + void cp_hint_hook(); + void timeset_hook(); + void timestop_hook(); + void cp__a(); + void cp__b(); + void cp__w(); + void cp__h(); + void cp__l(); + + // this interface (with subroutines in it) must to be defined right before + // the regular subroutines/functions - otherwise prettify.py will screw up. + INTERFACE + SUBROUTINE cp_abort_interface(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + END SUBROUTINE cp_abort_interface + + SUBROUTINE cp_warn_interface(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + END SUBROUTINE cp_warn_interface + + SUBROUTINE cp_hint_interface(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + END SUBROUTINE cp_hint_interface + + SUBROUTINE timeset_interface(routineN, handle) + CHARACTER(LEN=*), INTENT(IN) :: routineN + INTEGER, INTENT(OUT) :: handle + + END SUBROUTINE timeset_interface + + SUBROUTINE timestop_interface(handle) + INTEGER, INTENT(IN) :: handle + + END SUBROUTINE timestop_interface + END INTERFACE + + PROCEDURE(cp_abort_interface), POINTER :: cp_abort_hook => Null() + PROCEDURE(cp_warn_interface), POINTER :: cp_warn_hook => Null() + PROCEDURE(cp_hint_interface), POINTER :: cp_hint_hook => Null() + PROCEDURE(timeset_interface), POINTER :: timeset_hook => Null() + PROCEDURE(timestop_interface), POINTER :: timestop_hook => Null() + +} + +/*************************************************************************************************** + \brief Terminate the program + \param location ... + \param message ... + \author Ole Schuett +***************************************************************************************************/ + SUBROUTINE cp_abort(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + IF (ASSOCIATED(cp_abort_hook)) THEN + CALL cp_abort_hook(location, message) + ELSE + WRITE (default_output_unit, *) "ABORT in "//TRIM(location)//" "//TRIM(message) + CALL m_flush(default_output_unit) + CALL m_abort() + END IF + ! compiler hint + STOP "Never return from here" + END SUBROUTINE cp_abort + +! ************************************************************************************************** +!> \brief Issue a warning +!> \param location ... +!> \param message ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE cp_warn(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + IF (ASSOCIATED(cp_warn_hook)) THEN + CALL cp_warn_hook(location, message) + ELSE + WRITE (default_output_unit, *) "WARNING in "//TRIM(location)//" "//TRIM(message) + CALL m_flush(default_output_unit) + END IF + END SUBROUTINE cp_warn + +! ************************************************************************************************** +!> \brief Issue a hint +!> \param location ... +!> \param message ... +!> \author Hans Pabst +! ************************************************************************************************** + SUBROUTINE cp_hint(location, message) + CHARACTER(len=*), INTENT(in) :: location, message + + IF (ASSOCIATED(cp_hint_hook)) THEN + CALL cp_hint_hook(location, message) + ELSE + WRITE (default_output_unit, *) "HINT in "//TRIM(location)//" "//TRIM(message) + CALL m_flush(default_output_unit) + END IF + END SUBROUTINE cp_hint + +! ************************************************************************************************** +!> \brief Start timer +!> \param routineN ... +!> \param handle ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE timeset(routineN, handle) + CHARACTER(LEN=*), INTENT(IN) :: routineN + INTEGER, INTENT(OUT) :: handle + + IF (ASSOCIATED(timeset_hook)) THEN + CALL timeset_hook(routineN, handle) + ELSE + handle = -1 + END IF + END SUBROUTINE timeset + +! ************************************************************************************************** +!> \brief Stop timer +!> \param handle ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE timestop(handle) + INTEGER, INTENT(IN) :: handle + + IF (ASSOCIATED(timestop_hook)) THEN + CALL timestop_hook(handle) + ELSE + IF (handle /= -1) & + CALL cp_abort(cp__l("base_hooks.F", __LINE__), "Got wrong handle") + END IF + END SUBROUTINE timestop + +! ************************************************************************************************** +!> \brief CPASSERT handler +!> \param filename ... +!> \param lineNr ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE cp__a(filename, lineNr) + CHARACTER(len=*), INTENT(in) :: filename + INTEGER, INTENT(in) :: lineNr + + CALL cp_abort(location=cp__l(filename, lineNr), message="CPASSERT failed") + ! compiler hint + STOP "Never return from here" + END SUBROUTINE cp__a + +! ************************************************************************************************** +!> \brief CPABORT handler +!> \param filename ... +!> \param lineNr ... +!> \param message ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE cp__b(filename, lineNr, message) + CHARACTER(len=*), INTENT(in) :: filename + INTEGER, INTENT(in) :: lineNr + CHARACTER(len=*), INTENT(in) :: message + + CALL cp_abort(location=cp__l(filename, lineNr), message=message) + ! compiler hint + STOP "Never return from here" + END SUBROUTINE cp__b + +! ************************************************************************************************** +!> \brief CPWARN handler +!> \param filename ... +!> \param lineNr ... +!> \param message ... +!> \author Ole Schuett +! ************************************************************************************************** + SUBROUTINE cp__w(filename, lineNr, message) + CHARACTER(len=*), INTENT(in) :: filename + INTEGER, INTENT(in) :: lineNr + CHARACTER(len=*), INTENT(in) :: message + + CALL cp_warn(location=cp__l(filename, lineNr), message=message) + END SUBROUTINE cp__w + +! ************************************************************************************************** +!> \brief CPHINT handler +!> \param filename ... +!> \param lineNr ... +!> \param message ... +!> \author Hans Pabst +! ************************************************************************************************** + SUBROUTINE cp__h(filename, lineNr, message) + CHARACTER(len=*), INTENT(in) :: filename + INTEGER, INTENT(in) :: lineNr + CHARACTER(len=*), INTENT(in) :: message + + CALL cp_hint(location=cp__l(filename, lineNr), message=message) + END SUBROUTINE cp__h + +/*************************************************************************************************** + \brief Helper routine to assemble __LOCATION__ + \param filename ... + \param lineNr ... + \return ... + \author Ole Schuett +***************************************************************************************************/ +char* cp__l(char *filename, int lineNr) +{ + char *location; + + char lineNr_str[15]; + + WRITE (lineNr_str, FMT='(I10)') lineNr + location = TRIM(filename)//":"//TRIM(ADJUSTL(lineNr_str)) + + return location; +} diff --git a/src/base/base_uses.cpp b/src/base/base_uses.cpp new file mode 100644 index 0000000..e0b5a4c --- /dev/null +++ b/src/base/base_uses.cpp @@ -0,0 +1,59 @@ +#ifndef _BASE_USES_H +#define _BASE_USES_H + +// Basic use statements and preprocessor macros +// should be included in the use statements + +int cp__a, cp__b, cp__w, cp__h, cp__l, cp_abort, cp_warn, cp_hint, timeset, timestop; + +#if defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP) +#define __OFFLOAD +#endif + +// Check for OpenMP early on - ideally before the compiler fails with a cryptic message. +#if !defined(_OPENMP) + printf("OpenMP is required. Please add the corresponding flag (eg. -fopenmp for GFortran) to your Fortran compiler flags."); +#endif + +// Dangerous: Full path can be arbitrarily long and might overflow Fortran line. +#if !defined(__SHORT_FILE__) +#define __SHORT_FILE__ __FILE__ +#endif + +#define __LOCATION__ cp__l(__SHORT_FILE__,__LINE__) +#define CPWARN(msg) CALL cp__w(__SHORT_FILE__,__LINE__,msg) +#define CPABORT(msg) CALL cp__b(__SHORT_FILE__,__LINE__,msg) +// In contrast to CPWARN, the warning counter is not increased +#define CPHINT(msg) CALL cp__h(__SHORT_FILE__,__LINE__,msg) + +// CPASSERT can be elided if NDEBUG is defined. +#if defined(NDEBUG) +# define CPASSERT(cond) +#else +# define CPASSERT(cond) IF(.NOT.(cond))CALL cp__a(__SHORT_FILE__,__LINE__) +#endif + +// The MARK_USED macro can be used to mark an argument/variable as used. It is intended to make +// it possible to switch on -Werror=unused-dummy-argument, but deal elegantly with, e.g., +// library wrapper routines that take arguments only used if the library is linked in. +// This code should be valid for any Fortran variable, is always standard conforming, +// and will be optimized away completely by the compiler +#define MARK_USED(foo) IF(.FALSE.)THEN; DO ; IF(SIZE(SHAPE(foo))==-1) EXIT ; END DO ; ENDIF + +// Calculate version number from 2 or 3 components. Can be used for comparison, e.g., +// CPVERSION3(4, 9, 0) <= CPVERSION3(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__) +// CPVERSION(8, 0) <= CPVERSION(__GNUC__, __GNUC_MINOR__) +#define CPVERSION2(MAJOR, MINOR) ((MAJOR) * 10000 + (MINOR) * 100) +#define CPVERSION3(MAJOR, MINOR, UPDATE) (CPVERSION2(MAJOR, MINOR) + (UPDATE)) +#define CPVERSION CPVERSION2 + +// gfortran before 8.3 complains about internal symbols not being specified in +// any data clause when using DEFAULT(NONE) and OOP procedures are called from +// within the parallel region. +#if __GNUC__ < 8 || (__GNUC__ == 8 && (__GNUC_MINOR__ < 3)) +#define OMP_DEFAULT_NONE_WITH_OOP SHARED +#else +#define OMP_DEFAULT_NONE_WITH_OOP NONE +#endif + +#endif diff --git a/src/base/kinds.cpp b/src/base/kinds.cpp new file mode 100644 index 0000000..cb77ae3 --- /dev/null +++ b/src/base/kinds.cpp @@ -0,0 +1,51 @@ +/*------------------------------------------------------------------------------------------------ + CP2K: A general program to perform molecular dynamics simulations + Copyright 2000-2023 CP2K developers group + + SPDX-License-Identifier: GPL-2.0-or-later +-------------------------------------------------------------------------------------------------*/ + +class kinds { + +public: +/*************************************************************************************************** + \brief Print informations about the used data types. + \param iw ... + \par History + Adapted by JGH for Cp2k + \author Matthias Krack +***************************************************************************************************/ +void print_kind_info(int iw) { + + WRITE (iw, '( /, T2, A )') 'DATA TYPE INFORMATION:' + + WRITE (iw, '( /,T2,A,T79,A,2(/,T2,A,T75,I6),3(/,T2,A,T67,E14.8) )') & + 'REAL: Data type name:', 'dp', ' Kind value:', KIND(0.0_dp), & + ' Precision:', PRECISION(0.0_dp), & + ' Smallest non-negligible quantity relative to 1:', & + EPSILON(0.0_dp), & + ' Smallest positive number:', TINY(0.0_dp), & + ' Largest representable number:', HUGE(0.0_dp) + WRITE (iw, '( /,T2,A,T79,A,2(/,T2,A,T75,I6),3(/,T2,A,T67,E14.8) )') & + ' Data type name:', 'sp', ' Kind value:', KIND(0.0_sp), & + ' Precision:', PRECISION(0.0_sp), & + ' Smallest non-negligible quantity relative to 1:', & + EPSILON(0.0_sp), & + ' Smallest positive number:', TINY(0.0_sp), & + ' Largest representable number:', HUGE(0.0_sp) + WRITE (iw, '( /,T2,A,T72,A,4(/,T2,A,T61,I20) )') & + 'INTEGER: Data type name:', '(default)', ' Kind value:', & + KIND(0), & + ' Bit size:', BIT_SIZE(0), & + ' Largest representable number:', HUGE(0) + WRITE (iw, '( /,T2,A,T72,A,/,T2,A,T75,I6,/ )') & + 'LOGICAL: Data type name:', '(default)', & + ' Kind value:', KIND(.TRUE.) + WRITE (iw, '( /,T2,A,T72,A,/,T2,A,T75,I6,/ )') & + 'CHARACTER: Data type name:', '(default)', & + ' Kind value:', KIND('C') + +} + + +} diff --git a/src/base/machine.c b/src/base/machine.c new file mode 100644 index 0000000..4aae8e5 --- /dev/null +++ b/src/base/machine.c @@ -0,0 +1,117 @@ +//------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2022 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//------------------------------------------------------------------------------------------------// + +// ************************************************************************************************* +//> \brief Machine interface based on Fortran 2003 and POSIX +//> \par History +//> JGH (05.07.2001) : added G95 interface +//> - m_flush added (12.06.2002,MK) +//> - Missing print_memory added (24.09.2002,MK) +//> - Migrate to generic implementation based on F2003 + POSIX (2014, Ole Schuett) +//> \author APSI, JGH, Ole Schuett +// ************************************************************************************************* + +#include + +#include "machine_cpuid.h" + +#if defined(__LIBXSMM) + //do something +#endif + + // Except for some error handling code, all code should + // get a unit number from the print keys or from the logger, in order + // to guarantee correct output behavior, + // for example in farming or path integral runs + // default_input_unit should never be used + // but we need to know what it is, as we should not try to open it for output + + int default_output_unit = output_unit; + int default_input_unit = input_unit; + + // Enumerates the target architectures or instruction set extensions. + // A feature is present if within range for the respective architecture. + // For example, to check for MACHINE_X86_AVX the following is true: + // MACHINE_X86_AVX <= m_cpuid() and MACHINE_X86 >= m_cpuid(). + // For example, to check for MACHINE_ARM_SOME the following is true: + // MACHINE_ARM_SOME <= m_cpuid() and MACHINE_ARM >= m_cpuid(). + + int MACHINE_CPU_GENERIC = CP_MACHINE_CPU_GENERIC; + int MACHINE_X86_SSE4 = CP_MACHINE_X86_SSE4; + int MACHINE_X86_AVX = CP_MACHINE_X86_AVX; + int MACHINE_X86_AVX2 = CP_MACHINE_X86_AVX2; + int MACHINE_X86_AVX512 = CP_MACHINE_X86_AVX512; + int MACHINE_X86 = MACHINE_X86_AVX512; // marks end of range + + // other arch to be added as needed e.g., + //MACHINE_ARM_SOME = 2000 + //MACHINE_ARM_ELSE = 2001 + //MACHINE_ARM = MACHINE_ARM_ELSE + //MACHINE_PWR_???? = 3000 + + // Flushing is enabled by default because without it crash reports can get lost. + // For performance reasons it can be disabled via the input in &GLOBAL. + bool flush_should_flush = true; + + int m_memory_max = 0; + +// ************************************************************************************************* +//> \brief flushes units if the &GLOBAL flag is set accordingly +//> \param lunit ... +//> \par History +//> 10.2008 created [Joost VandeVondele] +//> \note +//> flushing might degrade performance significantly (30% and more) +// ************************************************************************************************* + +void m_flush(std::ofstream &lunit) { + if (flush_should_flush) {std::flush(lunit);} +} + +// ************************************************************************************************* +//> \brief returns time from a real-time clock, protected against rolling +//> early/easily +//> \return ... +//> \par History +//> 03.2006 created [Joost VandeVondele] +//> \note +//> same implementation for all machines. +//> might still roll, if not called multiple times per count_max/count_rate +// ************************************************************************************************* +double m_walltime() { + double wt; + + #ifdef __LIBXSMM + wt = libxsmm_timer_duration(0_int_8, libxsmm_timer_tick()); + #else + wt = omp_get_wtime(); + #endif + + return wt; +} + +// ************************************************************************************************* +//> \brief reads /proc/cpuinfo if it exists (i.e. Linux) to return relevant info +//> \param model_name as obtained from the 'model name' field, UNKNOWN otherwise +// ************************************************************************************************* +void m_cpuinfo(std::string model_name) { + int bufferlen = 2048; + + char buffer[bufferlen]; + int i, icol, iline, imod, stat; + + model_name = "UNKNOWN"; + buffer = ""; + + +} + + + + + + diff --git a/src/base/machine_cpuid.c b/src/base/machine_cpuid.c new file mode 100644 index 0000000..b4c04f7 --- /dev/null +++ b/src/base/machine_cpuid.c @@ -0,0 +1,42 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2021 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/* shared between C and Fortran */ +#include "machine_cpuid.h" + +#if defined(__cplusplus) +extern "C" { +#endif + +/******************************************************************************* + * \brief This routine determines the CPUID according to the given compiler + * flags (expected to be similar to Fortran). Similar to other Fortran + * compilers, "gfortran -E -dM -mavx - < /dev/null | grep AVX" defines a + * variety of predefined macros (also similar to C). However, with a + * Fortran translation unit only a subset of these definitions disappears + * ("gfortran -E -dM -mavx my.F | grep AVX") + * hence an implementation in C is used. + ******************************************************************************/ +int m_cpuid_static(void); /* avoid pedantic warning about missing prototype */ +int m_cpuid_static(void) { +#if (__AVX512F__ && __AVX512CD__ && __AVX2__ && __FMA__ && __AVX__ && \ + __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX512; +#elif (__AVX2__ && __FMA__ && __AVX__ && __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX2; +#elif (__AVX__ && __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX; +#elif (__SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_SSE4; +#else + return CP_MACHINE_CPU_GENERIC; +#endif +} + +#if defined(__cplusplus) +} +#endif diff --git a/src/base/machine_cpuid.h b/src/base/machine_cpuid.h new file mode 100644 index 0000000..968c41b --- /dev/null +++ b/src/base/machine_cpuid.h @@ -0,0 +1,12 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2021 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +#define CP_MACHINE_CPU_GENERIC 0 +#define CP_MACHINE_X86_SSE4 1000 +#define CP_MACHINE_X86_AVX 1001 +#define CP_MACHINE_X86_AVX2 1002 +#define CP_MACHINE_X86_AVX512 1003 diff --git a/src/cp2k_info.cpp b/src/cp2k_info.cpp new file mode 100644 index 0000000..7de5745 --- /dev/null +++ b/src/cp2k_info.cpp @@ -0,0 +1,357 @@ +/*-------------------------------------------------------------------------------------------------- + CP2K: A general program to perform molecular dynamics simulations + Copyright 2000-2023 CP2K developers group + + SPDX-License-Identifier: GPL-2.0-or-later +--------------------------------------------------------------------------------------------------*/ + +/*************************************************************************************************** + \brief some minimal info about CP2K, including its version and license + \par History + - created (2007-09, Joost VandeVondele) + - moved into this module information related to runtime:pid, user_name, + host_name, cwd, datx (2009-06, Teodoro Laino) + \author Joost VandeVondele +***************************************************************************************************/ + +#include +#include +#include + +#include "./base/base.h" + +class cp2k_info { + +private: + char *moduleN = "cp2k_info"; + +public: + void cp2k_flags(); + void print_cp2k_license(); + void get_runtime_info(); + void write_restart_header(); + int cp2k_version, cp2k_year, cp2k_home; + int compile_arch, compile_date, compile_host, compile_revision; +#ifdef __COMPILE_REVISION + std::string compile_revision = __COMPILE_REVISION; +#else + std::string compile_revision = "unknown"; +#endif + +std::string cp2k_version = "CP2K version 2022.1 (Development Version)"; +std::string cp2k_year = "2022"; +std::string cp2k_home = "https://www.cp2k.org/"; + +// compile time information +#ifdef __COMPILE_ARCH + std::string compile_arch = __COMPILE_ARCH; +#else + std::string compile_arch = "unknown: -D__COMPILE_ARCH=?"; +#endif + +#ifdef __COMPILE_DATE + std::string compile_date = __COMPILE_DATE; +#else + std::string compile_date = "unknown: -D__COMPILE_DATE=?"; +#endif + +#ifdef __COMPILE_HOST + std::string compile_host = __COMPILE_HOST; +#else + std::string compile_host = "unknown: -D__COMPILE_HOST=?"; +#endif + +// local runtime informations +std::string r_datx, r_cwd, r_host_name, r_user_name; +int r_pid; + +std::string moduleN = "cp2k_info"; + +}; + +/*************************************************************************************************** + \brief list all compile time options that influence the capabilities of cp2k. + All new flags should be added here (and be unique grep-able) + \return ... +***************************************************************************************************/ + +std::string cp2k_info::cp2k_flags() { + std::string flags, tmp_str; + + flags = "cp2kflags:"; + + // Ensure that tmp_str is used to silence compiler warnings + tmp_str = ""; + flags = rtrim(flags)+tmp_str; + + #ifdef NBEBUG + flags = rtrim(flags)+" ndebug"; + #endif + flags = rtrim(flags)+" omp" + #if defined(__LIBINT) + flags = rtrim(flags)+" libint" + #endif + #if defined(__FFTW3) + flags = rtrim(flags)+" fftw3" + #endif + #if defined(__FFTW3_MKL) + flags = rtrim(flags)+" fftw3_mkl" + #endif + #if defined(__LIBXC) + flags = rtrim(flags)+" libxc" + #endif + #if defined(__LIBPEXSI) + flags = rtrim(flags)+" pexsi" + #endif + #if defined(__ELPA) + flags = rtrim(flags)+" elpa" + #endif + #if defined(__ELPA_NVIDIA_GPU) + flags = rtrim(flags)+" elpa_nvidia_gpu" + #endif + #if defined(__ELPA_AMD_GPU) + flags = rtrim(flags)+" elpa_amd_gpu" + #endif + #if defined(__ELPA_INTEL_GPU) + flags = rtrim(flags)+" elpa_intel_gpu" + #endif + #if defined(__parallel) + flags = rtrim(flags)+" parallel" + #if !defined(__MPI_VERSION) || (__MPI_VERSION > 2) + flags = rtrim(flags)+" mpi3" + #else + flags = rtrim(flags)+" mpi2" + #endif + #endif + #if defined(__SCALAPACK) + flags = rtrim(flags)+" scalapack" + #endif + #if defined(__COSMA) + flags = rtrim(flags)+" cosma" + #endif + + #if defined(__QUIP) + flags = rtrim(flags)+" quip" + #endif + + #if defined(__HAS_PATCHED_CUFFT_70) + flags = rtrim(flags)+" patched_cufft_70" + #endif + + #if defined(__PW_FPGA) + flags = rtrim(flags)+" pw_fpga" + #endif + #if defined(__PW_FPGA_SP) + flags = rtrim(flags)+" pw_fpga_sp" + #endif + + #if defined(__LIBXSMM) + flags = rtrim(flags)+" xsmm" + #endif + + #if defined(__CRAY_PM_ACCEL_ENERGY) + flags = rtrim(flags)+" cray_pm_accel_energy" + #endif + #if defined(__CRAY_PM_ENERGY) + flags = rtrim(flags)+" cray_pm_energy" + #endif + #if defined(__CRAY_PM_FAKE_ENERGY) + flags = rtrim(flags)+" cray_pm_fake_energy" + #endif + #if defined(__DBCSR_ACC) + flags = rtrim(flags)+" dbcsr_acc" + #endif + #if defined(__MAX_CONTR) + CALL integer_to_string(__MAX_CONTR, tmp_str) + flags = rtrim(flags)+" max_contr="+rtrim(tmp_str) + #endif + #if defined(__NO_IPI_DRIVER) + flags = rtrim(flags)+" no_ipi_driver" + #endif + #if defined(__NO_MPI_THREAD_SUPPORT_CHECK) + flags = rtrim(flags)+" no_mpi_thread_support_check" + #endif + #if defined(__NO_STATM_ACCESS) + flags = rtrim(flags)+" no_statm_access" + #endif + #if defined(__MINGW) + flags = rtrim(flags)+" mingw" + #endif + #if defined(__PW_CUDA_NO_HOSTALLOC) + flags = rtrim(flags)+" pw_cuda_no_hostalloc" + #endif + #if defined(__STATM_RESIDENT) + flags = rtrim(flags)+" statm_resident" + #endif + #if defined(__STATM_TOTAL) + flags = rtrim(flags)+" statm_total" + #endif + #if defined(__PLUMED2) + flags = rtrim(flags)+" plumed2" + #endif + #if defined(__HAS_IEEE_EXCEPTIONS) + flags = rtrim(flags)+" has_ieee_exceptions" + #endif + #if defined(__NO_ABORT) + flags = rtrim(flags)+" no_abort" + #endif + #if defined(__SPGLIB) + flags = rtrim(flags)+" spglib" + #endif + #if defined(__ACCELERATE) + flags = rtrim(flags)+" accelerate" + #endif + #if defined(__MKL) + flags = rtrim(flags)+" mkl" + #endif + #if defined(__SIRIUS) + flags = rtrim(flags)+" sirius" + #endif + #if defined(__CHECK_DIAG) + flags = rtrim(flags)+" check_diag" + #endif + #if defined(__LIBVORI) + flags = rtrim(flags)+" libvori" + flags = rtrim(flags)+" libbqb" + #endif + #if defined(__LIBMAXWELL) + flags = rtrim(flags)+" libmaxwell" + #endif + #if defined(__LIBTORCH) + flags = rtrim(flags)+" libtorch" + #endif + #if defined(__OFFLOAD_CUDA) + flags = rtrim(flags)+" offload_cuda" + #endif + #if defined(__OFFLOAD_HIP) + flags = rtrim(flags)+" offload_hip" + #endif + #if defined(__NO_OFFLOAD_GRID) + flags = rtrim(flags)+" no_offload_grid" + #endif + #if defined(__NO_OFFLOAD_DBM) + flags = rtrim(flags)+" no_offload_dbm" + #endif + #if defined(__NO_OFFLOAD_PW) + flags = rtrim(flags)+" no_offload_pw" + #endif + #if defined(__OFFLOAD_PROFILING) + flags = rtrim(flags)+" offload_profiling" + #endif + + #if defined(__SPLA) && defined(__OFFLOAD_GEMM) + flags = rtrim(flags)+" spla_gemm_offloading" + #endif + + #if defined(__CUSOLVERMP) + flags = rtrim(flags)+" cusolvermp" + #endif + + #if defined(__LIBVDWXC) + flags = rtrim(flags)+" libvdwxc" + #endif + + #if defined(__HDF5) + flags = rtrim(flags)+" hdf5" + #endif + + return flags; +} + +/*************************************************************************************************** + \brief ... + \param iunit ... +***************************************************************************************************/ +void cp2k_info::print_cp2k_license(std::ofstream iunit) +{ + + iunit << "!-----------------------------------------------------------------------------!" << std::endl; + iunit << "! !" << std::endl; + iunit << "! CP2K: A general program to perform molecular dynamics simulations !" << std::endl; + iunit << "! Copyright (C) 2000, 2001, 2002, 2003 CP2K developers group !" << std::endl; + iunit << "! Copyright (C) 2004, 2005, 2006, 2007 CP2K developers group !" << std::endl; + iunit << "! Copyright (C) 2008, 2009, 2010, 2011 CP2K developers group !" << std::endl; + iunit << "! Copyright (C) 2012, 2013, 2014, 2015 CP2K developers group !" << std::endl; + iunit << "! Copyright (C) 2016, 2017, 2018, 2019 CP2K developers group !" << std::endl; + iunit << "! Copyright (C) 2020 CP2K developers group !" << std::endl; + iunit << "! !" << std::endl; + iunit << "! This program is free software; you can redistribute it and/or modify !" << std::endl; + iunit << "! it under the terms of the GNU General Public License as published by !" << std::endl; + iunit << "! the Free Software Foundation; either version 2 of the License, or !" << std::endl; + iunit << "! (at your option) any later version. !" << std::endl; + iunit << "! !" << std::endl; + iunit << "! This program is distributed in the hope that it will be useful, !" << std::endl; + iunit << "! but WITHOUT ANY WARRANTY; without even the implied warranty of !" << std::endl; + iunit << "! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !" << std::endl; + iunit << "! GNU General Public License for more details. !" << std::endl; + iunit << "! !" << std::endl; + iunit << "! You should have received a copy of the GNU General Public License !" << std::endl; + iunit << "! along with this program; if not, write to the Free Software !" << std::endl; + iunit << "! Foundation, Inc., 51 Franklin Street, Fifth Floor, !" << std::endl; + iunit << "! Boston, MA 02110-1301, USA. !" << std::endl; + iunit << "! !" << std::endl; + iunit << "! See also https://www.fsf.org/licensing/licenses/gpl.html !" << std::endl; + iunit << "! !" << std::endl; + iunit << "!-----------------------------------------------------------------------------!" << std::endl; + iunit << "! CP2K, including its sources and pointers to the authors !" << std::endl; + iunit << "! can be found at https://www.cp2k.org/ !" << std::endl; + iunit << "!-----------------------------------------------------------------------------!" << std::endl; +} + +/*************************************************************************************************** + \brief ... +***************************************************************************************************/ +void cp2k_info::get_runtime_info() +{ + r_datx = ""; + r_cwd = ""; + r_host_name = ""; + r_user_name = ""; + r_pid = -1; + + m_getpid(&r_pid); + m_getlog(&r_user_name); + m_hostnm(&r_host_name); + m_datum(&r_datx); + m_getcwd(&r_cwd); + +} + +/*************************************************************************************************** + \brief Writes the header for the restart file + \param iunit ... + \par History + 01.2008 [created] - Split from write_restart + \author Teodoro Laino - University of Zurich - 01.2008 +***************************************************************************************************/ +void cp2k_info::write_restart_header(std::ofstream &iunit) +{ + char cwd[255], datx[255]; + + m_datum(&datx); + m_getcwd(&cwd); + + iunit << " # Version information for this restart file " << std::endl; + iunit << " # current date "+rtrim(datx) << std::endl; + iunit << " # current working dir "+rtrim(cwd) << std::endl; + + iunit << " # Program compiled at "; + iunit.width(50); + if (compile_date.size() > 50) {compile_date.resize(50);} + iunit << std::right << compile_date; + + iunit << " # Program compiled on "; + iunit.width(50); + if (compile_host.size() > 50) {compile_host.resize(50);} + iunit << std::right << compile_host; + + iunit << " # Program compiled for "; + iunit.width(50); + if (compile_arch.size() > 50) {compile_arch.resize(50);} + iunit << std::right << compile_arch; + + iunit << " # Source code revision number"; + iunit.width(50); + iunit << std::right << compile_revision; + +} diff --git a/src/header.cpp b/src/header.cpp new file mode 100644 index 0000000..67c7535 --- /dev/null +++ b/src/header.cpp @@ -0,0 +1,14 @@ +//------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2022 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//------------------------------------------------------------------------------------------------// + +// ************************************************************************************************* +//> \par History +//> none +//> \author APSI & CJM & JGH +// ************************************************************************************************* + +#include "base/base_uses.cpp" diff --git a/src/motion/simpar_methods.c b/src/motion/simpar_methods.c new file mode 100644 index 0000000..1d966d8 --- /dev/null +++ b/src/motion/simpar_methods.c @@ -0,0 +1,62 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Methods for storing MD parameters type + * \author CJM + * \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 + * reorganization of the original routines/modules + ******************************************************************************/ + +#include "../base/base_uses.h" + +/******************************************************************************* + * \brief Reads the MD section and setup the simulation parameters type + * \param simpar ... + * \param motion_section ... + * \param md_section ... + * \author Teodoro Laino + ******************************************************************************/ + +void read_md_section(struct simpar_type *simpar, struct section_vals_type *motion_section, + struct section_vals_type *md_section) { + + char* filename[default_path_length]; + int iprint, iw; + double tmp_r1, tmp_r2, tmp_r3; + struct cp_logger_type *logger; + struct enumeration_type *enumer; + struct keyword_type *keyword; + struct section_type *section; + struct section_vals_type *print_key; + + // free(logger,print_key,enumer,keyword,section); + + logger = &cp_get_default_logger(); + iw = cp_print_key_uint_nr(logger, md_section, "PRINT.PROGRAM_RUN_INFO", ".log"); + + read_md_low(simpar, motion_section, md_section); + if (iw > 0) + printf("%i",iw); + + // Begin setup Langevin dynamics + if (simpar.ensemble == langevin_ensemble) { + cite_reference(Ricci2003); + if (simpar.noisy_gamma > 0.0) + cite_reference(Kuhne2007); + + if (simpar.shadow_gamma > 0.0) + cite_reference(Rengaraj2020); + + // Normalization factor using a normal Gaussian random number distribution + simpar.var_w = 2.0 * simpar.temp_ext * simpar.dt * (simpar.gamma + simpar.noisy_gamma); + if (iw > 0) { + + } + } + +} \ No newline at end of file diff --git a/src/mpiwrap/PACKAGE b/src/mpiwrap/PACKAGE new file mode 100644 index 0000000..7ddf13d --- /dev/null +++ b/src/mpiwrap/PACKAGE @@ -0,0 +1,5 @@ +{ + "description": "wrappers of the mpi routines", + "requires": ["../base"], + "implicit": "MPI_.*", +} \ No newline at end of file diff --git a/src/mpiwrap/message_passing.c b/src/mpiwrap/message_passing.c new file mode 100644 index 0000000..e69de29 diff --git a/src/mpiwrap/message_passing.h b/src/mpiwrap/message_passing.h new file mode 100644 index 0000000..46a7b29 --- /dev/null +++ b/src/mpiwrap/message_passing.h @@ -0,0 +1,60 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Interface to the message passing library MPI + * \par History + * JGH (02-Jan-2001): New error handling + * Performance tools + * JGH (14-Jan-2001): New routines mp_comm_compare, mp_cart_coords, + * mp_rank_compare, mp_alltoall + * JGH (06-Feb-2001): New routines mp_comm_free + * JGH (22-Mar-2001): New routines mp_comm_dup + * fawzi (04-NOV-2004): storable performance info (for f77 interface) + * Wrapper routine for mpi_gatherv added (22.12.2005,MK) + * JGH (13-Feb-2006): Flexible precision + * JGH (15-Feb-2006): single precision mp_alltoall + * \author JGH +*******************************************************************************/ + +#include + +#include "../base/base_uses.h" + +#if defined(__parallel) + USE mpi +//subroutines: unfortunately, mpi implementations do not provide interfaces for all subroutines (problems with types and ranks explosion), +// we do not quite know what is in the module, so we can not include any.... +// to nevertheless get checking for what is included, we use the mpi module without use clause, getting all there is +//USE mpi, ONLY: mpi_allgather, mpi_allgatherv, mpi_alloc_mem, mpi_allreduce, mpi_alltoall, mpi_alltoallv, mpi_bcast,& +// mpi_cart_coords, mpi_cart_create, mpi_cart_get, mpi_cart_rank, mpi_cart_sub, mpi_dims_create, mpi_file_close,& +// mpi_file_get_size, mpi_file_open, mpi_file_read_at_all, mpi_file_read_at, mpi_file_write_at_all,& +// mpi_file_write_at, mpi_free_mem, mpi_gather, mpi_gatherv, mpi_get_address, mpi_group_translate_ranks, mpi_irecv,& +// mpi_isend, mpi_recv, mpi_reduce, mpi_reduce_scatter, mpi_rget, mpi_scatter, mpi_send,& +// mpi_sendrecv, mpi_sendrecv_replace, mpi_testany, mpi_waitall, mpi_waitany, mpi_win_create +//functions +//USE mpi, ONLY: mpi_wtime +//constants +//USE mpi, ONLY: MPI_DOUBLE_PRECISION, MPI_DOUBLE_COMPLEX, MPI_REAL, MPI_COMPLEX, MPI_ANY_TAG,& +// MPI_ANY_SOURCE, MPI_COMM_NULL, MPI_REQUEST_NULL, MPI_WIN_NULL, MPI_STATUS_SIZE, MPI_STATUS_IGNORE, MPI_STATUSES_IGNORE, & +// MPI_ADDRESS_KIND, MPI_OFFSET_KIND, MPI_MODE_CREATE, MPI_MODE_RDONLY, MPI_MODE_WRONLY,& +// MPI_MODE_RDWR, MPI_MODE_EXCL, MPI_COMM_SELF, MPI_COMM_WORLD, MPI_THREAD_SERIALIZED,& +// MPI_ERRORS_RETURN, MPI_SUCCESS, MPI_MAX_PROCESSOR_NAME, MPI_MAX_ERROR_STRING, MPI_IDENT,& +// MPI_UNEQUAL, MPI_MAX, MPI_SUM, MPI_INFO_NULL, MPI_IN_PLACE, MPI_CONGRUENT, MPI_SIMILAR, MPI_MIN, MPI_SOURCE,& +// MPI_TAG, MPI_INTEGER8, MPI_INTEGER, MPI_MAXLOC, MPI_2INTEGER, MPI_MINLOC, MPI_LOGICAL, MPI_2DOUBLE_PRECISION,& +// MPI_LOR, MPI_CHARACTER, MPI_BOTTOM, MPI_MODE_NOCHECK, MPI_2REAL +#endif + + // parameters that might be needed +#if defined(__parallel) + const int MP_STD_REAL = MPI_DOUBLE_PRECISION; + const int MP_STD_COMPLEX = MPI_DOUBLE_COMPLEX; + const int MP_STD_HALF_REAL = MPI_REAL; + const int MP_STD_HALF_COMPLEX = MPI_COMPLEX; + + const bool cp2k_is_parallel = true; + diff --git a/src/qs_environment_types.cpp b/src/qs_environment_types.cpp new file mode 100644 index 0000000..deaad7a --- /dev/null +++ b/src/qs_environment_types.cpp @@ -0,0 +1,1675 @@ +/*-------------------------------------------------------------------------------------------------- + CP2K: A general program to perform molecular dynamics simulations + Copyright 2000-2023 CP2K developers group + + SPDX-License-Identifier: GPL-2.0-or-later +--------------------------------------------------------------------------------------------------*/ + +/*************************************************************************************************** + \par History + - mo_set_p_type added to qs_env (23.04.02,MK) + - qs_force_type added to qs_env (05.06.02,MK) + \author MK (23.01.2002) +***************************************************************************************************/ + +#include "./base/base_uses.h" + +class qs_environment_types { + +private: + char *moduleN = "qs_environment_types"; + +public: +// *** Public data types *** + + struct qs_environment_type; + +// *** Public subroutines *** + void get_qs_env(); + void qs_env_create(); + void qs_env_release(); + void qs_env_part_release(); + void set_qs_env(); + +/*************************************************************************************************** + \param local_rho_set contains the atomic, compensations and core densities + and the local parts of the xc terms + \param hartree_local contains the 1, 2 and 3 centers coulomb terms + \param requires_mo_derivs logical, true if dE/dC is required (e.g. OT) + \param has_unit_metric logical, true if the S matrix is considered unity for the SCF + \param mo_derivs the actual derivatives of the total energy wrt to MO coeffs (divided by +2*f_i) + \param xas_env temporary information for xas calculation + \param dftb_potential pair potentials for use with DFTB + \param dispersion_env environment for use with QS dispersion + + compatibility get (things that you should get from the subsys): + \param atomic_kind_set array with infos about the species (atomic_kinds) + present in the system + \param particle_set info on the atoms you simulate, pos,... + \param local_particles which particles ar local to this processor + new: + \param local_molecules which molecules are local to this processor + \param molecule_kind_set description of the molecule kinds + \param molecule_set all the molecule description + \param rtp all data needed for real time propagation + \param x contains data used in Hartree-Fock-Exchange calculations + \param task_list the list of tasks used in collocate and integrate + \param task_list_soft the list of tasks used in collocate and integrate in case of soft basis +functions + \param mo_loc_history if a history of localized wfn is kept, they are stored here. + \param molecular_scf_guess_env contains inforamation about and results of claculations + on separate molecules + \par History + 11.2002 added doc and attribute description [fawzi] + 08.2004 renamed some of the very short names (s,c,k,h) for easier grepping + 06.2018 polar_env added (MK) + \author Matthias Krack & fawzi +***************************************************************************************************/ + + struct qs_environment_type { + bool qmmm, qmmm_periodic; + bool requires_mo_derivs; + bool requires_matrix_vxc; + bool has_unit_metric; + bool run_rtp; + bool linres_run; + bool calc_image_preconditioner; + bool do_transport; + bool single_point_run; + bool given_embed_pot; + bool energy_correction; + double sim_time; + double start_time, target_time; + double *image_matrix; + double *image_coeff; + int *ipiv; + int sim_step; + TYPE(ls_scf_env_type), POINTER :: ls_scf_env + TYPE(almo_scf_env_type), POINTER :: almo_scf_env + TYPE(transport_env_type), POINTER :: transport_env + TYPE(cell_type), POINTER :: super_cell + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_loc_history + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs + TYPE(scf_control_type), POINTER :: scf_control + TYPE(rel_control_type), POINTER :: rel_control + // ZMP adding variables + TYPE(qs_rho_type), POINTER :: rho_external + TYPE(pw_type), POINTER :: external_vxc + TYPE(pw_type), POINTER :: mask + TYPE(qs_charges_type), POINTER :: qs_charges + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env + TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm + TYPE(qs_wf_history_type), POINTER :: wf_history + TYPE(qs_scf_env_type), POINTER :: scf_env + TYPE(qs_matrix_pools_type), POINTER :: mpools + TYPE(oce_matrix_type), POINTER :: oce + TYPE(local_rho_type), POINTER :: local_rho_set + TYPE(hartree_local_type), POINTER :: hartree_local + TYPE(section_vals_type), POINTER :: input + TYPE(linres_control_type), POINTER :: linres_control + TYPE(xas_environment_type), POINTER :: xas_env + TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env + TYPE(cp_ddapc_ewald_type), POINTER :: cp_ddapc_ewald + REAL(KIND=dp), DIMENSION(:, :), POINTER :: outer_scf_history + INTEGER :: outer_scf_ihistory + REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, & + variable_history + TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data + TYPE(et_coupling_type), POINTER :: et_coupling + TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), POINTER :: dftb_potential + TYPE(admm_type), POINTER :: admm_env + TYPE(active_space_type), POINTER :: active_space + // LRI + TYPE(lri_environment_type), POINTER :: lri_env + TYPE(lri_density_type), POINTER :: lri_density + // Energy correction + TYPE(energy_correction_type), POINTER :: ec_env + // Excited States + bool excited_state; + TYPE(excited_energy_type), POINTER :: exstate_env + // Empirical dispersion + TYPE(qs_dispersion_type), POINTER :: dispersion_env + // Empirical geometrical BSSE correction + TYPE(qs_gcp_type), POINTER :: gcp_env + // Semi-empirical and DFTB types + TYPE(ewald_environment_type), POINTER :: ewald_env + TYPE(ewald_pw_type), POINTER :: ewald_pw + // Semi-empirical types + TYPE(se_taper_type), POINTER :: se_taper + TYPE(semi_empirical_si_type), POINTER :: se_store_int_env + TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole + TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env + TYPE(rt_prop_type), POINTER :: rtp + TYPE(efield_berry_type), POINTER :: efield + // a history for the broyden ot + double broyden_adaptive_sigma + TYPE(mp2_type), POINTER :: mp2_env => NULL() + TYPE(post_scf_bandstructure_type), POINTER :: bs_env => +NULL() + TYPE(kg_environment_type), POINTER :: kg_env + TYPE(wannier_centres_type), POINTER, DIMENSION(:) :: WannierCentres => NULL() + TYPE(molecular_scf_guess_env_type), POINTER :: molecular_scf_guess_env => NULL() + // Subsystem densities + TYPE(qs_rho_p_type), DIMENSION(:), POINTER :: subsys_dens + // Embedding potential + TYPE(pw_type), POINTER :: embed_pot + TYPE(pw_type), POINTER :: spin_embed_pot + // Polarizability tensor + TYPE(polar_env_type), POINTER :: polar_env + // Resp charges + REAL(KIND=dp), DIMENSION(:), POINTER :: rhs => NULL() + double total_zeff_corr, surface_dipole_moment + bool surface_dipole_switch_off; + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_last_converged + } + +} + +/*************************************************************************************************** + \brief Get the QUICKSTEP environment. + \param qs_env ... + \param atomic_kind_set ... + \param qs_kind_set ... + \param cell ... + \param super_cell ... + \param cell_ref ... + \param use_ref_cell ... + \param kpoints ... + \param dft_control ... + \param mos ... + \param sab_orb ... + \param sab_all ... + \param qmmm ... + \param qmmm_periodic ... + \param sac_ae ... + \param sac_ppl ... + \param sac_lri ... + \param sap_ppnl ... + \param sab_vdw ... + \param sab_scp ... + \param sap_oce ... + \param sab_lrc ... + \param sab_se ... + \param sab_xtbe ... + \param sab_tbe ... + \param sab_core ... + \param sab_xb ... + \param sab_xtb_nonbond ... + \param sab_almo ... + \param sab_kp ... + \param sab_kp_nosym ... + \param particle_set ... + \param energy ... + \param force ... + \param matrix_h ... + \param matrix_h_im ... + \param matrix_ks ... + \param matrix_ks_im ... + \param matrix_vxc ... + \param run_rtp ... + \param rtp ... + \param matrix_h_kp ... + \param matrix_h_im_kp ... + \param matrix_ks_kp ... + \param matrix_ks_im_kp ... + \param matrix_vxc_kp ... + \param kinetic_kp ... + \param matrix_s_kp ... + \param matrix_w_kp ... + \param matrix_s_RI_aux_kp ... + \param matrix_s ... + \param matrix_s_RI_aux ... + \param matrix_w ... + \param matrix_p_mp2 ... + \param matrix_p_mp2_admm ... + \param rho ... + \param rho_xc ... + \param pw_env ... + \param ewald_env ... + \param ewald_pw ... + \param active_space ... + \param mpools ... + \param input ... + \param para_env ... + \param blacs_env ... + \param scf_control ... + \param rel_control ... + \param kinetic ... + \param qs_charges ... + \param vppl ... + \param rho_core ... + \param rho_nlcc ... + \param rho_nlcc_g ... + \param ks_env ... + \param ks_qmmm_env ... + \param wf_history ... + \param scf_env ... + \param local_particles ... + \param local_molecules ... + \param distribution_2d ... + \param dbcsr_dist ... + \param molecule_kind_set ... + \param molecule_set ... + \param subsys ... + \param cp_subsys ... + \param oce ... + \param local_rho_set ... + \param rho_atom_set ... + \param task_list ... + \param task_list_soft ... + \param rho0_atom_set ... + \param rho0_mpole ... + \param rhoz_set ... + \param ecoul_1c ... + \param rho0_s_rs ... + \param rho0_s_gs ... + \param do_kpoints ... + \param has_unit_metric ... + \param requires_mo_derivs ... + \param mo_derivs ... + \param mo_loc_history ... + \param nkind ... + \param natom ... + \param nelectron_total ... + \param nelectron_spin ... + \param efield ... + \param neighbor_list_id ... + \param linres_control ... + \param xas_env ... + \param virial ... + \param cp_ddapc_env ... + \param cp_ddapc_ewald ... + \param outer_scf_history ... + \param outer_scf_ihistory ... + \param x_data ... + \param et_coupling ... + \param dftb_potential ... + \param results ... + \param se_taper ... + \param se_store_int_env ... + \param se_nddo_mpole ... + \param se_nonbond_env ... + \param admm_env ... + \param lri_env ... + \param lri_density ... + \param exstate_env ... + \param ec_env ... + \param dispersion_env ... + \param gcp_env ... + \param vee ... + \param rho_external ... + \param external_vxc ... + \param mask ... + \param mp2_env ... + \param bs_env ... + \param kg_env ... + \param WannierCentres ... + \param atprop ... + \param ls_scf_env ... + \param do_transport ... + \param transport_env ... + \param v_hartree_rspace ... + \param s_mstruct_changed ... + \param rho_changed ... + \param potential_changed ... + \param forces_up_to_date ... + \param mscfg_env ... + \param almo_scf_env ... + \param gradient_history ... + \param variable_history ... + \param embed_pot ... + \param spin_embed_pot ... + \param polar_env ... + \param mos_last_converged ... [SGh] + \param rhs ... + \date 23.01.2002 + \author MK + \version 1.0 +***************************************************************************************************/ +void qs_environment_types::get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, +cell_ref, use_ref_cell, kpoints, & + dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, +sac_ppl, sac_lri, & + sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, +sab_core, & + sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, +energy, force, & + matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, +rtp, & + matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, +matrix_vxc_kp, kinetic_kp, matrix_s_kp, & + matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, & + matrix_p_mp2, matrix_p_mp2_admm, rho, & + rho_xc, pw_env, ewald_env, ewald_pw, active_space, & + mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, +qs_charges, & + vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, +scf_env, local_particles, & + local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, & + molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, & + task_list, & + task_list_soft, & + rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, & + rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, +mo_derivs, & + mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, & + neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, +cp_ddapc_ewald, & + outer_scf_history, outer_scf_ihistory, x_data, et_coupling, +dftb_potential, results, & + se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, & + lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, +& + rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, & + WannierCentres, atprop, ls_scf_env, do_transport, transport_env, +v_hartree_rspace, & + s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, +mscfg_env, almo_scf_env, & + gradient_history, variable_history, embed_pot, spin_embed_pot, +polar_env, mos_last_converged, rhs) + TYPE(qs_environment_type), INTENT(IN) :: qs_env + TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, & + POINTER :: atomic_kind_set + TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, & + POINTER :: qs_kind_set + TYPE(cell_type), OPTIONAL, POINTER :: cell, super_cell, cell_ref + LOGICAL, OPTIONAL :: use_ref_cell + TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints + TYPE(dft_control_type), OPTIONAL, POINTER :: dft_control + TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + OPTIONAL, POINTER :: sab_orb, sab_all + LOGICAL, OPTIONAL :: qmmm, qmmm_periodic + TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sac_ae, sac_ppl, +sac_lri, & + sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, & + sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym + TYPE(particle_type), DIMENSION(:), OPTIONAL, & + POINTER :: particle_set + TYPE(qs_energy_type), OPTIONAL, POINTER :: energy + TYPE(qs_force_type), DIMENSION(:), OPTIONAL, & + POINTER :: force + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: matrix_h, matrix_h_im, matrix_ks, & + matrix_ks_im, matrix_vxc + LOGICAL, OPTIONAL :: run_rtp + TYPE(rt_prop_type), OPTIONAL, POINTER :: rtp + TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, POINTER :: matrix_h_kp, matrix_h_im_kp, & + matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, & + matrix_s_RI_aux_kp + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: matrix_s, matrix_s_RI_aux, matrix_w, +& + matrix_p_mp2, matrix_p_mp2_admm + TYPE(qs_rho_type), OPTIONAL, POINTER :: rho, rho_xc + TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env + TYPE(ewald_environment_type), OPTIONAL, POINTER :: ewald_env + TYPE(ewald_pw_type), OPTIONAL, POINTER :: ewald_pw + TYPE(active_space_type), OPTIONAL, POINTER :: active_space + TYPE(qs_matrix_pools_type), OPTIONAL, POINTER :: mpools + TYPE(section_vals_type), OPTIONAL, POINTER :: input + TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env + TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env + TYPE(scf_control_type), OPTIONAL, POINTER :: scf_control + TYPE(rel_control_type), OPTIONAL, POINTER :: rel_control + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: kinetic + TYPE(qs_charges_type), OPTIONAL, POINTER :: qs_charges + TYPE(pw_type), OPTIONAL, POINTER :: vppl, rho_core, rho_nlcc, rho_nlcc_g + TYPE(qs_ks_env_type), OPTIONAL, POINTER :: ks_env + TYPE(qs_ks_qmmm_env_type), OPTIONAL, POINTER :: ks_qmmm_env + TYPE(qs_wf_history_type), OPTIONAL, POINTER :: wf_history + TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env + TYPE(distribution_1d_type), OPTIONAL, POINTER :: local_particles, local_molecules + TYPE(distribution_2d_type), OPTIONAL, POINTER :: distribution_2d + TYPE(dbcsr_distribution_type), OPTIONAL, POINTER :: dbcsr_dist + TYPE(molecule_kind_type), DIMENSION(:), OPTIONAL, & + POINTER :: molecule_kind_set + TYPE(molecule_type), DIMENSION(:), OPTIONAL, & + POINTER :: molecule_set + TYPE(qs_subsys_type), OPTIONAL, POINTER :: subsys + TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys + TYPE(oce_matrix_type), OPTIONAL, POINTER :: oce + TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set + TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, & + POINTER :: rho_atom_set + TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_soft + TYPE(rho0_atom_type), DIMENSION(:), OPTIONAL, & + POINTER :: rho0_atom_set + TYPE(rho0_mpole_type), OPTIONAL, POINTER :: rho0_mpole + TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER :: rhoz_set + TYPE(ecoul_1center_type), DIMENSION(:), OPTIONAL, & + POINTER :: ecoul_1c + TYPE(pw_type), OPTIONAL, POINTER :: rho0_s_rs, rho0_s_gs + LOGICAL, OPTIONAL :: do_kpoints, has_unit_metric, & + requires_mo_derivs + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: mo_derivs + TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history + INTEGER, OPTIONAL :: nkind, natom, nelectron_total + INTEGER, DIMENSION(2), OPTIONAL :: nelectron_spin + TYPE(efield_berry_type), OPTIONAL, POINTER :: efield + INTEGER, OPTIONAL :: neighbor_list_id + TYPE(linres_control_type), OPTIONAL, POINTER :: linres_control + TYPE(xas_environment_type), OPTIONAL, POINTER :: xas_env + TYPE(virial_type), OPTIONAL, POINTER :: virial + TYPE(cp_ddapc_type), OPTIONAL, POINTER :: cp_ddapc_env + TYPE(cp_ddapc_ewald_type), OPTIONAL, POINTER :: cp_ddapc_ewald + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: outer_scf_history + INTEGER, INTENT(out), OPTIONAL :: outer_scf_ihistory + TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, POINTER :: x_data + TYPE(et_coupling_type), OPTIONAL, POINTER :: et_coupling + TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & + OPTIONAL, POINTER :: dftb_potential + TYPE(cp_result_type), OPTIONAL, POINTER :: results + TYPE(se_taper_type), OPTIONAL, POINTER :: se_taper + TYPE(semi_empirical_si_type), OPTIONAL, POINTER :: se_store_int_env + TYPE(nddo_mpole_type), OPTIONAL, POINTER :: se_nddo_mpole + TYPE(fist_nonbond_env_type), OPTIONAL, POINTER :: se_nonbond_env + TYPE(admm_type), OPTIONAL, POINTER :: admm_env + TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env + TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density + TYPE(excited_energy_type), OPTIONAL, POINTER :: exstate_env + TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env + TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env + TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env + TYPE(pw_type), OPTIONAL, POINTER :: vee + TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_external + TYPE(pw_type), OPTIONAL, POINTER :: external_vxc, mask + TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env + TYPE(post_scf_bandstructure_type), OPTIONAL, & + POINTER :: bs_env + TYPE(kg_environment_type), OPTIONAL, POINTER :: kg_env + TYPE(wannier_centres_type), DIMENSION(:), & + OPTIONAL, POINTER :: WannierCentres + TYPE(atprop_type), OPTIONAL, POINTER :: atprop + TYPE(ls_scf_env_type), OPTIONAL, POINTER :: ls_scf_env + LOGICAL, OPTIONAL :: do_transport + TYPE(transport_env_type), OPTIONAL, POINTER :: transport_env + TYPE(pw_type), OPTIONAL, POINTER :: v_hartree_rspace + LOGICAL, OPTIONAL :: s_mstruct_changed, rho_changed, & + potential_changed, forces_up_to_date + TYPE(molecular_scf_guess_env_type), OPTIONAL, & + POINTER :: mscfg_env + TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history + TYPE(pw_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot + TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env + TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged + REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rhs + + TYPE(rho0_mpole_type), POINTER :: rho0_m + + NULLIFY (rho0_m) + CPASSERT(ASSOCIATED(qs_env%ks_env)) + + IF (PRESENT(outer_scf_history)) outer_scf_history => qs_env%outer_scf_history + IF (PRESENT(outer_scf_ihistory)) outer_scf_ihistory = qs_env%outer_scf_ihistory + IF (PRESENT(gradient_history)) gradient_history => qs_env%gradient_history + IF (PRESENT(variable_history)) variable_history => qs_env%variable_history + IF (PRESENT(mp2_env)) mp2_env => qs_env%mp2_env + IF (PRESENT(bs_env)) bs_env => qs_env%bs_env + IF (PRESENT(kg_env)) kg_env => qs_env%kg_env + IF (PRESENT(super_cell)) super_cell => qs_env%super_cell + IF (PRESENT(qmmm)) qmmm = qs_env%qmmm + IF (PRESENT(qmmm_periodic)) qmmm_periodic = qs_env%qmmm_periodic + IF (PRESENT(mos)) mos => qs_env%mos + IF (PRESENT(mos_last_converged)) mos_last_converged => qs_env%mos_last_converged + IF (PRESENT(ewald_env)) ewald_env => qs_env%ewald_env + IF (PRESENT(ewald_pw)) ewald_pw => qs_env%ewald_pw + IF (PRESENT(mpools)) mpools => qs_env%mpools + IF (PRESENT(scf_control)) scf_control => qs_env%scf_control + IF (PRESENT(rel_control)) rel_control => qs_env%rel_control + // ZMP pointing vectors + IF (PRESENT(rho_external)) rho_external => qs_env%rho_external + IF (PRESENT(external_vxc)) external_vxc => qs_env%external_vxc + IF (PRESENT(mask)) mask => qs_env%mask + IF (PRESENT(qs_charges)) qs_charges => qs_env%qs_charges + IF (PRESENT(ks_env)) ks_env => qs_env%ks_env + IF (PRESENT(ks_qmmm_env)) ks_qmmm_env => qs_env%ks_qmmm_env + IF (PRESENT(wf_history)) wf_history => qs_env%wf_history + IF (PRESENT(scf_env)) scf_env => qs_env%scf_env + IF (PRESENT(oce)) oce => qs_env%oce + IF (PRESENT(requires_mo_derivs)) requires_mo_derivs = qs_env%requires_mo_derivs + IF (PRESENT(has_unit_metric)) has_unit_metric = qs_env%has_unit_metric + IF (PRESENT(mo_derivs)) mo_derivs => qs_env%mo_derivs + IF (PRESENT(mo_loc_history)) mo_loc_history => qs_env%mo_loc_history + IF (PRESENT(linres_control)) linres_control => qs_env%linres_control + IF (PRESENT(se_taper)) se_taper => qs_env%se_taper + IF (PRESENT(se_store_int_env)) se_store_int_env => qs_env%se_store_int_env + IF (PRESENT(se_nddo_mpole)) se_nddo_mpole => qs_env%se_nddo_mpole + IF (PRESENT(se_nonbond_env)) se_nonbond_env => qs_env%se_nonbond_env + IF (PRESENT(lri_env)) lri_env => qs_env%lri_env + IF (PRESENT(lri_density)) lri_density => qs_env%lri_density + IF (PRESENT(ec_env)) ec_env => qs_env%ec_env + IF (PRESENT(exstate_env)) exstate_env => qs_env%exstate_env + IF (PRESENT(dispersion_env)) dispersion_env => qs_env%dispersion_env + IF (PRESENT(gcp_env)) gcp_env => qs_env%gcp_env + IF (PRESENT(run_rtp)) run_rtp = qs_env%run_rtp + IF (PRESENT(rtp)) rtp => qs_env%rtp + IF (PRESENT(ls_scf_env)) ls_scf_env => qs_env%ls_scf_env + IF (PRESENT(almo_scf_env)) almo_scf_env => qs_env%almo_scf_env + IF (PRESENT(do_transport)) do_transport = qs_env%do_transport + IF (PRESENT(transport_env)) transport_env => qs_env%transport_env + IF (PRESENT(mscfg_env)) mscfg_env => qs_env%molecular_scf_guess_env + IF (PRESENT(active_space)) active_space => qs_env%active_space + IF (PRESENT(admm_env)) admm_env => qs_env%admm_env + + // Embedding potential + IF (PRESENT(embed_pot)) embed_pot => qs_env%embed_pot + IF (PRESENT(spin_embed_pot)) spin_embed_pot => qs_env%spin_embed_pot + + // Polarisability tensor + IF (PRESENT(polar_env)) polar_env => qs_env%polar_env + + // Resp charges + IF (PRESENT(rhs)) rhs => qs_env%rhs + + IF (PRESENT(local_rho_set)) & + local_rho_set => qs_env%local_rho_set + IF (PRESENT(rho_atom_set)) & + CALL get_local_rho(qs_env%local_rho_set, rho_atom_set=rho_atom_set) + IF (PRESENT(rho0_atom_set)) & + CALL get_local_rho(qs_env%local_rho_set, rho0_atom_set=rho0_atom_set) + IF (PRESENT(rho0_mpole)) & + CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_mpole) + IF (PRESENT(rhoz_set)) & + CALL get_local_rho(qs_env%local_rho_set, rhoz_set=rhoz_set) + IF (PRESENT(ecoul_1c)) & + CALL get_hartree_local(qs_env%hartree_local, ecoul_1c=ecoul_1c) + IF (PRESENT(rho0_s_rs)) THEN + CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_m) + IF (ASSOCIATED(rho0_m)) THEN + rho0_s_rs => rho0_m%rho0_s_rs + END IF + END IF + IF (PRESENT(rho0_s_gs)) THEN + CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_m) + IF (ASSOCIATED(rho0_m)) THEN + rho0_s_gs => rho0_m%rho0_s_gs + END IF + END IF + + IF (PRESENT(xas_env)) xas_env => qs_env%xas_env + IF (PRESENT(input)) input => qs_env%input + IF (PRESENT(cp_ddapc_env)) cp_ddapc_env => qs_env%cp_ddapc_env + IF (PRESENT(cp_ddapc_ewald)) cp_ddapc_ewald => qs_env%cp_ddapc_ewald + IF (PRESENT(x_data)) x_data => qs_env%x_data + IF (PRESENT(et_coupling)) et_coupling => qs_env%et_coupling + IF (PRESENT(dftb_potential)) dftb_potential => qs_env%dftb_potential + IF (PRESENT(efield)) efield => qs_env%efield + IF (PRESENT(WannierCentres)) WannierCentres => qs_env%WannierCentres + + CALL get_ks_env(qs_env%ks_env, & + v_hartree_rspace=v_hartree_rspace, & + s_mstruct_changed=s_mstruct_changed, & + rho_changed=rho_changed, & + potential_changed=potential_changed, & + forces_up_to_date=forces_up_to_date, & + matrix_h=matrix_h, & + matrix_h_im=matrix_h_im, & + matrix_ks=matrix_ks, & + matrix_ks_im=matrix_ks_im, & + matrix_vxc=matrix_vxc, & + kinetic=kinetic, & + matrix_s=matrix_s, & + matrix_s_RI_aux=matrix_s_RI_aux, & + matrix_ks_im_kp=matrix_ks_im_kp, & + matrix_w=matrix_w, & + matrix_p_mp2=matrix_p_mp2, & + matrix_p_mp2_admm=matrix_p_mp2_admm, & + matrix_h_kp=matrix_h_kp, & + matrix_h_im_kp=matrix_h_im_kp, & + matrix_ks_kp=matrix_ks_kp, & + matrix_vxc_kp=matrix_vxc_kp, & + kinetic_kp=kinetic_kp, & + matrix_s_kp=matrix_s_kp, & + matrix_w_kp=matrix_w_kp, & + matrix_s_RI_aux_kp=matrix_s_RI_aux_kp, & + rho=rho, & + rho_xc=rho_xc, & + rho_core=rho_core, & + rho_nlcc=rho_nlcc, & + rho_nlcc_g=rho_nlcc_g, & + vppl=vppl, & + vee=vee, & + neighbor_list_id=neighbor_list_id, & + sab_orb=sab_orb, & + sab_all=sab_all, & + sab_scp=sab_scp, & + sab_vdw=sab_vdw, & + sac_ae=sac_ae, & + sac_ppl=sac_ppl, & + sac_lri=sac_lri, & + sap_ppnl=sap_ppnl, & + sap_oce=sap_oce, & + sab_se=sab_se, & + sab_lrc=sab_lrc, & + sab_tbe=sab_tbe, & + sab_xtbe=sab_xtbe, & + sab_core=sab_core, & + sab_xb=sab_xb, & + sab_xtb_nonbond=sab_xtb_nonbond, & + sab_almo=sab_almo, & + sab_kp=sab_kp, & + sab_kp_nosym=sab_kp_nosym, & + task_list=task_list, & + task_list_soft=task_list_soft, & + kpoints=kpoints, & + do_kpoints=do_kpoints, & + local_molecules=local_molecules, & + local_particles=local_particles, & + atprop=atprop, & + virial=virial, & + results=results, & + cell=cell, & + cell_ref=cell_ref, & + use_ref_cell=use_ref_cell, & + energy=energy, & + force=force, & + qs_kind_set=qs_kind_set, & + subsys=subsys, & + cp_subsys=cp_subsys, & + atomic_kind_set=atomic_kind_set, & + particle_set=particle_set, & + molecule_kind_set=molecule_kind_set, & + molecule_set=molecule_set, & + natom=natom, & + nkind=nkind, & + dft_control=dft_control, & + dbcsr_dist=dbcsr_dist, & + distribution_2d=distribution_2d, & + pw_env=pw_env, & + para_env=para_env, & + blacs_env=blacs_env, & + nelectron_total=nelectron_total, & + nelectron_spin=nelectron_spin) + + END SUBROUTINE get_qs_env + +/*************************************************************************************************** + \brief Initialise the QUICKSTEP environment. + \param qs_env ... + \param globenv ... + \date 25.01.2002 + \author MK + \version 1.0 +***************************************************************************************************/ +void init_qs_env(qs_env, globenv) + + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(global_environment_type), OPTIONAL, POINTER :: globenv + + NULLIFY (qs_env%ls_scf_env) + NULLIFY (qs_env%almo_scf_env) + NULLIFY (qs_env%transport_env) + NULLIFY (qs_env%image_matrix) + NULLIFY (qs_env%ipiv) + NULLIFY (qs_env%image_coeff) + NULLIFY (qs_env%super_cell) + NULLIFY (qs_env%mos) + NULLIFY (qs_env%mos_last_converged) + NULLIFY (qs_env%mpools) + NULLIFY (qs_env%ewald_env) + NULLIFY (qs_env%ewald_pw) + NULLIFY (qs_env%scf_control) + NULLIFY (qs_env%rel_control) + NULLIFY (qs_env%qs_charges) + // ZMP initializing arrays + NULLIFY (qs_env%rho_external) + NULLIFY (qs_env%external_vxc) + NULLIFY (qs_env%mask) + // Embedding potential + NULLIFY (qs_env%embed_pot) + NULLIFY (qs_env%spin_embed_pot) + + // Polarisability tensor + NULLIFY (qs_env%polar_env) + + NULLIFY (qs_env%ks_env) + NULLIFY (qs_env%ks_qmmm_env) + NULLIFY (qs_env%wf_history) + NULLIFY (qs_env%scf_env) + NULLIFY (qs_env%oce) + NULLIFY (qs_env%local_rho_set) + NULLIFY (qs_env%hartree_local) + NULLIFY (qs_env%input) + NULLIFY (qs_env%linres_control) + NULLIFY (qs_env%xas_env) + NULLIFY (qs_env%cp_ddapc_env) + NULLIFY (qs_env%cp_ddapc_ewald) + NULLIFY (qs_env%outer_scf_history) + NULLIFY (qs_env%gradient_history) + NULLIFY (qs_env%variable_history) + NULLIFY (qs_env%x_data) + NULLIFY (qs_env%et_coupling) + NULLIFY (qs_env%dftb_potential) + NULLIFY (qs_env%active_space) + + NULLIFY (qs_env%se_taper) + NULLIFY (qs_env%se_store_int_env) + NULLIFY (qs_env%se_nddo_mpole) + NULLIFY (qs_env%se_nonbond_env) + NULLIFY (qs_env%admm_env) + NULLIFY (qs_env%efield) + NULLIFY (qs_env%lri_env) + NULLIFY (qs_env%ec_env) + NULLIFY (qs_env%exstate_env) + NULLIFY (qs_env%lri_density) + NULLIFY (qs_env%gcp_env) + NULLIFY (qs_env%rtp) + NULLIFY (qs_env%mp2_env) + NULLIFY (qs_env%bs_env) + NULLIFY (qs_env%kg_env) + NULLIFY (qs_env%ec_env) + NULLIFY (qs_env%WannierCentres) + + qs_env%outer_scf_ihistory = 0 + qs_env%broyden_adaptive_sigma = -1.0_dp + + CALL local_rho_set_create(qs_env%local_rho_set) + CALL hartree_local_create(qs_env%hartree_local) + qs_env%run_rtp = .FALSE. + qs_env%linres_run = .FALSE. + qs_env%single_point_run = .FALSE. + qs_env%qmmm = .FALSE. + qs_env%qmmm_periodic = .FALSE. + qs_env%requires_mo_derivs = .FALSE. + qs_env%requires_matrix_vxc = .FALSE. + qs_env%has_unit_metric = .FALSE. + qs_env%calc_image_preconditioner = .TRUE. + qs_env%do_transport = .FALSE. + qs_env%given_embed_pot = .FALSE. + IF (PRESENT(globenv)) THEN + qs_env%target_time = globenv%cp2k_target_time + qs_env%start_time = globenv%cp2k_start_time + qs_env%single_point_run = (globenv%run_type_id == energy_run .OR. & + globenv%run_type_id == energy_force_run) + ELSE + qs_env%target_time = 0.0_dp + qs_env%start_time = 0.0_dp + END IF + + qs_env%sim_time = 0._dp + qs_env%sim_step = 0 + + qs_env%total_zeff_corr = 0.0_dp + qs_env%surface_dipole_moment = 0.0_dp + qs_env%surface_dipole_switch_off = .FALSE. + + // Zero all variables containing results + NULLIFY (qs_env%mo_derivs) + NULLIFY (qs_env%mo_loc_history) + + IF (.NOT. ASSOCIATED(qs_env%molecular_scf_guess_env)) ALLOCATE +(qs_env%molecular_scf_guess_env) + + END SUBROUTINE init_qs_env + +/*************************************************************************************************** + \brief Set the QUICKSTEP environment. + \param qs_env ... + \param super_cell ... + \param mos ... + \param qmmm ... + \param qmmm_periodic ... + \param ewald_env ... + \param ewald_pw ... + \param mpools ... + \param rho_external ... + \param external_vxc ... + \param mask ... + \param scf_control ... + \param rel_control ... + \param qs_charges ... + \param ks_env ... + \param ks_qmmm_env ... + \param wf_history ... + \param scf_env ... + \param active_space ... + \param input ... + \param oce ... + \param rho_atom_set ... + \param rho0_atom_set ... + \param rho0_mpole ... + \param run_rtp ... + \param rtp ... + \param rhoz_set ... + \param rhoz_tot ... + \param ecoul_1c ... + \param has_unit_metric ... + \param requires_mo_derivs ... + \param mo_derivs ... + \param mo_loc_history ... + \param efield ... + \param linres_control ... + \param xas_env ... + \param cp_ddapc_env ... + \param cp_ddapc_ewald ... + \param outer_scf_history ... + \param outer_scf_ihistory ... + \param x_data ... + \param et_coupling ... + \param dftb_potential ... + \param se_taper ... + \param se_store_int_env ... + \param se_nddo_mpole ... + \param se_nonbond_env ... + \param admm_env ... + \param ls_scf_env ... + \param do_transport ... + \param transport_env ... + \param lri_env ... + \param lri_density ... + \param exstate_env ... + \param ec_env ... + \param dispersion_env ... + \param gcp_env ... + \param mp2_env ... + \param bs_env ... + \param kg_env ... + \param force ... + \param kpoints ... + \param WannierCentres ... + \param almo_scf_env ... + \param gradient_history ... + \param variable_history ... + \param embed_pot ... + \param spin_embed_pot ... + \param polar_env ... + \param mos_last_converged ... [SGh] + \param rhs ... + \date 23.01.2002 + \author MK + \version 1.0 +***************************************************************************************************/ +void set_qs_env(qs_env, super_cell, & + mos, qmmm, qmmm_periodic, & + ewald_env, ewald_pw, mpools, & + rho_external, external_vxc, mask, & + scf_control, rel_control, qs_charges, ks_env, & + ks_qmmm_env, wf_history, scf_env, active_space, & + input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, & + rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, +mo_derivs, & + mo_loc_history, efield, & + linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, & + outer_scf_history, outer_scf_ihistory, x_data, et_coupling, +dftb_potential, & + se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, +ls_scf_env, & + do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, +dispersion_env, & + gcp_env, mp2_env, bs_env, kg_env, force, & + kpoints, WannierCentres, almo_scf_env, gradient_history, +variable_history, embed_pot, & + spin_embed_pot, polar_env, mos_last_converged, rhs) + + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(cell_type), OPTIONAL, POINTER :: super_cell + TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos + LOGICAL, OPTIONAL :: qmmm, qmmm_periodic + TYPE(ewald_environment_type), OPTIONAL, POINTER :: ewald_env + TYPE(ewald_pw_type), OPTIONAL, POINTER :: ewald_pw + TYPE(qs_matrix_pools_type), OPTIONAL, POINTER :: mpools + TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_external + TYPE(pw_type), OPTIONAL, POINTER :: external_vxc, mask + TYPE(scf_control_type), OPTIONAL, POINTER :: scf_control + TYPE(rel_control_type), OPTIONAL, POINTER :: rel_control + TYPE(qs_charges_type), OPTIONAL, POINTER :: qs_charges + TYPE(qs_ks_env_type), OPTIONAL, POINTER :: ks_env + TYPE(qs_ks_qmmm_env_type), OPTIONAL, POINTER :: ks_qmmm_env + TYPE(qs_wf_history_type), OPTIONAL, POINTER :: wf_history + TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env + TYPE(active_space_type), OPTIONAL, POINTER :: active_space + TYPE(section_vals_type), OPTIONAL, POINTER :: input + TYPE(oce_matrix_type), OPTIONAL, POINTER :: oce + TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, & + POINTER :: rho_atom_set + TYPE(rho0_atom_type), DIMENSION(:), OPTIONAL, & + POINTER :: rho0_atom_set + TYPE(rho0_mpole_type), OPTIONAL, POINTER :: rho0_mpole + LOGICAL, OPTIONAL :: run_rtp + TYPE(rt_prop_type), OPTIONAL, POINTER :: rtp + TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER :: rhoz_set + REAL(dp), OPTIONAL :: rhoz_tot + TYPE(ecoul_1center_type), DIMENSION(:), OPTIONAL, & + POINTER :: ecoul_1c + LOGICAL, OPTIONAL :: has_unit_metric, requires_mo_derivs + TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & + POINTER :: mo_derivs + TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history + TYPE(efield_berry_type), OPTIONAL, POINTER :: efield + TYPE(linres_control_type), OPTIONAL, POINTER :: linres_control + TYPE(xas_environment_type), OPTIONAL, POINTER :: xas_env + TYPE(cp_ddapc_type), OPTIONAL, POINTER :: cp_ddapc_env + TYPE(cp_ddapc_ewald_type), OPTIONAL, POINTER :: cp_ddapc_ewald + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: outer_scf_history + INTEGER, INTENT(IN), OPTIONAL :: outer_scf_ihistory + TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, POINTER :: x_data + TYPE(et_coupling_type), OPTIONAL, POINTER :: et_coupling + TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & + OPTIONAL, POINTER :: dftb_potential + TYPE(se_taper_type), OPTIONAL, POINTER :: se_taper + TYPE(semi_empirical_si_type), OPTIONAL, POINTER :: se_store_int_env + TYPE(nddo_mpole_type), OPTIONAL, POINTER :: se_nddo_mpole + TYPE(fist_nonbond_env_type), OPTIONAL, POINTER :: se_nonbond_env + TYPE(admm_type), OPTIONAL, POINTER :: admm_env + TYPE(ls_scf_env_type), OPTIONAL, POINTER :: ls_scf_env + LOGICAL, OPTIONAL :: do_transport + TYPE(transport_env_type), OPTIONAL, POINTER :: transport_env + TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env + TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density + TYPE(excited_energy_type), OPTIONAL, POINTER :: exstate_env + TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env + TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env + TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env + TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env + TYPE(post_scf_bandstructure_type), OPTIONAL, & + POINTER :: bs_env + TYPE(kg_environment_type), OPTIONAL, POINTER :: kg_env + TYPE(qs_force_type), DIMENSION(:), OPTIONAL, & + POINTER :: force + TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints + TYPE(wannier_centres_type), DIMENSION(:), & + OPTIONAL, POINTER :: WannierCentres + TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history + TYPE(pw_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot + TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env + TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged + REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rhs + + TYPE(qs_subsys_type), POINTER :: subsys + + IF (PRESENT(mp2_env)) qs_env%mp2_env => mp2_env + IF (PRESENT(bs_env)) qs_env%bs_env => bs_env + IF (PRESENT(kg_env)) qs_env%kg_env => kg_env + IF (PRESENT(super_cell)) THEN + CALL cell_retain(super_cell) + CALL cell_release(qs_env%super_cell) + qs_env%super_cell => super_cell + END IF + ! + IF (PRESENT(qmmm)) qs_env%qmmm = qmmm + IF (PRESENT(qmmm_periodic)) qs_env%qmmm_periodic = qmmm_periodic + IF (PRESENT(mos)) qs_env%mos => mos + IF (PRESENT(mos_last_converged)) qs_env%mos_last_converged => mos_last_converged + IF (PRESENT(ls_scf_env)) qs_env%ls_scf_env => ls_scf_env + IF (PRESENT(almo_scf_env)) qs_env%almo_scf_env => almo_scf_env + IF (PRESENT(do_transport)) qs_env%do_transport = do_transport + IF (PRESENT(transport_env)) qs_env%transport_env => transport_env + // if intels checking (-C) complains here, you have rediscovered a bug in the intel + // compiler (present in at least 10.0.025). A testcase has been submitted to intel. + IF (PRESENT(oce)) qs_env%oce => oce + IF (PRESENT(outer_scf_history)) qs_env%outer_scf_history => outer_scf_history + IF (PRESENT(gradient_history)) qs_env%gradient_history => gradient_history + IF (PRESENT(variable_history)) qs_env%variable_history => variable_history + IF (PRESENT(outer_scf_ihistory)) qs_env%outer_scf_ihistory = outer_scf_ihistory + IF (PRESENT(requires_mo_derivs)) qs_env%requires_mo_derivs = requires_mo_derivs + IF (PRESENT(has_unit_metric)) qs_env%has_unit_metric = has_unit_metric + IF (PRESENT(mo_derivs)) qs_env%mo_derivs => mo_derivs + IF (PRESENT(mo_loc_history)) qs_env%mo_loc_history => mo_loc_history + IF (PRESENT(run_rtp)) qs_env%run_rtp = run_rtp + IF (PRESENT(rtp)) qs_env%rtp => rtp + IF (PRESENT(efield)) qs_env%efield => efield + IF (PRESENT(active_space)) qs_env%active_space => active_space + + IF (PRESENT(ewald_env)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%ewald_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%ewald_env, ewald_env)) THEN + CALL ewald_env_release(qs_env%ewald_env) + DEALLOCATE (qs_env%ewald_env) + END IF + END IF + qs_env%ewald_env => ewald_env + END IF + IF (PRESENT(ewald_pw)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%ewald_pw)) THEN + IF (.NOT. ASSOCIATED(ewald_pw, qs_env%ewald_pw)) THEN + CALL ewald_pw_release(qs_env%ewald_pw) + DEALLOCATE (qs_env%ewald_pw) + END IF + END IF + qs_env%ewald_pw => ewald_pw + END IF + IF (PRESENT(scf_control)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%scf_control)) THEN + IF (.NOT. ASSOCIATED(qs_env%scf_control, scf_control)) THEN + CALL scf_c_release(qs_env%scf_control) + DEALLOCATE (qs_env%scf_control) + END IF + END IF + qs_env%scf_control => scf_control + END IF + IF (PRESENT(rel_control)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%rel_control)) THEN + IF (.NOT. ASSOCIATED(qs_env%rel_control, rel_control)) THEN + CALL rel_c_release(qs_env%rel_control) + DEALLOCATE (qs_env%rel_control) + END IF + END IF + qs_env%rel_control => rel_control + END IF + IF (PRESENT(linres_control)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%linres_control)) THEN + IF (.NOT. ASSOCIATED(qs_env%linres_control, linres_control)) THEN + CALL linres_control_release(qs_env%linres_control) + DEALLOCATE (qs_env%linres_control) + END IF + END IF + qs_env%linres_control => linres_control + END IF + // ZMP associating variables + IF (PRESENT(rho_external)) THEN + IF (ASSOCIATED(qs_env%rho_external)) THEN + IF (.NOT. ASSOCIATED(qs_env%rho_external, rho_external)) THEN + CALL qs_rho_release(qs_env%rho_external) + DEALLOCATE (qs_env%rho_external) + END IF + END IF + qs_env%rho_external => rho_external + END IF + IF (PRESENT(external_vxc)) qs_env%external_vxc => external_vxc + IF (PRESENT(mask)) qs_env%mask => mask + // Embedding potential + IF (PRESENT(embed_pot)) qs_env%embed_pot => embed_pot + IF (PRESENT(spin_embed_pot)) qs_env%spin_embed_pot => spin_embed_pot + + // Polarisability tensor + IF (PRESENT(polar_env)) qs_env%polar_env => polar_env + + IF (PRESENT(qs_charges)) THEN + IF (ASSOCIATED(qs_env%qs_charges)) THEN + IF (.NOT. ASSOCIATED(qs_env%qs_charges, qs_charges)) THEN + CALL qs_charges_release(qs_env%qs_charges) + DEALLOCATE (qs_env%qs_charges) + END IF + END IF + qs_env%qs_charges => qs_charges + END IF + IF (PRESENT(ks_qmmm_env)) THEN + IF (ASSOCIATED(qs_env%ks_qmmm_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%ks_qmmm_env, ks_qmmm_env)) THEN + CALL qs_ks_qmmm_release(qs_env%ks_qmmm_env) + DEALLOCATE (qs_env%ks_qmmm_env) + END IF + END IF + qs_env%ks_qmmm_env => ks_qmmm_env + END IF + IF (PRESENT(ks_env)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%ks_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%ks_env, ks_env)) THEN + CALL qs_ks_release(qs_env%ks_env) + DEALLOCATE (qs_env%ks_env) + END IF + END IF + qs_env%ks_env => ks_env + END IF + IF (PRESENT(wf_history)) THEN // accept also null pointers ? + CALL wfi_retain(wf_history) + CALL wfi_release(qs_env%wf_history) + qs_env%wf_history => wf_history + END IF + IF (PRESENT(scf_env)) THEN // accept also null pointers ? + IF (ASSOCIATED(qs_env%scf_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%scf_env, scf_env)) THEN + CALL scf_env_release(qs_env%scf_env) + DEALLOCATE (qs_env%scf_env) + END IF + END IF + qs_env%scf_env => scf_env + END IF + IF (PRESENT(xas_env)) THEN // accept also null pointers? + IF (ASSOCIATED(qs_env%xas_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%xas_env, xas_env)) THEN + CALL xas_env_release(qs_env%xas_env) + DEALLOCATE (qs_env%xas_env) + END IF + END IF + qs_env%xas_env => xas_env + END IF + IF (PRESENT(mpools)) THEN + CALL mpools_retain(mpools) + CALL mpools_release(qs_env%mpools) + qs_env%mpools => mpools + END IF + IF (PRESENT(rho_atom_set)) THEN + CALL set_local_rho(qs_env%local_rho_set, rho_atom_set=rho_atom_set) + END IF + IF (PRESENT(rho0_atom_set)) THEN + CALL set_local_rho(qs_env%local_rho_set, rho0_atom_set=rho0_atom_set) + END IF + IF (PRESENT(rho0_mpole)) THEN + CALL set_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_mpole) + END IF + IF (PRESENT(rhoz_set)) THEN + CALL set_local_rho(qs_env%local_rho_set, rhoz_set=rhoz_set) + END IF + IF (PRESENT(rhoz_tot)) qs_env%local_rho_set%rhoz_tot = rhoz_tot + IF (PRESENT(ecoul_1c)) THEN + CALL set_hartree_local(qs_env%hartree_local, ecoul_1c=ecoul_1c) + END IF + IF (PRESENT(input)) THEN + CALL section_vals_retain(input) + CALL section_vals_release(qs_env%input) + qs_env%input => input + END IF + IF (PRESENT(cp_ddapc_env)) THEN + IF (ASSOCIATED(qs_env%cp_ddapc_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%cp_ddapc_env, cp_ddapc_env)) THEN + CALL cp_ddapc_release(qs_env%cp_ddapc_env) + DEALLOCATE (qs_env%cp_ddapc_env) + END IF + END IF + qs_env%cp_ddapc_env => cp_ddapc_env + END IF + IF (PRESENT(cp_ddapc_ewald)) THEN + qs_env%cp_ddapc_ewald => cp_ddapc_ewald + END IF + IF (PRESENT(x_data)) qs_env%x_data => x_data + IF (PRESENT(et_coupling)) qs_env%et_coupling => et_coupling + IF (PRESENT(dftb_potential)) qs_env%dftb_potential => dftb_potential + IF (PRESENT(se_taper)) THEN + CALL se_taper_release(qs_env%se_taper) + qs_env%se_taper => se_taper + END IF + IF (PRESENT(se_store_int_env)) THEN + CALL semi_empirical_si_release(qs_env%se_store_int_env) + qs_env%se_store_int_env => se_store_int_env + END IF + IF (PRESENT(se_nddo_mpole)) THEN + CALL nddo_mpole_release(qs_env%se_nddo_mpole) + qs_env%se_nddo_mpole => se_nddo_mpole + END IF + IF (PRESENT(se_nonbond_env)) THEN + IF (ASSOCIATED(qs_env%se_nonbond_env)) THEN + IF (.NOT. ASSOCIATED(qs_env%se_nonbond_env, se_nonbond_env)) THEN + CALL fist_nonbond_env_release(qs_env%se_nonbond_env) + DEALLOCATE (qs_env%se_nonbond_env) + END IF + END IF + qs_env%se_nonbond_env => se_nonbond_env + END IF + IF (PRESENT(admm_env)) qs_env%admm_env => admm_env + IF (PRESENT(lri_env)) qs_env%lri_env => lri_env + IF (PRESENT(lri_density)) qs_env%lri_density => lri_density + IF (PRESENT(ec_env)) qs_env%ec_env => ec_env + IF (PRESENT(exstate_env)) qs_env%exstate_env => exstate_env + IF (PRESENT(dispersion_env)) qs_env%dispersion_env => dispersion_env + IF (PRESENT(gcp_env)) qs_env%gcp_env => gcp_env + IF (PRESENT(WannierCentres)) qs_env%WannierCentres => WannierCentres + IF (PRESENT(kpoints)) CALL set_ks_env(qs_env%ks_env, kpoints=kpoints) + + // Resp charges + IF (PRESENT(rhs)) qs_env%rhs => rhs + + IF (PRESENT(force)) THEN + CALL get_qs_env(qs_env, subsys=subsys) + CALL qs_subsys_set(subsys, force=force) + END IF + + END SUBROUTINE set_qs_env + +/*************************************************************************************************** + \brief allocates and intitializes a qs_env + \param qs_env the object to create + \param globenv ... + \par History + 12.2002 created [fawzi] + \author Fawzi Mohamed +***************************************************************************************************/ +void qs_env_create(qs_env, globenv) + TYPE(qs_environment_type), INTENT(OUT) :: qs_env + TYPE(global_environment_type), OPTIONAL, POINTER :: globenv + + CALL init_qs_env(qs_env, globenv=globenv) + END SUBROUTINE qs_env_create + +/*************************************************************************************************** + \brief releases the given qs_env (see doc/ReferenceCounting.html) + \param qs_env the object to release + \par History + 12.2002 created [fawzi] + 06.2018 polar_env added (MK) + \author Fawzi Mohamed +***************************************************************************************************/ +void qs_env_release(qs_env) + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + + INTEGER :: i + + CALL cell_release(qs_env%super_cell) + IF (ASSOCIATED(qs_env%mos)) THEN + DO i = 1, SIZE(qs_env%mos) + CALL deallocate_mo_set(qs_env%mos(i)) + END DO + DEALLOCATE (qs_env%mos) + END IF + IF (ASSOCIATED(qs_env%mos_last_converged)) THEN + DO i = 1, SIZE(qs_env%mos_last_converged) + CALL deallocate_mo_set(qs_env%mos_last_converged(i)) + END DO + DEALLOCATE (qs_env%mos_last_converged) + END IF + + IF (ASSOCIATED(qs_env%mo_derivs)) THEN + DO I = 1, SIZE(qs_env%mo_derivs) + CALL dbcsr_release_p(qs_env%mo_derivs(I)%matrix) + END DO + DEALLOCATE (qs_env%mo_derivs) + END IF + + CALL cp_fm_release(qs_env%mo_loc_history) + + IF (ASSOCIATED(qs_env%rtp)) THEN + CALL rt_prop_release(qs_env%rtp) + DEALLOCATE (qs_env%rtp) + END IF + IF (ASSOCIATED(qs_env%outer_scf_history)) THEN + DEALLOCATE (qs_env%outer_scf_history) + qs_env%outer_scf_ihistory = 0 + END IF + IF (ASSOCIATED(qs_env%gradient_history)) & + DEALLOCATE (qs_env%gradient_history) + IF (ASSOCIATED(qs_env%variable_history)) & + DEALLOCATE (qs_env%variable_history) + IF (ASSOCIATED(qs_env%oce)) CALL deallocate_oce_set(qs_env%oce) + IF (ASSOCIATED(qs_env%local_rho_set)) THEN + CALL local_rho_set_release(qs_env%local_rho_set) + END IF + IF (ASSOCIATED(qs_env%hartree_local)) THEN + CALL hartree_local_release(qs_env%hartree_local) + END IF + IF (ASSOCIATED(qs_env%scf_control)) THEN + CALL scf_c_release(qs_env%scf_control) + DEALLOCATE (qs_env%scf_control) + END IF + IF (ASSOCIATED(qs_env%rel_control)) THEN + CALL rel_c_release(qs_env%rel_control) + DEALLOCATE (qs_env%rel_control) + END IF + + IF (ASSOCIATED(qs_env%linres_control)) THEN + CALL linres_control_release(qs_env%linres_control) + DEALLOCATE (qs_env%linres_control) + END IF + + IF (ASSOCIATED(qs_env%almo_scf_env)) THEN + CALL almo_scf_env_release(qs_env%almo_scf_env) + END IF + + IF (ASSOCIATED(qs_env%ls_scf_env)) THEN + CALL ls_scf_release(qs_env%ls_scf_env) + END IF + IF (ASSOCIATED(qs_env%molecular_scf_guess_env)) THEN + CALL molecular_scf_guess_env_destroy(qs_env%molecular_scf_guess_env) + DEALLOCATE (qs_env%molecular_scf_guess_env) + END IF + + IF (ASSOCIATED(qs_env%transport_env)) THEN + CALL transport_env_release(qs_env%transport_env) + END IF + + !Only if do_xas_calculation + IF (ASSOCIATED(qs_env%xas_env)) THEN + CALL xas_env_release(qs_env%xas_env) + DEALLOCATE (qs_env%xas_env) + END IF + IF (ASSOCIATED(qs_env%ewald_env)) THEN + CALL ewald_env_release(qs_env%ewald_env) + DEALLOCATE (qs_env%ewald_env) + END IF + IF (ASSOCIATED(qs_env%ewald_pw)) THEN + CALL ewald_pw_release(qs_env%ewald_pw) + DEALLOCATE (qs_env%ewald_pw) + END IF + IF (ASSOCIATED(qs_env%image_matrix)) THEN + DEALLOCATE (qs_env%image_matrix) + END IF + IF (ASSOCIATED(qs_env%ipiv)) THEN + DEALLOCATE (qs_env%ipiv) + END IF + IF (ASSOCIATED(qs_env%image_coeff)) THEN + DEALLOCATE (qs_env%image_coeff) + END IF + // ZMP + IF (ASSOCIATED(qs_env%rho_external)) THEN + CALL qs_rho_release(qs_env%rho_external) + DEALLOCATE (qs_env%rho_external) + END IF + IF (ASSOCIATED(qs_env%external_vxc)) THEN + CALL qs_env%external_vxc%release() + DEALLOCATE (qs_env%external_vxc) + END IF + IF (ASSOCIATED(qs_env%mask)) THEN + CALL qs_env%mask%release() + DEALLOCATE (qs_env%mask) + END IF + IF (ASSOCIATED(qs_env%active_space)) THEN + CALL release_active_space_type(qs_env%active_space) + END IF + // Embedding potentials if provided as input + IF (qs_env%given_embed_pot) THEN + CALL qs_env%embed_pot%release() + DEALLOCATE (qs_env%embed_pot) + IF (ASSOCIATED(qs_env%spin_embed_pot)) THEN + CALL qs_env%spin_embed_pot%release() + DEALLOCATE (qs_env%spin_embed_pot) + END IF + END IF + + // Polarisability tensor + CALL polar_env_release(qs_env%polar_env) + + IF (ASSOCIATED(qs_env%qs_charges)) THEN + CALL qs_charges_release(qs_env%qs_charges) + DEALLOCATE (qs_env%qs_charges) + END IF + IF (ASSOCIATED(qs_env%ks_env)) THEN + CALL qs_ks_release(qs_env%ks_env) + DEALLOCATE (qs_env%ks_env) + END IF + IF (ASSOCIATED(qs_env%ks_qmmm_env)) THEN + CALL qs_ks_qmmm_release(qs_env%ks_qmmm_env) + DEALLOCATE (qs_env%ks_qmmm_env) + END IF + CALL wfi_release(qs_env%wf_history) + IF (ASSOCIATED(qs_env%scf_env)) THEN + CALL scf_env_release(qs_env%scf_env) + DEALLOCATE (qs_env%scf_env) + END IF + CALL mpools_release(qs_env%mpools) + CALL section_vals_release(qs_env%input) + IF (ASSOCIATED(qs_env%cp_ddapc_env)) THEN + CALL cp_ddapc_release(qs_env%cp_ddapc_env) + DEALLOCATE (qs_env%cp_ddapc_env) + END IF + CALL cp_ddapc_ewald_release(qs_env%cp_ddapc_ewald) + CALL efield_berry_release(qs_env%efield) + IF (ASSOCIATED(qs_env%x_data)) THEN + CALL hfx_release(qs_env%x_data) + END IF + IF (ASSOCIATED(qs_env%et_coupling)) THEN + CALL et_coupling_release(qs_env%et_coupling) + END IF + IF (ASSOCIATED(qs_env%dftb_potential)) THEN + CALL qs_dftb_pairpot_release(qs_env%dftb_potential) + END IF + IF (ASSOCIATED(qs_env%se_taper)) THEN + CALL se_taper_release(qs_env%se_taper) + END IF + IF (ASSOCIATED(qs_env%se_store_int_env)) THEN + CALL semi_empirical_si_release(qs_env%se_store_int_env) + END IF + IF (ASSOCIATED(qs_env%se_nddo_mpole)) THEN + CALL nddo_mpole_release(qs_env%se_nddo_mpole) + END IF + IF (ASSOCIATED(qs_env%se_nonbond_env)) THEN + CALL fist_nonbond_env_release(qs_env%se_nonbond_env) + DEALLOCATE (qs_env%se_nonbond_env) + END IF + IF (ASSOCIATED(qs_env%admm_env)) THEN + CALL admm_env_release(qs_env%admm_env) + END IF + IF (ASSOCIATED(qs_env%lri_env)) THEN + CALL lri_env_release(qs_env%lri_env) + DEALLOCATE (qs_env%lri_env) + END IF + IF (ASSOCIATED(qs_env%lri_density)) THEN + CALL lri_density_release(qs_env%lri_density) + DEALLOCATE (qs_env%lri_density) + END IF + IF (ASSOCIATED(qs_env%ec_env)) THEN + CALL ec_env_release(qs_env%ec_env) + END IF + IF (ASSOCIATED(qs_env%exstate_env)) THEN + CALL exstate_release(qs_env%exstate_env) + END IF + IF (ASSOCIATED(qs_env%mp2_env)) THEN + CALL mp2_env_release(qs_env%mp2_env) + DEALLOCATE (qs_env%mp2_env) + NULLIFY (qs_env%mp2_env) + END IF + IF (ASSOCIATED(qs_env%bs_env)) THEN + CALL bs_env_release(qs_env%bs_env) + END IF + IF (ASSOCIATED(qs_env%kg_env)) THEN + CALL kg_env_release(qs_env%kg_env) + END IF + + // dispersion + CALL qs_dispersion_release(qs_env%dispersion_env) + // gCP + IF (ASSOCIATED(qs_env%gcp_env)) THEN + CALL qs_gcp_release(qs_env%gcp_env) + END IF + + IF (ASSOCIATED(qs_env%WannierCentres)) THEN + DO i = 1, SIZE(qs_env%WannierCentres) + DEALLOCATE (qs_env%WannierCentres(i)%WannierHamDiag) + DEALLOCATE (qs_env%WannierCentres(i)%centres) + END DO + DEALLOCATE (qs_env%WannierCentres) + END IF + // Resp charges + IF (ASSOCIATED(qs_env%rhs)) DEALLOCATE (qs_env%rhs) + + END SUBROUTINE qs_env_release + +/*************************************************************************************************** + \brief releases part of the given qs_env in order to save memory + \param qs_env the object to release + \par History + 04.2022 created [JGH] +***************************************************************************************************/ +void qs_env_part_release(qs_env) + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + + INTEGER :: i + + IF (ASSOCIATED(qs_env%mos_last_converged)) THEN + DO i = 1, SIZE(qs_env%mos_last_converged) + CALL deallocate_mo_set(qs_env%mos_last_converged(i)) + END DO + DEALLOCATE (qs_env%mos_last_converged) + END IF + + IF (ASSOCIATED(qs_env%mo_derivs)) THEN + DO I = 1, SIZE(qs_env%mo_derivs) + CALL dbcsr_release_p(qs_env%mo_derivs(I)%matrix) + END DO + DEALLOCATE (qs_env%mo_derivs) + END IF + + CALL cp_fm_release(qs_env%mo_loc_history) + + IF (ASSOCIATED(qs_env%rtp)) THEN + CALL rt_prop_release(qs_env%rtp) + DEALLOCATE (qs_env%rtp) + END IF + IF (ASSOCIATED(qs_env%outer_scf_history)) THEN + DEALLOCATE (qs_env%outer_scf_history) + qs_env%outer_scf_ihistory = 0 + END IF + IF (ASSOCIATED(qs_env%gradient_history)) & + DEALLOCATE (qs_env%gradient_history) + IF (ASSOCIATED(qs_env%variable_history)) & + DEALLOCATE (qs_env%variable_history) + IF (ASSOCIATED(qs_env%oce)) CALL deallocate_oce_set(qs_env%oce) + IF (ASSOCIATED(qs_env%local_rho_set)) THEN + CALL local_rho_set_release(qs_env%local_rho_set) + END IF + IF (ASSOCIATED(qs_env%hartree_local)) THEN + CALL hartree_local_release(qs_env%hartree_local) + END IF + IF (ASSOCIATED(qs_env%scf_control)) THEN + CALL scf_c_release(qs_env%scf_control) + DEALLOCATE (qs_env%scf_control) + END IF + IF (ASSOCIATED(qs_env%rel_control)) THEN + CALL rel_c_release(qs_env%rel_control) + DEALLOCATE (qs_env%rel_control) + END IF + + IF (ASSOCIATED(qs_env%linres_control)) THEN + CALL linres_control_release(qs_env%linres_control) + DEALLOCATE (qs_env%linres_control) + END IF + + IF (ASSOCIATED(qs_env%almo_scf_env)) THEN + CALL almo_scf_env_release(qs_env%almo_scf_env) + END IF + + IF (ASSOCIATED(qs_env%ls_scf_env)) THEN + CALL ls_scf_release(qs_env%ls_scf_env) + END IF + IF (ASSOCIATED(qs_env%molecular_scf_guess_env)) THEN + CALL molecular_scf_guess_env_destroy(qs_env%molecular_scf_guess_env) + DEALLOCATE (qs_env%molecular_scf_guess_env) + END IF + + IF (ASSOCIATED(qs_env%transport_env)) THEN + CALL transport_env_release(qs_env%transport_env) + END IF + + !Only if do_xas_calculation + IF (ASSOCIATED(qs_env%xas_env)) THEN + CALL xas_env_release(qs_env%xas_env) + DEALLOCATE (qs_env%xas_env) + END IF + IF (ASSOCIATED(qs_env%ewald_env)) THEN + CALL ewald_env_release(qs_env%ewald_env) + DEALLOCATE (qs_env%ewald_env) + END IF + IF (ASSOCIATED(qs_env%ewald_pw)) THEN + CALL ewald_pw_release(qs_env%ewald_pw) + DEALLOCATE (qs_env%ewald_pw) + END IF + IF (ASSOCIATED(qs_env%image_matrix)) THEN + DEALLOCATE (qs_env%image_matrix) + END IF + IF (ASSOCIATED(qs_env%ipiv)) THEN + DEALLOCATE (qs_env%ipiv) + END IF + IF (ASSOCIATED(qs_env%image_coeff)) THEN + DEALLOCATE (qs_env%image_coeff) + END IF + // ZMP + IF (ASSOCIATED(qs_env%rho_external)) THEN + CALL qs_rho_release(qs_env%rho_external) + DEALLOCATE (qs_env%rho_external) + END IF + IF (ASSOCIATED(qs_env%external_vxc)) THEN + CALL qs_env%external_vxc%release() + DEALLOCATE (qs_env%external_vxc) + END IF + IF (ASSOCIATED(qs_env%mask)) THEN + CALL qs_env%mask%release() + DEALLOCATE (qs_env%mask) + END IF + IF (ASSOCIATED(qs_env%active_space)) THEN + CALL release_active_space_type(qs_env%active_space) + END IF + // Embedding potentials if provided as input + IF (qs_env%given_embed_pot) THEN + CALL qs_env%embed_pot%release() + DEALLOCATE (qs_env%embed_pot) + IF (ASSOCIATED(qs_env%spin_embed_pot)) THEN + CALL qs_env%spin_embed_pot%release() + DEALLOCATE (qs_env%spin_embed_pot) + END IF + END IF + + // Polarisability tensor + CALL polar_env_release(qs_env%polar_env) + + IF (ASSOCIATED(qs_env%qs_charges)) THEN + CALL qs_charges_release(qs_env%qs_charges) + DEALLOCATE (qs_env%qs_charges) + END IF + CALL qs_ks_part_release(qs_env%ks_env) + IF (ASSOCIATED(qs_env%ks_qmmm_env)) THEN + CALL qs_ks_qmmm_release(qs_env%ks_qmmm_env) + DEALLOCATE (qs_env%ks_qmmm_env) + END IF + CALL wfi_release(qs_env%wf_history) + IF (ASSOCIATED(qs_env%scf_env)) THEN + CALL scf_env_release(qs_env%scf_env) + DEALLOCATE (qs_env%scf_env) + END IF + IF (ASSOCIATED(qs_env%cp_ddapc_env)) THEN + CALL cp_ddapc_release(qs_env%cp_ddapc_env) + DEALLOCATE (qs_env%cp_ddapc_env) + END IF + CALL cp_ddapc_ewald_release(qs_env%cp_ddapc_ewald) + CALL efield_berry_release(qs_env%efield) + IF (ASSOCIATED(qs_env%x_data)) THEN + CALL hfx_release(qs_env%x_data) + END IF + IF (ASSOCIATED(qs_env%et_coupling)) THEN + CALL et_coupling_release(qs_env%et_coupling) + END IF + IF (ASSOCIATED(qs_env%dftb_potential)) THEN + CALL qs_dftb_pairpot_release(qs_env%dftb_potential) + END IF + IF (ASSOCIATED(qs_env%se_taper)) THEN + CALL se_taper_release(qs_env%se_taper) + END IF + IF (ASSOCIATED(qs_env%se_store_int_env)) THEN + CALL semi_empirical_si_release(qs_env%se_store_int_env) + END IF + IF (ASSOCIATED(qs_env%se_nddo_mpole)) THEN + CALL nddo_mpole_release(qs_env%se_nddo_mpole) + END IF + IF (ASSOCIATED(qs_env%se_nonbond_env)) THEN + CALL fist_nonbond_env_release(qs_env%se_nonbond_env) + DEALLOCATE (qs_env%se_nonbond_env) + END IF + IF (ASSOCIATED(qs_env%admm_env)) THEN + CALL admm_env_release(qs_env%admm_env) + END IF + IF (ASSOCIATED(qs_env%lri_env)) THEN + CALL lri_env_release(qs_env%lri_env) + DEALLOCATE (qs_env%lri_env) + END IF + IF (ASSOCIATED(qs_env%lri_density)) THEN + CALL lri_density_release(qs_env%lri_density) + DEALLOCATE (qs_env%lri_density) + END IF + IF (ASSOCIATED(qs_env%ec_env)) THEN + CALL ec_env_release(qs_env%ec_env) + END IF + IF (ASSOCIATED(qs_env%exstate_env)) THEN + CALL exstate_release(qs_env%exstate_env) + END IF + IF (ASSOCIATED(qs_env%mp2_env)) THEN + CALL mp2_env_release(qs_env%mp2_env) + DEALLOCATE (qs_env%mp2_env) + NULLIFY (qs_env%mp2_env) + END IF + IF (ASSOCIATED(qs_env%kg_env)) THEN + CALL kg_env_release(qs_env%kg_env) + END IF + + // dispersion + CALL qs_dispersion_release(qs_env%dispersion_env) + // gCP + IF (ASSOCIATED(qs_env%gcp_env)) THEN + CALL qs_gcp_release(qs_env%gcp_env) + END IF + + IF (ASSOCIATED(qs_env%WannierCentres)) THEN + DO i = 1, SIZE(qs_env%WannierCentres) + DEALLOCATE (qs_env%WannierCentres(i)%WannierHamDiag) + DEALLOCATE (qs_env%WannierCentres(i)%centres) + END DO + DEALLOCATE (qs_env%WannierCentres) + END IF + // Resp charges + IF (ASSOCIATED(qs_env%rhs)) DEALLOCATE (qs_env%rhs) + + END SUBROUTINE qs_env_part_release diff --git a/src/qs_ot_eigensolver.c b/src/qs_ot_eigensolver.c new file mode 100644 index 0000000..fc7877e --- /dev/null +++ b/src/qs_ot_eigensolver.c @@ -0,0 +1,77 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + + + +#include + +#include "./base/base_uses.h" + +/******************************************************************************* + * \brief an eigen-space solver for the generalised symmetric eigenvalue problem + * for sparse matrices, needing only multiplications + * \author Joost VandeVondele (25.08.2002) +*******************************************************************************/ + +const char* moduleN = "qs_ot_eigensolver"; + +// on input c contains the initial guess (should not be zero !) +// on output c spans the subspace +/******************************************************************************* + * \brief ... + * \param matrix_h ... + * \param matrix_s ... + * \param matrix_orthogonal_space_fm ... + * \param matrix_c_fm ... + * \param preconditioner ... + * \param eps_gradient ... + * \param iter_max ... + * \param size_ortho_space ... + * \param silent ... + * \param ot_settings ... +*******************************************************************************/ + +void ot_eigensolver(struct dbcsr_type *matrix_h, struct dbcsr_type *matrix_s, + struct cp_fm_type matrix_orthogonal_space_fm, struct cp_fm_type matrix_c_fm, + struct preconditioner_type *preconditioner, double eps_gradient, int iter_max, + int size_ortho_space, bool silent, struct qs_ot_settings_type ot_settings) { + + const char* routineN = "ot_eigensolver"; + const int max_iter_inner_loop = 40; + double rone = 1.0; + double rzero = 0.0; + + int handle, ieigensolver, iter_total, k, n, ortho_k, ortho_space_k, output_unit; + + bool energy_only, my_silent, ortho; + double delta, energy; + + struct dbcsr_p_type *matrix_hc[]; + struct dbcsr_type *matrix_buf1_ortho, *matrix_buf2_ortho, *matrix_c, *matrix_orthogonal_space, + *matrix_os_ortho, *matrix_s_ortho; + struct qs_ot_type *qs_ot_env[]; + + timeset(routineN, handle); + + output_unit = cp_logger_get_default_io_unit(); + + if (PRESENT(silent)) + my_silent = silent; + else + my_silent = false; + + matrix_c = NULL; // fm->dbcsr + + cp_fm_get_info(matrix_c_fm, n, k); // fm->dbcsr + cp_fm_to_dbcsr_row_template(matrix_c, matrix_c_fm, matrix_h); + + iter_total = 0; + + + + +} \ No newline at end of file diff --git a/src/qs_tenosrs.c b/src/qs_tenosrs.c new file mode 100644 index 0000000..113fb4a --- /dev/null +++ b/src/qs_tenosrs.c @@ -0,0 +1,112 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Utility methods to build 3-center integral tensors of various types. +*******************************************************************************/ + +#include + +#include "./base/base_uses.h" + +void build_2c_neighbor_lists(int** ij_list, struct gto_basis_set_p_type* basis_i, struct gto_basis_set_p_type* basis_j, + int potential_parameter, char* name, int* qs_env, bool sym_ij, bool molecular, + float* dist_2d, int pot_to_rad) { + + int ikind, nkind, pot_to_rad_prv; + bool* i_present, j_present; + double* pair_radius; + double subcells; + double* i_radius, j_radius; + int *atomic_kind_set; + int *cell; + int *local_particles; + int *dist_2d_prv; + int *atom2d; + int *molecule_set; + int *particle_set; + + if (PRESENT(pot_to_rad)) { + pot_to_rad_prv = pot_to_rad; + } else { + pot_to_rad_prv = 1; + } + + get_qs_env(qs_env, nkind, cell, particle_set, atomic_kind_set, + local_particles, dist_2d_prv, molecule_set); + + section_vals_val_get(qs_env.input, "DFT.SUBCELLS", subcells); + + i_present = calloc(sizeof(bool) * nkind); + j_present = calloc(sizeof(bool) * nkind); + i_radius = calloc(sizeof(double) * nkind); + j_radius = calloc(sizeof(double) * nkind); + + if (PRESENT(dist_2d)) + dist_2d_prv = &dist2d; + + // Set up the radii, depending on the operator type + if (potential_parameter.potential_type == do_potential_id) { + + //overlap => use the kind radius for both i and j + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind)); + } + if (ASSOCIATED(basis_j(ikind)%gto_basis_set)) { + j_present(ikind) = true; + get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind)); + } + } + } else if (potential_parameter.potential_type == do_potential_coulomb) { + + //Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + if (pot_to_rad_prv == 1) + i_radius(ikind) = 1000000.0; + } + if (ASSOCIATED(basis_j(ikind).gto_basis_set)) { + j_present(ikind) = true; + if (pot_to_rad_prv == 2) + j_radius(ikind) = 1000000.0; + } + } //ikind + + } else if (potential_parameter.potential_type == do_potential_truncated || + potential_parameter.potential_type == do_potential_short) { + + //Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + get_gto_basis_set(basis_i(ikind).gto_basis_set, kind_radius=i_radius(ikind)); + if (pot_to_rad_prv == 1) + i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter.cutoff_radius; + } + if (ASSOCIATED(basis_j(ikind).gto_basis_set)) { + j_present(ikind) = true; + get_gto_basis_set(basis_j(ikind).gto_basis_set, kind_radius=j_radius(ikind)); + if (pot_to_rad_prv == 2) + j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter.cutoff_radius; + } + } + + } else { + CPABORT("Operator not implemented."); + } + + pair_radius = calloc(sizeof(double)*nkind*nkind); + + + + + +} + diff --git a/src/sockets.c b/src/sockets.c new file mode 100644 index 0000000..f9ef530 --- /dev/null +++ b/src/sockets.c @@ -0,0 +1,177 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/*----------------------------------------------------------------------------*/ +/* Copyright (C) 2013, Joshua More and Michele Ceriotti */ +/* */ +/* Permission is hereby granted, free of charge, to any person obtaining */ +/* a copy of this software and associated documentation files (the */ +/* "Software"), to deal in the Software without restriction, including */ +/* without limitation the rights to use, copy, modify, merge, publish, */ +/* distribute, sublicense, and/or sell copies of the Software, and to */ +/* permit persons to whom the Software is furnished to do so, subject to */ +/* the following conditions: */ +/* */ +/* The above copyright notice and this permission notice shall be included */ +/* in all copies or substantial portions of the Software. */ +/* */ +/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, */ +/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF */ +/* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. */ +/* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY */ +/* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, */ +/* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE */ +/* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief A minimal wrapper for socket communication. + * Contains both the functions that transmit data to the socket and read + * the data back out again once finished, and the function which opens + * the socket initially. Can be linked to a FORTRAN code that does not + * support sockets natively. + * \author Joshua More and Michele Ceriotti + ******************************************************************************/ +#ifndef __NO_IPI_DRIVER + +#define _POSIX_C_SOURCE 200809L + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/******************************************************************************* + * \brief Opens a socket. + * \param psockfd The id of the socket that will be created. + * \param inet An integer that determines whether the socket will be an inet + * or unix domain socket. Gives unix if 0, inet otherwise. + * \param port The port number for the socket to be created. Low numbers are + * often reserved for important channels, so use of numbers of 4 + * or more digits is recommended. + * \param host The name of the host server. + * \note Fortran passes an extra argument for the string length, but this is + * ignored here for C compatibility. + ******************************************************************************/ +void open_socket(int *psockfd, int *inet, int *port, char *host) { + int sockfd, ai_err; + + if (*inet > 0) { // creates an internet socket + + // fetches information on the host + struct addrinfo hints, *res; + char service[256]; + + memset(&hints, 0, sizeof(hints)); + hints.ai_socktype = SOCK_STREAM; + hints.ai_family = AF_INET; + hints.ai_flags = AI_PASSIVE; + + sprintf(service, "%d", *port); // convert the port number to a string + ai_err = getaddrinfo(host, service, &hints, &res); + if (ai_err != 0) { + perror("Error fetching host data. Wrong host name?"); + exit(-1); + } + + // creates socket + sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol); + if (sockfd < 0) { + perror("Error opening socket"); + exit(-1); + } + + // makes connection + if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0) { + perror("Error opening INET socket: wrong port or server unreachable"); + exit(-1); + } + freeaddrinfo(res); + } else { // creates a unix socket + struct sockaddr_un serv_addr; + + // fills up details of the socket address + memset(&serv_addr, 0, sizeof(serv_addr)); + serv_addr.sun_family = AF_UNIX; + strcpy(serv_addr.sun_path, "/tmp/ipi_"); + strcpy(serv_addr.sun_path + 9, host); + + // creates the socket + sockfd = socket(AF_UNIX, SOCK_STREAM, 0); + + // connects + if (connect(sockfd, (struct sockaddr *)&serv_addr, sizeof(serv_addr)) < 0) { + perror( + "Error opening UNIX socket: path unavailable, or already existing"); + exit(-1); + } + } + + *psockfd = sockfd; +} + +/******************************************************************************* + * \brief Writes to a socket. + * \param psockfd The id of the socket that will be written to. + * \param data The data to be written to the socket. + * \param plen The length of the data in bytes. + ******************************************************************************/ +void writebuffer(int *psockfd, char *data, int *plen) { + int n; + int sockfd = *psockfd; + int len = *plen; + + n = write(sockfd, data, len); + if (n < 0) { + perror("Error writing to socket: server has quit or connection broke"); + exit(-1); + } +} + +/******************************************************************************* + * \brief Reads from a socket. + * \param psockfd The id of the socket that will be read from. + * \param data The storage array for data read from the socket. + * \param plen The length of the data in bytes. + ******************************************************************************/ +void readbuffer(int *psockfd, char *data, int *plen) { + int n, nr; + int sockfd = *psockfd; + int len = *plen; + + n = nr = read(sockfd, data, len); + + while (nr > 0 && n < len) { + nr = read(sockfd, &data[n], len - n); + n += nr; + } + + if (n == 0) { + perror("Error reading from socket: server has quit or connection broke"); + exit(-1); + } + } + +/******************************************************************************* + * \brief Mini-wrapper to nanosleep + * \param dsec number of seconds to wait (float values accepted) + ******************************************************************************/ +void uwait(double *dsec) { + struct timespec wt, rem; + wt.tv_sec = floor(*dsec); + wt.tv_nsec = (*dsec - wt.tv_sec) * 1000000000; + nanosleep(&wt, &rem); +} + +#endif \ No newline at end of file diff --git a/src/start/cp2k.c b/src/start/cp2k.c new file mode 100644 index 0000000..a50a19e --- /dev/null +++ b/src/start/cp2k.c @@ -0,0 +1,103 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Main program of CP2K +//> \par Copyright +//> CP2K: A general program to perform molecular dynamics simulations +//> Copyright (C) 2000, 2001, 2002, 2003 CP2K developers group +//> Copyright (C) 2004, 2005, 2006, 2007 CP2K developers group +//> Copyright (C) 2008, 2009, 2010, 2011 CP2K developers group +//> Copyright (C) 2012, 2013, 2014, 2015 CP2K developers group +//> Copyright (C) 2016 CP2K developers group +//> \par +//> This program is free software; you can redistribute it and/or modify +//> it under the terms of the GNU General Public License as published by +//> the Free Software Foundation; either version 2 of the License, or +//> (at your option) any later version. +//> \par +//> This program is distributed in the hope that it will be useful, +//> but WITHOUT ANY WARRANTY; without even the implied warranty of +//> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//> GNU General Public License for more details. +//> \par +//> You should have received a copy of the GNU General Public License +//> along with this program; if not, write to the Free Software +//> Foundation, Inc., 51 Franklin Street, Fifth Floor, +//> Boston, MA 02110-1301, USA. +//> \par +//> See also https://www.fsf.org/licensing/licenses/gpl.html +//> \par +//> CP2K, including its sources and pointers to the authors +//> can be found at https://www.cp2k.org/ +//> \note +//> should be kept as lean as possible. +//> see cp2k_run for more comments +//> \author Joost VandeVondele +// ************************************************************************************************** +#include + +#include "../base/kinds.h" +#include "../base/base_uses.h" + +int main(int argc, char** argv) { + + char input_file_name[default_path_length]; + + int output_unit, l, i, var_set_sep, inp_var_idx, ierr, i_arg; + + bool check, usage, echo_input, command_line_error,run_it; + bool force_run, has_input, xml, print_version, print_license, shell_mode; + int* input_declaration; + + // output goes to the screen by default + output_unit = default_output_unit; + + // set default behaviour for the command line switches + check = false; + usage = false; + echo_input = false; + has_input = false; + run_it = true; + shell_mode = false; + force_run = false; + print_version = false; + print_license = false; + + command_line_error = false; + xml = false; + input_file_name = "Missing input file name" // no default + output_file_name = "__STD_OUT__" // by default we go to std_out + ALLOCATE (initial_variables(2, 1:0)) + + // Get command and strip path + GET_COMMAND_ARGUMENT(NUMBER=0, VALUE=command, STATUS=ierr) + CPASSERT(ierr == 0) + l = LEN_TRIM(command) + DO i = l, 1, -1 + IF (command(i:i) == "/" || command(i:i) == "\\") EXIT + END DO + command = command(i + 1:l) + + // Consider output redirection + i_arg = 0; + + + + if (i == 0) { + else { + section_release(input_declaration); + } + else { + printf("initial setup (MPI ?) error"); + } + + // and the final cleanup + finalize_cp2k(finalize_mpi=true, ierr=ierr); + delete initial_variables; + CPASSERT(ierr == 0); +} diff --git a/src/start/cp2k.cpp b/src/start/cp2k.cpp new file mode 100644 index 0000000..f970da7 --- /dev/null +++ b/src/start/cp2k.cpp @@ -0,0 +1,337 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Main program of CP2K +//> \par Copyright +//> CP2K: A general program to perform molecular dynamics simulations +//> Copyright (C) 2000, 2001, 2002, 2003 CP2K developers group +//> Copyright (C) 2004, 2005, 2006, 2007 CP2K developers group +//> Copyright (C) 2008, 2009, 2010, 2011 CP2K developers group +//> Copyright (C) 2012, 2013, 2014, 2015 CP2K developers group +//> Copyright (C) 2016 CP2K developers group +//> \par +//> This program is free software; you can redistribute it and/or modify +//> it under the terms of the GNU General Public License as published by +//> the Free Software Foundation; either version 2 of the License, or +//> (at your option) any later version. +//> \par +//> This program is distributed in the hope that it will be useful, +//> but WITHOUT ANY WARRANTY; without even the implied warranty of +//> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//> GNU General Public License for more details. +//> \par +//> You should have received a copy of the GNU General Public License +//> along with this program; if not, write to the Free Software +//> Foundation, Inc., 51 Franklin Street, Fifth Floor, +//> Boston, MA 02110-1301, USA. +//> \par +//> See also https://www.fsf.org/licensing/licenses/gpl.html +//> \par +//> CP2K, including its sources and pointers to the authors +//> can be found at https://www.cp2k.org/ +//> \note +//> should be kept as lean as possible. +//> see cp2k_run for more comments +//> \author Joost VandeVondele +// ************************************************************************************************** + +#include +#include +#include +#include +#include +#include + +#include "base_uses.h" + +int main(int argc, char** argv) { + + std::string input_file_name, output_file_name, arg_att, command; + std::vector initial_variables, initial_variables_tmp; + std::string compiler_options_string; + + int output_unit, l, i, var_set_sep, inp_var_idx, ierr, i_arg; + + bool check, usage, echo_input, command_line_error, run_it; + bool force_run, has_input, xml, print_version, print_license, shell_mode; + int *input_declaration; + + // output goes to the screen by default + output_unit = default_output_unit; + + // set default behaviour for the command line switches + check = false; + usage = false; + echo_input = false; + has_input = false; + run_it = true; + shell_mode = false; + force_run = false; + print_version = false; + print_license = false; + command_line_error = false; + xml = false; + input_file_name = "Missing input file name"; // no default + output_file_name = "__STD_OUT__"; // by default we go to std_out + + // Get command and strip path + get_command_argument(number=0, value=command, status=ierr); + assert(ierr == 0); + l = len_trim(command); + for (i = l-1; i >= 0; i--) { + if (command[i] == '/' or command[i] == '\\') exit(1); + } + command = command.substr(i); + + // Consider output redirection + i_arg = 0; + while (i_arg < command_argument_count()) { + i_arg++; + get_command_argument(number=&i_arg, value=&arg_att, status=&ierr); + assert(ierr == 0); + switch (arg_att) { + case "-o": + if (output_file_name == "__STD_OUT__") { + // Consider only the first -o flag + i_arg++; + get_command_argument(number=&i_arg, value=&arg_att, status=&ierr); + assert(ierr == 0); + if (arg_att[0] == '-') { + std::cout << "ERROR: The output file name " << arg_att << " starts with -" << std::endl; + command_line_error = true; + } else { + output_file_name = arg_att; + open_file(file_name=output_file_name, + file_status="UNKNOWN", + file_action="WRITE", + file_position="APPEND", + skip_get_unit_number=true, + unit_number=output_unit); + } + } else { + i_arg++; + std::cout << "ERROR: The command line flag -o has been specified multiple times" << std::endl; + command_line_error = true; + } + break; + } + } + + // Check if binary was invoked as cp2k_shell + if (command.compare("cp2k_shell")) { + shell_mode = true; + run_it = false; + } else if (command_argument_count() < 1) { + std::cout << "ERROR: At least one command line argument must be specified" << std::endl; + command_line_error = true; + } + + // Check if binary was invoked as sopt or popt alias + l = len_trim(command); + if (command[l-4:l] == ".sopt" or command[l-4:l] == ".popt") { + omp_set_num_threads(1); + } + +#ifdef __ACCELERATE + if (omp_get_max_threads() > 1) { + std::string _block_env_var; + int _block_veclib_max_threads, _block_ierr; + get_environment_variable("VECLIB_MAXIMUM_THREADS", _block_env_var, status=_block_ierr); + _block_veclib_max_threads = 0; + if (_block_ierr == 0) { + _block_env_var = _block_veclib_max_threads; + } + if (_block_ierr == 1 or (_block_ierr == 0 and _block_veclib_max_threads > 1)) { + cp_warn(__LOCATON__, "macOS' Accelerate framework has its own threading enabled which may interfere" + " with the OpenMP threading. You can disable the Accelerate threading by setting" + " the environment variable VECLIB_MAXIMUM_THREADS=1") + } + } +#endif + + i_arg = 0; + while (i_arg < command_argument_count()) { + i_arg++; + get_command_argument(&i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + switch (arg_att) { + case "--check": case "-c": + check = true; + run_it = false; + echo_input = false; + break; + case "--echo": case "-e": + check = true; + run_it = false; + echo_input = true; + break; + case "-v": case "--version": + print_version = true; + run_it = false; + break; + case "--license": + print_license = true; + run_it = false; + break; + case "--run": case "-r": + force_run = true; + break; + case "--shell": case "-s": + shell_mode = true; + run_it = false; + break; + case "-help": case "--help": case "-h": + usage = true; + run_it = false; + break; + case "-i": + i_arg++; + get_command_argument(i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + // argument does not start with a - it is a filename + if (!(arg_att[0] == '-')) { + input_file_name = arg_att; + has_input = true; + } else { + std::cout << "ERROR: The input file name " << arg_att << " starts with -" << std::endl; + command_line_error = true; + exit(1); + } + break; + case "-E": case "--set": + i_arg++; + get_command_argument(i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + + var_set_sep = arg_att.find("="); + + if (var_set_sep < 1) { + std::cout << "ERROR: Invalid initializer for preprocessor variable: " << arg_att << std::endl; + command_line_error = false; + exit(1); + } + + for (inp_var_idx = 0; inp_var_idx < sizeof(initial_variables, std::string); inp_var_idx++) { + // check whether the variable was already set, in this case, overwrite + if (initial_variables[0:inp_var_idx] == arg_att[0:var_set_sep - 1]) exit(1); + } + + break; + } + } + + if ((run_it || force_run || check || echo_input) && !has_input && !command_line_error) { + std::cout << "\n ERROR: An input file name is required" << std::endl; + command_line_error = true; + } + + init_cp2k(init_mpi=true, ierr=&ierr); + + if (ierr == 0){ + // some first info concerning how to run CP2K + if (usage || command_line_error) { + if (default_para_env.is_source()) { + l = len_trim(command); + std::cout << "\n " << command << " [-c|--check] [-e|--echo] [-h|--help]\n"; + std::cout << std::string(l, ' ') << " [-i] \n"; + std::cout << std::string(l, ' ') << " [-mpi-mapping|--mpi-mapping] \n"; + std::cout << std::string(l, ' ') << " [-o] \n"; + std::cout << std::string(l, ' ') << " [-r|-run] [-s|--shell] [--xml]\n"; + + std::cout << "\n starts the CP2K program, see \n"; + std::cout << "\n The easiest way is " << command << " \n"; + std::cout << "\n The following options can be used:\n"; + std::cout << "\n -i : provides an input file name, if it is the last\n"; + std::cout << std::string(24, ' ') << " argument, the -i flag is not needed\n"; + std::cout << " -o : provides an output file name [default: screen]\n"; + + std::cout << "\n These switches skip the simulation, unless [-r|-run] is specified:\n"; + std::cout << "\n --check, -c : performs a syntax check of the \n"; + std::cout << " --echo, -e : echoes the , and make all defaults explicit\n"; + std::cout << " The input is also checked, but only a failure is reported\n"; + std::cout << " --help, -h : writes this message\n"; + std::cout << " --license : prints the CP2K license\n"; + std::cout << " --mpi-mapping : applies a given MPI reordering to CP2K\n"; + std::cout << " --run, -r : forces a CP2K run regardless of other specified flags\n"; + std::cout << " --shell, -s : start interactive shell mode\n"; + std::cout << " --version, -v : prints the CP2K version and the revision number\n"; + std::cout << " --xml : dumps the whole CP2K input structure as a XML file\n"; + std::cout << " xml2htm generates a HTML manual from this XML file\n"; + std::cout << " --set, -E name=value : set the initial value of a preprocessor value\n " << std::endl; + } + } + + if (!command_line_error) { + // write the version string + if (print_version) { + if (default_para_env.is_source()) { + std::cout << " " << cp2k_version << "\n Source code revision " << compile_revision << "\n " << cp2k_flags() << "\n"; + compiler_options_string = compiler_options(); + std::cout << " compiler: " << compiler_version() << "\n"; + std::cout << " compiler options:"; + for (i = 0; i < (len(compiler_options_string) - 1)/68; i++) { + std::cout << compiler_options_string[i*68:fmin(len(compiler_options_string), (i+1)*68)] << "\n"; + } + std::cout << std::endl; + // delete compiler_options_string; + } + } + + // write the license + if (print_license) { + if (default_para_env.is_source()) { + print_cp2k_license(output_unit); + } + } + + if (xml) { + if (default_para_env.is_source()) { + write_xml_file(); + } + } + + create_cp2k_root_section(input_declaration); + + if (check) { + check_input(input_declaration, input_file_name, output_file_name, echo_input=echo_input, + ierr=&ierr, initial_variables=initial_variables); + + if (default_para_env.is_source()) { + if (ierr == 0) { + if (!echo_input) { + std::cout << "SUCCESS, the input could be parsed correctly.\n"; + std::cout << " This does not guarantee that this input is meaningful\n"; + std::cout << " or will run successfully" << std::endl; + } + } else { + std::cout << "ERROR, the input could *NOT* be parsed correctly."; + std::cout << " Please, check and correct it" << std::endl; + } + } + } + + if (shell_mode) { + launch_cp2k_shell(input_declaration) + } + + if (run_it || force_run) { + run_input(input_declaration, input_file_name, output_file_name, initial_variables) + } + + section_release(input_declaration); + } + } else { + std::cout << "initial setup (MPI ?) error" << std::endl; + } + + // and final cleanup + finalize_cp2k(finalize_mpi=true, ierr=&ierr); + // delete initial_variables; + assert(ierr == 0); + return 0; +} \ No newline at end of file diff --git a/src/start/python/README.md b/src/start/python/README.md new file mode 100644 index 0000000..da12fe7 --- /dev/null +++ b/src/start/python/README.md @@ -0,0 +1,30 @@ +# CP2K Python Bindings + +## Installation + +There is a target `py-cython-bindings` in the global `Makefile` to build the +Python bindings. The shared object can be found in: +`/lib///python` + +Only the Python headers and a NumPy installation are required. + +## Development + +To regenerate the C file from the `cp2k.pyx`, `Cython` is required and should +be called as follows: + +```sh +cd /src/start/python +cython cp2k.pyx +``` + +Unittests can be found in the `test/` directory. They must be run in separate +Python interpreter instances due to side effects in the library. + +## Known Issues + +* If libcp2k is built with MPI support, you may get an MPI initialization error + depending on your MPI implementation/configuration. In that case MPI must be + initialized first by using Mpi4py and the Fortran MPI communicator handler + must be passed down the CP2K via the respective `...comm` functions. + The reason for this is documented here: diff --git a/src/start/python/cp2k.pyx b/src/start/python/cp2k.pyx new file mode 100644 index 0000000..e06fe5c --- /dev/null +++ b/src/start/python/cp2k.pyx @@ -0,0 +1,161 @@ +# cython: language_level=2 +# vim: set ts=4 sw=4 tw=0 : + +import numpy as np + +from cpython.mem cimport PyMem_Malloc, PyMem_Free + +cdef extern from "../libcp2k.h": + ctypedef int force_env_t + + void cp2k_get_version(char* version_str, int str_length) + void cp2k_init() + void cp2k_init_without_mpi() + void cp2k_finalize() + void cp2k_finalize_without_mpi() + void cp2k_create_force_env(force_env_t* new_force_env, const char* input_file_path, const char* output_file_path) + void cp2k_create_force_env_comm(force_env_t* new_force_env, const char* input_file_path, const char* output_file_path, int mpi_comm) + void cp2k_destroy_force_env(force_env_t force_env) + void cp2k_set_positions(force_env_t force_env, const double* new_pos, int n_el) + void cp2k_set_velocities(force_env_t force_env, const double* new_vel, int n_el) + void cp2k_get_result(force_env_t force_env, const char* description, double* result, int n_el) + void cp2k_get_natom(force_env_t force_env, int* natom) + void cp2k_get_nparticle(force_env_t force_env, int* nparticle) + void cp2k_get_positions(force_env_t force_env, double* pos, int n_el) + void cp2k_get_forces(force_env_t force_env, double* force, int n_el) + void cp2k_get_potential_energy(force_env_t force_env, double* e_pot) + void cp2k_calc_energy_force(force_env_t force_env) + void cp2k_calc_energy(force_env_t force_env) + void cp2k_run_input(const char* input_file_path, const char* output_file_path) + void cp2k_run_input_comm(const char* input_file_path, const char* output_file_path, int mpi_comm) + +def get_version_string(): + n = 255 * sizeof(char) + + data = PyMem_Malloc(n) + if not data: + raise MemoryError() + + versionstr = '' + try: + cp2k_get_version(data, n) + versionstr = data.decode('UTF-8') + finally: + PyMem_Free(data) + + return versionstr + +def init(manage_mpi = True): + if manage_mpi: + cp2k_init() + else: + cp2k_init_without_mpi() + +def finalize(manage_mpi = True): + if manage_mpi: + cp2k_finalize() + else: + cp2k_finalize_without_mpi() + +def run_input(input_file_path, output_file_path = None, mpi_comm = None): + input_file_path = input_file_path.encode('UTF-8') + + if output_file_path is None: + output_file_path = u'__STD_OUT__'.encode('UTF-8') + else: + output_file_path = output_file_path.encode('UTF-8') + + if mpi_comm: + cp2k_run_input_comm(input_file_path, output_file_path, mpi_comm) + else: + cp2k_run_input(input_file_path, output_file_path) + +def create_force_env(input_file_path, output_file_path, mpi_comm = None): + cdef force_env_t fenv + if mpi_comm: + cp2k_create_force_env_comm(&fenv, input_file_path, output_file_path, mpi_comm) + else: + cp2k_create_force_env(&fenv, input_file_path, output_file_path) + +cdef class ForceEnvironment(object): + cdef force_env_t _force_env + + def __init__(self, input_file_path not None, output_file_path = None, mpi_comm = None): + input_file_path = input_file_path.encode('UTF-8') + + if output_file_path is None: + output_file_path = u'__STD_OUT__'.encode('UTF-8') + else: + output_file_path = output_file_path.encode('UTF-8') + + if mpi_comm: + cp2k_create_force_env_comm(&self._force_env, input_file_path, output_file_path, mpi_comm) + else: + cp2k_create_force_env(&self._force_env, input_file_path, output_file_path) + + def destroy(self): + cp2k_destroy_force_env(self._force_env) + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_value, traceback): + self.destroy() + + + property positions: + + def __get__(self): + positions = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] positions_view = positions + cp2k_get_positions(self._force_env, &positions_view[0], positions_view.shape[0]) + return positions + + def __set__(self, double[::1] positions not None): + if positions.shape[0] != 3*self.nparticle: + raise ValueError('the positions array must have exactly {} (3*nparticle) elements'.format(3*self.nparticle)) + + cp2k_set_positions(self._force_env, &positions[0], positions.shape[0]) + + def set_velocities(self, double[::1] velocities not None): + if velocities.shape[0] != 3*self.nparticle: + raise ValueError('the velocities array must have exactly {} (3*nparticle) elements'.format(3*self.nparticle)) + + cp2k_set_velocities(self._force_env, &velocities[0], velocities.shape[0]) + + property natom: + def __get__(self): + cdef int natom + cp2k_get_natom(self._force_env, &natom) + return natom + + property nparticle: + def __get__(self): + cdef int nparticle + cp2k_get_nparticle(self._force_env, &nparticle) + return nparticle + + def get_result(self, description): + results = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] results_view = results + cp2k_get_result(self._force_env, description.encode('UTF-8'), &results_view[0], results_view.shape[0]) + return results + + property forces: + def __get__(self): + forces = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] forces_view = forces + cp2k_get_forces(self._force_env, &forces_view[0], forces_view.shape[0]) + return forces + + def calc_energy_force(self): + cp2k_calc_energy_force(self._force_env) + + def calc_energy(self): + cp2k_calc_energy(self._force_env) + + property potential_energy: + def __get__(self): + cdef double e_pot + cp2k_get_potential_energy(self._force_env, &e_pot) + return e_pot diff --git a/src/subcell_types.c b/src/subcell_types.c new file mode 100644 index 0000000..120e046 --- /dev/null +++ b/src/subcell_types.c @@ -0,0 +1,211 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief subcell types and allocation routines +//> \par History +//> - Separated from qs_neighbor_lists (25.07.2010,jhu) +//> \author Matthias Krack +// ************************************************************************************************** +#include "./base/base_uses.h" + +void subcell_types() { + + USE cell_types, ONLY: cell_type,& + real_to_scaled,& + scaled_to_real + USE kinds, ONLY: dp + USE util, ONLY: sort + + IMPLICIT NONE + + PRIVATE + +// ************************************************************************************************** + TYPE subcell_type + INTEGER :: natom + REAL(KIND=dp), DIMENSION(3) :: s_max, s_min + INTEGER, DIMENSION(:), POINTER :: atom_list + REAL(KIND=dp), DIMENSION(3, 8) :: corners + END TYPE subcell_type + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'subcell_types' + + PUBLIC :: subcell_type, allocate_subcell, deallocate_subcell + PUBLIC :: reorder_atoms_subcell, give_ijk_subcell + +// ************************************************************************************************** + +CONTAINS + +// ************************************************************************************************** +//> \brief Allocate and initialize a subcell grid structure for the atomic neighbor search. +//> \param subcell ... +//> \param nsubcell ... +//> \param maxatom ... +//> \param cell ... +//> \date 12.06.2003 +//> \author MK +//> \version 1.0 +// ************************************************************************************************** + SUBROUTINE allocate_subcell(subcell, nsubcell, maxatom, cell) + + TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell + INTEGER, DIMENSION(3), INTENT(IN) :: nsubcell + INTEGER, INTENT(IN), OPTIONAL :: maxatom + TYPE(cell_type), OPTIONAL, POINTER :: cell + + INTEGER :: i, j, k, na, nb, nc + REAL(dp) :: a_max, a_min, b_max, b_min, c_max, & + c_min, delta_a, delta_b, delta_c + + na = nsubcell(1) + nb = nsubcell(2) + nc = nsubcell(3) + + ALLOCATE (subcell(na, nb, nc)) + + delta_a = 1.0_dp/REAL(na, dp) + delta_b = 1.0_dp/REAL(nb, dp) + delta_c = 1.0_dp/REAL(nc, dp) + + c_min = -0.5_dp + + DO k = 1, nc + c_max = c_min + delta_c + b_min = -0.5_dp + DO j = 1, nb + b_max = b_min + delta_b + a_min = -0.5_dp + DO i = 1, na + a_max = a_min + delta_a + subcell(i, j, k)%s_min(1) = a_min + subcell(i, j, k)%s_min(2) = b_min + subcell(i, j, k)%s_min(3) = c_min + subcell(i, j, k)%s_max(1) = a_max + subcell(i, j, k)%s_max(2) = b_max + subcell(i, j, k)%s_max(3) = c_max + subcell(i, j, k)%natom = 0 + IF (PRESENT(cell)) THEN + CALL scaled_to_real(subcell(i, j, k)%corners(:, 1), (/a_min, b_min, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 2), (/a_max, b_min, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 3), (/a_min, b_max, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 4), (/a_max, b_max, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 5), (/a_min, b_min, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 6), (/a_max, b_min, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 7), (/a_min, b_max, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 8), (/a_max, b_max, c_max/), cell) + END IF + IF (PRESENT(maxatom)) THEN + ALLOCATE (subcell(i, j, k)%atom_list(maxatom)) + END IF + a_min = a_max + END DO + b_min = b_max + END DO + c_min = c_max + END DO + + END SUBROUTINE allocate_subcell + +// ************************************************************************************************** +//> \brief Deallocate a subcell grid structure. +//> \param subcell ... +//> \date 16.06.2003 +//> \author MK +//> \version 1.0 +// ************************************************************************************************** + SUBROUTINE deallocate_subcell(subcell) + + TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell + + INTEGER :: i, j, k + + IF (ASSOCIATED(subcell)) THEN + + DO k = 1, SIZE(subcell, 3) + DO j = 1, SIZE(subcell, 2) + DO i = 1, SIZE(subcell, 1) + DEALLOCATE (subcell(i, j, k)%atom_list) + END DO + END DO + END DO + + DEALLOCATE (subcell) + ELSE + CPABORT("") + END IF + + END SUBROUTINE deallocate_subcell + +// ************************************************************************************************** +//> \brief ... +//> \param atom_list ... +//> \param kind_of ... +//> \param work ... +//> \par History +//> 08.2006 created [tlaino] +//> \author Teodoro Laino +// ************************************************************************************************** + SUBROUTINE reorder_atoms_subcell(atom_list, kind_of, work) + // work needs to be dimensioned 3xSIZE(atom_list) + INTEGER, DIMENSION(:), POINTER :: atom_list + INTEGER, DIMENSION(:), INTENT(IN) :: kind_of + INTEGER, DIMENSION(:) :: work + + INTEGER :: i, i0, i1, i2, j0, j1, j2 + + i0 = 1 + j0 = SIZE(atom_list) + i1 = j0 + 1 + j1 = 2*j0 + i2 = j1 + 1 + j2 = 3*j0 + // Sort kind + DO i = 1, SIZE(atom_list) + work(i0 + i - 1) = kind_of(atom_list(i)) + END DO + CALL sort(work(i0:j0), SIZE(atom_list), work(i1:j1)) + work(i2:j2) = atom_list + DO i = 1, SIZE(atom_list) + atom_list(i) = work(i2 + work(i1 + i - 1) - 1) + END DO + END SUBROUTINE reorder_atoms_subcell + +// ************************************************************************************************** +//> \brief ... +//> \param r ... +//> \param i ... +//> \param j ... +//> \param k ... +//> \param cell ... +//> \param nsubcell ... +//> \par History +//> 08.2006 created [tlaino] +//> \author Teodoro Laino +// ************************************************************************************************** + SUBROUTINE give_ijk_subcell(r, i, j, k, cell, nsubcell) + REAL(KIND=dp) :: r(3) + INTEGER, INTENT(OUT) :: i, j, k + TYPE(cell_type), POINTER :: cell + INTEGER, DIMENSION(3), INTENT(IN) :: nsubcell + + REAL(KIND=dp) :: r_pbc(3), s(3), s_pbc(3) + + r_pbc = r + CALL real_to_scaled(s_pbc, r_pbc, cell) + s(:) = s_pbc + 0.5_dp + i = INT(s(1)*REAL(nsubcell(1), KIND=dp)) + 1 + j = INT(s(2)*REAL(nsubcell(2), KIND=dp)) + 1 + k = INT(s(3)*REAL(nsubcell(3), KIND=dp)) + 1 + i = MIN(MAX(i, 1), nsubcell(1)) + j = MIN(MAX(j, 1), nsubcell(2)) + k = MIN(MAX(k, 1), nsubcell(3)) + + END SUBROUTINE give_ijk_subcell + +}