diff --git a/src/admm_utils.c b/src/admm_utils.c new file mode 100644 index 0000000..56af759 --- /dev/null +++ b/src/admm_utils.c @@ -0,0 +1,179 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Contains methods used in the context of density fitting +//> \par History +//> 04.2008 created [Manuel Guidon] +//> 02.2013 moved from admm_methods +//> \author Manuel Guidon +// ************************************************************************************************** +#include "./base/base_uses.h" + +void admm_utils() { + USE admm_types, ONLY: admm_type + USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr + USE cp_gemm_interface, ONLY: cp_gemm + USE dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_create,& + dbcsr_deallocate_matrix,& + dbcsr_set,& + dbcsr_type,& + dbcsr_type_symmetric + USE input_constants, ONLY: do_admm_purify_cauchy,& + do_admm_purify_cauchy_subspace,& + do_admm_purify_mo_diag,& + do_admm_purify_mo_no_diag,& + do_admm_purify_none + USE kinds, ONLY: dp + + IMPLICIT NONE + PRIVATE + + PUBLIC :: admm_correct_for_eigenvalues, & + admm_uncorrect_for_eigenvalues + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_utils' + +//*** + +CONTAINS + +// ************************************************************************************************** +//> \brief ... +//> \param ispin ... +//> \param admm_env ... +//> \param ks_matrix ... +// ************************************************************************************************** + SUBROUTINE admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix) + INTEGER, INTENT(IN) :: ispin + TYPE(admm_type), POINTER :: admm_env + TYPE(dbcsr_type), POINTER :: ks_matrix + + INTEGER :: nao_aux_fit, nao_orb + TYPE(dbcsr_type), POINTER :: work + + nao_aux_fit = admm_env%nao_aux_fit + nao_orb = admm_env%nao_orb + + IF (.NOT. admm_env%block_dm) THEN + SELECT CASE (admm_env%purification_method) + CASE (do_admm_purify_cauchy_subspace) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + + // ** calculate A^T*H_tilde*A + CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%K(ispin)%matrix, admm_env%A, 0.0_dp, & + admm_env%work_aux_orb) + CALL cp_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, & + admm_env%H_corr(ispin)%matrix) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_diag) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + // ** calculate A^T*H_tilde*A + CALL cp_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%K(ispin)%matrix, admm_env%A, 0.0_dp, & + admm_env%work_aux_orb) + CALL cp_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, & + 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, & + admm_env%H_corr(ispin)%matrix) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_cauchy) + // do nothing + END SELECT + END IF + + END SUBROUTINE admm_correct_for_eigenvalues + +// ************************************************************************************************** +//> \brief ... +//> \param ispin ... +//> \param admm_env ... +//> \param ks_matrix ... +// ************************************************************************************************** + SUBROUTINE admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix) + INTEGER, INTENT(IN) :: ispin + TYPE(admm_type), POINTER :: admm_env + TYPE(dbcsr_type), POINTER :: ks_matrix + + INTEGER :: nao_aux_fit, nao_orb + TYPE(dbcsr_type), POINTER :: work + + nao_aux_fit = admm_env%nao_aux_fit + nao_orb = admm_env%nao_orb + + IF (.NOT. admm_env%block_dm) THEN + SELECT CASE (admm_env%purification_method) + CASE (do_admm_purify_cauchy_subspace) + //* remove what has been added and add the correction + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_set(work, 0.0_dp) + CALL copy_fm_to_dbcsr(admm_env%ks_to_be_merged(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, 1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_diag) + NULLIFY (work) + ALLOCATE (work) + CALL dbcsr_create(work, template=ks_matrix, name='work', matrix_type=dbcsr_type_symmetric) + + CALL dbcsr_copy(work, ks_matrix) + CALL dbcsr_set(work, 0.0_dp) + + CALL copy_fm_to_dbcsr(admm_env%H_corr(ispin)%matrix, work, keep_sparsity=.TRUE.) + + CALL dbcsr_add(ks_matrix, work, 1.0_dp, -1.0_dp) + CALL dbcsr_deallocate_matrix(work) + + CASE (do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_cauchy) + // do nothing + END SELECT + END IF + END SUBROUTINE admm_uncorrect_for_eigenvalues + +} diff --git a/src/base/PACKAGE b/src/base/PACKAGE new file mode 100644 index 0000000..397ef51 --- /dev/null +++ b/src/base/PACKAGE @@ -0,0 +1,4 @@ +{ + "description": "base routines needed to abstract away some machine/compiler dependent functionality", + "requires": [], +} diff --git a/src/base/base_hooks.cpp b/src/base/base_hooks.cpp index 602befe..8b099f0 100644 --- a/src/base/base_hooks.cpp +++ b/src/base/base_hooks.cpp @@ -69,14 +69,12 @@ public: } -! -************************************************************************************************** -!> \brief Terminate the program -!> \param location ... -!> \param message ... -!> \author Ole Schuett -! -************************************************************************************************** +/*************************************************************************************************** + \brief Terminate the program + \param location ... + \param message ... + \author Ole Schuett +***************************************************************************************************/ SUBROUTINE cp_abort(location, message) CHARACTER(len=*), INTENT(in) :: location, message @@ -91,14 +89,12 @@ public: 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 @@ -110,14 +106,12 @@ public: 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 @@ -129,14 +123,12 @@ public: 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 @@ -148,13 +140,11 @@ public: END IF END SUBROUTINE timeset -! -************************************************************************************************** +! ************************************************************************************************** !> \brief Stop timer !> \param handle ... !> \author Ole Schuett -! -************************************************************************************************** +! ************************************************************************************************** SUBROUTINE timestop(handle) INTEGER, INTENT(IN) :: handle @@ -166,14 +156,12 @@ public: 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 @@ -183,15 +171,13 @@ public: 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 @@ -202,15 +188,13 @@ public: 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 @@ -219,15 +203,13 @@ public: 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 diff --git a/src/base/base_uses.cpp b/src/base/base_uses.cpp index 8f02b98..e0b5a4c 100644 --- a/src/base/base_uses.cpp +++ b/src/base/base_uses.cpp @@ -12,8 +12,7 @@ int cp__a, cp__b, cp__w, cp__h, cp__l, cp_abort, cp_warn, cp_hint, timeset, time // Check for OpenMP early on - ideally before the compiler fails with a cryptic message. #if !defined(_OPENMP) - "OpenMP is required. Please add the corresponding flag (eg. -fopenmp for GFortran) to your -Fortran compiler flags." + 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. diff --git a/src/base/kinds.cpp b/src/base/kinds.cpp index 8dbbd49..cb77ae3 100644 --- a/src/base/kinds.cpp +++ b/src/base/kinds.cpp @@ -15,8 +15,7 @@ public: Adapted by JGH for Cp2k \author Matthias Krack ***************************************************************************************************/ -void print_kind_info(int iw) -{ +void print_kind_info(int iw) { WRITE (iw, '( /, T2, A )') 'DATA TYPE INFORMATION:' diff --git a/src/base/machine.c b/src/base/machine.c new file mode 100644 index 0000000..4aae8e5 --- /dev/null +++ b/src/base/machine.c @@ -0,0 +1,117 @@ +//------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2022 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//------------------------------------------------------------------------------------------------// + +// ************************************************************************************************* +//> \brief Machine interface based on Fortran 2003 and POSIX +//> \par History +//> JGH (05.07.2001) : added G95 interface +//> - m_flush added (12.06.2002,MK) +//> - Missing print_memory added (24.09.2002,MK) +//> - Migrate to generic implementation based on F2003 + POSIX (2014, Ole Schuett) +//> \author APSI, JGH, Ole Schuett +// ************************************************************************************************* + +#include + +#include "machine_cpuid.h" + +#if defined(__LIBXSMM) + //do something +#endif + + // Except for some error handling code, all code should + // get a unit number from the print keys or from the logger, in order + // to guarantee correct output behavior, + // for example in farming or path integral runs + // default_input_unit should never be used + // but we need to know what it is, as we should not try to open it for output + + int default_output_unit = output_unit; + int default_input_unit = input_unit; + + // Enumerates the target architectures or instruction set extensions. + // A feature is present if within range for the respective architecture. + // For example, to check for MACHINE_X86_AVX the following is true: + // MACHINE_X86_AVX <= m_cpuid() and MACHINE_X86 >= m_cpuid(). + // For example, to check for MACHINE_ARM_SOME the following is true: + // MACHINE_ARM_SOME <= m_cpuid() and MACHINE_ARM >= m_cpuid(). + + int MACHINE_CPU_GENERIC = CP_MACHINE_CPU_GENERIC; + int MACHINE_X86_SSE4 = CP_MACHINE_X86_SSE4; + int MACHINE_X86_AVX = CP_MACHINE_X86_AVX; + int MACHINE_X86_AVX2 = CP_MACHINE_X86_AVX2; + int MACHINE_X86_AVX512 = CP_MACHINE_X86_AVX512; + int MACHINE_X86 = MACHINE_X86_AVX512; // marks end of range + + // other arch to be added as needed e.g., + //MACHINE_ARM_SOME = 2000 + //MACHINE_ARM_ELSE = 2001 + //MACHINE_ARM = MACHINE_ARM_ELSE + //MACHINE_PWR_???? = 3000 + + // Flushing is enabled by default because without it crash reports can get lost. + // For performance reasons it can be disabled via the input in &GLOBAL. + bool flush_should_flush = true; + + int m_memory_max = 0; + +// ************************************************************************************************* +//> \brief flushes units if the &GLOBAL flag is set accordingly +//> \param lunit ... +//> \par History +//> 10.2008 created [Joost VandeVondele] +//> \note +//> flushing might degrade performance significantly (30% and more) +// ************************************************************************************************* + +void m_flush(std::ofstream &lunit) { + if (flush_should_flush) {std::flush(lunit);} +} + +// ************************************************************************************************* +//> \brief returns time from a real-time clock, protected against rolling +//> early/easily +//> \return ... +//> \par History +//> 03.2006 created [Joost VandeVondele] +//> \note +//> same implementation for all machines. +//> might still roll, if not called multiple times per count_max/count_rate +// ************************************************************************************************* +double m_walltime() { + double wt; + + #ifdef __LIBXSMM + wt = libxsmm_timer_duration(0_int_8, libxsmm_timer_tick()); + #else + wt = omp_get_wtime(); + #endif + + return wt; +} + +// ************************************************************************************************* +//> \brief reads /proc/cpuinfo if it exists (i.e. Linux) to return relevant info +//> \param model_name as obtained from the 'model name' field, UNKNOWN otherwise +// ************************************************************************************************* +void m_cpuinfo(std::string model_name) { + int bufferlen = 2048; + + char buffer[bufferlen]; + int i, icol, iline, imod, stat; + + model_name = "UNKNOWN"; + buffer = ""; + + +} + + + + + + diff --git a/src/base/machine_cpuid.c b/src/base/machine_cpuid.c new file mode 100644 index 0000000..b4c04f7 --- /dev/null +++ b/src/base/machine_cpuid.c @@ -0,0 +1,42 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2021 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/* shared between C and Fortran */ +#include "machine_cpuid.h" + +#if defined(__cplusplus) +extern "C" { +#endif + +/******************************************************************************* + * \brief This routine determines the CPUID according to the given compiler + * flags (expected to be similar to Fortran). Similar to other Fortran + * compilers, "gfortran -E -dM -mavx - < /dev/null | grep AVX" defines a + * variety of predefined macros (also similar to C). However, with a + * Fortran translation unit only a subset of these definitions disappears + * ("gfortran -E -dM -mavx my.F | grep AVX") + * hence an implementation in C is used. + ******************************************************************************/ +int m_cpuid_static(void); /* avoid pedantic warning about missing prototype */ +int m_cpuid_static(void) { +#if (__AVX512F__ && __AVX512CD__ && __AVX2__ && __FMA__ && __AVX__ && \ + __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX512; +#elif (__AVX2__ && __FMA__ && __AVX__ && __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX2; +#elif (__AVX__ && __SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_AVX; +#elif (__SSE4_2__ && __SSE4_1__ && __SSE3__) + return CP_MACHINE_X86_SSE4; +#else + return CP_MACHINE_CPU_GENERIC; +#endif +} + +#if defined(__cplusplus) +} +#endif diff --git a/src/base/machine_cpuid.h b/src/base/machine_cpuid.h new file mode 100644 index 0000000..968c41b --- /dev/null +++ b/src/base/machine_cpuid.h @@ -0,0 +1,12 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2021 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +#define CP_MACHINE_CPU_GENERIC 0 +#define CP_MACHINE_X86_SSE4 1000 +#define CP_MACHINE_X86_AVX 1001 +#define CP_MACHINE_X86_AVX2 1002 +#define CP_MACHINE_X86_AVX512 1003 diff --git a/src/cp2k_info.cpp b/src/cp2k_info.cpp index 9173900..7de5745 100644 --- a/src/cp2k_info.cpp +++ b/src/cp2k_info.cpp @@ -15,6 +15,8 @@ ***************************************************************************************************/ #include +#include +#include #include "./base/base.h" @@ -30,42 +32,40 @@ public: void write_restart_header(); int cp2k_version, cp2k_year, cp2k_home; int compile_arch, compile_date, compile_host, compile_revision; - -#if defined(__COMPILE_REVISION) - char *compile_revision = __COMPILE_REVISION; +#ifdef __COMPILE_REVISION + std::string compile_revision = __COMPILE_REVISION; #else - char *compile_revision = "unknown"; + std::string compile_revision = "unknown"; #endif - char *cp2k_version = "CP2K version 2023.2 (Development Version)"; - char *cp2k_year = "2023"; - char *cp2k_home = "https://www.cp2k.org/"; +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 -#if defined(__COMPILE_ARCH) - char *compile_arch = __COMPILE_ARCH; +// compile time information +#ifdef __COMPILE_ARCH + std::string compile_arch = __COMPILE_ARCH; #else - char *compile_arch = "unknown: -D__COMPILE_ARCH=?"; + std::string compile_arch = "unknown: -D__COMPILE_ARCH=?"; #endif -#if defined(__COMPILE_DATE) - char *compile_date = __COMPILE_DATE; +#ifdef __COMPILE_DATE + std::string compile_date = __COMPILE_DATE; #else - char *compile_date = "unknown: -D__COMPILE_DATE=?"; + std::string compile_date = "unknown: -D__COMPILE_DATE=?"; #endif -#if defined(__COMPILE_HOST) - char *compile_host = __COMPILE_HOST; +#ifdef __COMPILE_HOST + std::string compile_host = __COMPILE_HOST; #else - char *compile_host = "unknown: -D__COMPILE_HOST=?"; + std::string compile_host = "unknown: -D__COMPILE_HOST=?"; #endif - // Local runtime informations - char r_datx[26]; - char r_cwd[default_path_length]; - char r_host_name[default_string_length], r_user_name[default_string_length]; - int r_pid; +// local runtime informations +std::string r_datx, r_cwd, r_host_name, r_user_name; +int r_pid; +std::string moduleN = "cp2k_info"; }; @@ -74,186 +74,186 @@ public: All new flags should be added here (and be unique grep-able) \return ... ***************************************************************************************************/ -char* cp2k_info::cp2k_flags() -{ - char flags[10*default_string_length]; - char tmp_str[default_string_length]; +std::string cp2k_info::cp2k_flags() { + std::string flags, tmp_str; - std::strncpy(flags,"cp2kflags:", sizeof(flags)); + flags = "cp2kflags:"; // Ensure that tmp_str is used to silence compiler warnings - std::strncpy(tmp_str, "", sizeof(tmp_str); - flags = TRIM(flags)//TRIM(tmp_str) + tmp_str = ""; + flags = rtrim(flags)+tmp_str; -#if defined(NDEBUG); - flags = TRIM(flags)//" ndebug" -#endif - flags = TRIM(flags)//" omp" -#if defined(__LIBINT) - flags = TRIM(flags)//" libint" -#endif -#if defined(__FFTW3) - flags = TRIM(flags)//" fftw3" -#endif -#if defined(__FFTW3_MKL) - flags = TRIM(flags)//" fftw3_mkl" -#endif -#if defined(__LIBXC) - flags = TRIM(flags)//" libxc" -#endif -#if defined(__LIBPEXSI) - flags = TRIM(flags)//" pexsi" -#endif -#if defined(__ELPA) - flags = TRIM(flags)//" elpa" -#endif -#if defined(__ELPA_NVIDIA_GPU) - flags = TRIM(flags)//" elpa_nvidia_gpu" -#endif -#if defined(__ELPA_AMD_GPU) - flags = TRIM(flags)//" elpa_amd_gpu" -#endif -#if defined(__ELPA_INTEL_GPU) - flags = TRIM(flags)//" elpa_intel_gpu" -#endif -#if defined(__parallel) - flags = TRIM(flags)//" parallel" -#endif -#if defined(__MPI_F08) - flags = TRIM(flags)//" mpi_f08" -#endif -#if defined(__SCALAPACK) - flags = TRIM(flags)//" scalapack" -#endif -#if defined(__COSMA) - flags = TRIM(flags)//" cosma" -#endif + #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 = TRIM(flags)//" quip" -#endif + #if defined(__QUIP) + flags = rtrim(flags)+" quip" + #endif -#if defined(__HAS_PATCHED_CUFFT_70) - flags = TRIM(flags)//" patched_cufft_70" -#endif + #if defined(__HAS_PATCHED_CUFFT_70) + flags = rtrim(flags)+" patched_cufft_70" + #endif -#if defined(__PW_FPGA) - flags = TRIM(flags)//" pw_fpga" -#endif -#if defined(__PW_FPGA_SP) - flags = TRIM(flags)//" pw_fpga_sp" -#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 = TRIM(flags)//" xsmm" -#endif + #if defined(__LIBXSMM) + flags = rtrim(flags)+" xsmm" + #endif -#if defined(__CRAY_PM_ACCEL_ENERGY) - flags = TRIM(flags)//" cray_pm_accel_energy" -#endif -#if defined(__CRAY_PM_ENERGY) - flags = TRIM(flags)//" cray_pm_energy" -#endif -#if defined(__CRAY_PM_FAKE_ENERGY) - flags = TRIM(flags)//" cray_pm_fake_energy" -#endif -#if defined(__DBCSR_ACC) - flags = TRIM(flags)//" dbcsr_acc" -#endif -#if defined(__MAX_CONTR) - CALL integer_to_string(__MAX_CONTR, tmp_str) - flags = TRIM(flags)//" max_contr="//TRIM(tmp_str) -#endif -#if defined(__NO_SOCKETS) - flags = TRIM(flags)//" no_sockets" -#endif -#if defined(__NO_MPI_THREAD_SUPPORT_CHECK) - flags = TRIM(flags)//" no_mpi_thread_support_check" -#endif -#if defined(__NO_STATM_ACCESS) - flags = TRIM(flags)//" no_statm_access" -#endif -#if defined(__MINGW) - flags = TRIM(flags)//" mingw" -#endif -#if defined(__PW_CUDA_NO_HOSTALLOC) - flags = TRIM(flags)//" pw_cuda_no_hostalloc" -#endif -#if defined(__STATM_RESIDENT) - flags = TRIM(flags)//" statm_resident" -#endif -#if defined(__STATM_TOTAL) - flags = TRIM(flags)//" statm_total" -#endif -#if defined(__PLUMED2) - flags = TRIM(flags)//" plumed2" -#endif -#if defined(__HAS_IEEE_EXCEPTIONS) - flags = TRIM(flags)//" has_ieee_exceptions" -#endif -#if defined(__NO_ABORT) - flags = TRIM(flags)//" no_abort" -#endif -#if defined(__SPGLIB) - flags = TRIM(flags)//" spglib" -#endif -#if defined(__ACCELERATE) - flags = TRIM(flags)//" accelerate" -#endif -#if defined(__MKL) - flags = TRIM(flags)//" mkl" -#endif -#if defined(__SIRIUS) - flags = TRIM(flags)//" sirius" -#endif -#if defined(__CHECK_DIAG) - flags = TRIM(flags)//" check_diag" -#endif -#if defined(__LIBVORI) - flags = TRIM(flags)//" libvori" - flags = TRIM(flags)//" libbqb" -#endif -#if defined(__LIBMAXWELL) - flags = TRIM(flags)//" libmaxwell" -#endif -#if defined(__LIBTORCH) - flags = TRIM(flags)//" libtorch" -#endif -#if defined(__OFFLOAD_CUDA) - flags = TRIM(flags)//" offload_cuda" -#endif -#if defined(__OFFLOAD_HIP) - flags = TRIM(flags)//" offload_hip" -#endif -#if defined(__NO_OFFLOAD_GRID) - flags = TRIM(flags)//" no_offload_grid" -#endif -#if defined(__NO_OFFLOAD_DBM) - flags = TRIM(flags)//" no_offload_dbm" -#endif -#if defined(__NO_OFFLOAD_PW) - flags = TRIM(flags)//" no_offload_pw" -#endif -#if defined(__OFFLOAD_PROFILING) - flags = TRIM(flags)//" offload_profiling" -#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 = TRIM(flags)//" spla_gemm_offloading" -#endif + #if defined(__SPLA) && defined(__OFFLOAD_GEMM) + flags = rtrim(flags)+" spla_gemm_offloading" + #endif -#if defined(__CUSOLVERMP) - flags = TRIM(flags)//" cusolvermp" -#endif + #if defined(__CUSOLVERMP) + flags = rtrim(flags)+" cusolvermp" + #endif -#if defined(__LIBVDWXC) - flags = TRIM(flags)//" libvdwxc" -#endif + #if defined(__LIBVDWXC) + flags = rtrim(flags)+" libvdwxc" + #endif -#if defined(__HDF5) - flags = TRIM(flags)//" hdf5" -#endif + #if defined(__HDF5) + flags = rtrim(flags)+" hdf5" + #endif return flags; } @@ -262,30 +262,40 @@ char* cp2k_info::cp2k_flags() \brief ... \param iunit ... ***************************************************************************************************/ -void cp2k_info::print_cp2k_license(int iunit) +void cp2k_info::print_cp2k_license(std::ofstream iunit) { - WRITE (UNIT=iunit, FMT="(T2,A)") & - "******************************************************************************", & - "* *", & - "* CP2K: A general program to perform molecular dynamics simulations *", & - "* Copyright (C) 2000-2023 CP2K developer group *", & - "* *", & - "* 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. *", & - "* *", & - "* 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. *", & - "* *", & - "* You should have received a copy of the GNU General Public License *", & - "* along with this program. If not, see . *", & - "* *", & - "******************************************************************************"; - + 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; } /*************************************************************************************************** @@ -299,11 +309,11 @@ void cp2k_info::get_runtime_info() r_user_name = ""; r_pid = -1; - m_getpid(r_pid) - CALL m_getlog(r_user_name) - CALL m_hostnm(r_host_name) - CALL m_datum(r_datx) - CALL m_getcwd(r_cwd) + m_getpid(&r_pid); + m_getlog(&r_user_name); + m_hostnm(&r_host_name); + m_datum(&r_datx); + m_getcwd(&r_cwd); } @@ -314,28 +324,34 @@ void cp2k_info::get_runtime_info() 01.2008 [created] - Split from write_restart \author Teodoro Laino - University of Zurich - 01.2008 ***************************************************************************************************/ -void cp2k_info::write_restart_header(int iunit) +void cp2k_info::write_restart_header(std::ofstream &iunit) { - char cwd[256], datx[256]; + char cwd[255], datx[255]; - m_datum(datx); - m_getcwd(cwd); + m_datum(&datx); + m_getcwd(&cwd); - WRITE (UNIT=iunit, FMT="(T2,A)") "# Version information for this restart file " - WRITE (UNIT=iunit, FMT="(T2,A)") "# current date "//TRIM(datx) - WRITE (UNIT=iunit, FMT="(T2,A)") "# current working dir "//TRIM(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; - WRITE (UNIT=iunit, FMT="(T2,A,T31,A50)") & - "# Program compiled at", & - ADJUSTR(compile_date(1:MIN(50, LEN(compile_date)))) - WRITE (UNIT=iunit, FMT="(T2,A,T31,A50)") & - "# Program compiled on", & - ADJUSTR(compile_host(1:MIN(50, LEN(compile_host)))) - WRITE (UNIT=iunit, FMT="(T2,A,T31,A50)") & - "# Program compiled for", & - ADJUSTR(compile_arch(1:MIN(50, LEN(compile_arch)))) - WRITE (UNIT=iunit, FMT="(T2,A,T31,A50)") & - "# Source code revision number", & - ADJUSTR(compile_revision) + iunit << " # Program compiled at "; + iunit.width(50); + if (compile_date.size() > 50) {compile_date.resize(50);} + iunit << std::right << compile_date; + + iunit << " # Program compiled on "; + iunit.width(50); + if (compile_host.size() > 50) {compile_host.resize(50);} + iunit << std::right << compile_host; + + iunit << " # Program compiled for "; + iunit.width(50); + if (compile_arch.size() > 50) {compile_arch.resize(50);} + iunit << std::right << compile_arch; + + iunit << " # Source code revision number"; + iunit.width(50); + iunit << std::right << compile_revision; } diff --git a/src/header.cpp b/src/header.cpp new file mode 100644 index 0000000..67c7535 --- /dev/null +++ b/src/header.cpp @@ -0,0 +1,14 @@ +//------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2022 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//------------------------------------------------------------------------------------------------// + +// ************************************************************************************************* +//> \par History +//> none +//> \author APSI & CJM & JGH +// ************************************************************************************************* + +#include "base/base_uses.cpp" diff --git a/src/motion/simpar_methods.c b/src/motion/simpar_methods.c new file mode 100644 index 0000000..1d966d8 --- /dev/null +++ b/src/motion/simpar_methods.c @@ -0,0 +1,62 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Methods for storing MD parameters type + * \author CJM + * \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 + * reorganization of the original routines/modules + ******************************************************************************/ + +#include "../base/base_uses.h" + +/******************************************************************************* + * \brief Reads the MD section and setup the simulation parameters type + * \param simpar ... + * \param motion_section ... + * \param md_section ... + * \author Teodoro Laino + ******************************************************************************/ + +void read_md_section(struct simpar_type *simpar, struct section_vals_type *motion_section, + struct section_vals_type *md_section) { + + char* filename[default_path_length]; + int iprint, iw; + double tmp_r1, tmp_r2, tmp_r3; + struct cp_logger_type *logger; + struct enumeration_type *enumer; + struct keyword_type *keyword; + struct section_type *section; + struct section_vals_type *print_key; + + // free(logger,print_key,enumer,keyword,section); + + logger = &cp_get_default_logger(); + iw = cp_print_key_uint_nr(logger, md_section, "PRINT.PROGRAM_RUN_INFO", ".log"); + + read_md_low(simpar, motion_section, md_section); + if (iw > 0) + printf("%i",iw); + + // Begin setup Langevin dynamics + if (simpar.ensemble == langevin_ensemble) { + cite_reference(Ricci2003); + if (simpar.noisy_gamma > 0.0) + cite_reference(Kuhne2007); + + if (simpar.shadow_gamma > 0.0) + cite_reference(Rengaraj2020); + + // Normalization factor using a normal Gaussian random number distribution + simpar.var_w = 2.0 * simpar.temp_ext * simpar.dt * (simpar.gamma + simpar.noisy_gamma); + if (iw > 0) { + + } + } + +} \ No newline at end of file diff --git a/src/mpiwrap/PACKAGE b/src/mpiwrap/PACKAGE new file mode 100644 index 0000000..7ddf13d --- /dev/null +++ b/src/mpiwrap/PACKAGE @@ -0,0 +1,5 @@ +{ + "description": "wrappers of the mpi routines", + "requires": ["../base"], + "implicit": "MPI_.*", +} \ No newline at end of file diff --git a/src/mpiwrap/message_passing.c b/src/mpiwrap/message_passing.c new file mode 100644 index 0000000..e69de29 diff --git a/src/mpiwrap/message_passing.h b/src/mpiwrap/message_passing.h new file mode 100644 index 0000000..46a7b29 --- /dev/null +++ b/src/mpiwrap/message_passing.h @@ -0,0 +1,60 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Interface to the message passing library MPI + * \par History + * JGH (02-Jan-2001): New error handling + * Performance tools + * JGH (14-Jan-2001): New routines mp_comm_compare, mp_cart_coords, + * mp_rank_compare, mp_alltoall + * JGH (06-Feb-2001): New routines mp_comm_free + * JGH (22-Mar-2001): New routines mp_comm_dup + * fawzi (04-NOV-2004): storable performance info (for f77 interface) + * Wrapper routine for mpi_gatherv added (22.12.2005,MK) + * JGH (13-Feb-2006): Flexible precision + * JGH (15-Feb-2006): single precision mp_alltoall + * \author JGH +*******************************************************************************/ + +#include + +#include "../base/base_uses.h" + +#if defined(__parallel) + USE mpi +//subroutines: unfortunately, mpi implementations do not provide interfaces for all subroutines (problems with types and ranks explosion), +// we do not quite know what is in the module, so we can not include any.... +// to nevertheless get checking for what is included, we use the mpi module without use clause, getting all there is +//USE mpi, ONLY: mpi_allgather, mpi_allgatherv, mpi_alloc_mem, mpi_allreduce, mpi_alltoall, mpi_alltoallv, mpi_bcast,& +// mpi_cart_coords, mpi_cart_create, mpi_cart_get, mpi_cart_rank, mpi_cart_sub, mpi_dims_create, mpi_file_close,& +// mpi_file_get_size, mpi_file_open, mpi_file_read_at_all, mpi_file_read_at, mpi_file_write_at_all,& +// mpi_file_write_at, mpi_free_mem, mpi_gather, mpi_gatherv, mpi_get_address, mpi_group_translate_ranks, mpi_irecv,& +// mpi_isend, mpi_recv, mpi_reduce, mpi_reduce_scatter, mpi_rget, mpi_scatter, mpi_send,& +// mpi_sendrecv, mpi_sendrecv_replace, mpi_testany, mpi_waitall, mpi_waitany, mpi_win_create +//functions +//USE mpi, ONLY: mpi_wtime +//constants +//USE mpi, ONLY: MPI_DOUBLE_PRECISION, MPI_DOUBLE_COMPLEX, MPI_REAL, MPI_COMPLEX, MPI_ANY_TAG,& +// MPI_ANY_SOURCE, MPI_COMM_NULL, MPI_REQUEST_NULL, MPI_WIN_NULL, MPI_STATUS_SIZE, MPI_STATUS_IGNORE, MPI_STATUSES_IGNORE, & +// MPI_ADDRESS_KIND, MPI_OFFSET_KIND, MPI_MODE_CREATE, MPI_MODE_RDONLY, MPI_MODE_WRONLY,& +// MPI_MODE_RDWR, MPI_MODE_EXCL, MPI_COMM_SELF, MPI_COMM_WORLD, MPI_THREAD_SERIALIZED,& +// MPI_ERRORS_RETURN, MPI_SUCCESS, MPI_MAX_PROCESSOR_NAME, MPI_MAX_ERROR_STRING, MPI_IDENT,& +// MPI_UNEQUAL, MPI_MAX, MPI_SUM, MPI_INFO_NULL, MPI_IN_PLACE, MPI_CONGRUENT, MPI_SIMILAR, MPI_MIN, MPI_SOURCE,& +// MPI_TAG, MPI_INTEGER8, MPI_INTEGER, MPI_MAXLOC, MPI_2INTEGER, MPI_MINLOC, MPI_LOGICAL, MPI_2DOUBLE_PRECISION,& +// MPI_LOR, MPI_CHARACTER, MPI_BOTTOM, MPI_MODE_NOCHECK, MPI_2REAL +#endif + + // parameters that might be needed +#if defined(__parallel) + const int MP_STD_REAL = MPI_DOUBLE_PRECISION; + const int MP_STD_COMPLEX = MPI_DOUBLE_COMPLEX; + const int MP_STD_HALF_REAL = MPI_REAL; + const int MP_STD_HALF_COMPLEX = MPI_COMPLEX; + + const bool cp2k_is_parallel = true; + diff --git a/src/qs_ot_eigensolver.c b/src/qs_ot_eigensolver.c new file mode 100644 index 0000000..fc7877e --- /dev/null +++ b/src/qs_ot_eigensolver.c @@ -0,0 +1,77 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + + + +#include + +#include "./base/base_uses.h" + +/******************************************************************************* + * \brief an eigen-space solver for the generalised symmetric eigenvalue problem + * for sparse matrices, needing only multiplications + * \author Joost VandeVondele (25.08.2002) +*******************************************************************************/ + +const char* moduleN = "qs_ot_eigensolver"; + +// on input c contains the initial guess (should not be zero !) +// on output c spans the subspace +/******************************************************************************* + * \brief ... + * \param matrix_h ... + * \param matrix_s ... + * \param matrix_orthogonal_space_fm ... + * \param matrix_c_fm ... + * \param preconditioner ... + * \param eps_gradient ... + * \param iter_max ... + * \param size_ortho_space ... + * \param silent ... + * \param ot_settings ... +*******************************************************************************/ + +void ot_eigensolver(struct dbcsr_type *matrix_h, struct dbcsr_type *matrix_s, + struct cp_fm_type matrix_orthogonal_space_fm, struct cp_fm_type matrix_c_fm, + struct preconditioner_type *preconditioner, double eps_gradient, int iter_max, + int size_ortho_space, bool silent, struct qs_ot_settings_type ot_settings) { + + const char* routineN = "ot_eigensolver"; + const int max_iter_inner_loop = 40; + double rone = 1.0; + double rzero = 0.0; + + int handle, ieigensolver, iter_total, k, n, ortho_k, ortho_space_k, output_unit; + + bool energy_only, my_silent, ortho; + double delta, energy; + + struct dbcsr_p_type *matrix_hc[]; + struct dbcsr_type *matrix_buf1_ortho, *matrix_buf2_ortho, *matrix_c, *matrix_orthogonal_space, + *matrix_os_ortho, *matrix_s_ortho; + struct qs_ot_type *qs_ot_env[]; + + timeset(routineN, handle); + + output_unit = cp_logger_get_default_io_unit(); + + if (PRESENT(silent)) + my_silent = silent; + else + my_silent = false; + + matrix_c = NULL; // fm->dbcsr + + cp_fm_get_info(matrix_c_fm, n, k); // fm->dbcsr + cp_fm_to_dbcsr_row_template(matrix_c, matrix_c_fm, matrix_h); + + iter_total = 0; + + + + +} \ No newline at end of file diff --git a/src/qs_tenosrs.c b/src/qs_tenosrs.c new file mode 100644 index 0000000..113fb4a --- /dev/null +++ b/src/qs_tenosrs.c @@ -0,0 +1,112 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief Utility methods to build 3-center integral tensors of various types. +*******************************************************************************/ + +#include + +#include "./base/base_uses.h" + +void build_2c_neighbor_lists(int** ij_list, struct gto_basis_set_p_type* basis_i, struct gto_basis_set_p_type* basis_j, + int potential_parameter, char* name, int* qs_env, bool sym_ij, bool molecular, + float* dist_2d, int pot_to_rad) { + + int ikind, nkind, pot_to_rad_prv; + bool* i_present, j_present; + double* pair_radius; + double subcells; + double* i_radius, j_radius; + int *atomic_kind_set; + int *cell; + int *local_particles; + int *dist_2d_prv; + int *atom2d; + int *molecule_set; + int *particle_set; + + if (PRESENT(pot_to_rad)) { + pot_to_rad_prv = pot_to_rad; + } else { + pot_to_rad_prv = 1; + } + + get_qs_env(qs_env, nkind, cell, particle_set, atomic_kind_set, + local_particles, dist_2d_prv, molecule_set); + + section_vals_val_get(qs_env.input, "DFT.SUBCELLS", subcells); + + i_present = calloc(sizeof(bool) * nkind); + j_present = calloc(sizeof(bool) * nkind); + i_radius = calloc(sizeof(double) * nkind); + j_radius = calloc(sizeof(double) * nkind); + + if (PRESENT(dist_2d)) + dist_2d_prv = &dist2d; + + // Set up the radii, depending on the operator type + if (potential_parameter.potential_type == do_potential_id) { + + //overlap => use the kind radius for both i and j + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind)); + } + if (ASSOCIATED(basis_j(ikind)%gto_basis_set)) { + j_present(ikind) = true; + get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind)); + } + } + } else if (potential_parameter.potential_type == do_potential_coulomb) { + + //Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + if (pot_to_rad_prv == 1) + i_radius(ikind) = 1000000.0; + } + if (ASSOCIATED(basis_j(ikind).gto_basis_set)) { + j_present(ikind) = true; + if (pot_to_rad_prv == 2) + j_radius(ikind) = 1000000.0; + } + } //ikind + + } else if (potential_parameter.potential_type == do_potential_truncated || + potential_parameter.potential_type == do_potential_short) { + + //Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii + for (ikind = 0; ikind < nkind; ikind++) { + if (ASSOCIATED(basis_i(ikind).gto_basis_set)) { + i_present(ikind) = true; + get_gto_basis_set(basis_i(ikind).gto_basis_set, kind_radius=i_radius(ikind)); + if (pot_to_rad_prv == 1) + i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter.cutoff_radius; + } + if (ASSOCIATED(basis_j(ikind).gto_basis_set)) { + j_present(ikind) = true; + get_gto_basis_set(basis_j(ikind).gto_basis_set, kind_radius=j_radius(ikind)); + if (pot_to_rad_prv == 2) + j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter.cutoff_radius; + } + } + + } else { + CPABORT("Operator not implemented."); + } + + pair_radius = calloc(sizeof(double)*nkind*nkind); + + + + + +} + diff --git a/src/sockets.c b/src/sockets.c new file mode 100644 index 0000000..f9ef530 --- /dev/null +++ b/src/sockets.c @@ -0,0 +1,177 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2023 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +/*----------------------------------------------------------------------------*/ +/* Copyright (C) 2013, Joshua More and Michele Ceriotti */ +/* */ +/* Permission is hereby granted, free of charge, to any person obtaining */ +/* a copy of this software and associated documentation files (the */ +/* "Software"), to deal in the Software without restriction, including */ +/* without limitation the rights to use, copy, modify, merge, publish, */ +/* distribute, sublicense, and/or sell copies of the Software, and to */ +/* permit persons to whom the Software is furnished to do so, subject to */ +/* the following conditions: */ +/* */ +/* The above copyright notice and this permission notice shall be included */ +/* in all copies or substantial portions of the Software. */ +/* */ +/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, */ +/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF */ +/* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. */ +/* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY */ +/* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, */ +/* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE */ +/* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ +/*----------------------------------------------------------------------------*/ + +/******************************************************************************* + * \brief A minimal wrapper for socket communication. + * Contains both the functions that transmit data to the socket and read + * the data back out again once finished, and the function which opens + * the socket initially. Can be linked to a FORTRAN code that does not + * support sockets natively. + * \author Joshua More and Michele Ceriotti + ******************************************************************************/ +#ifndef __NO_IPI_DRIVER + +#define _POSIX_C_SOURCE 200809L + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/******************************************************************************* + * \brief Opens a socket. + * \param psockfd The id of the socket that will be created. + * \param inet An integer that determines whether the socket will be an inet + * or unix domain socket. Gives unix if 0, inet otherwise. + * \param port The port number for the socket to be created. Low numbers are + * often reserved for important channels, so use of numbers of 4 + * or more digits is recommended. + * \param host The name of the host server. + * \note Fortran passes an extra argument for the string length, but this is + * ignored here for C compatibility. + ******************************************************************************/ +void open_socket(int *psockfd, int *inet, int *port, char *host) { + int sockfd, ai_err; + + if (*inet > 0) { // creates an internet socket + + // fetches information on the host + struct addrinfo hints, *res; + char service[256]; + + memset(&hints, 0, sizeof(hints)); + hints.ai_socktype = SOCK_STREAM; + hints.ai_family = AF_INET; + hints.ai_flags = AI_PASSIVE; + + sprintf(service, "%d", *port); // convert the port number to a string + ai_err = getaddrinfo(host, service, &hints, &res); + if (ai_err != 0) { + perror("Error fetching host data. Wrong host name?"); + exit(-1); + } + + // creates socket + sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol); + if (sockfd < 0) { + perror("Error opening socket"); + exit(-1); + } + + // makes connection + if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0) { + perror("Error opening INET socket: wrong port or server unreachable"); + exit(-1); + } + freeaddrinfo(res); + } else { // creates a unix socket + struct sockaddr_un serv_addr; + + // fills up details of the socket address + memset(&serv_addr, 0, sizeof(serv_addr)); + serv_addr.sun_family = AF_UNIX; + strcpy(serv_addr.sun_path, "/tmp/ipi_"); + strcpy(serv_addr.sun_path + 9, host); + + // creates the socket + sockfd = socket(AF_UNIX, SOCK_STREAM, 0); + + // connects + if (connect(sockfd, (struct sockaddr *)&serv_addr, sizeof(serv_addr)) < 0) { + perror( + "Error opening UNIX socket: path unavailable, or already existing"); + exit(-1); + } + } + + *psockfd = sockfd; +} + +/******************************************************************************* + * \brief Writes to a socket. + * \param psockfd The id of the socket that will be written to. + * \param data The data to be written to the socket. + * \param plen The length of the data in bytes. + ******************************************************************************/ +void writebuffer(int *psockfd, char *data, int *plen) { + int n; + int sockfd = *psockfd; + int len = *plen; + + n = write(sockfd, data, len); + if (n < 0) { + perror("Error writing to socket: server has quit or connection broke"); + exit(-1); + } +} + +/******************************************************************************* + * \brief Reads from a socket. + * \param psockfd The id of the socket that will be read from. + * \param data The storage array for data read from the socket. + * \param plen The length of the data in bytes. + ******************************************************************************/ +void readbuffer(int *psockfd, char *data, int *plen) { + int n, nr; + int sockfd = *psockfd; + int len = *plen; + + n = nr = read(sockfd, data, len); + + while (nr > 0 && n < len) { + nr = read(sockfd, &data[n], len - n); + n += nr; + } + + if (n == 0) { + perror("Error reading from socket: server has quit or connection broke"); + exit(-1); + } + } + +/******************************************************************************* + * \brief Mini-wrapper to nanosleep + * \param dsec number of seconds to wait (float values accepted) + ******************************************************************************/ +void uwait(double *dsec) { + struct timespec wt, rem; + wt.tv_sec = floor(*dsec); + wt.tv_nsec = (*dsec - wt.tv_sec) * 1000000000; + nanosleep(&wt, &rem); +} + +#endif \ No newline at end of file diff --git a/src/start/cp2k.c b/src/start/cp2k.c new file mode 100644 index 0000000..a50a19e --- /dev/null +++ b/src/start/cp2k.c @@ -0,0 +1,103 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Main program of CP2K +//> \par Copyright +//> CP2K: A general program to perform molecular dynamics simulations +//> Copyright (C) 2000, 2001, 2002, 2003 CP2K developers group +//> Copyright (C) 2004, 2005, 2006, 2007 CP2K developers group +//> Copyright (C) 2008, 2009, 2010, 2011 CP2K developers group +//> Copyright (C) 2012, 2013, 2014, 2015 CP2K developers group +//> Copyright (C) 2016 CP2K developers group +//> \par +//> This program is free software; you can redistribute it and/or modify +//> it under the terms of the GNU General Public License as published by +//> the Free Software Foundation; either version 2 of the License, or +//> (at your option) any later version. +//> \par +//> This program is distributed in the hope that it will be useful, +//> but WITHOUT ANY WARRANTY; without even the implied warranty of +//> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//> GNU General Public License for more details. +//> \par +//> You should have received a copy of the GNU General Public License +//> along with this program; if not, write to the Free Software +//> Foundation, Inc., 51 Franklin Street, Fifth Floor, +//> Boston, MA 02110-1301, USA. +//> \par +//> See also https://www.fsf.org/licensing/licenses/gpl.html +//> \par +//> CP2K, including its sources and pointers to the authors +//> can be found at https://www.cp2k.org/ +//> \note +//> should be kept as lean as possible. +//> see cp2k_run for more comments +//> \author Joost VandeVondele +// ************************************************************************************************** +#include + +#include "../base/kinds.h" +#include "../base/base_uses.h" + +int main(int argc, char** argv) { + + char input_file_name[default_path_length]; + + int output_unit, l, i, var_set_sep, inp_var_idx, ierr, i_arg; + + bool check, usage, echo_input, command_line_error,run_it; + bool force_run, has_input, xml, print_version, print_license, shell_mode; + int* input_declaration; + + // output goes to the screen by default + output_unit = default_output_unit; + + // set default behaviour for the command line switches + check = false; + usage = false; + echo_input = false; + has_input = false; + run_it = true; + shell_mode = false; + force_run = false; + print_version = false; + print_license = false; + + command_line_error = false; + xml = false; + input_file_name = "Missing input file name" // no default + output_file_name = "__STD_OUT__" // by default we go to std_out + ALLOCATE (initial_variables(2, 1:0)) + + // Get command and strip path + GET_COMMAND_ARGUMENT(NUMBER=0, VALUE=command, STATUS=ierr) + CPASSERT(ierr == 0) + l = LEN_TRIM(command) + DO i = l, 1, -1 + IF (command(i:i) == "/" || command(i:i) == "\\") EXIT + END DO + command = command(i + 1:l) + + // Consider output redirection + i_arg = 0; + + + + if (i == 0) { + else { + section_release(input_declaration); + } + else { + printf("initial setup (MPI ?) error"); + } + + // and the final cleanup + finalize_cp2k(finalize_mpi=true, ierr=ierr); + delete initial_variables; + CPASSERT(ierr == 0); +} diff --git a/src/start/cp2k.cpp b/src/start/cp2k.cpp new file mode 100644 index 0000000..f970da7 --- /dev/null +++ b/src/start/cp2k.cpp @@ -0,0 +1,337 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief Main program of CP2K +//> \par Copyright +//> CP2K: A general program to perform molecular dynamics simulations +//> Copyright (C) 2000, 2001, 2002, 2003 CP2K developers group +//> Copyright (C) 2004, 2005, 2006, 2007 CP2K developers group +//> Copyright (C) 2008, 2009, 2010, 2011 CP2K developers group +//> Copyright (C) 2012, 2013, 2014, 2015 CP2K developers group +//> Copyright (C) 2016 CP2K developers group +//> \par +//> This program is free software; you can redistribute it and/or modify +//> it under the terms of the GNU General Public License as published by +//> the Free Software Foundation; either version 2 of the License, or +//> (at your option) any later version. +//> \par +//> This program is distributed in the hope that it will be useful, +//> but WITHOUT ANY WARRANTY; without even the implied warranty of +//> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//> GNU General Public License for more details. +//> \par +//> You should have received a copy of the GNU General Public License +//> along with this program; if not, write to the Free Software +//> Foundation, Inc., 51 Franklin Street, Fifth Floor, +//> Boston, MA 02110-1301, USA. +//> \par +//> See also https://www.fsf.org/licensing/licenses/gpl.html +//> \par +//> CP2K, including its sources and pointers to the authors +//> can be found at https://www.cp2k.org/ +//> \note +//> should be kept as lean as possible. +//> see cp2k_run for more comments +//> \author Joost VandeVondele +// ************************************************************************************************** + +#include +#include +#include +#include +#include +#include + +#include "base_uses.h" + +int main(int argc, char** argv) { + + std::string input_file_name, output_file_name, arg_att, command; + std::vector initial_variables, initial_variables_tmp; + std::string compiler_options_string; + + int output_unit, l, i, var_set_sep, inp_var_idx, ierr, i_arg; + + bool check, usage, echo_input, command_line_error, run_it; + bool force_run, has_input, xml, print_version, print_license, shell_mode; + int *input_declaration; + + // output goes to the screen by default + output_unit = default_output_unit; + + // set default behaviour for the command line switches + check = false; + usage = false; + echo_input = false; + has_input = false; + run_it = true; + shell_mode = false; + force_run = false; + print_version = false; + print_license = false; + command_line_error = false; + xml = false; + input_file_name = "Missing input file name"; // no default + output_file_name = "__STD_OUT__"; // by default we go to std_out + + // Get command and strip path + get_command_argument(number=0, value=command, status=ierr); + assert(ierr == 0); + l = len_trim(command); + for (i = l-1; i >= 0; i--) { + if (command[i] == '/' or command[i] == '\\') exit(1); + } + command = command.substr(i); + + // Consider output redirection + i_arg = 0; + while (i_arg < command_argument_count()) { + i_arg++; + get_command_argument(number=&i_arg, value=&arg_att, status=&ierr); + assert(ierr == 0); + switch (arg_att) { + case "-o": + if (output_file_name == "__STD_OUT__") { + // Consider only the first -o flag + i_arg++; + get_command_argument(number=&i_arg, value=&arg_att, status=&ierr); + assert(ierr == 0); + if (arg_att[0] == '-') { + std::cout << "ERROR: The output file name " << arg_att << " starts with -" << std::endl; + command_line_error = true; + } else { + output_file_name = arg_att; + open_file(file_name=output_file_name, + file_status="UNKNOWN", + file_action="WRITE", + file_position="APPEND", + skip_get_unit_number=true, + unit_number=output_unit); + } + } else { + i_arg++; + std::cout << "ERROR: The command line flag -o has been specified multiple times" << std::endl; + command_line_error = true; + } + break; + } + } + + // Check if binary was invoked as cp2k_shell + if (command.compare("cp2k_shell")) { + shell_mode = true; + run_it = false; + } else if (command_argument_count() < 1) { + std::cout << "ERROR: At least one command line argument must be specified" << std::endl; + command_line_error = true; + } + + // Check if binary was invoked as sopt or popt alias + l = len_trim(command); + if (command[l-4:l] == ".sopt" or command[l-4:l] == ".popt") { + omp_set_num_threads(1); + } + +#ifdef __ACCELERATE + if (omp_get_max_threads() > 1) { + std::string _block_env_var; + int _block_veclib_max_threads, _block_ierr; + get_environment_variable("VECLIB_MAXIMUM_THREADS", _block_env_var, status=_block_ierr); + _block_veclib_max_threads = 0; + if (_block_ierr == 0) { + _block_env_var = _block_veclib_max_threads; + } + if (_block_ierr == 1 or (_block_ierr == 0 and _block_veclib_max_threads > 1)) { + cp_warn(__LOCATON__, "macOS' Accelerate framework has its own threading enabled which may interfere" + " with the OpenMP threading. You can disable the Accelerate threading by setting" + " the environment variable VECLIB_MAXIMUM_THREADS=1") + } + } +#endif + + i_arg = 0; + while (i_arg < command_argument_count()) { + i_arg++; + get_command_argument(&i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + switch (arg_att) { + case "--check": case "-c": + check = true; + run_it = false; + echo_input = false; + break; + case "--echo": case "-e": + check = true; + run_it = false; + echo_input = true; + break; + case "-v": case "--version": + print_version = true; + run_it = false; + break; + case "--license": + print_license = true; + run_it = false; + break; + case "--run": case "-r": + force_run = true; + break; + case "--shell": case "-s": + shell_mode = true; + run_it = false; + break; + case "-help": case "--help": case "-h": + usage = true; + run_it = false; + break; + case "-i": + i_arg++; + get_command_argument(i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + // argument does not start with a - it is a filename + if (!(arg_att[0] == '-')) { + input_file_name = arg_att; + has_input = true; + } else { + std::cout << "ERROR: The input file name " << arg_att << " starts with -" << std::endl; + command_line_error = true; + exit(1); + } + break; + case "-E": case "--set": + i_arg++; + get_command_argument(i_arg, &arg_att, status=&ierr); + assert(ierr == 0); + + var_set_sep = arg_att.find("="); + + if (var_set_sep < 1) { + std::cout << "ERROR: Invalid initializer for preprocessor variable: " << arg_att << std::endl; + command_line_error = false; + exit(1); + } + + for (inp_var_idx = 0; inp_var_idx < sizeof(initial_variables, std::string); inp_var_idx++) { + // check whether the variable was already set, in this case, overwrite + if (initial_variables[0:inp_var_idx] == arg_att[0:var_set_sep - 1]) exit(1); + } + + break; + } + } + + if ((run_it || force_run || check || echo_input) && !has_input && !command_line_error) { + std::cout << "\n ERROR: An input file name is required" << std::endl; + command_line_error = true; + } + + init_cp2k(init_mpi=true, ierr=&ierr); + + if (ierr == 0){ + // some first info concerning how to run CP2K + if (usage || command_line_error) { + if (default_para_env.is_source()) { + l = len_trim(command); + std::cout << "\n " << command << " [-c|--check] [-e|--echo] [-h|--help]\n"; + std::cout << std::string(l, ' ') << " [-i] \n"; + std::cout << std::string(l, ' ') << " [-mpi-mapping|--mpi-mapping] \n"; + std::cout << std::string(l, ' ') << " [-o] \n"; + std::cout << std::string(l, ' ') << " [-r|-run] [-s|--shell] [--xml]\n"; + + std::cout << "\n starts the CP2K program, see \n"; + std::cout << "\n The easiest way is " << command << " \n"; + std::cout << "\n The following options can be used:\n"; + std::cout << "\n -i : provides an input file name, if it is the last\n"; + std::cout << std::string(24, ' ') << " argument, the -i flag is not needed\n"; + std::cout << " -o : provides an output file name [default: screen]\n"; + + std::cout << "\n These switches skip the simulation, unless [-r|-run] is specified:\n"; + std::cout << "\n --check, -c : performs a syntax check of the \n"; + std::cout << " --echo, -e : echoes the , and make all defaults explicit\n"; + std::cout << " The input is also checked, but only a failure is reported\n"; + std::cout << " --help, -h : writes this message\n"; + std::cout << " --license : prints the CP2K license\n"; + std::cout << " --mpi-mapping : applies a given MPI reordering to CP2K\n"; + std::cout << " --run, -r : forces a CP2K run regardless of other specified flags\n"; + std::cout << " --shell, -s : start interactive shell mode\n"; + std::cout << " --version, -v : prints the CP2K version and the revision number\n"; + std::cout << " --xml : dumps the whole CP2K input structure as a XML file\n"; + std::cout << " xml2htm generates a HTML manual from this XML file\n"; + std::cout << " --set, -E name=value : set the initial value of a preprocessor value\n " << std::endl; + } + } + + if (!command_line_error) { + // write the version string + if (print_version) { + if (default_para_env.is_source()) { + std::cout << " " << cp2k_version << "\n Source code revision " << compile_revision << "\n " << cp2k_flags() << "\n"; + compiler_options_string = compiler_options(); + std::cout << " compiler: " << compiler_version() << "\n"; + std::cout << " compiler options:"; + for (i = 0; i < (len(compiler_options_string) - 1)/68; i++) { + std::cout << compiler_options_string[i*68:fmin(len(compiler_options_string), (i+1)*68)] << "\n"; + } + std::cout << std::endl; + // delete compiler_options_string; + } + } + + // write the license + if (print_license) { + if (default_para_env.is_source()) { + print_cp2k_license(output_unit); + } + } + + if (xml) { + if (default_para_env.is_source()) { + write_xml_file(); + } + } + + create_cp2k_root_section(input_declaration); + + if (check) { + check_input(input_declaration, input_file_name, output_file_name, echo_input=echo_input, + ierr=&ierr, initial_variables=initial_variables); + + if (default_para_env.is_source()) { + if (ierr == 0) { + if (!echo_input) { + std::cout << "SUCCESS, the input could be parsed correctly.\n"; + std::cout << " This does not guarantee that this input is meaningful\n"; + std::cout << " or will run successfully" << std::endl; + } + } else { + std::cout << "ERROR, the input could *NOT* be parsed correctly."; + std::cout << " Please, check and correct it" << std::endl; + } + } + } + + if (shell_mode) { + launch_cp2k_shell(input_declaration) + } + + if (run_it || force_run) { + run_input(input_declaration, input_file_name, output_file_name, initial_variables) + } + + section_release(input_declaration); + } + } else { + std::cout << "initial setup (MPI ?) error" << std::endl; + } + + // and final cleanup + finalize_cp2k(finalize_mpi=true, ierr=&ierr); + // delete initial_variables; + assert(ierr == 0); + return 0; +} \ No newline at end of file diff --git a/src/start/python/README.md b/src/start/python/README.md new file mode 100644 index 0000000..da12fe7 --- /dev/null +++ b/src/start/python/README.md @@ -0,0 +1,30 @@ +# CP2K Python Bindings + +## Installation + +There is a target `py-cython-bindings` in the global `Makefile` to build the +Python bindings. The shared object can be found in: +`/lib///python` + +Only the Python headers and a NumPy installation are required. + +## Development + +To regenerate the C file from the `cp2k.pyx`, `Cython` is required and should +be called as follows: + +```sh +cd /src/start/python +cython cp2k.pyx +``` + +Unittests can be found in the `test/` directory. They must be run in separate +Python interpreter instances due to side effects in the library. + +## Known Issues + +* If libcp2k is built with MPI support, you may get an MPI initialization error + depending on your MPI implementation/configuration. In that case MPI must be + initialized first by using Mpi4py and the Fortran MPI communicator handler + must be passed down the CP2K via the respective `...comm` functions. + The reason for this is documented here: diff --git a/src/start/python/cp2k.pyx b/src/start/python/cp2k.pyx new file mode 100644 index 0000000..e06fe5c --- /dev/null +++ b/src/start/python/cp2k.pyx @@ -0,0 +1,161 @@ +# cython: language_level=2 +# vim: set ts=4 sw=4 tw=0 : + +import numpy as np + +from cpython.mem cimport PyMem_Malloc, PyMem_Free + +cdef extern from "../libcp2k.h": + ctypedef int force_env_t + + void cp2k_get_version(char* version_str, int str_length) + void cp2k_init() + void cp2k_init_without_mpi() + void cp2k_finalize() + void cp2k_finalize_without_mpi() + void cp2k_create_force_env(force_env_t* new_force_env, const char* input_file_path, const char* output_file_path) + void cp2k_create_force_env_comm(force_env_t* new_force_env, const char* input_file_path, const char* output_file_path, int mpi_comm) + void cp2k_destroy_force_env(force_env_t force_env) + void cp2k_set_positions(force_env_t force_env, const double* new_pos, int n_el) + void cp2k_set_velocities(force_env_t force_env, const double* new_vel, int n_el) + void cp2k_get_result(force_env_t force_env, const char* description, double* result, int n_el) + void cp2k_get_natom(force_env_t force_env, int* natom) + void cp2k_get_nparticle(force_env_t force_env, int* nparticle) + void cp2k_get_positions(force_env_t force_env, double* pos, int n_el) + void cp2k_get_forces(force_env_t force_env, double* force, int n_el) + void cp2k_get_potential_energy(force_env_t force_env, double* e_pot) + void cp2k_calc_energy_force(force_env_t force_env) + void cp2k_calc_energy(force_env_t force_env) + void cp2k_run_input(const char* input_file_path, const char* output_file_path) + void cp2k_run_input_comm(const char* input_file_path, const char* output_file_path, int mpi_comm) + +def get_version_string(): + n = 255 * sizeof(char) + + data = PyMem_Malloc(n) + if not data: + raise MemoryError() + + versionstr = '' + try: + cp2k_get_version(data, n) + versionstr = data.decode('UTF-8') + finally: + PyMem_Free(data) + + return versionstr + +def init(manage_mpi = True): + if manage_mpi: + cp2k_init() + else: + cp2k_init_without_mpi() + +def finalize(manage_mpi = True): + if manage_mpi: + cp2k_finalize() + else: + cp2k_finalize_without_mpi() + +def run_input(input_file_path, output_file_path = None, mpi_comm = None): + input_file_path = input_file_path.encode('UTF-8') + + if output_file_path is None: + output_file_path = u'__STD_OUT__'.encode('UTF-8') + else: + output_file_path = output_file_path.encode('UTF-8') + + if mpi_comm: + cp2k_run_input_comm(input_file_path, output_file_path, mpi_comm) + else: + cp2k_run_input(input_file_path, output_file_path) + +def create_force_env(input_file_path, output_file_path, mpi_comm = None): + cdef force_env_t fenv + if mpi_comm: + cp2k_create_force_env_comm(&fenv, input_file_path, output_file_path, mpi_comm) + else: + cp2k_create_force_env(&fenv, input_file_path, output_file_path) + +cdef class ForceEnvironment(object): + cdef force_env_t _force_env + + def __init__(self, input_file_path not None, output_file_path = None, mpi_comm = None): + input_file_path = input_file_path.encode('UTF-8') + + if output_file_path is None: + output_file_path = u'__STD_OUT__'.encode('UTF-8') + else: + output_file_path = output_file_path.encode('UTF-8') + + if mpi_comm: + cp2k_create_force_env_comm(&self._force_env, input_file_path, output_file_path, mpi_comm) + else: + cp2k_create_force_env(&self._force_env, input_file_path, output_file_path) + + def destroy(self): + cp2k_destroy_force_env(self._force_env) + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_value, traceback): + self.destroy() + + + property positions: + + def __get__(self): + positions = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] positions_view = positions + cp2k_get_positions(self._force_env, &positions_view[0], positions_view.shape[0]) + return positions + + def __set__(self, double[::1] positions not None): + if positions.shape[0] != 3*self.nparticle: + raise ValueError('the positions array must have exactly {} (3*nparticle) elements'.format(3*self.nparticle)) + + cp2k_set_positions(self._force_env, &positions[0], positions.shape[0]) + + def set_velocities(self, double[::1] velocities not None): + if velocities.shape[0] != 3*self.nparticle: + raise ValueError('the velocities array must have exactly {} (3*nparticle) elements'.format(3*self.nparticle)) + + cp2k_set_velocities(self._force_env, &velocities[0], velocities.shape[0]) + + property natom: + def __get__(self): + cdef int natom + cp2k_get_natom(self._force_env, &natom) + return natom + + property nparticle: + def __get__(self): + cdef int nparticle + cp2k_get_nparticle(self._force_env, &nparticle) + return nparticle + + def get_result(self, description): + results = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] results_view = results + cp2k_get_result(self._force_env, description.encode('UTF-8'), &results_view[0], results_view.shape[0]) + return results + + property forces: + def __get__(self): + forces = np.zeros([3*self.nparticle], dtype=np.double) + cdef double [::1] forces_view = forces + cp2k_get_forces(self._force_env, &forces_view[0], forces_view.shape[0]) + return forces + + def calc_energy_force(self): + cp2k_calc_energy_force(self._force_env) + + def calc_energy(self): + cp2k_calc_energy(self._force_env) + + property potential_energy: + def __get__(self): + cdef double e_pot + cp2k_get_potential_energy(self._force_env, &e_pot) + return e_pot diff --git a/src/subcell_types.c b/src/subcell_types.c new file mode 100644 index 0000000..120e046 --- /dev/null +++ b/src/subcell_types.c @@ -0,0 +1,211 @@ +//--------------------------------------------------------------------------------------------------// +// CP2K: A general program to perform molecular dynamics simulations // +// Copyright 2000-2021 CP2K developers group // +// // +// SPDX-License-Identifier: GPL-2.0-or-later // +//--------------------------------------------------------------------------------------------------// + +// ************************************************************************************************** +//> \brief subcell types and allocation routines +//> \par History +//> - Separated from qs_neighbor_lists (25.07.2010,jhu) +//> \author Matthias Krack +// ************************************************************************************************** +#include "./base/base_uses.h" + +void subcell_types() { + + USE cell_types, ONLY: cell_type,& + real_to_scaled,& + scaled_to_real + USE kinds, ONLY: dp + USE util, ONLY: sort + + IMPLICIT NONE + + PRIVATE + +// ************************************************************************************************** + TYPE subcell_type + INTEGER :: natom + REAL(KIND=dp), DIMENSION(3) :: s_max, s_min + INTEGER, DIMENSION(:), POINTER :: atom_list + REAL(KIND=dp), DIMENSION(3, 8) :: corners + END TYPE subcell_type + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'subcell_types' + + PUBLIC :: subcell_type, allocate_subcell, deallocate_subcell + PUBLIC :: reorder_atoms_subcell, give_ijk_subcell + +// ************************************************************************************************** + +CONTAINS + +// ************************************************************************************************** +//> \brief Allocate and initialize a subcell grid structure for the atomic neighbor search. +//> \param subcell ... +//> \param nsubcell ... +//> \param maxatom ... +//> \param cell ... +//> \date 12.06.2003 +//> \author MK +//> \version 1.0 +// ************************************************************************************************** + SUBROUTINE allocate_subcell(subcell, nsubcell, maxatom, cell) + + TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell + INTEGER, DIMENSION(3), INTENT(IN) :: nsubcell + INTEGER, INTENT(IN), OPTIONAL :: maxatom + TYPE(cell_type), OPTIONAL, POINTER :: cell + + INTEGER :: i, j, k, na, nb, nc + REAL(dp) :: a_max, a_min, b_max, b_min, c_max, & + c_min, delta_a, delta_b, delta_c + + na = nsubcell(1) + nb = nsubcell(2) + nc = nsubcell(3) + + ALLOCATE (subcell(na, nb, nc)) + + delta_a = 1.0_dp/REAL(na, dp) + delta_b = 1.0_dp/REAL(nb, dp) + delta_c = 1.0_dp/REAL(nc, dp) + + c_min = -0.5_dp + + DO k = 1, nc + c_max = c_min + delta_c + b_min = -0.5_dp + DO j = 1, nb + b_max = b_min + delta_b + a_min = -0.5_dp + DO i = 1, na + a_max = a_min + delta_a + subcell(i, j, k)%s_min(1) = a_min + subcell(i, j, k)%s_min(2) = b_min + subcell(i, j, k)%s_min(3) = c_min + subcell(i, j, k)%s_max(1) = a_max + subcell(i, j, k)%s_max(2) = b_max + subcell(i, j, k)%s_max(3) = c_max + subcell(i, j, k)%natom = 0 + IF (PRESENT(cell)) THEN + CALL scaled_to_real(subcell(i, j, k)%corners(:, 1), (/a_min, b_min, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 2), (/a_max, b_min, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 3), (/a_min, b_max, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 4), (/a_max, b_max, c_min/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 5), (/a_min, b_min, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 6), (/a_max, b_min, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 7), (/a_min, b_max, c_max/), cell) + CALL scaled_to_real(subcell(i, j, k)%corners(:, 8), (/a_max, b_max, c_max/), cell) + END IF + IF (PRESENT(maxatom)) THEN + ALLOCATE (subcell(i, j, k)%atom_list(maxatom)) + END IF + a_min = a_max + END DO + b_min = b_max + END DO + c_min = c_max + END DO + + END SUBROUTINE allocate_subcell + +// ************************************************************************************************** +//> \brief Deallocate a subcell grid structure. +//> \param subcell ... +//> \date 16.06.2003 +//> \author MK +//> \version 1.0 +// ************************************************************************************************** + SUBROUTINE deallocate_subcell(subcell) + + TYPE(subcell_type), DIMENSION(:, :, :), POINTER :: subcell + + INTEGER :: i, j, k + + IF (ASSOCIATED(subcell)) THEN + + DO k = 1, SIZE(subcell, 3) + DO j = 1, SIZE(subcell, 2) + DO i = 1, SIZE(subcell, 1) + DEALLOCATE (subcell(i, j, k)%atom_list) + END DO + END DO + END DO + + DEALLOCATE (subcell) + ELSE + CPABORT("") + END IF + + END SUBROUTINE deallocate_subcell + +// ************************************************************************************************** +//> \brief ... +//> \param atom_list ... +//> \param kind_of ... +//> \param work ... +//> \par History +//> 08.2006 created [tlaino] +//> \author Teodoro Laino +// ************************************************************************************************** + SUBROUTINE reorder_atoms_subcell(atom_list, kind_of, work) + // work needs to be dimensioned 3xSIZE(atom_list) + INTEGER, DIMENSION(:), POINTER :: atom_list + INTEGER, DIMENSION(:), INTENT(IN) :: kind_of + INTEGER, DIMENSION(:) :: work + + INTEGER :: i, i0, i1, i2, j0, j1, j2 + + i0 = 1 + j0 = SIZE(atom_list) + i1 = j0 + 1 + j1 = 2*j0 + i2 = j1 + 1 + j2 = 3*j0 + // Sort kind + DO i = 1, SIZE(atom_list) + work(i0 + i - 1) = kind_of(atom_list(i)) + END DO + CALL sort(work(i0:j0), SIZE(atom_list), work(i1:j1)) + work(i2:j2) = atom_list + DO i = 1, SIZE(atom_list) + atom_list(i) = work(i2 + work(i1 + i - 1) - 1) + END DO + END SUBROUTINE reorder_atoms_subcell + +// ************************************************************************************************** +//> \brief ... +//> \param r ... +//> \param i ... +//> \param j ... +//> \param k ... +//> \param cell ... +//> \param nsubcell ... +//> \par History +//> 08.2006 created [tlaino] +//> \author Teodoro Laino +// ************************************************************************************************** + SUBROUTINE give_ijk_subcell(r, i, j, k, cell, nsubcell) + REAL(KIND=dp) :: r(3) + INTEGER, INTENT(OUT) :: i, j, k + TYPE(cell_type), POINTER :: cell + INTEGER, DIMENSION(3), INTENT(IN) :: nsubcell + + REAL(KIND=dp) :: r_pbc(3), s(3), s_pbc(3) + + r_pbc = r + CALL real_to_scaled(s_pbc, r_pbc, cell) + s(:) = s_pbc + 0.5_dp + i = INT(s(1)*REAL(nsubcell(1), KIND=dp)) + 1 + j = INT(s(2)*REAL(nsubcell(2), KIND=dp)) + 1 + k = INT(s(3)*REAL(nsubcell(3), KIND=dp)) + 1 + i = MIN(MAX(i, 1), nsubcell(1)) + j = MIN(MAX(j, 1), nsubcell(2)) + k = MIN(MAX(k, 1), nsubcell(3)) + + END SUBROUTINE give_ijk_subcell + +}