Added new files found in filesystem

.gitignore

@@ -110,6 +110,9 @@ auto-save-list
 tramp
 .\#*
 
+### VS Code ###
+.vscode
+
 # Org-mode
 .org-id-locations
 *_archive

src/admm_dm_methods.cpp

@@ -0,0 +1,533 @@
+/*--------------------------------------------------------------------------------------------------
+    CP2K: A general program to perform molecular dynamics simulations
+    Copyright 2000-2023 CP2K developers group <https://cp2k.org>
+
+    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);
+
+}

src/admm_utils.c

@@ -0,0 +1,179 @@
+//--------------------------------------------------------------------------------------------------//
+//   CP2K: A general program to perform molecular dynamics simulations                              //
+//   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   //
+//                                                                                                  //
+//   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
+
+}

src/base/PACKAGE

@@ -0,0 +1,4 @@
+{
+    "description": "base routines needed to abstract away some machine/compiler dependent functionality",
+    "requires": [],
+}

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

src/base/base_hooks.cpp

@@ -0,0 +1,238 @@
+/*--------------------------------------------------------------------------------------------------
+    CP2K: A general program to perform molecular dynamics simulations
+    Copyright 2000-2023 CP2K developers group <https://cp2k.org>
+
+    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;
+}

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

src/base/kinds.cpp

@@ -0,0 +1,51 @@
+/*------------------------------------------------------------------------------------------------
+   CP2K: A general program to perform molecular dynamics simulations
+   Copyright 2000-2023 CP2K developers group <https://cp2k.org>
+
+   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')
+
+}
+
+
+}

src/base/machine.c

@@ -0,0 +1,117 @@
+//------------------------------------------------------------------------------------------------//
+//  CP2K: A general program to perform molecular dynamics simulations                             //
+//  Copyright 2000-2022 CP2K developers group <https://cp2k.org>                                  //
+//                                                                                                //
+//  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 <omp.h>
+
+#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 = "";
+
+
+}
+
+
+
+
+
+

src/base/machine_cpuid.c

@@ -0,0 +1,42 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2021 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  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

src/base/machine_cpuid.h

@@ -0,0 +1,12 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2021 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  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

src/cp2k_info.cpp

@@ -0,0 +1,357 @@
+/*--------------------------------------------------------------------------------------------------
+    CP2K: A general program to perform molecular dynamics simulations
+    Copyright 2000-2023 CP2K developers group <https://cp2k.org>
+
+    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 <string>
+#include <fstream>
+#include <algorithm>
+
+#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;
+
+}

src/header.cpp

@@ -0,0 +1,14 @@
+//------------------------------------------------------------------------------------------------//
+//  CP2K: A general program to perform molecular dynamics simulations                             //
+//  Copyright 2000-2022 CP2K developers group <https://cp2k.org>                                  //
+//                                                                                                //
+//  SPDX-License-Identifier: GPL-2.0-or-later                                                     //
+//------------------------------------------------------------------------------------------------//
+
+// *************************************************************************************************
+//> \par History
+//>      none
+//> \author APSI & CJM & JGH
+// *************************************************************************************************
+
+#include "base/base_uses.cpp"

src/motion/simpar_methods.c

@@ -0,0 +1,62 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2023 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  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) {
+
+        }
+    }
+
+}

src/mpiwrap/PACKAGE

@@ -0,0 +1,5 @@
+{
+    "description": "wrappers of the mpi routines",
+    "requires": ["../base"],
+    "implicit": "MPI_.*",
+}

src/mpiwrap/message_passing.h

@@ -0,0 +1,60 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2023 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  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 <stdbool.h>
+
+#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;
+    

src/qs_ot_eigensolver.c

@@ -0,0 +1,77 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2023 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+
+
+#include <stdbool.h>
+
+#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;
+
+
+
+
+}

src/qs_tenosrs.c

@@ -0,0 +1,112 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2023 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+/*******************************************************************************
+ * \brief Utility methods to build 3-center integral tensors of various types.
+*******************************************************************************/
+
+#include <stdbool.h>
+
+#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);
+
+
+
+
+
+}
+

src/sockets.c

@@ -0,0 +1,177 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2023 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  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 <math.h>
+#include <netdb.h>
+#include <netinet/in.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/socket.h>
+#include <sys/types.h>
+#include <sys/un.h>
+#include <time.h>
+#include <unistd.h>
+
+/*******************************************************************************
+ * \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

src/start/cp2k.c

@@ -0,0 +1,103 @@
+//--------------------------------------------------------------------------------------------------//
+//   CP2K: A general program to perform molecular dynamics simulations                              //
+//   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   //
+//                                                                                                  //
+//   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 <stdbool.h>
+
+#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);
+}

src/start/cp2k.cpp

@@ -0,0 +1,337 @@
+//--------------------------------------------------------------------------------------------------//
+//   CP2K: A general program to perform molecular dynamics simulations                              //
+//   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   //
+//                                                                                                  //
+//   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 <omp.h>
+#include <cmath>
+#include <string>
+#include <vector>
+#include <cassert>
+#include <iostream>
+
+#include "base_uses.h"
+
+int main(int argc, char** argv) {
+
+    std::string input_file_name, output_file_name, arg_att, command;
+    std::vector<std::string> 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] <input_file>\n";
+                std::cout << std::string(l, ' ') << " [-mpi-mapping|--mpi-mapping] <method>\n";
+                std::cout << std::string(l, ' ') << " [-o] <output_file>\n";
+                std::cout << std::string(l, ' ') << " [-r|-run] [-s|--shell] [--xml]\n";
+
+                std::cout << "\n starts the CP2K program, see <https://www.cp2k.org/>\n";
+                std::cout << "\n The easiest way is " << command << " <input_file>\n";
+                std::cout << "\n The following options can be used:\n";
+                std::cout << "\n  -i <input_file>      : 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 <output_file>     : 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 <input_file>\n";
+                std::cout << "  --echo, -e           : echoes the <input_file>, 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;
+}

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:
+`<CP2K_SOURCE_DIR>/lib/<ARCH>/<VERSION>/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 <CP2K_SOURCE_DIR>/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: <https://github.com/jhedev/mpi_python>

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 = <char *>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

src/subcell_types.c

@@ -0,0 +1,211 @@
+//--------------------------------------------------------------------------------------------------//
+//   CP2K: A general program to perform molecular dynamics simulations                              //
+//   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   //
+//                                                                                                  //
+//   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
+
+}