diff --git a/INSTALL.md b/INSTALL.md index a291b16afa..9ba785ddf2 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -116,7 +116,7 @@ If using this framework, `-D__ACCELERATE` must be defined to account for some interface incompatibilities between Accelerate and reference BLAS/LAPACK. When building on/for Windows using the Minimalist GNU for Windows (MinGW) environment, -you must set `-D__MINGW`, `-D__NO_STATM_ACCESS` and `-D__NO_IPI_DRIVER` to avoid +you must set `-D__MINGW`, `-D__NO_STATM_ACCESS` and `-D__NO_SOCKETS` to avoid undefined references during linking, respectively errors while printing the statistics. ### 2e. MPI and SCALAPACK (optional, required for MPI parallel builds) @@ -547,12 +547,15 @@ partially depending on installed libraries (see 2.) switch on energy profiling on Cray systems - `-D__NO_ABORT` to avoid calling abort, but STOP instead (useful for coverage testing, and to avoid core dumps on some systems) +- `-D__HDF5` enables hdf5 support. This is a hard dependency for SIRIUS, but can also + be used by itself to allow read/write functionalities of QCSchema files in the + active space module. Features useful to deal with legacy systems - `-D__NO_MPI_THREAD_SUPPORT_CHECK` - Workaround for MPI libraries that do not declare they are thread safe (serialized). -- `-D__NO_IPI_DRIVER` disables the socket interface in case of troubles compiling +- `-D__NO_SOCKETS` disables the socket interface in case of troubles compiling on systems that do not support POSIX sockets. - `-D__HAS_IEEE_EXCEPTIONS` disables trapping temporarily for libraries like scalapack. - The Makefile automatically compiles in the path to the data directory via the diff --git a/arch/Linux-x86-64-mingw64-minimal.sopt b/arch/Linux-x86-64-mingw64-minimal.sopt index 915f005375..16b150b801 100644 --- a/arch/Linux-x86-64-mingw64-minimal.sopt +++ b/arch/Linux-x86-64-mingw64-minimal.sopt @@ -14,7 +14,7 @@ CXX = x86_64-w64-mingw32-g++ FC = x86_64-w64-mingw32-gfortran LD = x86_64-w64-mingw32-gfortran AR = x86_64-w64-mingw32-ar -r -DFLAGS = -D__NO_STATM_ACCESS -D__NO_IPI_DRIVER -D__MINGW +DFLAGS = -D__NO_STATM_ACCESS -D__NO_SOCKETS -D__MINGW CFLAGS = $(DFLAGS) -O2 FCFLAGS = $(DFLAGS) -O2 -ffree-form -ffree-line-length-none \ -ftree-vectorize -funroll-loops -std=f2008 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0c81c36794..6b6135ff27 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -467,6 +467,10 @@ list( qs_2nd_kernel_ao.F qs_active_space_methods.F qs_active_space_types.F + qs_active_space_utils.F + qcschema.F + hdf5_wrapper.F + sockets_interface.F qs_atomic_block.F qs_band_structure.F qs_basis_gradient.F diff --git a/src/cp2k_info.F b/src/cp2k_info.F index ff0424b8c6..2b6542303b 100644 --- a/src/cp2k_info.F +++ b/src/cp2k_info.F @@ -161,8 +161,8 @@ CONTAINS CALL integer_to_string(__MAX_CONTR, tmp_str) flags = TRIM(flags)//" max_contr="//TRIM(tmp_str) #endif -#if defined(__NO_IPI_DRIVER) - flags = TRIM(flags)//" no_ipi_driver" +#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" @@ -247,6 +247,10 @@ CONTAINS flags = TRIM(flags)//" libvdwxc" #endif +#if defined(__HDF5) + flags = TRIM(flags)//" hdf5" +#endif + END FUNCTION cp2k_flags ! ************************************************************************************************** diff --git a/src/hdf5_wrapper.F b/src/hdf5_wrapper.F new file mode 100644 index 0000000000..13fbd15c92 --- /dev/null +++ b/src/hdf5_wrapper.F @@ -0,0 +1,762 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2023 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief A wrapper around the HDF5 Fortran API +!> \par History +!> 04.2023 created [SB] +!> \author Stefano Battaglia +! ************************************************************************************************** +MODULE hdf5_wrapper + + USE hdf5, ONLY: & + h5aclose_f, h5acreate_f, h5aopen_f, h5aread_f, h5awrite_f, h5close_f, h5dclose_f, & + h5dcreate_f, h5dget_space_f, h5dopen_f, h5dread_f, h5dwrite_f, h5f_acc_rdonly_f, & + h5f_acc_trunc_f, h5fclose_f, h5fcreate_f, h5fopen_f, h5gclose_f, h5gcreate_f, h5gopen_f, & + h5open_f, h5s_scalar_f, h5sclose_f, h5screate_f, h5screate_simple_f, & + h5sget_simple_extent_npoints_f, h5t_c_s1, h5t_cset_utf8_f, h5t_enum_f, h5t_native_double, & + h5t_native_integer, h5t_str_nullpad_f, h5t_string, h5tclose_f, h5tcopy_f, h5tcreate_f, & + h5tenum_insert_f, h5tset_cset_f, h5tset_size_f, h5tset_strpad_f, hid_t, hsize_t, size_t + USE iso_c_binding, ONLY: C_LOC,& + c_ptr + USE kinds, ONLY: dp +#include "./base/base_uses.f90" + + IMPLICIT NONE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hdf5_wrapper' + INTEGER, PARAMETER, PUBLIC :: hdf5_id = hid_t + +CONTAINS + +! ************************************************************************************************** +!> \brief Initialize the HDF5 fortran API +! ************************************************************************************************** + SUBROUTINE h5open() + INTEGER :: error + + CALL h5open_f(error) + IF (error < 0) CPABORT('ERROR: failed to initialize HDF5 interface') + + END SUBROUTINE h5open + +! ************************************************************************************************** +!> \brief Close the HDF5 fortran API +! ************************************************************************************************** + SUBROUTINE h5close() + INTEGER :: error + + CALL h5close_f(error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 interface') + + END SUBROUTINE h5close + +! ************************************************************************************************** +!> \brief Create a HDF5 file +!> \param filename the name of the hdf5 file +!> \param file_id the file id of the hdf5 file +! ************************************************************************************************** + SUBROUTINE h5fcreate(filename, file_id) + CHARACTER(LEN=*), INTENT(IN) :: filename + INTEGER(KIND=hid_t), INTENT(OUT) :: file_id + + INTEGER :: error + + CALL h5fcreate_f(filename, h5f_acc_trunc_f, file_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 file') + + END SUBROUTINE h5fcreate + +! ************************************************************************************************** +!> \brief Open a HDF5 file +!> \param filename the name of the hdf5 file +!> \param file_id the file id of the hdf5 file +! ************************************************************************************************** + SUBROUTINE h5fopen(filename, file_id) + CHARACTER(LEN=*), INTENT(IN) :: filename + INTEGER(KIND=hid_t), INTENT(OUT) :: file_id + + INTEGER :: error + + CALL h5fopen_f(TRIM(filename), h5f_acc_rdonly_f, file_id, error) + IF (error < 0) CPABORT('ERROR: failed to open HDF5 file') + + END SUBROUTINE h5fopen + +! ************************************************************************************************** +!> \brief Close a HDF5 file +!> \param file_id the file id of the hdf5 file +! ************************************************************************************************** + SUBROUTINE h5fclose(file_id) + INTEGER(KIND=hid_t), INTENT(IN) :: file_id + + INTEGER :: error + + CALL h5fclose_f(file_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 file') + + END SUBROUTINE h5fclose + +! ************************************************************************************************** +!> \brief Create a HDF5 group +!> \param loc_id file or group identifier +!> \param name name of the group +!> \param grp_id group identifier +! ************************************************************************************************** + SUBROUTINE h5gcreate(loc_id, name, grp_id) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(KIND=hid_t), INTENT(OUT) :: grp_id + + INTEGER :: error + + CALL h5gcreate_f(loc_id, name, grp_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 group') + + END SUBROUTINE h5gcreate + +! ************************************************************************************************** +!> \brief Open a HDF5 group +!> \param loc_id file or group identifier +!> \param name name of the group +!> \param grp_id group identifier +! ************************************************************************************************** + SUBROUTINE h5gopen(loc_id, name, grp_id) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(KIND=hid_t), INTENT(OUT) :: grp_id + + INTEGER :: error + + CALL h5gopen_f(loc_id, name, grp_id, error) + IF (error < 0) CPABORT('ERROR: failed to open HDF5 group') + + END SUBROUTINE h5gopen + +! ************************************************************************************************** +!> \brief Close a HDF5 group +!> \param grp_id group identifier +! ************************************************************************************************** + SUBROUTINE h5gclose(grp_id) + INTEGER(KIND=hid_t), INTENT(IN) :: grp_id + + INTEGER :: error + + CALL h5gclose_f(grp_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 group') + + END SUBROUTINE h5gclose + +! ************************************************************************************************** +!> \brief Write a variable-length string attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the string to write +! ************************************************************************************************** + SUBROUTINE h5awrite_varlen_string(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + CHARACTER(LEN=*), INTENT(IN), TARGET :: attr_data + + INTEGER :: error, output_unit + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + TYPE(c_ptr) :: buffer + TYPE(c_ptr), TARGET :: in_between_ptr + + ! create a scalar dataspace + CALL h5screate_f(h5s_scalar_f, space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 dataspace' + RETURN + END IF + + ! create a variable-length string type + CALL h5tcopy_f(h5t_string, type_id, error) + CALL h5tset_cset_f(type_id, h5t_cset_utf8_f, error) + CALL h5tset_strpad_f(type_id, h5t_str_nullpad_f, error) + + ! create the attribute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 attribute' + RETURN + END IF + + ! weird in-between pointer needed for variable-length + ! string to a scalar dataspace + in_between_ptr = C_LOC(attr_data) + ! the actual pointer to be passed + buffer = C_LOC(in_between_ptr) + + ! write the string attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to write HDF5 attribute' + RETURN + END IF + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 attribute' + RETURN + END IF + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 dataspace' + RETURN + END IF + + ! close datatype + CALL h5tclose_f(type_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 datatype' + RETURN + END IF + + END SUBROUTINE h5awrite_varlen_string + +! ************************************************************************************************** +!> \brief Write a fixed-length string attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the string to write +! ************************************************************************************************** + SUBROUTINE h5awrite_fixlen_string(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + CHARACTER(LEN=*), INTENT(IN), TARGET :: attr_data + + INTEGER :: error, output_unit + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + TYPE(c_ptr) :: buffer + + ! create a scalar dataspace + CALL h5screate_f(h5s_scalar_f, space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 dataspace' + RETURN + END IF + + ! create a fixed-length string datatype + CALL h5tcopy_f(h5t_c_s1, type_id, error) + CALL h5tset_cset_f(type_id, h5t_cset_utf8_f, error) + CALL h5tset_size_f(type_id, LEN(attr_data, size_t), error) + + ! create the attribute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 attribute' + RETURN + END IF + + ! the actual pointer to be passed + buffer = C_LOC(attr_data) + + ! write the string attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to write HDF5 attribute' + RETURN + END IF + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 attribute' + RETURN + END IF + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 dataspace' + RETURN + END IF + + ! close datatype + CALL h5tclose_f(type_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 datatype' + RETURN + END IF + + END SUBROUTINE h5awrite_fixlen_string + +! ************************************************************************************************** +!> \brief Write a boolean attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the logical to write (.true. or .false.) +! ************************************************************************************************** + SUBROUTINE h5awrite_boolean(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + LOGICAL, INTENT(IN) :: attr_data + + INTEGER :: error, output_unit + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + INTEGER(KIND=size_t) :: enum_size = 1 + INTEGER, TARGET :: attr_data_to_int + TYPE(c_ptr) :: buffer + + ! 8-bit integers in enum bool_type + + ! create a scalar dataspace + CALL h5screate_f(h5s_scalar_f, space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 dataspace' + RETURN + END IF + + ! create the datatype + CALL h5tcreate_f(h5t_enum_f, enum_size, type_id, error) + CALL h5tenum_insert_f(type_id, "FALSE", 0, error) + CALL h5tenum_insert_f(type_id, "TRUE", 1, error) + + IF (attr_data) THEN + attr_data_to_int = 1 + ELSE + attr_data_to_int = 0 + END IF + ! the C pointer to the actual data + buffer = C_LOC(attr_data_to_int) + + ! create the attribute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 attribute' + RETURN + END IF + + ! write the string attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to write HDF5 attribute' + RETURN + END IF + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 attribute' + RETURN + END IF + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 dataspace' + RETURN + END IF + + ! close datatype + CALL h5tclose_f(type_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 datatype' + RETURN + END IF + + END SUBROUTINE h5awrite_boolean + +! ************************************************************************************************** +!> \brief Write a (scalar) integer attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the integer to write +! ************************************************************************************************** + SUBROUTINE h5awrite_integer_scalar(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + INTEGER, INTENT(IN), TARGET :: attr_data + + INTEGER :: error, output_unit + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + TYPE(c_ptr) :: buffer + + ! create a scalar dataspace + CALL h5screate_f(h5s_scalar_f, space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 dataspace' + RETURN + END IF + + ! the C pointer to the actual data + buffer = C_LOC(attr_data) + + ! set the type of data + type_id = h5t_native_integer + + ! create the attribute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to create HDF5 attribute' + RETURN + END IF + + ! write the string attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to write HDF5 attribute' + RETURN + END IF + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 attribute' + RETURN + END IF + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & + ' ERROR: failed to close HDF5 dataspace' + RETURN + END IF + + END SUBROUTINE h5awrite_integer_scalar + +! ************************************************************************************************** +!> \brief Write a (scalar) double precision attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the double to write +! ************************************************************************************************** + SUBROUTINE h5awrite_double_scalar(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + REAL(KIND=dp), INTENT(IN), TARGET :: attr_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + TYPE(c_ptr) :: buffer + + ! create a scalar dataspace + CALL h5screate_f(h5s_scalar_f, space_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataspace') + + ! the C pointer to the actual data + buffer = C_LOC(attr_data) + + ! set the type of data + type_id = h5t_native_double + + ! create the attribute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 attribute') + + ! write the string attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to write HDF5 attribute') + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 attribute') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + END SUBROUTINE h5awrite_double_scalar + +! ************************************************************************************************** +!> \brief Write an array of fixed-length string attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the array of strings +! ************************************************************************************************** + SUBROUTINE h5awrite_string_simple(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + CHARACTER(LEN=*), DIMENSION(:), INTENT(IN), TARGET :: attr_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + INTEGER(KIND=hsize_t), DIMENSION(2) :: dims + TYPE(c_ptr) :: buffer + + dims(1) = LEN(attr_data(1), kind=hsize_t) ! length of a string entry + dims(2) = SIZE(attr_data, kind=hsize_t) ! length of array of strings + + ! create a fixed-length string datatype + CALL h5tcopy_f(h5t_c_s1, type_id, error) + CALL h5tset_cset_f(type_id, h5t_cset_utf8_f, error) + CALL h5tset_size_f(type_id, dims(1), error) + + ! create a simple dataspace + CALL h5screate_simple_f(1, dims(2:2), space_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataspace') + + ! create the atrtibute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 attribute') + + ! the actual pointer to be passed + buffer = C_LOC(attr_data(1)) + + ! write the string array attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to write HDF5 attribute') + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 attribute') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + ! close datatype + CALL h5tclose_f(type_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 datatype') + + END SUBROUTINE h5awrite_string_simple + +! ************************************************************************************************** +!> \brief Write an array of doubles attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the array of doubles +! ************************************************************************************************** + SUBROUTINE h5awrite_double_simple(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + REAL(KIND=dp), DIMENSION(:), INTENT(IN), TARGET :: attr_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + INTEGER(KIND=hsize_t), DIMENSION(1) :: dims + TYPE(c_ptr) :: buffer + + dims(1) = SIZE(attr_data, kind=hsize_t) ! length of array of strings + + ! set the type of data + type_id = h5t_native_double + + ! create a simple dataspace + CALL h5screate_simple_f(1, dims, space_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataspace') + + ! create the atrtibute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 attribute') + + ! the actual pointer to be passed + buffer = C_LOC(attr_data(1)) + + ! write the string array attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to write HDF5 attribute') + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 attribute') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + END SUBROUTINE h5awrite_double_simple + +! ************************************************************************************************** +!> \brief Write an array of integers attribute +!> \param loc_id either file id or group id +!> \param attr_name the name of the attribute +!> \param attr_data the attribute data, i.e. the array of integers +! ************************************************************************************************** + SUBROUTINE h5awrite_integer_simple(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + INTEGER, DIMENSION(:), INTENT(IN), TARGET :: attr_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: attr_id, space_id, type_id + INTEGER(KIND=hsize_t), DIMENSION(1) :: dims + TYPE(c_ptr) :: buffer + + dims(1) = SIZE(attr_data, kind=hsize_t) ! length of array of strings + + ! set the type of data + type_id = h5t_native_integer + + ! create a simple dataspace + CALL h5screate_simple_f(1, dims, space_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataspace') + + ! create the atrtibute + CALL h5acreate_f(loc_id, attr_name, type_id, space_id, attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 attribute') + + ! the actual pointer to be passed + buffer = C_LOC(attr_data(1)) + + ! write the string array attribute to file + CALL h5awrite_f(attr_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to write HDF5 attribute') + + ! close attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 attribute') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + END SUBROUTINE h5awrite_integer_simple + +! ************************************************************************************************** +!> \brief Write a dataset containing an array of doubles +!> \param loc_id either file id or group id +!> \param dset_name the name of the dataset +!> \param dset_data the dataset data, i.e. the array of doubles +! ************************************************************************************************** + SUBROUTINE h5dwrite_double_simple(loc_id, dset_name, dset_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: dset_name + REAL(KIND=dp), DIMENSION(:), INTENT(IN), TARGET :: dset_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: dset_id, space_id, type_id + INTEGER(KIND=hsize_t), DIMENSION(1) :: dims + TYPE(c_ptr) :: buffer + + dims(1) = SIZE(dset_data, kind=hsize_t) ! length of array + + ! set the type of data + type_id = h5t_native_double + + ! create a simple dataspace + CALL h5screate_simple_f(1, dims, space_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataspace') + + ! create the dataset + CALL h5dcreate_f(loc_id, dset_name, type_id, space_id, dset_id, error) + IF (error < 0) CPABORT('ERROR: failed to create HDF5 dataset') + + ! the actual pointer to be passed + buffer = C_LOC(dset_data(1)) + + ! write the string array attribute to file + CALL h5dwrite_f(dset_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to write HDF5 dataset') + + ! close dataset + CALL h5dclose_f(dset_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataset') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + END SUBROUTINE h5dwrite_double_simple + +! ************************************************************************************************** +!> \brief Read a dataset containing an array of doubles +!> \param loc_id either file id or group id +!> \param dset_name the name of the dataset +!> \param dset_data where the read dataset data will be written +! ************************************************************************************************** + SUBROUTINE h5dread_double_simple(loc_id, dset_name, dset_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: dset_name + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: dset_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: dset_id, npoints, space_id, type_id + INTEGER(KIND=hsize_t), DIMENSION(1) :: dims + + dims(1) = SIZE(dset_data, kind=hsize_t) ! length of array + + ! set the type of data + type_id = h5t_native_double + + ! open the dataset + CALL h5dopen_f(loc_id, dset_name, dset_id, error) + IF (error < 0) CPABORT('ERROR: failed to open HDF5 dataset') + + ! get information on the dataspace + CALL h5dget_space_f(dset_id, space_id, error) + IF (error < 0) CPABORT('ERROR: failed to fetch HDF5 dataspace info') + + ! get dataspace dims + CALL h5sget_simple_extent_npoints_f(space_id, npoints, error) + IF (error < 0) CPABORT('ERROR: failed to fetch HDF5 dataspace dimension') + + ! read the data + CALL h5dread_f(dset_id, type_id, dset_data, dims, error) + IF (error < 0) CPABORT('ERROR: failed to read HDF5 dataset') + + ! close dataset + CALL h5dclose_f(dset_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataset') + + ! close dataspace + CALL h5sclose_f(space_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 dataspace') + + END SUBROUTINE h5dread_double_simple + +! ************************************************************************************************** +!> \brief Read an attribute containing a scalar double +!> \param loc_id either file id or group id +!> \param attr_name ... +!> \param attr_data ... +! ************************************************************************************************** + SUBROUTINE h5aread_double_scalar(loc_id, attr_name, attr_data) + INTEGER(KIND=hid_t), INTENT(IN) :: loc_id + CHARACTER(LEN=*), INTENT(IN) :: attr_name + REAL(KIND=dp), INTENT(OUT), TARGET :: attr_data + + INTEGER :: error + INTEGER(KIND=hid_t) :: attr_id, type_id + TYPE(c_ptr) :: buffer + + ! set the type of data + type_id = h5t_native_double + + ! open the attribute + CALL h5aopen_f(loc_id, attr_name, attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to open HDF5 attribute') + + buffer = C_LOC(attr_data) + ! read the data + CALL h5aread_f(attr_id, type_id, buffer, error) + IF (error < 0) CPABORT('ERROR: failed to read HDF5 attribute') + + ! close the attribute + CALL h5aclose_f(attr_id, error) + IF (error < 0) CPABORT('ERROR: failed to close HDF5 attribute') + + END SUBROUTINE h5aread_double_scalar + +END MODULE hdf5_wrapper diff --git a/src/input_constants.F b/src/input_constants.F index 2466aadfc2..c6b822a118 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1097,11 +1097,15 @@ MODULE input_constants do_lri_opt_coeff = 1, & do_lri_opt_exps = 2 - ! Active smape model parameters + ! Active space model parameters INTEGER, PARAMETER, PUBLIC :: hf_model = 100, & rsdft_model = 101, & dmft_model = 102 + ! Active space external solvers + INTEGER, PARAMETER, PUBLIC :: no_solver = 200, & + qiskit_solver = 201 + ! callgraph parameters INTEGER, PARAMETER, PUBLIC :: callgraph_none = 0, & callgraph_master = 1, & diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 62acc54508..a3de6befba 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -74,11 +74,11 @@ MODULE input_cp2k_dft jacobian_fd1_central, jacobian_fd2, jacobian_fd2_backward, kg_color_dsatur, & kg_color_greedy, kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_embed_ri, kg_tnadd_none, & ls_2pnt, ls_3pnt, ls_gold, ls_none, manual_selection, mao_basis_ext, mao_basis_orb, & - mao_basis_prim, mao_projection, mopac_guess, no_excitations, no_guess, numerical, oe_gllb, & - oe_lb, oe_none, oe_saop, oe_sic, orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, & - orb_px, orb_py, orb_pz, orb_s, ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, & - ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, ot_poly_irac, & - ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, & + mao_basis_prim, mao_projection, mopac_guess, no_excitations, no_guess, no_solver, & + numerical, oe_gllb, oe_lb, oe_none, oe_saop, oe_sic, orb_dx2, orb_dxy, orb_dy2, orb_dyz, & + orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s, ot_algo_irac, ot_algo_taylor_or_diag, & + ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, & + ot_poly_irac, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, & ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, & ot_precond_solver_default, ot_precond_solver_direct, ot_precond_solver_inv_chol, & ot_precond_solver_update, outer_scf_basis_center_opt, outer_scf_becke_constraint, & @@ -87,9 +87,9 @@ MODULE input_cp2k_dft outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, & outer_scf_optimizer_none, outer_scf_optimizer_sd, outer_scf_optimizer_secant, & outer_scf_s2_constraint, plus_u_lowdin, plus_u_mulliken, plus_u_mulliken_charges, & - radius_covalent, radius_default, radius_single, radius_user, radius_vdw, random_guess, & - real_time_propagation, ref_charge_atomic, ref_charge_mulliken, rel_dkh, rel_none, & - rel_pot_erfc, rel_pot_full, rel_sczora_mp, rel_trans_atom, rel_trans_full, & + qiskit_solver, radius_covalent, radius_default, radius_single, radius_user, radius_vdw, & + random_guess, real_time_propagation, ref_charge_atomic, ref_charge_mulliken, rel_dkh, & + rel_none, rel_pot_erfc, rel_pot_full, rel_sczora_mp, rel_trans_atom, rel_trans_full, & rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, restart_guess, rsdft_model, & sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, & sccs_derivative_fft, sccs_fattebert_gygi, shape_function_density, shape_function_gaussian, & @@ -510,6 +510,10 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_active_space_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_dft_section ! ************************************************************************************************** @@ -1078,10 +1082,6 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) - CALL create_active_space_section(subsection) - CALL section_add_subsection(section, subsection) - CALL section_release(subsection) - CALL section_create(subsection, __LOCATION__, name="GAPW", & description="Controls the printing of some gapw related information (debug).", & n_keywords=0, n_subsections=1, repeats=.FALSE.) @@ -9054,24 +9054,9 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="INACTIVE_ELECTRONS", & - description="The number of inactive electrons in the CAS space", & - usage="INACTIVE_ELECTRONS 2 {2}", n_var=-1, default_i_vals=(/-1, -1/), & - type_of_var=integer_t) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="ACTIVE_ORBITALS", & description="The number of active orbitals defining the CAS space.", & - usage="ACTIVE_ORBITALS 2 {2}", n_var=-1, default_i_vals=(/-1, -1/), & - type_of_var=integer_t) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="INACTIVE_ORBITALS", & - description="The number of nactive orbitals.", & - usage="INACTIVE_ORBITALS 2 {2}", n_var=-1, default_i_vals=(/-1, -1/), & - type_of_var=integer_t) + usage="ACTIVE_ORBITALS 2", n_var=1, default_i_val=-1, type_of_var=integer_t) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -9122,16 +9107,42 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="READ_P_ACTIVE", & - description="Active density matrix must be read", & + CALL keyword_create(keyword, __LOCATION__, name="SCF_EMBEDDING", & + description="Whether to turn on the self-consistent embedding scheme", & default_l_val=.FALSE., lone_keyword_l_val=.FALSE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="P_ACTIVE_FILE_NAME", & - description="Name of the active DM file, may include a path", & - usage="P_ACTIVE_FILE_NAME ", & - type_of_var=lchar_t, repeats=.FALSE.) + CALL keyword_create(keyword, __LOCATION__, name="QCSCHEMA", & + description="Name of the QCSchema file, may include a path", & + usage="QCSCHEMA ", & + type_of_var=lchar_t, repeats=.FALSE., & + default_lc_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="AS_SOLVER", & + description="The external active space solver for the embedding approach", & + usage="AS_SOLVER QISKIT", & + default_i_val=no_solver, & + enum_c_vals=s2a("NONE", "QISKIT"), & + enum_i_vals=(/no_solver, qiskit_solver/), & + enum_desc=s2a("NO solver, used to produce FCIDUMP/QCSchema files", & + "QISKIT active space solver")) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_ITER", & + description="Energy convergence threshold of the DFT embedding scheme.", & + usage="EPS_ITER 1.0E-6 ", type_of_var=real_t, & + default_r_val=1.0E-6_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER", & + description="Max number of iterations for the DFT embedding scheme.", & + usage="MAX_ITER 50", type_of_var=integer_t, & + default_i_val=50) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -9151,8 +9162,50 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_socket_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_active_space_section +! ************************************************************************************************** +!> \brief ... +!> \param section ... +! ************************************************************************************************** + SUBROUTINE create_socket_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="SOCKET", & + description="Parameters to set up the socket communicating to the external active space solver.", & + n_keywords=3, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="UNIX", & + description="Use a UNIX socket rather than an INET socket.", & + usage="unix LOGICAL", & + default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="PORT", & + description="Port number for the socket client.", & + usage="port ", & + default_i_val=12345) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="HOST", & + description="Host name for the socket client.", & + usage="host ", & + default_c_val="embedding_socket") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_socket_section + ! ************************************************************************************************** !> \brief ... !> \param section ... @@ -9228,7 +9281,7 @@ CONTAINS enum_i_vals=(/eri_method_full_gpw, eri_method_gpw_ht/), & enum_desc=s2a("Use the GPW approach with MOs", & "Use the GPW approach for half-transformed MO ERIs"), & - default_i_val=eri_method_gpw_ht) + default_i_val=eri_method_full_gpw) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -9251,7 +9304,7 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="OPERATOR_PARAMETER", & description="Range parameter for ERI operator.", & usage="OPERATOR_PARAMETER 4.0", type_of_var=real_t, & - default_r_val=8.0_dp) + default_r_val=0.5_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/ipi_driver.F b/src/ipi_driver.F index e797b67458..f78e9611bf 100644 --- a/src/ipi_driver.F +++ b/src/ipi_driver.F @@ -12,39 +12,45 @@ !> \author Michele Ceriotti 03.2012 ! ************************************************************************************************** MODULE ipi_driver - USE ISO_C_BINDING, ONLY: C_CHAR,& - C_DOUBLE,& - C_INT,& - C_LOC,& - C_NULL_CHAR,& - C_PTR - USE bibliography, ONLY: Ceriotti2014,& - Kapil2016,& - cite_reference - USE cell_methods, ONLY: cell_create,& - init_cell - USE cell_types, ONLY: cell_release,& - cell_type - USE cp_external_control, ONLY: external_control - USE cp_log_handling, ONLY: cp_logger_get_default_io_unit - USE cp_subsys_types, ONLY: cp_subsys_get,& - cp_subsys_set,& - cp_subsys_type - USE force_env_methods, ONLY: force_env_calc_energy_force - USE force_env_types, ONLY: force_env_get,& - force_env_type - USE global_types, ONLY: global_environment_type - USE input_section_types, ONLY: section_vals_get_subs_vals,& - section_vals_type,& - section_vals_val_get - USE kinds, ONLY: default_path_length,& - default_string_length,& - dp,& - int_4 - USE message_passing, ONLY: mp_para_env_type,& - mp_request_type,& - mp_testany - USE virial_types, ONLY: virial_type + USE ISO_C_BINDING, ONLY: C_CHAR, & + C_DOUBLE, & + C_INT, & + C_LOC, & + C_NULL_CHAR, & + C_PTR + USE bibliography, ONLY: Ceriotti2014, & + Kapil2016, & + cite_reference + USE cell_methods, ONLY: cell_create, & + init_cell + USE cell_types, ONLY: cell_release, & + cell_type + USE cp_external_control, ONLY: external_control + USE cp_log_handling, ONLY: cp_logger_get_default_io_unit + USE cp_subsys_types, ONLY: cp_subsys_get, & + cp_subsys_set, & + cp_subsys_type + USE force_env_methods, ONLY: force_env_calc_energy_force + USE force_env_types, ONLY: force_env_get, & + force_env_type + USE global_types, ONLY: global_environment_type + USE input_section_types, ONLY: section_vals_get_subs_vals, & + section_vals_type, & + section_vals_val_get + USE kinds, ONLY: default_path_length, & + default_string_length, & + dp, & + int_4 + USE message_passing, ONLY: mp_para_env_type, & + mp_request_type, & + mp_testany +#ifndef __NO_SOCKETS + USE sockets_interface, ONLY: writebuffer, & + readbuffer, & + open_connect_socket, & + uwait +#endif + USE virial_types, ONLY: virial_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -55,244 +61,8 @@ MODULE ipi_driver PUBLIC :: run_driver -#ifndef __NO_IPI_DRIVER - INTERFACE writebuffer - MODULE PROCEDURE writebuffer_s, & - writebuffer_d, writebuffer_dv, & - writebuffer_i - - END INTERFACE - - INTERFACE readbuffer - MODULE PROCEDURE readbuffer_s, & - readbuffer_dv, readbuffer_d, & - readbuffer_i - - END INTERFACE - - INTERFACE - SUBROUTINE uwait(sec) BIND(C, NAME="uwait") - USE ISO_C_BINDING, ONLY: C_DOUBLE - REAL(C_DOUBLE) :: sec - - END SUBROUTINE - END INTERFACE - - INTERFACE - SUBROUTINE open_socket(psockfd, inet, port, host) BIND(C) - IMPORT - INTEGER(KIND=C_INT) :: psockfd, inet, port - CHARACTER(KIND=C_CHAR), DIMENSION(*) :: host - - END SUBROUTINE open_socket - - SUBROUTINE writebuffer_csocket(psockfd, pdata, plen) BIND(C, name="writebuffer") - IMPORT - INTEGER(KIND=C_INT) :: psockfd - TYPE(C_PTR), VALUE :: pdata - INTEGER(KIND=C_INT) :: plen - - END SUBROUTINE writebuffer_csocket - - SUBROUTINE readbuffer_csocket(psockfd, pdata, plen) BIND(C, name="readbuffer") - IMPORT - INTEGER(KIND=C_INT) :: psockfd - TYPE(C_PTR), VALUE :: pdata - INTEGER(KIND=C_INT) :: plen - - END SUBROUTINE readbuffer_csocket - END INTERFACE -#endif - CONTAINS -#ifndef __NO_IPI_DRIVER -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -! ************************************************************************************************** - SUBROUTINE writebuffer_d(psockfd, fdata) - INTEGER, INTENT(IN) :: psockfd - REAL(KIND=dp), INTENT(IN) :: fdata - - CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_d' - - INTEGER :: handle - REAL(KIND=C_DOUBLE), TARGET :: cdata - - CALL timeset(routineN, handle) - - cdata = fdata - CALL writebuffer_csocket(psockfd, c_loc(cdata), 8) - - CALL timestop(handle) - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -! ************************************************************************************************** - SUBROUTINE writebuffer_i(psockfd, fdata) - INTEGER, INTENT(IN) :: psockfd, fdata - - CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_i' - - INTEGER :: handle - INTEGER(KIND=C_INT), TARGET :: cdata - - CALL timeset(routineN, handle) - - cdata = fdata - CALL writebuffer_csocket(psockfd, c_loc(cdata), 4) - - CALL timestop(handle) - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fstring ... -!> \param plen ... -! ************************************************************************************************** - SUBROUTINE writebuffer_s(psockfd, fstring, plen) - INTEGER, INTENT(IN) :: psockfd - CHARACTER(LEN=*), INTENT(IN) :: fstring - INTEGER, INTENT(IN) :: plen - - CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_s' - - INTEGER :: handle, i - CHARACTER(LEN=1, KIND=C_CHAR), TARGET :: cstring(plen) - - CALL timeset(routineN, handle) - - DO i = 1, plen - cstring(i) = fstring(i:i) - END DO - CALL writebuffer_csocket(psockfd, c_loc(cstring(1)), plen) - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -!> \param plen ... -! ************************************************************************************************** - SUBROUTINE writebuffer_dv(psockfd, fdata, plen) - INTEGER, INTENT(IN) :: psockfd, plen - REAL(KIND=dp), INTENT(IN), TARGET :: fdata(plen) - - CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_dv' - - INTEGER :: handle - - CALL timeset(routineN, handle) - - CALL writebuffer_csocket(psockfd, c_loc(fdata(1)), 8*plen) - - CALL timestop(handle) - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -! ************************************************************************************************** - SUBROUTINE readbuffer_d(psockfd, fdata) - INTEGER, INTENT(IN) :: psockfd - REAL(KIND=dp), INTENT(OUT) :: fdata - - CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_d' - - INTEGER :: handle - REAL(KIND=C_DOUBLE), TARGET :: cdata - - CALL timeset(routineN, handle) - - CALL readbuffer_csocket(psockfd, c_loc(cdata), 8) - fdata = cdata - - CALL timestop(handle) - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -! ************************************************************************************************** - SUBROUTINE readbuffer_i(psockfd, fdata) - INTEGER, INTENT(IN) :: psockfd - INTEGER, INTENT(OUT) :: fdata - - CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_i' - - INTEGER :: handle - INTEGER(KIND=C_INT), TARGET :: cdata - - CALL timeset(routineN, handle) - - CALL readbuffer_csocket(psockfd, c_loc(cdata), 4) - fdata = cdata - - CALL timestop(handle) - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fstring ... -!> \param plen ... -! ************************************************************************************************** - SUBROUTINE readbuffer_s(psockfd, fstring, plen) - INTEGER, INTENT(IN) :: psockfd - CHARACTER(LEN=*), INTENT(OUT) :: fstring - INTEGER, INTENT(IN) :: plen - - CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_s' - - INTEGER :: handle, i - CHARACTER(LEN=1, KIND=C_CHAR), TARGET :: cstring(plen) - - CALL timeset(routineN, handle) - - CALL readbuffer_csocket(psockfd, c_loc(cstring(1)), plen) - fstring = "" - DO i = 1, plen - fstring(i:i) = cstring(i) - END DO - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param psockfd ... -!> \param fdata ... -!> \param plen ... -! ************************************************************************************************** - SUBROUTINE readbuffer_dv(psockfd, fdata, plen) - INTEGER, INTENT(IN) :: psockfd, plen - REAL(KIND=dp), INTENT(OUT), TARGET :: fdata(plen) - - CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_dv' - - INTEGER :: handle - - CALL timeset(routineN, handle) - - CALL readbuffer_csocket(psockfd, c_loc(fdata(1)), 8*plen) - - CALL timestop(handle) - - END SUBROUTINE -#endif - ! ************************************************************************************************** !> \brief ... !> \param force_env ... @@ -308,10 +78,10 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'run_driver' -#ifdef __NO_IPI_DRIVER +#ifdef __NO_SOCKETS INTEGER :: handle CALL timeset(routineN, handle) - CPABORT("CP2K was compiled with the __NO_IPI_DRIVER option!") + CPABORT("CP2K was compiled with the __NO_SOCKETS option!") MARK_USED(globenv) MARK_USED(force_env) #else @@ -374,7 +144,7 @@ CONTAINS END IF c_hostname = TRIM(drv_hostname)//C_NULL_CHAR - IF (ionode) CALL open_socket(socket, i_drv_unix, drv_port, c_hostname) + IF (ionode) CALL open_connect_socket(socket, i_drv_unix, drv_port, c_hostname) NULLIFY (wait_msg) ALLOCATE (wait_msg(1)) diff --git a/src/qcschema.F b/src/qcschema.F new file mode 100644 index 0000000000..85ac01847f --- /dev/null +++ b/src/qcschema.F @@ -0,0 +1,1036 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2023 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief The module to read/write QCSchema HDF5 files for interfacing CP2K with other programs +!> \par History +!> 10.2022 created [SB] +!> \author Stefano Battaglia +! ************************************************************************************************** +MODULE qcschema + + USE atomic_kind_types, ONLY: get_atomic_kind + USE basis_set_types, ONLY: gto_basis_set_type + USE cp2k_info, ONLY: cp2k_version + USE cp_control_types, ONLY: dft_control_type + USE cp_log_handling, ONLY: cp_get_default_logger, & + cp_logger_get_default_io_unit, & + cp_logger_type +#ifdef __HDF5 + USE hdf5_wrapper, ONLY: & + h5aread_double_scalar, h5awrite_boolean, h5awrite_double_scalar, h5awrite_double_simple, & + h5awrite_fixlen_string, h5awrite_integer_scalar, h5awrite_integer_simple, & + h5awrite_string_simple, h5close, h5dread_double_simple, h5dwrite_double_simple, h5fclose, & + h5fcreate, h5fopen, h5gclose, h5gcreate, h5gopen, h5open, hdf5_id +#endif + USE input_section_types, ONLY: section_vals_get, & + section_vals_get_subs_vals, & + section_vals_type + USE kinds, ONLY: default_path_length, & + default_string_length, & + dp, & + int_8 + USE mp2_types, ONLY: mp2_type + USE particle_types, ONLY: particle_type + USE periodic_table, ONLY: get_ptable_info + USE qs_active_space_types, ONLY: active_space_type + USE qs_active_space_utils, ONLY: eri_to_array, & + subspace_matrix_to_array + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env, & + qs_environment_type + USE qs_force_types, ONLY: qs_force_type + USE qs_kind_types, ONLY: get_qs_kind, & + qs_kind_type + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_scf_types, ONLY: qs_scf_env_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qcschema' + + PUBLIC :: qcschema_type + PUBLIC :: qcschema_env_create, qcschema_env_release, qcschema_to_hdf5 + +! ************************************************************************************************** +!> \brief A derived type to store the program information that generated the QCSchema file. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#provenance +! ************************************************************************************************** + TYPE qcschema_provenance + CHARACTER(LEN=default_string_length) :: creator ! The name of the creator of this object + CHARACTER(LEN=default_string_length) :: version ! The version of the creator of this object + CHARACTER(LEN=default_string_length) :: routine ! The routine that was used to create this object + END TYPE qcschema_provenance + +! ************************************************************************************************** +!> \brief A derived type to store the topological information of the physical system. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#topology +! ************************************************************************************************** + TYPE qcschema_topology + CHARACTER(LEN=default_string_length) :: name ! of the molecule + CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: symbols ! of the atoms + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: geometry ! row major, in bohr + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: masses + INTEGER, DIMENSION(:), ALLOCATABLE :: atomic_numbers + INTEGER :: molecular_charge + INTEGER :: molecular_multiplicity + CHARACTER(LEN=default_string_length) :: schema_name + INTEGER :: schema_version + TYPE(qcschema_provenance) :: provenance + END TYPE qcschema_topology + +! ************************************************************************************************** +!> \brief A derived type to store the information of a single electron shell in a basis set. +!> For more information refer to: +!> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L43 +! ************************************************************************************************** + TYPE qcschema_electron_shell + ! The angular momenta of this electron shell as a list of integers + INTEGER, DIMENSION(:), POINTER :: angular_momentum => Null() + ! The type of this shell: spherical or cartesian + CHARACTER(LEN=9) :: harmonic_type + ! The exponents of this contracted shell. The official spec stores these values as strings + REAL(KIND=dp), DIMENSION(:), POINTER :: exponents => Null() + ! The general contraction coefficients of this contracted shell + REAL(KIND=dp), DIMENSION(:, :), POINTER :: coefficients => Null() + END TYPE qcschema_electron_shell + +! ************************************************************************************************** +!> \brief A derived type to store the information of an ECP in a basis set. +!> For more information refer to: +!> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L90 +! ************************************************************************************************** + TYPE qcschema_ecp + ! The type of this potential + CHARACTER(LEN=default_string_length) :: ecp_type + ! The angular momenta of this potential as a list of integers + INTEGER, DIMENSION(:), POINTER :: angular_momentum => Null() + ! The exponents of the r terms + INTEGER, DIMENSION(:), POINTER :: r_exponents => Null() + ! The exponents of the Gaussian terms + REAL(KIND=dp), DIMENSION(:), POINTER :: gaussian_exponents => Null() + ! The general contraction coefficients of this potential + REAL(KIND=dp), DIMENSION(:, :), POINTER :: coefficients => Null() + END TYPE qcschema_ecp + +! ************************************************************************************************** +!> \brief A derived type to store the information of a single atom/center in the basis. +!> For more information refer to: +!> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L146 +! ************************************************************************************************** + TYPE qcschema_center_basis + ! The list of electronic shells for this element + TYPE(qcschema_electron_shell), DIMENSION(:), POINTER :: electron_shells => Null() + ! The list of effective core potentials for this element + TYPE(qcschema_ecp), DIMENSION(:), POINTER :: ecp_potentials => Null() + ! The number of electrons replaced by an ECP + INTEGER :: ecp_electrons + END TYPE qcschema_center_basis + +! ************************************************************************************************** +!> \brief A derived type to store the information of the basis set used in the calculation. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/auto_basis.html#basis-set-schema +! ************************************************************************************************** + TYPE qcschema_basis_set + ! The name of the basis set + CHARACTER(LEN=default_string_length) :: name + ! A dictionary mapping the keys provided by `atom_map` to their basis center data + TYPE(qcschema_center_basis), DIMENSION(:), POINTER :: center_data => Null() + ! The list of atomic kinds, indicating the keys used to store the basis in `center_data` + ! Not clear if this will be of the length of the basis set size, or rather just one + ! entry for atomic kind. E.g. only one entry for hydrogen even though there might be + ! many hydrogen atoms in the molecule. If this is the case, then we really need a + ! hash table for `center_data` + CHARACTER(LEN=2), DIMENSION(:), POINTER :: atom_map => Null() + ! The version of this specific schema + INTEGER :: schema_version + ! The name of this schema. This value is expected to be `qcschema_basis` + CHARACTER(LEN=default_string_length) :: schema_name + ! A description of this basis set + CHARACTER(LEN=default_string_length) :: description + END TYPE qcschema_basis_set + +! ************************************************************************************************** +!> \brief A derived type to store any additional computed wavefunction properties. +!> Matrix quantities are stored as flat, column-major arrays. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/auto_wf.html#wavefunction-schema +! ************************************************************************************************** + TYPE qcschema_wavefunction + + ! The name of the method used to obtain the wf + CHARACTER(LEN=default_string_length) :: method + + ! The basis set used during the computation + TYPE(qcschema_basis_set) :: basis_set + + ! SCF quantities in AO or MO basis + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_orbitals_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_orbitals_b + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eigenvalues_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eigenvalues_b + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_occupations_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_occupations_b + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_density_mo_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_density_mo_b + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_fock_mo_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_fock_mo_b + + ! Electron repulsion integrals + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_aa + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_ab + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_bb + + ! Quantities with localized orbitals. All `nmo` orbitals are included, + ! even if only a subset were localized + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_orbitals_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_orbitals_b + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_fock_a + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_fock_b + + ! Whether the computation used restricted spin orbitals + LOGICAL :: restricted + + END TYPE qcschema_wavefunction + +! ************************************************************************************************** +!> \brief A derived type to store the computed properties of the original calculation. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/auto_props.html#properties-schema +! ************************************************************************************************** + TYPE qcschema_properties + + REAL(KIND=dp) :: return_energy = 0.0_dp + + INTEGER :: calcinfo_nbasis = 0 ! AO basis size + INTEGER :: calcinfo_nmo = 0 ! MO basis size + INTEGER :: calcinfo_nalpha = 0 ! # of alpha electrons + INTEGER :: calcinfo_nbeta = 0 ! # of beta electrons + INTEGER :: calcinfo_natom = 0 + + ! SCF results + INTEGER :: scf_iterations = 0 + REAL(KIND=dp) :: scf_one_electron_energy = 0.0_dp + REAL(KIND=dp) :: scf_two_electron_energy = 0.0_dp + REAL(KIND=dp) :: nuclear_repulsion_energy = 0.0_dp + REAL(KIND=dp) :: scf_vv10_energy = 0.0_dp + REAL(KIND=dp) :: scf_xc_energy = 0.0_dp + REAL(KIND=dp) :: scf_dispersion_correction_energy = 0.0_dp + REAL(KIND=dp) :: scf_total_energy = 0.0_dp + ! the dipole moment is calculated on the fly and not stored + REAL(KIND=dp), DIMENSION(3) :: scf_dipole_moment = 0.0_dp + + ! MP2 results + REAL(KIND=dp) :: mp2_same_spin_correlation_energy = 0.0_dp + REAL(KIND=dp) :: mp2_opposite_spin_correlation_energy = 0.0_dp + REAL(KIND=dp) :: mp2_singles_energy = 0.0_dp + REAL(KIND=dp) :: mp2_doubles_energy = 0.0_dp + ! these are the only two that are saved + REAL(KIND=dp) :: mp2_correlation_energy = 0.0_dp + REAL(KIND=dp) :: mp2_total_energy = 0.0_dp + + ! internal flags to know the type of calculation + LOGICAL :: mp2 = .FALSE. + + END TYPE qcschema_properties + +! ************************************************************************************************** +!> \brief The full QCSchema output type. +!> For more information refer to: +!> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#output-components +! ************************************************************************************************** + TYPE qcschema_type + TYPE(qcschema_topology) :: topology + TYPE(qcschema_provenance) :: provenance + TYPE(qcschema_properties) :: properties + TYPE(qcschema_wavefunction) :: wavefunction + TYPE(qcschema_basis_set) :: basis + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: return_result + CHARACTER(LEN=default_string_length) :: driver + LOGICAL :: success + END TYPE qcschema_type + +CONTAINS + +! ************************************************************************************************** +!> \brief Create and initialize a qcschema object from a quickstep environment +!> \param qcschema_env the qcschema environment to populate +!> \param qs_env the qs environment with all the info of the computation +! ************************************************************************************************** + SUBROUTINE qcschema_env_create(qcschema_env, qs_env) + TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env + TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'qcschema_env_create' + + CHARACTER(LEN=2) :: atomic_symbol + CHARACTER(LEN=default_string_length) :: basis_set_name, method + INTEGER :: atomic_number, handle, i, i_glb, iatom, & + ikind, nalpha, nao, natoms, nbeta, & + nel, nmo, nspins, output_unit + LOGICAL :: do_hfx + REAL(KIND=dp) :: dispersion, mass, one_el_en, two_el_en + TYPE(active_space_type), POINTER :: active_space_env + TYPE(cp_logger_type), POINTER :: logger + TYPE(dft_control_type), POINTER :: dft_control + TYPE(gto_basis_set_type), POINTER :: basis_set + TYPE(mp2_type), POINTER :: mp2_env + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(qs_scf_env_type), POINTER :: scf_env + TYPE(section_vals_type), POINTER :: hfx_sections, input + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + output_unit = cp_logger_get_default_io_unit(logger) + + ! reset everything + CALL qcschema_env_release(qcschema_env) + + ! collect environment info + IF (ASSOCIATED(qs_env)) THEN + CALL get_qs_env(qs_env, ks_env=ks_env, energy=energy, & + dft_control=dft_control, force=force, & + particle_set=particle_set, & + scf_env=scf_env, mp2_env=mp2_env, & + input=input, qs_kind_set=kind_set, & + active_space=active_space_env) + ELSE + CPABORT("QS environment not associated, QCSchema interface quitting") + END IF + + ! we need the AS environemnt to get all the SCF data + IF (.NOT. ASSOCIATED(active_space_env)) THEN + CPABORT("Active space environment not associated, QCSchema interface quitting") + END IF + + !========================================================================================! + ! *** QCSchema provenance *** + !========================================================================================! + + qcschema_env%provenance%creator = 'CP2K' + qcschema_env%provenance%version = cp2k_version + qcschema_env%provenance%routine = routineN + + !========================================================================================! + ! *** QCSchema topology *** + !========================================================================================! + + qcschema_env%topology%schema_name = 'qcschema' + qcschema_env%topology%schema_version = 3 + + natoms = SIZE(particle_set) + + ALLOCATE (qcschema_env%topology%geometry(3*natoms)) + ALLOCATE (qcschema_env%topology%symbols(natoms)) + ALLOCATE (qcschema_env%topology%atomic_numbers(natoms)) + ALLOCATE (qcschema_env%topology%masses(natoms)) + + DO iatom = 1, natoms + ! set the geometry as a flat array + qcschema_env%topology%geometry((iatom - 1)*3 + 1:(iatom)*3) = particle_set(iatom)%r(1:3) + + ! set the atomic symbols + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=atomic_symbol) + qcschema_env%topology%symbols(iatom) = atomic_symbol + + ! set the atomic numbers and masses + CALL get_ptable_info(atomic_symbol, number=atomic_number, amass=mass) + qcschema_env%topology%atomic_numbers(iatom) = atomic_number + qcschema_env%topology%masses(iatom) = mass + END DO + + qcschema_env%topology%molecular_charge = dft_control%charge + qcschema_env%topology%molecular_multiplicity = dft_control%multiplicity + + !========================================================================================! + ! *** QCSchema properties *** + !========================================================================================! + + nspins = active_space_env%nspins + + nao = active_space_env%mos_active(1)%nao + nmo = active_space_env%nmo_active + nel = active_space_env%nelec_active + + IF (nspins == 1) THEN + nalpha = active_space_env%nelec_active/2 + nbeta = nalpha + ELSE + nalpha = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2 + nbeta = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2 + END IF + + qcschema_env%properties%calcinfo_natom = natoms + qcschema_env%properties%calcinfo_nbasis = nao + qcschema_env%properties%calcinfo_nmo = nmo + qcschema_env%properties%calcinfo_nalpha = nalpha + qcschema_env%properties%calcinfo_nbeta = nbeta + + ! energy results + qcschema_env%properties%return_energy = energy%total + qcschema_env%properties%scf_total_energy = energy%total + ! here we abuse the nuclear repulsion energy to store the inactive energy + qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive + ! SCF info + qcschema_env%properties%scf_iterations = scf_env%iter_count + ! one-electron energy is the sum of all core terms + one_el_en = energy%core_overlap + energy%core_self + energy%core + qcschema_env%properties%scf_two_electron_energy = one_el_en + ! two-electron energy is the sum of hartree and exact exchange (if there) + two_el_en = energy%hartree + energy%ex + energy%hartree_1c + qcschema_env%properties%scf_one_electron_energy = two_el_en + ! xc energy + qcschema_env%properties%scf_xc_energy = & + energy%exc + energy%exc_aux_fit + energy%exc1 + energy%exc1_aux_fit + ! dispersion energy + dispersion = energy%dispersion + energy%gcp + qcschema_env%properties%scf_dispersion_correction_energy = dispersion + + ! Some methods of CP2K are not supported by QCSchema, let's warn the user + IF (dft_control%smear) CPABORT('WARNING: smearing not supported in QCSchema') + IF (dft_control%dft_plus_u) CPABORT('WARNING: DFT+U not supported in QCSchema') + IF (dft_control%do_sccs) CPABORT('WARNING: SCCS not supported in QCSchema') + IF (qs_env%qmmm) CPABORT('WARNING: QM/MM not supported in QCSchema') + IF (dft_control%qs_control%mulliken_restraint) & + CPABORT('WARNING: Mulliken restrains not supported in QCSchema') + IF (dft_control%qs_control%semi_empirical) & + CPABORT('WARNING: semi_empirical methods not supported in QCSchema') + IF (dft_control%qs_control%dftb) CPABORT('WARNING: DFTB not supported in QCSchema') + IF (dft_control%qs_control%xtb) CPABORT('WARNING: xTB not supported in QCSchema') + + ! MP2 info + IF (ASSOCIATED(qs_env%mp2_env)) THEN + qcschema_env%properties%mp2 = .TRUE. + ! this info is computed on the fly, but not stored! + ! qcschema_env%properties%mp2_same_spin_correlation_energy + ! qcschema_env%properties%mp2_opposite_spin_correlation_energy + + qcschema_env%properties%mp2_correlation_energy = energy%mp2 + qcschema_env%properties%mp2_total_energy = energy%total + + ! update the scf energy + qcschema_env%properties%scf_total_energy = energy%total - energy%mp2 + END IF + + !========================================================================================! + ! *** QCSchema wavefunction *** + !========================================================================================! + + IF (nspins == 1) THEN + qcschema_env%wavefunction%restricted = .TRUE. + ELSE + qcschema_env%wavefunction%restricted = .FALSE. + END IF + + ! alpha MO energies + ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a(nmo)) + DO i = 1, nmo + i_glb = active_space_env%active_orbitals(i, 1) + qcschema_env%wavefunction%scf_eigenvalues_a(i) = & + active_space_env%mos_active(1)%eigenvalues(i_glb) + END DO + + ! alpha MO occupations + ALLOCATE (qcschema_env%wavefunction%scf_occupations_a(nmo)) + DO i = 1, nmo + i_glb = active_space_env%active_orbitals(i, 1) + qcschema_env%wavefunction%scf_occupations_a(i) = & + active_space_env%mos_active(1)%occupation_numbers(i_glb) + END DO + + ! alpha Fock matrix + ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a(nmo*nmo)) + CALL subspace_matrix_to_array(active_space_env%fock_sub(1), & + qcschema_env%wavefunction%scf_fock_mo_a, & + active_space_env%active_orbitals(:, 1), & + active_space_env%active_orbitals(:, 1)) + + ! alpha density matrix + ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo)) + CALL subspace_matrix_to_array(active_space_env%p_active(1), & + qcschema_env%wavefunction%scf_density_mo_a, & + active_space_env%active_orbitals(:, 1), & + active_space_env%active_orbitals(:, 1)) + + ! alpha MOs coefficients + ALLOCATE (qcschema_env%wavefunction%scf_orbitals_a(nao*nmo)) + CALL subspace_matrix_to_array(active_space_env%mos_active(1)%mo_coeff, & + qcschema_env%wavefunction%scf_orbitals_a, & + (/(i, i=1, nao)/), active_space_env%active_orbitals(:, 1)) + + IF (nspins == 2) THEN + ! beta MO energies + ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b(nmo)) + DO i = 1, nmo + i_glb = active_space_env%active_orbitals(i, 2) + qcschema_env%wavefunction%scf_eigenvalues_b(i) = & + active_space_env%mos_active(2)%eigenvalues(i_glb) + END DO + + ! beta MO occupations + ALLOCATE (qcschema_env%wavefunction%scf_occupations_b(nmo)) + DO i = 1, nmo + i_glb = active_space_env%active_orbitals(i, 2) + qcschema_env%wavefunction%scf_occupations_b(i) = & + active_space_env%mos_active(2)%occupation_numbers(i_glb) + END DO + + ! beta Fock matrix + ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b(nmo*nmo)) + CALL subspace_matrix_to_array(active_space_env%fock_sub(2), & + qcschema_env%wavefunction%scf_fock_mo_b, & + active_space_env%active_orbitals(:, 2), & + active_space_env%active_orbitals(:, 2)) + + ! beta density matrix + ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo)) + CALL subspace_matrix_to_array(active_space_env%p_active(2), & + qcschema_env%wavefunction%scf_density_mo_b, & + active_space_env%active_orbitals(:, 2), & + active_space_env%active_orbitals(:, 2)) + + ! beta MOs coefficients + ALLOCATE (qcschema_env%wavefunction%scf_orbitals_b(nao*nmo)) + CALL subspace_matrix_to_array(active_space_env%mos_active(2)%mo_coeff, & + qcschema_env%wavefunction%scf_orbitals_b, & + (/(i, i=1, nao)/), active_space_env%active_orbitals(:, 2)) + END IF + + ! get the alpha-alpha eri + ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa(nmo**4)) + CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_aa, & + active_space_env%active_orbitals, 1, 1) + + IF (nspins == 2) THEN + ! get the alpha-beta eri + ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab(nmo**4)) + CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_ab, & + active_space_env%active_orbitals, 1, 2) + + ! get the beta-beta eri + ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb(nmo**4)) + CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_bb, & + active_space_env%active_orbitals, 2, 2) + END IF + + !========================================================================================! + ! *** QCSchema model *** + !========================================================================================! + + DO iatom = 1, natoms + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) + CALL get_qs_kind(kind_set(ikind), basis_set=basis_set) + + basis_set_name = basis_set%name + + ! make sure that we do not run a mixed basis set + IF (iatom > 1) THEN + CPASSERT(basis_set_name == basis_set%name) + END IF + END DO + qcschema_env%wavefunction%basis_set%name = basis_set_name + + ! figure out which method was used for the calculation + IF (dft_control%uks) THEN + method = 'U' + ELSE IF (dft_control%roks) THEN + method = 'RO' + ELSE + method = 'R' + END IF + + hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") + CALL section_vals_get(hfx_sections, explicit=do_hfx) + + IF (do_hfx) THEN + method = TRIM(method)//'HF' + ELSE IF (qcschema_env%properties%mp2) THEN + method = TRIM(method)//'MP2' + ELSE + method = TRIM(method)//'KS' + END IF + + qcschema_env%wavefunction%method = TRIM(method) + + !========================================================================================! + ! *** QCSchema root *** + !========================================================================================! + + ! driver + IF (ASSOCIATED(force)) THEN + qcschema_env%driver = 'gradient' + ELSE + qcschema_env%driver = 'energy' + END IF + + ! success + ! TODO: how to check if the calculation was succesful? + qcschema_env%success = .TRUE. + + ! return result + IF (qcschema_env%success) THEN + IF (qcschema_env%driver == 'energy') THEN + ALLOCATE (qcschema_env%return_result(1)) + qcschema_env%return_result(1) = energy%total + ELSE + ALLOCATE (qcschema_env%return_result(3*SIZE(particle_set))) + ! TODO: populate with forces!! + qcschema_env%return_result = 0.0_dp + END IF + ELSE + CPABORT("The calculation to build the AS is unsuccessful") + END IF + + CALL timestop(handle) + + END SUBROUTINE qcschema_env_create + +! ************************************************************************************************** +!> \brief Releases the allocated memory of a qcschema environment +!> \param qcschema_env the qcschema environment to release +! ************************************************************************************************** + SUBROUTINE qcschema_env_release(qcschema_env) + TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env + + IF (ALLOCATED(qcschema_env%return_result)) THEN + DEALLOCATE (qcschema_env%return_result) + END IF + + IF (ALLOCATED(qcschema_env%topology%atomic_numbers)) THEN + DEALLOCATE (qcschema_env%topology%atomic_numbers) + END IF + + IF (ALLOCATED(qcschema_env%topology%masses)) THEN + DEALLOCATE (qcschema_env%topology%masses) + END IF + + IF (ALLOCATED(qcschema_env%topology%geometry)) THEN + DEALLOCATE (qcschema_env%topology%geometry) + END IF + + IF (ALLOCATED(qcschema_env%topology%symbols)) THEN + DEALLOCATE (qcschema_env%topology%symbols) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_occupations_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_occupations_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eri)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eri) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_aa)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_bb)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_ab)) THEN + DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_b) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_a)) THEN + DEALLOCATE (qcschema_env%wavefunction%localized_fock_a) + END IF + + IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_b)) THEN + DEALLOCATE (qcschema_env%wavefunction%localized_fock_b) + END IF + + END SUBROUTINE qcschema_env_release + +! ************************************************************************************************** +!> \brief Updates the Fock matrix and the inactive energy in a qcschema object +!> \param qcschema_env the qcschema environment +!> \param active_space_env the active space environment with the updated data +! ************************************************************************************************** + SUBROUTINE qcschema_update_fock(qcschema_env, active_space_env) + TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env + TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env + + ! alpha Fock matrix + CALL subspace_matrix_to_array(active_space_env%fock_sub(1), & + qcschema_env%wavefunction%scf_fock_mo_a, & + active_space_env%active_orbitals(:, 1), & + active_space_env%active_orbitals(:, 1)) + + ! beta Fock matrix + IF (active_space_env%nspins == 2) THEN + CALL subspace_matrix_to_array(active_space_env%fock_sub(2), & + qcschema_env%wavefunction%scf_fock_mo_b, & + active_space_env%active_orbitals(:, 2), & + active_space_env%active_orbitals(:, 2)) + END IF + + ! update inactive energy + qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive + + END SUBROUTINE qcschema_update_fock + +! ************************************************************************************************** +!> \brief Writes a qcschema object to an hdf5 file +!> \param qcschema_env the qcschema environment to write to file +!> \param filename ... +! ************************************************************************************************** + SUBROUTINE qcschema_to_hdf5(qcschema_env, filename) + TYPE(qcschema_type), INTENT(IN) :: qcschema_env + CHARACTER(LEN=default_path_length), INTENT(IN) :: filename +#ifndef __HDF5 + CPABORT("CP2K was compiled without the HDF5 library") + MARK_USED(filename) + MARK_USED(qcschema_env) +#else + INTEGER :: output_unit + INTEGER(KIND=hdf5_id) :: file_id, group_id + INTEGER(KIND=int_8) :: nresult + TYPE(cp_logger_type), POINTER :: logger + + logger => cp_get_default_logger() + output_unit = cp_logger_get_default_io_unit(logger) + + ! initialize HDF5 Fortran API + CALL h5open() + + ! create qcschema hdf5 file + ! filename = TRIM(logger%iter_info%project_name) // 'hdf5' + CALL h5fcreate(TRIM(filename), file_id) + + ! !===========================================================================! + ! *** Root group *** + ! !===========================================================================! + ! driver + CALL h5awrite_fixlen_string(file_id, 'driver', TRIM(qcschema_env%driver)) + ! return result + nresult = SIZE(qcschema_env%return_result) + IF (SIZE(qcschema_env%return_result) == 1) THEN + CALL h5awrite_double_scalar(file_id, 'return_result', qcschema_env%return_result(1)) + ELSE + CALL h5awrite_double_simple(file_id, 'return_result', qcschema_env%return_result) + END IF + ! schema name + CALL h5awrite_fixlen_string(file_id, 'schema_name', TRIM(qcschema_env%topology%schema_name)) + ! schema version + CALL h5awrite_integer_scalar(file_id, 'schema_version', qcschema_env%topology%schema_version) + ! success + CALL h5awrite_boolean(file_id, 'success', qcschema_env%success) + + !========================================================================================! + ! *** QCSchema provenance *** + !========================================================================================! + ! create the provenance group + CALL h5gcreate(file_id, 'provenance', group_id) + ! populate provenance + CALL h5awrite_fixlen_string(group_id, 'creator', TRIM(qcschema_env%provenance%creator)) + CALL h5awrite_fixlen_string(group_id, 'routine', TRIM(qcschema_env%provenance%routine)) + CALL h5awrite_fixlen_string(group_id, 'version', TRIM(qcschema_env%provenance%version)) + ! close provenance group + CALL h5gclose(group_id) + + !========================================================================================! + ! *** QCSchema molecule *** + !========================================================================================! + ! create the molecule group + CALL h5gcreate(file_id, 'molecule', group_id) + ! populate molecule + CALL h5awrite_double_simple(group_id, 'geometry', qcschema_env%topology%geometry) + CALL h5awrite_integer_simple(group_id, 'atomic_numbers', qcschema_env%topology%atomic_numbers) + CALL h5awrite_double_simple(group_id, 'masses', qcschema_env%topology%masses) + CALL h5awrite_integer_scalar(group_id, 'molecular_charge', qcschema_env%topology%molecular_charge) + CALL h5awrite_integer_scalar(group_id, 'molecular_multiplicity', qcschema_env%topology%molecular_multiplicity) + CALL h5awrite_string_simple(group_id, 'symbols', qcschema_env%topology%symbols) + + CALL h5awrite_fixlen_string(group_id, 'schema_name', 'qcschema_molecule') + CALL h5awrite_integer_scalar(group_id, 'schema_version', 2) + ! close molecule group + CALL h5gclose(group_id) + + !========================================================================================! + ! *** QCSchema properties *** + !========================================================================================! + ! create the properties group + CALL h5gcreate(file_id, 'properties', group_id) + ! populate properties + CALL h5awrite_integer_scalar(group_id, 'calcinfo_natom', qcschema_env%properties%calcinfo_natom) + CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbasis', qcschema_env%properties%calcinfo_nbasis) + CALL h5awrite_integer_scalar(group_id, 'calcinfo_nmo', qcschema_env%properties%calcinfo_nmo) + CALL h5awrite_integer_scalar(group_id, 'calcinfo_nalpha', qcschema_env%properties%calcinfo_nalpha) + CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbeta', qcschema_env%properties%calcinfo_nbeta) + + ! CALL h5dwrite_double_simple(group_id, 'scf_dipole_moment', & + ! qcschema_env%properties%scf_dipole_moment) + + ! energies, scf, mp2, ... + CALL h5awrite_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy) + CALL h5awrite_double_scalar(group_id, 'scf_total_energy', qcschema_env%properties%scf_total_energy) + CALL h5awrite_double_scalar(group_id, 'nuclear_repulsion_energy', & + qcschema_env%properties%nuclear_repulsion_energy) + + IF (qcschema_env%properties%scf_iterations /= 0) THEN + CALL h5awrite_integer_scalar(group_id, 'scf_iterations', qcschema_env%properties%scf_iterations) + END IF + + IF (qcschema_env%properties%scf_one_electron_energy /= 0.0_dp) THEN + CALL h5awrite_double_scalar(group_id, 'scf_one_electron_energy', & + qcschema_env%properties%scf_one_electron_energy) + END IF + + IF (qcschema_env%properties%scf_two_electron_energy /= 0.0_dp) THEN + CALL h5awrite_double_scalar(group_id, 'scf_two_electron_energy', & + qcschema_env%properties%scf_two_electron_energy) + END IF + + IF (qcschema_env%properties%scf_xc_energy /= 0.0_dp) THEN + CALL h5awrite_double_scalar(group_id, 'scf_xc_energy', & + qcschema_env%properties%scf_xc_energy) + END IF + + IF (qcschema_env%properties%scf_dispersion_correction_energy /= 0.0_dp) THEN + CALL h5awrite_double_scalar(group_id, 'scf_dispersion_correction_energy', & + qcschema_env%properties%scf_dispersion_correction_energy) + END IF + + IF (qcschema_env%properties%mp2) THEN + CALL h5awrite_double_scalar(group_id, 'mp2_correlation_energy', & + qcschema_env%properties%mp2_correlation_energy) + END IF + + ! close properties group + CALL h5gclose(group_id) + + !========================================================================================! + ! *** QCSchema wavefunction *** + !========================================================================================! + ! create the wavefunction group + CALL h5gcreate(file_id, 'wavefunction', group_id) + + CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name)) + + CALL h5dwrite_double_simple(group_id, 'scf_orbitals_a', & + qcschema_env%wavefunction%scf_orbitals_a) + + CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_a', & + qcschema_env%wavefunction%scf_eigenvalues_a) + + CALL h5dwrite_double_simple(group_id, 'scf_occupations_a', & + qcschema_env%wavefunction%scf_occupations_a) + + CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_a', & + qcschema_env%wavefunction%scf_fock_mo_a) + + CALL h5dwrite_double_simple(group_id, 'scf_density_mo_a', & + qcschema_env%wavefunction%scf_density_mo_a) + + CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_aa', & + qcschema_env%wavefunction%scf_eri_mo_aa) + + IF (.NOT. qcschema_env%wavefunction%restricted) THEN + CALL h5dwrite_double_simple(group_id, 'scf_orbitals_b', & + qcschema_env%wavefunction%scf_orbitals_b) + + CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_b', & + qcschema_env%wavefunction%scf_eigenvalues_b) + + CALL h5dwrite_double_simple(group_id, 'scf_occupations_b', & + qcschema_env%wavefunction%scf_occupations_b) + + CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_b', & + qcschema_env%wavefunction%scf_fock_mo_b) + + CALL h5dwrite_double_simple(group_id, 'scf_density_mo_b', & + qcschema_env%wavefunction%scf_density_mo_b) + + CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_bb', & + qcschema_env%wavefunction%scf_eri_mo_bb) + + CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_ba', & + qcschema_env%wavefunction%scf_eri_mo_ab) + + END IF + + ! close wavefunction group + CALL h5gclose(group_id) + + !========================================================================================! + ! *** QCSchema model *** + !========================================================================================! + ! create the model group + CALL h5gcreate(file_id, 'model', group_id) + CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name)) + CALL h5awrite_fixlen_string(group_id, 'method', TRIM(qcschema_env%wavefunction%method)) + ! close model group + CALL h5gclose(group_id) + + ! create the keywords group + CALL h5gcreate(file_id, 'keywords', group_id) + ! close keywords group + CALL h5gclose(group_id) + + CALL h5fclose(file_id) + CALL h5close() +#endif + + END SUBROUTINE qcschema_to_hdf5 + +#ifdef __HDF5 +! ************************************************************************************************** +!> \brief Reads the electron density from a qcschema hdf5 file +!> \param filename the path to the qcschema hdf5 file +!> \param qcschema_env the qcschema environment onto which it writes the density +! ************************************************************************************************** + SUBROUTINE read_pmat_from_hdf5(filename, qcschema_env) + CHARACTER(LEN=default_path_length), INTENT(IN) :: filename + TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env + + INTEGER :: nmo + INTEGER(KIND=hdf5_id) :: file_id, group_id + + ! initialize HDF5 Fortran API + CALL h5open() + + ! open qcschema hdf5 file + CALL h5fopen(TRIM(filename), file_id) + + ! open the wave function group + CALL h5gopen(file_id, 'wavefunction', group_id) + + ! allocate the space for the array containing the density + nmo = qcschema_env%properties%calcinfo_nmo + IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN + ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo)) + END IF + + ! read the alpha density + CALL h5dread_double_simple(group_id, 'scf_density_mo_a', qcschema_env%wavefunction%scf_density_mo_a) + + IF (.NOT. qcschema_env%wavefunction%restricted) THEN + IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN + ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo)) + END IF + ! read the beta density + CALL h5dread_double_simple(group_id, 'scf_density_mo_b', qcschema_env%wavefunction%scf_density_mo_b) + END IF + + ! close everything + CALL h5gclose(group_id) + CALL h5fclose(file_id) + CALL h5close() + + END SUBROUTINE read_pmat_from_hdf5 + +! ************************************************************************************************** +!> \brief Reads the return energy from a qcschema hdf5 file +!> \param filename the path to the qcschema hdf5 file +!> \param qcschema_env the qcschema environment onto which it writes the energy +! ************************************************************************************************** + SUBROUTINE read_return_energy_from_hdf5(filename, qcschema_env) + CHARACTER(LEN=default_path_length), INTENT(IN) :: filename + TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env + + INTEGER(KIND=hdf5_id) :: file_id, group_id + + ! initialize HDF5 Fortran API + CALL h5open() + + ! open qcschema hdf5 file + CALL h5fopen(TRIM(filename), file_id) + + ! open the properties group + CALL h5gopen(file_id, 'properties', group_id) + + ! read the return energy + CALL h5aread_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy) + + ! close everything + CALL h5gclose(group_id) + CALL h5fclose(file_id) + CALL h5close() + + END SUBROUTINE read_return_energy_from_hdf5 + +! ************************************************************************************************** +!> \brief Reads the active space energy from a qcschema file and stores it in active_space_env +!> \param active_space_env ... +!> \param qcschema_env ... +!> \author Stefano Battaglia +! ************************************************************************************************** + SUBROUTINE read_active_energy_from_hdf5(active_space_env, qcschema_env) + TYPE(active_space_type), POINTER :: active_space_env + TYPE(qcschema_type) :: qcschema_env + + CHARACTER(LEN=default_path_length) :: qcschema_filename + + ! File name + qcschema_filename = active_space_env%qcschema_filename + ! read active space energy + CALL read_return_energy_from_hdf5(qcschema_filename, qcschema_env) + + active_space_env%energy_active = qcschema_env%properties%return_energy + active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active + + END SUBROUTINE read_active_energy_from_hdf5 +#endif + +END MODULE qcschema diff --git a/src/qs_active_space_methods.F b/src/qs_active_space_methods.F index 02a0a739a8..e21bfba65e 100644 --- a/src/qs_active_space_methods.F +++ b/src/qs_active_space_methods.F @@ -12,125 +12,140 @@ !> \author JGH ! ************************************************************************************************** MODULE qs_active_space_methods - USE atomic_kind_types, ONLY: atomic_kind_type - USE basis_set_types, ONLY: allocate_sto_basis_set,& - create_gto_from_sto_basis,& - deallocate_sto_basis_set,& - gto_basis_set_type,& - init_orb_basis_set,& - set_sto_basis_set,& - srules,& - sto_basis_set_type - USE cell_types, ONLY: cell_type - USE cp_blacs_env, ONLY: cp_blacs_env_type - USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,& - cp_dbcsr_sm_fm_multiply,& - dbcsr_allocate_matrix_set - USE cp_external_control, ONLY: external_control - USE cp_files, ONLY: close_file,& - open_file - USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale - USE cp_fm_diag, ONLY: cp_fm_syevd - USE cp_fm_struct, ONLY: cp_fm_struct_create,& - cp_fm_struct_release,& - cp_fm_struct_type - USE cp_fm_types, ONLY: & - cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_init_random, cp_fm_release, & - cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm, cp_fm_type - USE cp_log_handling, ONLY: cp_get_default_logger,& - cp_logger_get_default_io_unit,& - cp_logger_type - USE cp_output_handling, ONLY: & - cp_p_file, cp_print_key_finished_output, cp_print_key_should_output, cp_print_key_unit_nr, & - debug_print_level, high_print_level, low_print_level, medium_print_level, & - silent_print_level - USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube - USE dbcsr_api, ONLY: dbcsr_copy,& - dbcsr_csr_create,& - dbcsr_csr_type,& - dbcsr_p_type,& - dbcsr_type - USE input_constants, ONLY: & - casci_canonical, dmft_model, eri_method_full_gpw, eri_method_gpw_ht, eri_operator_coulomb, & - eri_operator_erf, eri_operator_erfc, eri_operator_gaussian, eri_operator_yukawa, hf_model, & - manual_selection, mao_projection, rsdft_model, wannier_projection - USE input_section_types, ONLY: section_vals_get,& - section_vals_get_subs_vals,& - section_vals_type,& - section_vals_val_get - USE kinds, ONLY: default_path_length,& - dp,& - int_8 - USE machine, ONLY: m_walltime - USE mathconstants, ONLY: fourpi - USE memory_utilities, ONLY: reallocate - USE message_passing, ONLY: mp_comm_type,& - mp_para_env_type - USE parallel_gemm_api, ONLY: parallel_gemm - USE particle_list_types, ONLY: particle_list_type - USE particle_types, ONLY: particle_type - USE periodic_table, ONLY: ptable - USE preconditioner_types, ONLY: preconditioner_type - USE pw_env_methods, ONLY: pw_env_create,& - pw_env_rebuild - USE pw_env_types, ONLY: pw_env_get,& - pw_env_release,& - pw_env_type - USE pw_methods, ONLY: pw_integrate_function,& - pw_transfer - USE pw_poisson_methods, ONLY: pw_poisson_rebuild,& - pw_poisson_solve - USE pw_poisson_types, ONLY: ANALYTIC0D,& - PERIODIC3D,& - greens_fn_type,& - pw_poisson_analytic,& - pw_poisson_periodic,& - pw_poisson_type - USE pw_pool_types, ONLY: pw_pool_create_pw,& - pw_pool_give_back_pw,& - pw_pool_type - USE pw_types, ONLY: COMPLEXDATA1D,& - REALDATA3D,& - REALSPACE,& - RECIPROCALSPACE,& - pw_release,& - pw_type - USE qs_active_space_types, ONLY: active_space_type,& - create_active_space_type,& - csr_idx_from_combined,& - csr_idx_to_combined,& - eri_type,& - eri_type_eri_element_func,& - get_irange_csr - USE qs_collocate_density, ONLY: calculate_wavefunction - USE qs_density_matrices, ONLY: calculate_density_matrix - USE qs_energy_types, ONLY: qs_energy_type - USE qs_environment_types, ONLY: get_qs_env,& - qs_environment_type,& - set_qs_env - USE qs_kind_types, ONLY: qs_kind_type - USE qs_ks_methods, ONLY: qs_ks_update_qs_env - USE qs_ks_types, ONLY: qs_ks_did_change,& - qs_ks_env_type - USE qs_mo_io, ONLY: write_mo_set_to_output_unit - USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues - USE qs_mo_types, ONLY: allocate_mo_set,& - get_mo_set,& - init_mo_set,& - mo_set_type - USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type - USE qs_ot_eigensolver, ONLY: ot_eigensolver - USE qs_rho_methods, ONLY: qs_rho_update_rho - USE qs_rho_types, ONLY: qs_rho_get,& - qs_rho_type - USE qs_subsys_types, ONLY: qs_subsys_get,& - qs_subsys_type - USE scf_control_types, ONLY: scf_control_type - USE task_list_methods, ONLY: generate_qs_task_list - USE task_list_types, ONLY: allocate_task_list,& - deallocate_task_list,& - task_list_type + USE atomic_kind_types, ONLY: atomic_kind_type + USE basis_set_types, ONLY: allocate_sto_basis_set, & + create_gto_from_sto_basis, & + deallocate_sto_basis_set, & + gto_basis_set_type, & + init_orb_basis_set, & + set_sto_basis_set, & + srules, & + sto_basis_set_type + USE cell_types, ONLY: cell_type + USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t, & + cp_dbcsr_sm_fm_multiply, & + dbcsr_allocate_matrix_set + USE cp_files, ONLY: close_file, & + file_exists, & + open_file + USE cp_fm_struct, ONLY: cp_fm_struct_create, & + cp_fm_struct_release, & + cp_fm_struct_type + USE cp_fm_types, ONLY: & + cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_init_random, cp_fm_release, & + cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm, cp_fm_type + USE cp_log_handling, ONLY: cp_get_default_logger, & + cp_logger_get_default_io_unit, & + cp_logger_type + USE cp_output_handling, ONLY: & + cp_p_file, cp_print_key_finished_output, cp_print_key_should_output, cp_print_key_unit_nr, & + debug_print_level, high_print_level, low_print_level, medium_print_level, & + silent_print_level + USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube + USE dbcsr_api, ONLY: dbcsr_copy, & + dbcsr_csr_create, & + dbcsr_csr_type, & + dbcsr_p_type, & + dbcsr_type + USE input_constants, ONLY: & + casci_canonical, dmft_model, eri_method_full_gpw, eri_method_gpw_ht, eri_operator_coulomb, & + eri_operator_erf, eri_operator_erfc, eri_operator_gaussian, eri_operator_yukawa, hf_model, & + manual_selection, mao_projection, no_solver, qiskit_solver, rsdft_model, wannier_projection + USE input_section_types, ONLY: section_vals_get, & + section_vals_get_subs_vals, & + section_vals_type, & + section_vals_val_get + USE iso_c_binding, ONLY: c_null_char + USE kinds, ONLY: default_path_length, & + default_string_length, & + dp, & + int_8 + USE machine, ONLY: m_walltime + USE mathconstants, ONLY: fourpi + USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_comm_type, & + mp_para_env_type + USE parallel_gemm_api, ONLY: parallel_gemm + USE particle_list_types, ONLY: particle_list_type + USE particle_types, ONLY: particle_type + USE periodic_table, ONLY: ptable + USE preconditioner_types, ONLY: preconditioner_type + USE pw_env_methods, ONLY: pw_env_create, & + pw_env_rebuild + USE pw_env_types, ONLY: pw_env_get, & + pw_env_release, & + pw_env_type + USE pw_methods, ONLY: pw_integrate_function, & + pw_transfer + USE pw_poisson_methods, ONLY: pw_poisson_rebuild, & + pw_poisson_solve + USE pw_poisson_types, ONLY: ANALYTIC0D, & + PERIODIC3D, & + greens_fn_type, & + pw_poisson_analytic, & + pw_poisson_periodic, & + pw_poisson_type + USE pw_pool_types, ONLY: pw_pool_create_pw, & + pw_pool_give_back_pw, & + pw_pool_type + USE pw_types, ONLY: COMPLEXDATA1D, & + REALDATA3D, & + REALSPACE, & + RECIPROCALSPACE, & + pw_release, & + pw_type + USE qcschema, ONLY: qcschema_env_create, & + qcschema_env_release, & + qcschema_to_hdf5, & + qcschema_type + USE qs_active_space_types, ONLY: active_space_type, & + create_active_space_type, & + csr_idx_from_combined, & + csr_idx_to_combined, & + eri_type, & + eri_type_eri_element_func, & + get_irange_csr + USE qs_active_space_utils, ONLY: eri_to_array, & + subspace_matrix_to_array + USE qs_collocate_density, ONLY: calculate_wavefunction + USE qs_density_matrices, ONLY: calculate_density_matrix + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env, & + qs_environment_type, & + set_qs_env + USE qs_kind_types, ONLY: qs_kind_type + USE qs_ks_methods, ONLY: qs_ks_update_qs_env + USE qs_ks_types, ONLY: qs_ks_did_change, & + qs_ks_env_type + USE qs_mo_io, ONLY: write_mo_set_to_output_unit + USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues + USE qs_mo_types, ONLY: allocate_mo_set, & + get_mo_set, & + init_mo_set, & + mo_set_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_ot_eigensolver, ONLY: ot_eigensolver + USE qs_rho_methods, ONLY: qs_rho_update_rho + USE qs_rho_types, ONLY: qs_rho_get, & + qs_rho_type + USE qs_subsys_types, ONLY: qs_subsys_get, & + qs_subsys_type + USE scf_control_types, ONLY: scf_control_type +#ifndef __NO_SOCKETS + USE sockets_interface, ONLY: accept_socket, & + close_socket, & + listen_socket, & + open_bind_socket, & + readbuffer, & + remove_socket_file, & + writebuffer +#endif + USE task_list_methods, ONLY: generate_qs_task_list + USE task_list_types, ONLY: allocate_task_list, & + deallocate_task_list, & + task_list_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -154,8 +169,6 @@ MODULE qs_active_space_methods PROCEDURE :: func => eri_fcidump_checksum_func END TYPE eri_fcidump_checksum -! ************************************************************************************************** - CONTAINS ! ************************************************************************************************** @@ -173,41 +186,36 @@ CONTAINS ! ************************************************************************************************** !> \brief Main method for determining the active space Hamiltonian -!> \param input ... -!> \param logger ... !> \param qs_env ... ! ************************************************************************************************** - SUBROUTINE active_space_main(input, logger, qs_env) - TYPE(section_vals_type), POINTER :: input - TYPE(cp_logger_type), POINTER :: logger + SUBROUTINE active_space_main(qs_env) TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'active_space_main' CHARACTER(len=10) :: cshell, lnam(5) - CHARACTER(LEN=default_path_length) :: p_act_filename - INTEGER :: eri_method, eri_operator, eri_print, handle, i, iatom, ishell, isp, ispin, iw, j, & - jm, l, m, max_orb_ind, mselect, n1, n2, nao, natom, ncol, nel, nelec_active, & - nelec_inactive, nelec_inactive_alpha, nelec_inactive_beta, nelec_total, nmo, nn1, nn2, & - nrow_global, nspins, nvirt - INTEGER, ALLOCATABLE, DIMENSION(:, :) :: active_orbitals - INTEGER, DIMENSION(2) :: nactive_orb, nepol, ninactive_orb + CHARACTER(len=default_path_length) :: qcschema_filename + INTEGER :: as_solver, eri_method, eri_operator, eri_print, handle, i, iatom, ishell, isp, & + ispin, iw, j, jm, m, max_orb_ind, mselect, n1, n2, nao, natom, nel, nelec_active, & + nelec_inactive, nelec_total, nmo, nmo_active, nmo_available, nmo_inactive, & + nmo_inactive_remaining, nmo_occ, nmo_virtual, nn1, nn2, nrow_global, nspins INTEGER, DIMENSION(5) :: nshell INTEGER, DIMENSION(:), POINTER :: invals LOGICAL :: do_kpoints, ex_operator, ex_perd, & - explicit, isolated, read_p_act, & - stop_after_print, store_wfn + explicit, isolated, stop_after_print, & + store_wfn REAL(KIND=dp) :: eri_eps_grid, eri_eps_int, eri_gpw_cutoff, eri_op_param, eri_rcut, & eri_rel_cutoff, fel, focc, maxocc, nze_percentage REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues - REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virt + REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virtual TYPE(active_space_type), POINTER :: active_space_env TYPE(cell_type), POINTER :: cell TYPE(cp_blacs_env_type), POINTER :: context TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp - TYPE(cp_fm_type) :: fm_dummy, mo_virt - TYPE(cp_fm_type), POINTER :: fm_ref, fm_target, fm_target_active, & - fm_target_inactive, fmat, mo_coeff + TYPE(cp_fm_type) :: fm_dummy, mo_virtual + TYPE(cp_fm_type), POINTER :: fm_target_active, fm_target_inactive, & + fmat, mo_coeff, mo_ref, mo_target + TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_csr_type), POINTER :: eri_mat TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, rho_ao, s_matrix TYPE(dbcsr_type), POINTER :: denmat @@ -216,18 +224,22 @@ CONTAINS TYPE(mo_set_type), POINTER :: mo_set, mo_set_active, mo_set_inactive TYPE(mp_para_env_type), POINTER :: para_env TYPE(preconditioner_type), POINTER :: local_preconditioner + TYPE(qcschema_type) :: qcschema_env + TYPE(qs_energy_type), POINTER :: energy TYPE(qs_rho_type), POINTER :: rho TYPE(scf_control_type), POINTER :: scf_control - TYPE(section_vals_type), POINTER :: as_input, loc_print, loc_section, & + TYPE(section_vals_type), POINTER :: as_input, input, loc_print, loc_section, & print_orb !--------------------------------------------------------------------------------------------! - as_input => section_vals_get_subs_vals(input, "DFT%PRINT%ACTIVE_SPACE") + CALL get_qs_env(qs_env, input=input) + as_input => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE") CALL section_vals_get(as_input, explicit=explicit) IF (.NOT. explicit) RETURN CALL timeset(routineN, handle) + logger => cp_get_default_logger() iw = cp_logger_get_default_io_unit(logger) IF (iw > 0) THEN @@ -254,16 +266,27 @@ CONTAINS ! input options + ! figure out what needs to be printed/stored + IF (BTEST(cp_print_key_should_output(logger%iter_info, as_input, "FCIDUMP"), cp_p_file)) THEN + active_space_env%fcidump = .TRUE. + END IF + + CALL section_vals_val_get(as_input, "QCSCHEMA", c_val=qcschema_filename, explicit=explicit) + IF (explicit) THEN + active_space_env%qcschema = .TRUE. + active_space_env%qcschema_filename = qcschema_filename + END IF + ! model CALL section_vals_val_get(as_input, "MODEL", i_val=active_space_env%model) IF (iw > 0) THEN SELECT CASE (active_space_env%model) CASE (hf_model) - WRITE (iw, '(T4,A)') "Hartree-Fock model for interaction Hamiltonian" + WRITE (iw, '(T3,A)') "Hartree-Fock model for interaction Hamiltonian" CASE (rsdft_model) - WRITE (iw, '(T4,A)') "Range-separated DFT model for interaction Hamiltonian" + WRITE (iw, '(T3,A)') "Range-separated DFT model for interaction Hamiltonian" CASE (dmft_model) - WRITE (iw, '(T4,A)') "DMFT model for interaction Hamiltonian" + WRITE (iw, '(T3,A)') "DMFT model for interaction Hamiltonian" CASE DEFAULT CPABORT("Unknown Model") END SELECT @@ -274,132 +297,46 @@ CONTAINS active_space_env%molecule = isolated IF (iw > 0) THEN IF (active_space_env%molecule) THEN - WRITE (iw, '(T4,A)') "System is treated without periodicity" + WRITE (iw, '(T3,A)') "System is treated without periodicity" END IF END IF CALL section_vals_val_get(as_input, "ACTIVE_ELECTRONS", i_val=nelec_active) + ! actually nelec_spin tells me the number of electrons per spin channel from qs_env + ! CALL get_qs_env(qs_env, nelectron_total=nelec_total, nelectron_spin=nelec_spin) + CALL get_qs_env(qs_env, nelectron_total=nelec_total) + IF (nelec_active <= 0) CPABORT("Specify a positive number of active electrons.") - CALL get_qs_env(qs_env, nelectron_total=nelec_total, nelectron_spin=nepol) + IF (nelec_active > nelec_total) CPABORT("More active electrons than total electrons.") + + nelec_inactive = nelec_total - nelec_active + IF (MOD(nelec_inactive, 2) /= 0) THEN + CPABORT("The remaining number of inactive electrons has to be even.") + END IF + CALL get_qs_env(qs_env, dft_control=dft_control) nspins = dft_control%nspins - CALL section_vals_val_get(as_input, "INACTIVE_ELECTRONS", explicit=explicit) - CALL section_vals_val_get(as_input, "INACTIVE_ELECTRONS", i_vals=invals) - IF (nspins > 1) THEN - IF (.NOT. explicit) THEN - CALL cp_abort(__LOCATION__, "Number of Inactive Electrons has to be specified"// & - " in spin polarised case.") - END IF - nelec_inactive_alpha = invals(1) - nelec_inactive_beta = invals(2) - nelec_inactive = nelec_inactive_alpha + nelec_inactive_beta - ELSE - IF (explicit) THEN - nelec_inactive = invals(1) - ELSE - nelec_inactive = nelec_total - nelec_active - END IF - nelec_inactive_alpha = nelec_inactive - nelec_inactive_beta = 0 - END IF - IF (iw > 0) THEN - IF (nspins < 2) THEN - WRITE (iw, '(T4,A,T69,I10)') "Total number of electrons", nelec_total - WRITE (iw, '(T4,A,T69,I10)') "Number of active electrons", nelec_active - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive electrons", nelec_inactive - ELSE - WRITE (iw, '(T4,A,T69,I10)') "Total number of electrons", nelec_total - WRITE (iw, '(T4,A,T69,I10)') "Number of active electrons", nelec_active - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive electrons", nelec_inactive - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive electrons (alpha)", nelec_inactive_alpha - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive electrons (beta)", nelec_inactive_beta - END IF - END IF - - ! Read active density from file (generated by external program) - CALL section_vals_val_get(as_input, "READ_P_ACTIVE", l_val=read_p_act) - active_space_env%read_p_act = read_p_act - IF (read_p_act) THEN - CALL section_vals_val_get(as_input, "P_ACTIVE_FILE_NAME", c_val=p_act_filename) - active_space_env%p_act_filename = p_act_filename - END IF - - CPASSERT(nelec_inactive >= 0) - CPASSERT(nelec_inactive_alpha >= 0) - CPASSERT(nelec_inactive_beta >= 0) - CPASSERT(nelec_total == nelec_inactive + nelec_active) - IF (nspins > 1) THEN - CPASSERT(nepol(1) >= nelec_inactive_alpha) - CPASSERT(nepol(2) >= nelec_inactive_beta) + WRITE (iw, '(T3,A,T70,I10)') "Total number of electrons", nelec_total + WRITE (iw, '(T3,A,T70,I10)') "Number of inactive electrons", nelec_inactive + WRITE (iw, '(T3,A,T70,I10)') "Number of active electrons", nelec_active END IF active_space_env%nelec_active = nelec_active active_space_env%nelec_inactive = nelec_inactive - IF (nspins == 1) THEN - active_space_env%nelec_total(1) = nelec_total - active_space_env%nelec_total(2) = 0 - ELSE - active_space_env%nelec_total(1) = nepol(1) - active_space_env%nelec_total(2) = nepol(2) - END IF - active_space_env%nelec_inactive_spinwise(1) = nelec_inactive_alpha - active_space_env%nelec_inactive_spinwise(2) = nelec_inactive_beta - active_space_env%multiplicity = dft_control%multiplicity + active_space_env%nelec_total = nelec_total active_space_env%nspins = nspins + active_space_env%multiplicity = dft_control%multiplicity ! define the active/inactive space orbitals - CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit) - CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", i_vals=invals) + CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit, i_val=nmo_active) IF (.NOT. explicit) THEN CALL cp_abort(__LOCATION__, "Number of Active Orbitals has to be specified.") END IF - nactive_orb = 0 - CPASSERT(SIZE(invals) >= nspins) - DO ispin = 1, nspins - nactive_orb(ispin) = invals(ispin) - END DO - ninactive_orb = 0 - CALL section_vals_val_get(as_input, "INACTIVE_ORBITALS", explicit=explicit) - CALL section_vals_val_get(as_input, "INACTIVE_ORBITALS", i_vals=invals) - IF (.NOT. explicit) THEN - IF (nspins == 1) THEN - ninactive_orb(1) = active_space_env%nelec_inactive_spinwise(1)/2 - ELSE - DO ispin = 1, nspins - ninactive_orb(ispin) = active_space_env%nelec_inactive_spinwise(ispin) - END DO - END IF - ELSE - CPASSERT(SIZE(invals) >= nspins) - DO ispin = 1, nspins - ninactive_orb(ispin) = invals(ispin) - END DO - END IF - IF (nspins == 1) THEN - CPASSERT(MOD(nelec_inactive, 2) == 0) - CPASSERT(ninactive_orb(1) >= nelec_inactive/2) - ELSE - CPASSERT(ninactive_orb(1) >= nelec_inactive_alpha) - CPASSERT(ninactive_orb(2) >= nelec_inactive_beta) - END IF - - ! CALL get_qs_env(qs_env, mos=mos) - ! - ! CALL get_mo_set(mos(1), mo_coeff=fm_ref, nao=nao) - ! maxocc = 2.0_dp - ! IF (nspins > 1) maxocc = 1.0_dp - ! ALLOCATE (active_space_env%mos_active(nspins)) - ! ALLOCATE (active_space_env%mos_inactive(nspins)) - ! DO ispin = 1, nspins - ! nmo = nao - ! CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo, 0, 0.0_dp, maxocc, 0.0_dp) - ! CALL init_mo_set(active_space_env%mos_active(ispin), fm_ref=fm_ref, name="Active Space MO") - ! nmo = nao - ! nel = active_space_env%nelec_inactive_spinwise(ispin) - ! CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo, nel, REAL(nel, KIND=dp), maxocc, 0.0_dp) - ! CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_ref=fm_ref, name="Inactive Space MO") - ! END DO + active_space_env%nmo_active = nmo_active + ! this is safe because nelec_inactive is always even + nmo_inactive = nelec_inactive/2 + active_space_env%nmo_inactive = nmo_inactive CALL section_vals_val_get(as_input, "ORBITAL_SELECTION", i_val=mselect) IF (iw > 0) THEN @@ -407,29 +344,21 @@ CONTAINS CASE DEFAULT CPABORT("Unknown orbital selection method") CASE (casci_canonical) - WRITE (iw, '(/,T4,A)') & + WRITE (iw, '(/,T3,A)') & "Active space orbitals selected using energy ordered canonical orbitals" CASE (wannier_projection) - WRITE (iw, '(/,T4,A)') & + WRITE (iw, '(/,T3,A)') & "Active space orbitals selected using projected Wannier orbitals" CASE (mao_projection) - WRITE (iw, '(/,T4,A)') & + WRITE (iw, '(/,T3,A)') & "Active space orbitals selected using modified atomic orbitals (MAO)" CASE (manual_selection) - WRITE (iw, '(/,T4,A)') & + WRITE (iw, '(/,T3,A)') & "Active space orbitals selected manually" END SELECT - IF (nspins < 2) THEN - WRITE (iw, '(T4,A,T69,I10)') "Number of active orbitals", nactive_orb(1) - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive orbitals", ninactive_orb(1) - ELSE - WRITE (iw, '(T4,A,T69,I10)') "Number of active orbitals", SUM(nactive_orb) - WRITE (iw, '(T4,A,T69,I10)') "Number of active orbitals (alpha)", nactive_orb(1) - WRITE (iw, '(T4,A,T69,I10)') "Number of active orbitals (beta)", nactive_orb(2) - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive orbitals", SUM(ninactive_orb) - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive orbitals (alpha)", ninactive_orb(1) - WRITE (iw, '(T4,A,T69,I10)') "Number of inactive orbitals (beta)", ninactive_orb(2) - END IF + + WRITE (iw, '(T3,A,T70,I10)') "Number of inactive orbitals", nmo_inactive + WRITE (iw, '(T3,A,T70,I10)') "Number of active orbitals", nmo_active END IF ! get projection spaces @@ -438,7 +367,7 @@ CONTAINS CALL get_qs_env(qs_env, natom=natom) IF (iatom <= 0 .OR. iatom > natom) THEN IF (iw > 0) THEN - WRITE (iw, '(/,T4,A,I3)') "ERROR: SUBSPACE_ATOM number is not valid", iatom + WRITE (iw, '(/,T3,A,I3)') "ERROR: SUBSPACE_ATOM number is not valid", iatom END IF CPABORT("Select a valid SUBSPACE_ATOM") END IF @@ -467,153 +396,173 @@ CONTAINS CASE (casci_canonical) CALL get_qs_env(qs_env, mos=mos) - nmo = 0 - DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), nmo=m) - ! IF (m < ninactive_orb(ispin) + nactive_orb(ispin)) THEN - ! CPABORT("Not enough canonical orbitals available.") - ! END IF - nmo = MAX(m, nmo) - END DO + ! total number of occupied orbitals, i.e. inactive plus active MOs + nmo_occ = nmo_inactive + nmo_active - ! set inactive orbital indices - ALLOCATE (active_space_env%inactive_orbitals(MAXVAL(ninactive_orb), nspins)) + ! set inactive orbital indices, these are trivially 1...nmo_inactive + ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins)) DO ispin = 1, nspins - i = 1 - DO m = 1, ninactive_orb(ispin) - active_space_env%inactive_orbitals(i, ispin) = m - i = i + 1 + DO i = 1, nmo_inactive + active_space_env%inactive_orbitals(i, ispin) = i END DO END DO - max_orb_ind = 0 - ! set active orbital indices - ALLOCATE (active_space_env%active_orbitals(MAXVAL(nactive_orb), nspins)) + ! set active orbital indices, these are shifted by nmo_inactive + ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins)) DO ispin = 1, nspins - i = 1 - DO m = ninactive_orb(ispin) + 1, ninactive_orb(ispin) + nactive_orb(ispin) - active_space_env%active_orbitals(i, ispin) = m - max_orb_ind = MAX(m, max_orb_ind) - i = i + 1 + DO i = 1, nmo_active + active_space_env%active_orbitals(i, ispin) = nmo_inactive + i END DO END DO - CALL get_mo_set(mos(1), mo_coeff=fm_ref, nao=nao, nmo=nmo) - CALL cp_fm_get_info(fm_ref, context=context, para_env=para_env, & - nrow_global=nrow_global) - nvirt = max_orb_ind - nmo - nvirt = MAX(nvirt, 0) - nmo = max_orb_ind - + ! allocate and initialize inactive and active mo coefficients. + ! These are stored in a data structure for the full occupied space: + ! for inactive mos, the active subset is set to zero, vice versa for the active mos + ! TODO: allocate data structures only for the eaxct number MOs maxocc = 2.0_dp IF (nspins > 1) maxocc = 1.0_dp ALLOCATE (active_space_env%mos_active(nspins)) ALLOCATE (active_space_env%mos_inactive(nspins)) DO ispin = 1, nspins - CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo, 0, 0.0_dp, maxocc, 0.0_dp) + CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao) + CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global) + ! the right number of active electrons per spin channel is initialized further down + CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo_occ, 0, 0.0_dp, maxocc, 0.0_dp) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & - nrow_global=nrow_global, ncol_global=nmo) + nrow_global=nrow_global, ncol_global=nmo_occ) CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO") CALL cp_fm_struct_release(fm_struct_tmp) - nel = active_space_env%nelec_inactive_spinwise(ispin) - CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo, nel, REAL(nel, KIND=dp), maxocc, 0.0_dp) + IF (nspins == 2) THEN + nel = nelec_inactive/2 + ELSE + nel = nelec_inactive + END IF + CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo_occ, nel, & + REAL(nel, KIND=dp), maxocc, 0.0_dp) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & - nrow_global=nrow_global, ncol_global=nmo) + nrow_global=nrow_global, ncol_global=nmo_occ) CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO") CALL cp_fm_struct_release(fm_struct_tmp) END DO ! create canonical orbitals IF (dft_control%restricted) THEN - CPABORT(" Unclear how we define MOs in the restricted case ... stopping") + CPABORT("Unclear how we define MOs in the restricted case ... stopping") ELSE IF (dft_control%do_admm) THEN CPABORT("ADMM currently not possible for canonical orbital options") END IF - ALLOCATE (eigenvalues(nmo, nspins)) - eigenvalues = 100000._dp + ALLOCATE (eigenvalues(nmo_occ, nspins)) + eigenvalues = 0.0_dp CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control) + ! calculate virtual MOs and copy inactive and active orbitals + IF (iw > 0) THEN + WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..." + END IF DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=fm_ref, nmo=nmo) + ! nmo_available is the number of MOs available from the SCF calculation: + ! this is at least the number of occupied orbitals in the SCF, plus + ! any number of added MOs (virtuals) requested in the SCF section + CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available) - NULLIFY (evals_virt) - ALLOCATE (evals_virt(nvirt)) + ! calculate how many extra MOs we still have to compute + nmo_virtual = nmo_occ - nmo_available + nmo_virtual = MAX(nmo_virtual, 0) - CALL cp_fm_get_info(fm_ref, context=context, para_env=para_env, & + NULLIFY (evals_virtual) + ALLOCATE (evals_virtual(nmo_virtual)) + + CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, & nrow_global=nrow_global) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & - nrow_global=nrow_global, ncol_global=nvirt) - CALL cp_fm_create(mo_virt, fm_struct_tmp, name="virtual") + nrow_global=nrow_global, ncol_global=nmo_virtual) + CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual") CALL cp_fm_struct_release(fm_struct_tmp) - CALL cp_fm_init_random(mo_virt, nvirt) + CALL cp_fm_init_random(mo_virtual, nmo_virtual) NULLIFY (local_preconditioner) + ! compute missing virtual MOs CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, & - matrix_c_fm=mo_virt, matrix_orthogonal_space_fm=fm_ref, & + matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, & eps_gradient=scf_control%eps_lumos, & preconditioner=local_preconditioner, & iter_max=scf_control%max_iter_lumos, & - size_ortho_space=nmo) + size_ortho_space=nmo_available) - CALL calculate_subspace_eigenvalues(mo_virt, ks_matrix(ispin)%matrix, & - evals_virt) - - ! We need to send the copy of MOs to preserve the sign - CALL cp_fm_create(fm_dummy, fm_ref%matrix_struct) - CALL cp_fm_to_fm(fm_ref, fm_dummy) + ! get the eigenvalues + CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, evals_virtual) + ! TODO: double check that we really need this. + ! we need to send the copy of MOs to preserve the sign + CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct) + CALL cp_fm_to_fm(mo_ref, fm_dummy) CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, & evals_arg=eigenvalues(:, ispin), do_rotation=.TRUE.) + ! copy inactive orbitals mo_set => active_space_env%mos_inactive(ispin) - CALL get_mo_set(mo_set, mo_coeff=fm_target) + CALL get_mo_set(mo_set, mo_coeff=mo_target) DO i = 1, SIZE(active_space_env%inactive_orbitals, 1) m = active_space_env%inactive_orbitals(i, ispin) - CALL cp_fm_to_fm(fm_ref, fm_target, 1, m, m) + CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m) mo_set%eigenvalues(m) = eigenvalues(m, ispin) - mo_set%occupation_numbers(m) = 1.0 - END DO - ! copy active orbitals - mo_set => active_space_env%mos_active(ispin) - CALL get_mo_set(mo_set, mo_coeff=fm_target) - ncol = nactive_orb(ispin) - DO i = 1, SIZE(active_space_env%active_orbitals, 1) - m = active_space_env%active_orbitals(i, ispin) - IF (i > nmo) THEN - CALL cp_fm_to_fm(mo_virt, fm_target, 1, m - nmo, m) - mo_set%eigenvalues(m) = evals_virt(m - nmo) - mo_set%occupation_numbers(m) = 0.0 - ELSE - CALL cp_fm_to_fm(fm_ref, fm_target, 1, m, m) - mo_set%eigenvalues(m) = eigenvalues(m, ispin) + IF (nspins > 1) THEN mo_set%occupation_numbers(m) = 1.0 + ELSE + mo_set%occupation_numbers(m) = 2.0 END IF END DO - ! Release - DEALLOCATE (evals_virt) - CALL cp_fm_release(fm_dummy) - CALL cp_fm_release(mo_virt) - END DO - CALL get_qs_env(qs_env, mos=mos) + ! copy active orbitals + mo_set => active_space_env%mos_active(ispin) + CALL get_mo_set(mo_set, mo_coeff=mo_target) + ! for mult > 1, put the polarized electrons in the alpha channel + IF (nspins == 2) THEN + IF (ispin == 1) THEN + nel = (nelec_active + active_space_env%multiplicity - 1)/2 + ELSE + nel = (nelec_active - active_space_env%multiplicity + 1)/2 + END IF + ELSE + nel = nelec_active + END IF + mo_set%nelectron = nel + mo_set%n_el_f = REAL(nel, KIND=dp) + DO i = 1, nmo_active + m = active_space_env%active_orbitals(i, ispin) + IF (m > nmo_available) THEN + CALL cp_fm_to_fm(mo_virtual, mo_target, 1, m - nmo_available, m) + eigenvalues(m, ispin) = evals_virtual(m - nmo_available) + mo_set%occupation_numbers(m) = 0.0 + ELSE + CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m) + mo_set%occupation_numbers(m) = mos(ispin)%occupation_numbers(m) + END IF + mo_set%eigenvalues(m) = eigenvalues(m, ispin) + END DO + ! Release + DEALLOCATE (evals_virtual) + CALL cp_fm_release(fm_dummy) + CALL cp_fm_release(mo_virtual) + END DO IF (iw > 0) THEN DO ispin = 1, nspins - WRITE (iw, '(/,T4,A,I3,T65,A)') "Canonical Orbital Selection for spin", ispin, & + WRITE (iw, '(/,T3,A,I3,T66,A)') "Canonical Orbital Selection for spin", ispin, & "[atomic units]" - DO i = 1, ninactive_orb(ispin), 4 - jm = MIN(3, ninactive_orb(ispin) - i) + DO i = 1, nmo_inactive, 4 + jm = MIN(3, nmo_inactive - i) WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [I]", j=0, jm) END DO - DO i = ninactive_orb(ispin) + 1, ninactive_orb(ispin) + nactive_orb(ispin), 4 - jm = MIN(3, ninactive_orb(ispin) + nactive_orb(ispin) - i) + DO i = nmo_inactive + 1, nmo_inactive + nmo_active, 4 + jm = MIN(3, nmo_inactive + nmo_active - i) WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [A]", j=0, jm) END DO - WRITE (iw, '(/,T4,A,I3)') "Active Orbital Indices for spin", ispin + WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4 jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i) WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm) @@ -626,120 +575,107 @@ CONTAINS CASE (manual_selection) ! create canonical orbitals IF (dft_control%restricted) THEN - CPABORT(" Unclear how we define MOs in the restricted case ... stopping") + CPABORT("Unclear how we define MOs in the restricted case ... stopping") ELSE IF (dft_control%do_admm) THEN - CPABORT("ADMM currently not possible for canonical orbital options") + CPABORT("ADMM currently not possible for manual orbitals selection") END IF - IF (nspins > 1 .AND. .NOT. nactive_orb(1) == nactive_orb(2)) THEN - CALL cp_abort(__LOCATION__, "The manual orbital selection "// & - "requires identical orbital indices for each spin!") - END IF - - CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", explicit=explicit) - CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", i_vals=invals) + CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", explicit=explicit, i_vals=invals) IF (.NOT. explicit) THEN - CALL cp_abort(__LOCATION__, "When using the manual "// & - "orbital selection you MUST explicitly "// & + CALL cp_abort(__LOCATION__, "Manual orbital"// & + "selection requires to explicitly "// & "set the active orbital indices via "// & - "ACTIVE_ORBITAL_INDICES!") + "ACTIVE_ORBITAL_INDICES") END IF - ! TODO: assert size of invals - ALLOCATE (active_orbitals(MAXVAL(nactive_orb), nspins)) - max_orb_ind = 0 - ispin = 1 - DO i = 1, SIZE(invals) - IF (i > nactive_orb(1)) THEN - ispin = 2 - END IF - active_orbitals(i - nactive_orb(1)*(ispin - 1), ispin) = invals(i) - max_orb_ind = MAX(invals(i), max_orb_ind) + IF (nspins == 1) THEN + CPASSERT(SIZE(invals) == nmo_active) + ELSE + CPASSERT(SIZE(invals) == 2*nmo_active) + END IF + ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins)) + ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins)) + + DO ispin = 1, nspins + DO i = 1, nmo_active + active_space_env%active_orbitals(i, ispin) = invals(i + (ispin - 1)*nmo_active) + END DO END DO - IF (nspins == 2) THEN - IF (i < 2*SIZE(active_orbitals, 1)) THEN - ! reuse first spin active orbital indices for second spin - DO i = 1, SIZE(active_orbitals, 1) - active_orbitals(i, 2) = active_orbitals(i, 1) - END DO - END IF - END IF - CALL get_qs_env(qs_env, mos=mos) + ! include MOs up to the largest index in the list + max_orb_ind = MAXVAL(invals) maxocc = 2.0_dp IF (nspins > 1) maxocc = 1.0_dp ALLOCATE (active_space_env%mos_active(nspins)) ALLOCATE (active_space_env%mos_inactive(nspins)) DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=fm_ref, nao=nao, nmo=nmo) - CALL cp_fm_get_info(fm_ref, context=context, para_env=para_env, & - nrow_global=nrow_global) + ! init active orbitals + CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao) + CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global) CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, max_orb_ind, 0, 0.0_dp, maxocc, 0.0_dp) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & nrow_global=nrow_global, ncol_global=max_orb_ind) CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO") CALL cp_fm_struct_release(fm_struct_tmp) - nel = active_space_env%nelec_inactive_spinwise(ispin) + + ! init inactive orbitals + IF (nspins == 2) THEN + nel = nelec_inactive/2 + ELSE + nel = nelec_inactive + END IF CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, max_orb_ind, nel, REAL(nel, KIND=dp), maxocc, 0.0_dp) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & nrow_global=nrow_global, ncol_global=max_orb_ind) CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO") + ! small hack: set the correct inactive occupations down below + active_space_env%mos_inactive(ispin)%occupation_numbers = 0.0_dp CALL cp_fm_struct_release(fm_struct_tmp) END DO - nmo = 0 - DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), nmo=m) - ! IF (m < max_orb_ind) THEN - ! CPABORT("Not enough canonical orbitals available.") - ! END IF - nmo = MAX(m, nmo) - END DO - ALLOCATE (eigenvalues(nmo, nspins)) - eigenvalues = 100000._dp + ALLOCATE (eigenvalues(max_orb_ind, nspins)) + eigenvalues = 0.0_dp CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control) - ALLOCATE (active_space_env%inactive_orbitals(MAXVAL(ninactive_orb), nspins)) - ALLOCATE (active_space_env%active_orbitals(MAXVAL(nactive_orb), nspins)) - focc = 2 - IF (nspins == 2) THEN - focc = 1 + ! calculate virtual MOs and copy inactive and active orbitals + IF (iw > 0) THEN + WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..." END IF DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=fm_ref, nmo=nmo) - nvirt = max_orb_ind - nmo - nvirt = MAX(nvirt, 0) + CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available) + nmo_virtual = max_orb_ind - nmo_available + nmo_virtual = MAX(nmo_virtual, 0) - NULLIFY (evals_virt) - ALLOCATE (evals_virt(nvirt)) + NULLIFY (evals_virtual) + ALLOCATE (evals_virtual(nmo_virtual)) - CALL cp_fm_get_info(fm_ref, context=context, para_env=para_env, & + CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, & nrow_global=nrow_global) CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & - nrow_global=nrow_global, ncol_global=nvirt) - CALL cp_fm_create(mo_virt, fm_struct_tmp, name="virtual") + nrow_global=nrow_global, ncol_global=nmo_virtual) + CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual") CALL cp_fm_struct_release(fm_struct_tmp) - CALL cp_fm_init_random(mo_virt, nvirt) + CALL cp_fm_init_random(mo_virtual, nmo_virtual) NULLIFY (local_preconditioner) CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, & - matrix_c_fm=mo_virt, matrix_orthogonal_space_fm=fm_ref, & + matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, & eps_gradient=scf_control%eps_lumos, & preconditioner=local_preconditioner, & iter_max=scf_control%max_iter_lumos, & - size_ortho_space=nmo) + size_ortho_space=nmo_available) - CALL calculate_subspace_eigenvalues(mo_virt, ks_matrix(ispin)%matrix, & - evals_virt) + CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, & + evals_virtual) ! We need to send the copy of MOs to preserve the sign - CALL cp_fm_create(fm_dummy, fm_ref%matrix_struct) - CALL cp_fm_to_fm(fm_ref, fm_dummy) + CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct) + CALL cp_fm_to_fm(mo_ref, fm_dummy) CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, & evals_arg=eigenvalues(:, ispin), do_rotation=.TRUE.) @@ -750,47 +686,62 @@ CONTAINS CALL get_mo_set(mo_set_inactive, mo_coeff=fm_target_inactive) ! copy orbitals - nel = active_space_env%nelec_total(ispin) - l = 1 - DO i = 1, nmo + nvirt - nel = MAX(NINT(nel - focc), 0) - DO j = 1, SIZE(active_orbitals, 1) - m = active_orbitals(j, ispin) - IF (m == i) THEN - IF (i > nmo) THEN - CALL cp_fm_to_fm(mo_virt, fm_target_active, 1, i - nmo, m) - mo_set_active%eigenvalues(m) = evals_virt(i - nmo) - mo_set_active%occupation_numbers(m) = 0.0 - ELSE - CALL cp_fm_to_fm(fm_dummy, fm_target_active, 1, i, m) - mo_set_active%eigenvalues(m) = eigenvalues(i, ispin) - mo_set_active%occupation_numbers(m) = 1.0 - END IF - active_space_env%active_orbitals(j, ispin) = i - EXIT + nmo_inactive_remaining = nmo_inactive + DO i = 1, max_orb_ind + ! case for i being an active orbital + IF (ANY(active_space_env%active_orbitals(:, ispin) == i)) THEN + IF (i > nmo_available) THEN + CALL cp_fm_to_fm(mo_virtual, fm_target_active, 1, i - nmo_available, i) + eigenvalues(i, ispin) = evals_virtual(i - nmo_available) + mo_set_active%occupation_numbers(i) = 0.0 + ELSE + CALL cp_fm_to_fm(fm_dummy, fm_target_active, 1, i, i) + mo_set_active%occupation_numbers(i) = mos(ispin)%occupation_numbers(i) END IF - END DO - IF (j > SIZE(active_orbitals, 1)) THEN - ! inactive orbital - CALL cp_fm_to_fm(fm_dummy, fm_target_inactive, 1, i, l) - IF (nel > 0) THEN - active_space_env%inactive_orbitals(l, ispin) = i - mo_set_inactive%occupation_numbers(m) = 1.0 + mo_set_active%eigenvalues(i) = eigenvalues(i, ispin) + ! if it was not an active orbital, check whether it is an inactive orbital + ELSEIF (nmo_inactive_remaining > 0) THEN + CALL cp_fm_to_fm(fm_dummy, fm_target_inactive, 1, i, i) + ! store on the fly the mapping of inactive orbitals + active_space_env%inactive_orbitals(nmo_inactive - nmo_inactive_remaining + 1, ispin) = i + mo_set_inactive%eigenvalues(i) = eigenvalues(i, ispin) + mo_set_inactive%occupation_numbers(i) = mos(ispin)%occupation_numbers(i) + ! hack: set homo and lumo manually + IF (nmo_inactive_remaining == 1) THEN + mo_set_inactive%homo = i + mo_set_inactive%lfomo = i + 1 END IF - mo_set_inactive%eigenvalues(m) = eigenvalues(m, ispin) - l = l + 1 + nmo_inactive_remaining = nmo_inactive_remaining - 1 + ELSE + CYCLE END IF END DO ! Release - DEALLOCATE (evals_virt) + DEALLOCATE (evals_virtual) CALL cp_fm_release(fm_dummy) - CALL cp_fm_release(mo_virt) + CALL cp_fm_release(mo_virtual) END DO IF (iw > 0) THEN DO ispin = 1, nspins - WRITE (iw, '(/,T4,A,I3)') "Active Orbital Indices for spin", ispin + WRITE (iw, '(/,T3,A,I3,T66,A)') "Orbital Energies and Selection for spin", ispin, "[atomic units]" + + DO i = 1, max_orb_ind, 4 + jm = MIN(3, max_orb_ind - i) + WRITE (iw, '(T4)', advance="no") + DO j = 0, jm + IF (ANY(active_space_env%active_orbitals(:, ispin) == i + j)) THEN + WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [A]" + ELSEIF (ANY(active_space_env%inactive_orbitals(:, ispin) == i + j)) THEN + WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [I]" + ELSE + WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [V]" + END IF + END DO + WRITE (iw, *) + END DO + WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4 jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i) WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm) @@ -885,45 +836,44 @@ CONTAINS active_space_env%eri%periodicity(1:3) = cell%perd(1:3) END IF IF (iw > 0) THEN - WRITE (iw, '(/,T4,A)') "Calculation of Electron Repulsion Integrals" + WRITE (iw, '(/,T3,A)') "Calculation of Electron Repulsion Integrals" SELECT CASE (eri_method) CASE (eri_method_full_gpw) - WRITE (iw, '(T4,A,T48,A)') "Integration method", "GPW Fourier transform over MOs" + WRITE (iw, '(T3,A,T50,A)') "Integration method", "GPW Fourier transform over MOs" CASE (eri_method_gpw_ht) - WRITE (iw, '(T4,A,T43,A)') "Integration method", "Half transformed integrals from GPW" + WRITE (iw, '(T3,A,T44,A)') "Integration method", "Half transformed integrals from GPW" CASE DEFAULT CPABORT("Unknown ERI method") END SELECT SELECT CASE (eri_operator) CASE (eri_operator_coulomb) - WRITE (iw, '(T4,A,T40,A,T73,A)') "ERI operator", "Coulomb", "<1/R>" + WRITE (iw, '(T3,A,T66,A)') "ERI operator", "Coulomb <1/R>" CASE (eri_operator_yukawa) - WRITE (iw, '(T4,A,T40,A,T67,A)') "ERI operator", "Yukawa", "" - WRITE (iw, '(T4,A,T63,F14.6)') "ERI operator parameter", eri_op_param + WRITE (iw, '(T3,A,T59,A)') "ERI operator", "Yukawa " + WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param CASE (eri_operator_erf) - WRITE (iw, '(T4,A,T40,A,T67,A)') "ERI operator", "Error function", "" - WRITE (iw, '(T4,A,T63,F14.6)') "ERI operator parameter", eri_op_param + WRITE (iw, '(T3,A,T53,A)') "ERI operator", "Error function " + WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param CASE (eri_operator_erfc) - WRITE (iw, '(T4,A,T40,A,T66,A)') "ERI operator", "Compl. error function", "" - WRITE (iw, '(T4,A,T63,F14.6)') "ERI operator parameter", eri_op_param + WRITE (iw, '(T3,A,T45,A)') "ERI operator", "Compl. error function " + WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param CASE (eri_operator_gaussian) - WRITE (iw, '(T4,A,T40,A,T66,A)') "ERI operator", "Gaussian attenuated", "" - WRITE (iw, '(T4,A,T63,F14.6)') "ERI operator parameter", eri_op_param + WRITE (iw, '(T3,A,T41,A)') "ERI operator", "Gaussian attenuated " + WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param CASE DEFAULT CPABORT("Unknown ERI operator") END SELECT - WRITE (iw, '(T4,A,T66,E12.4)') "Accuracy of ERI", eri_eps_int - WRITE (iw, '(T4,A,T69,3I3)') "Periodicity", active_space_env%eri%periodicity(1:3) + WRITE (iw, '(T3,A,T68,E12.4)') "Accuracy of ERI", eri_eps_int + WRITE (iw, '(T3,A,T71,3I3)') "Periodicity", active_space_env%eri%periodicity(1:3) IF (PRODUCT(active_space_env%eri%periodicity(1:3)) == 0) THEN - IF (eri_rcut > 0.0_dp) WRITE (iw, '(T4,A,T63,F14.6)') "Periodicity (Cutoff)", eri_rcut + IF (eri_rcut > 0.0_dp) WRITE (iw, '(T3,A,T65,F14.6)') "Periodicity (Cutoff)", eri_rcut END IF IF (nspins < 2) THEN - WRITE (iw, '(T4,A,T66,I12)') "Total Number of ERI", (nactive_orb(1)**4)/8 + WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI", (nmo_active**4)/8 ELSE - WRITE (iw, '(T4,A,T66,I12)') "Total Number of ERI (aa|aa)", (nactive_orb(1)**4)/8 - WRITE (iw, '(T4,A,T66,I12)') "Total Number of ERI (bb|bb)", (nactive_orb(2)**4)/8 - WRITE (iw, '(T4,A,T66,I12)') "Total Number of ERI (aa|bb)", & - (nactive_orb(1)**2)*(nactive_orb(2)**2)/4 + WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|aa)", (nmo_active**4)/8 + WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (bb|bb)", (nmo_active**4)/8 + WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|bb)", (nmo_active**4)/4 END IF END IF @@ -963,10 +913,12 @@ CONTAINS active_space_env%eri%eri_gpw%print_level = eri_print CALL section_vals_val_get(as_input, "ERI_GPW%STORE_WFN", l_val=store_wfn) active_space_env%eri%eri_gpw%store_wfn = store_wfn - active_space_env%eri%eri_gpw%redo_poisson = (ex_operator .OR. ex_perd) + ! Always redo Poisson solver for now + active_space_env%eri%eri_gpw%redo_poisson = .TRUE. + ! active_space_env%eri%eri_gpw%redo_poisson = (ex_operator .OR. ex_perd) IF (iw > 0) THEN - WRITE (iw, '(T4,A,T68,F10.4)') "ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff - WRITE (iw, '(T4,A,T68,F10.4)') "ERI_GPW| Relative energy cutoff [Ry]", eri_rel_cutoff + WRITE (iw, '(/,T2,A,T71,F10.1)') "ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff + WRITE (iw, '(T2,A,T71,F10.1)') "ERI_GPW| Relative energy cutoff [Ry]", eri_rel_cutoff END IF ! CALL calculate_eri_gpw(active_space_env%mos_active, active_space_env%active_orbitals, active_space_env%eri, qs_env, iw) @@ -979,15 +931,15 @@ CONTAINS eri_mat => active_space_env%eri%eri(isp)%csr_mat nze_percentage = 100.0_dp*(REAL(eri_mat%nze_total, KIND=dp) & /REAL(eri_mat%nrows_total, KIND=dp))/REAL(eri_mat%ncols_total, KIND=dp) - WRITE (iw, '(T4,A,I2,T30,A,T66,I12)') "ERI_GPW| Spinmatrix:", isp, & + WRITE (iw, '(/,T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, & "Number of CSR non-zero elements:", eri_mat%nze_total - WRITE (iw, '(T4,A,I2,T30,A,T66,F12.4)') "ERI_GPW| Spinmatrix:", isp, & + WRITE (iw, '(T2,A,I2,T30,A,T68,F12.4)') "ERI_GPW| Spinmatrix:", isp, & "Percentage CSR non-zero elements:", nze_percentage - WRITE (iw, '(T4,A,I2,T30,A,T66,I12)') "ERI_GPW| Spinmatrix:", isp, & + WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, & "nrows_total", eri_mat%nrows_total - WRITE (iw, '(T4,A,I2,T30,A,T66,I12)') "ERI_GPW| Spinmatrix:", isp, & + WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, & "ncols_total", eri_mat%ncols_total - WRITE (iw, '(T4,A,I2,T30,A,T66,I12)') "ERI_GPW| Spinmatrix:", isp, & + WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, & "nrows_local", eri_mat%nrows_local END DO END IF @@ -1003,34 +955,27 @@ CONTAINS SELECT CASE (mselect) CASE DEFAULT CPABORT("Unknown orbital selection method") - CASE (casci_canonical) + CASE (casci_canonical, manual_selection) focc = 2.0_dp IF (nspins == 2) focc = 1.0_dp DO isp = 1, nspins fmat => active_space_env%p_active(isp) CALL cp_fm_set_all(fmat, alpha=0.0_dp) - n1 = active_space_env%nelec_total(isp) - active_space_env%nelec_inactive_spinwise(isp) - DO i = 1, nactive_orb(isp) + IF (nspins == 2) THEN + IF (isp == 1) THEN + nel = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2 + ELSE + nel = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2 + END IF + ELSE + nel = active_space_env%nelec_active + END IF + DO i = 1, nmo_active m = active_space_env%active_orbitals(i, isp) - fel = MIN(focc, REAL(n1, KIND=dp)) + fel = MIN(focc, REAL(nel, KIND=dp)) CALL cp_fm_set_element(fmat, m, m, fel) - n1 = n1 - NINT(fel) - n1 = MAX(n1, 0) - END DO - END DO - CASE (manual_selection) - focc = 2.0_dp - IF (nspins == 2) focc = 1.0_dp - DO isp = 1, nspins - fmat => active_space_env%p_active(isp) - CALL cp_fm_set_all(fmat, alpha=0.0_dp) - n1 = active_space_env%nelec_total(isp) - active_space_env%nelec_inactive_spinwise(isp) - DO i = 1, nactive_orb(isp) - m = active_space_env%active_orbitals(i, isp) - fel = MIN(focc, REAL(n1, KIND=dp)) - CALL cp_fm_set_element(fmat, m, m, fel) - n1 = n1 - NINT(fel) - n1 = MAX(n1, 0) + nel = nel - NINT(fel) + nel = MAX(nel, 0) END DO END DO CASE (wannier_projection) @@ -1039,17 +984,40 @@ CONTAINS CPABORT("NOT IMPLEMENTED") END SELECT - ! Read active space density matix if requested - !IF (active_space_env%read_p_act) CALL read_active_density(active_space_env, para_env) - - ! calculate one-electron operators in the subspace - CALL update_active_space(qs_env, active_space_env, as_input) - - ! Output a FCIDUMP file - ! CALL fcidump(active_space_env, as_input) + ! transform KS/Fock, Vxc and Hcore to AS MO basis + CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env) + ! set the reference energy in the active space + CALL get_qs_env(qs_env, energy=energy) + active_space_env%energy_ref = energy%total + ! calculate inactive energy and embedding potential + CALL subspace_fock_matrix(active_space_env) + ! associate the active space environment with the qs environment CALL set_qs_env(qs_env, active_space=active_space_env) + ! Perform the embedding calculation only if qiskit is specified + CALL section_vals_val_get(as_input, "AS_SOLVER", i_val=as_solver) + SELECT CASE (as_solver) + CASE (no_solver) + IF (iw > 0) THEN + WRITE (iw, '(/,T3,A)') "No active space solver specified, skipping embedding calculation" + END IF + CASE (qiskit_solver) + CALL rsdft_embedding(qs_env, active_space_env, as_input) + CASE DEFAULT + CPABORT("Unknown active space solver") + END SELECT + + ! Output a FCIDUMP file if requested + IF (active_space_env%fcidump) CALL fcidump(active_space_env, as_input) + + ! Output a QCSchema file if requested + IF (active_space_env%qcschema) THEN + CALL qcschema_env_create(qcschema_env, qs_env) + CALL qcschema_to_hdf5(qcschema_env, active_space_env%qcschema_filename) + CALL qcschema_env_release(qcschema_env) + END IF + IF (iw > 0) THEN WRITE (iw, '(/,T2,A)') & '!-------------------- End of Active Space Interface --------------------------!' @@ -1086,7 +1054,6 @@ CONTAINS ! Kohn-Sham / Fock operator CALL cp_fm_release(active_space_env%ks_sub) - CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix) IF (SIZE(ks_matrix, 2) > 1) THEN CPABORT("No k-points allowed at this point") @@ -1298,9 +1265,11 @@ CONTAINS CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, & poisson_env=poisson_env) IF (eri_env%eri_gpw%redo_poisson) THEN - IF (.NOT. (poisson_env%parameters%solver == pw_poisson_analytic .OR. & - poisson_env%parameters%solver == pw_poisson_periodic)) THEN - CPABORT("Only use simple analytic Poisson solvers with ERI calculation") + ! We need to rebuild the Poisson solver on the fly + IF (SUM(eri_env%periodicity) /= 0) THEN + poisson_env%parameters%solver = pw_poisson_periodic + ELSE + poisson_env%parameters%solver = pw_poisson_analytic END IF poisson_env%parameters%periodic = eri_env%periodicity CALL pw_poisson_rebuild(poisson_env) @@ -1314,9 +1283,31 @@ CONTAINS END DO poisson_env%green_fft%radius = rc END IF + CALL pw_eri_green_create(poisson_env%green_fft, eri_env) - IF (print1 .AND. iw > 0) THEN - WRITE (iw, "(T4,'ERI_GPW|',' Redefine Poisson Greens function ')") + + IF (iw > 0) THEN + CALL get_qs_env(qs_env, cell=cell) + IF (SUM(cell%perd) /= SUM(eri_env%periodicity)) THEN + IF (SUM(eri_env%periodicity) /= 0) THEN + WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") & + "ERI_GPW| Switching Poisson solver to", "PERIODIC" + ELSE + WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") & + "ERI_GPW| Switching Poisson solver to", "ANALYTIC" + END IF + END IF + ! print out the Greens function to check it matches the Poisson solver + SELECT CASE (poisson_env%green_fft%method) + CASE (PERIODIC3D) + WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") & + "ERI_GPW| Poisson Greens function", "PERIODIC" + CASE (ANALYTIC0D) + WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") & + "ERI_GPW| Poisson Greens function", "ANALYTIC" + CASE DEFAULT + CPABORT("Wrong Greens function setup") + END SELECT END IF END IF @@ -1407,6 +1398,7 @@ CONTAINS ! calculate charge distribution and potential rho_r%cr3d = wfn1%cr3d*wfn2%cr3d IF (print2 .AND. iw > 0) THEN + ! TODO: this call hangs in parallel erint = pw_integrate_function(rho_r)/dvol WRITE (iw, "(T4,'ERI_GPW| Integral rho_ab ',T32,2I4,' [',I1,']',T58,G20.14)") & iwa1, iwa2, isp1, erint @@ -1421,13 +1413,13 @@ CONTAINS CPABORT("Not available") ELSEIF (eri_env%method == eri_method_full_gpw) THEN DO isp2 = isp1, nspins - CALL get_mo_set(mo_set=mos(isp1), nmo=nmo2) + CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2) nx = (nmo2*(nmo2 + 1))/2 ALLOCATE (eri(nx), eri_index(nx)) icount2 = 0 iwbs = 1 IF (isp1 == isp2) iwbs = i1 - isp = (isp1 - 1)*isp2 - ((isp1 - 1)*(isp1 - 2))/2 + (isp2 - isp1 + 1) + isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1) DO i3 = iwbs, SIZE(orbitals, 1) iwb1 = orbitals(i3, isp2) IF (eri_env%eri_gpw%store_wfn) THEN @@ -1810,135 +1802,130 @@ CONTAINS checksum = 0.0_dp logger => cp_get_default_logger() - IF (BTEST(cp_print_key_should_output(logger%iter_info, & - as_input, "FCIDUMP"), cp_p_file)) THEN - iw = cp_print_key_unit_nr(logger, as_input, "FCIDUMP", & - extension=".fcidump", file_status="REPLACE", file_action="WRITE", file_form="FORMATTED") - ! - nspins = active_space_env%nspins - norb = SIZE(active_space_env%active_orbitals, 1) - IF (nspins == 1) THEN - ASSOCIATE (ms2 => active_space_env%multiplicity, & - nelec => active_space_env%nelec_active) + iw = cp_print_key_unit_nr(logger, as_input, "FCIDUMP", & + extension=".fcidump", file_status="REPLACE", file_action="WRITE", file_form="FORMATTED") + ! + nspins = active_space_env%nspins + norb = SIZE(active_space_env%active_orbitals, 1) + IF (nspins == 1) THEN + ASSOCIATE (ms2 => active_space_env%multiplicity, & + nelec => active_space_env%nelec_active) - IF (iw > 0) THEN - WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, "," - isym = 1 - WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb) - isym = 0 - WRITE (iw, "(A,I1,A)") " ISYM=", isym, "," - WRITE (iw, "(A)") " /" - END IF - ! - ! Print integrals: ERI - CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, & - eri_fcidump_print(iw, 1, 1), 1, 1) - CALL eri_checksum%set(1, 1) - CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1) + IF (iw > 0) THEN + WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, "," + isym = 1 + WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb) + isym = 0 + WRITE (iw, "(A,I1,A)") " ISYM=", isym, "," + WRITE (iw, "(A)") " /" + END IF + ! + ! Print integrals: ERI + CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, & + eri_fcidump_print(iw, 1, 1), 1, 1) + CALL eri_checksum%set(1, 1) + CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1) - ! Print integrals: Fij - ! replicate Fock matrix - nmo = active_space_env%eri%norb - ALLOCATE (fmat(nmo, nmo)) - ! TODO: extend to arbitrary active orbitals - CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat) - IF (iw > 0) THEN - i3 = 0; i4 = 0 - DO m1 = 1, SIZE(active_space_env%active_orbitals, 1) - i1 = active_space_env%active_orbitals(m1, 1) - DO m2 = m1, SIZE(active_space_env%active_orbitals, 1) - i2 = active_space_env%active_orbitals(m2, 1) - checksum = checksum + ABS(fmat(i1, i2)) - WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4 - END DO + ! Print integrals: Fij + ! replicate Fock matrix + nmo = active_space_env%eri%norb + ALLOCATE (fmat(nmo, nmo)) + CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat) + IF (iw > 0) THEN + i3 = 0; i4 = 0 + DO m1 = 1, SIZE(active_space_env%active_orbitals, 1) + i1 = active_space_env%active_orbitals(m1, 1) + DO m2 = m1, SIZE(active_space_env%active_orbitals, 1) + i2 = active_space_env%active_orbitals(m2, 1) + checksum = checksum + ABS(fmat(i1, i2)) + WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4 END DO - END IF - DEALLOCATE (fmat) - ! Print energy - esub = active_space_env%energy_inactive - i1 = 0; i2 = 0; i3 = 0; i4 = 0 - checksum = checksum + ABS(esub) - IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4 - END ASSOCIATE + END DO + END IF + DEALLOCATE (fmat) + ! Print energy + esub = active_space_env%energy_inactive + i1 = 0; i2 = 0; i3 = 0; i4 = 0 + checksum = checksum + ABS(esub) + IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4 + END ASSOCIATE - ELSE - ASSOCIATE (ms2 => active_space_env%multiplicity, & - nelec => active_space_env%nelec_active) + ELSE + ASSOCIATE (ms2 => active_space_env%multiplicity, & + nelec => active_space_env%nelec_active) - IF (iw > 0) THEN - WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, "," - isym = 1 - WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb) - isym = 0 - WRITE (iw, "(A,I1,A)") " ISYM=", isym, "," - WRITE (iw, "(A,I1,A)") " UHF=", 1, "," - WRITE (iw, "(A)") " /" - END IF - ! - ! Print integrals: ERI - ! alpha-alpha - CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, & - eri_fcidump_print(iw, 1, 1), 1, 1) - CALL eri_checksum%set(1, 1) - CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1) - ! alpha-beta - CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, & - eri_fcidump_print(iw, 1, norb + 1), 1, 2) - CALL eri_checksum%set(1, norb + 1) - CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, eri_checksum, 1, 2) - ! beta-beta - CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, & - eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2) - CALL eri_checksum%set(norb + 1, norb + 1) - CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2) - ! Print integrals: Fij - ! alpha - nmo = active_space_env%eri%norb - ALLOCATE (fmat(nmo, nmo)) - CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat) - IF (iw > 0) THEN - i3 = 0; i4 = 0 - DO m1 = 1, norb - i1 = active_space_env%active_orbitals(m1, 1) - DO m2 = m1, norb - i2 = active_space_env%active_orbitals(m2, 1) - checksum = checksum + ABS(fmat(i1, i2)) - WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4 - END DO + IF (iw > 0) THEN + WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, "," + isym = 1 + WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb) + isym = 0 + WRITE (iw, "(A,I1,A)") " ISYM=", isym, "," + WRITE (iw, "(A,I1,A)") " UHF=", 1, "," + WRITE (iw, "(A)") " /" + END IF + ! + ! Print integrals: ERI + ! alpha-alpha + CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, & + eri_fcidump_print(iw, 1, 1), 1, 1) + CALL eri_checksum%set(1, 1) + CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1) + ! alpha-beta + CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, & + eri_fcidump_print(iw, 1, norb + 1), 1, 2) + CALL eri_checksum%set(1, norb + 1) + CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, eri_checksum, 1, 2) + ! beta-beta + CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, & + eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2) + CALL eri_checksum%set(norb + 1, norb + 1) + CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2) + ! Print integrals: Fij + ! alpha + nmo = active_space_env%eri%norb + ALLOCATE (fmat(nmo, nmo)) + CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat) + IF (iw > 0) THEN + i3 = 0; i4 = 0 + DO m1 = 1, norb + i1 = active_space_env%active_orbitals(m1, 1) + DO m2 = m1, norb + i2 = active_space_env%active_orbitals(m2, 1) + checksum = checksum + ABS(fmat(i1, i2)) + WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4 END DO - END IF - DEALLOCATE (fmat) - ! beta - ALLOCATE (fmat(nmo, nmo)) - CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat) - IF (iw > 0) THEN - i3 = 0; i4 = 0 - DO m1 = 1, SIZE(active_space_env%active_orbitals, 1) - i1 = active_space_env%active_orbitals(m1, 2) - DO m2 = m1, SIZE(active_space_env%active_orbitals, 1) - i2 = active_space_env%active_orbitals(m2, 2) - checksum = checksum + ABS(fmat(i1, i2)) - WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4 - END DO + END DO + END IF + DEALLOCATE (fmat) + ! beta + ALLOCATE (fmat(nmo, nmo)) + CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat) + IF (iw > 0) THEN + i3 = 0; i4 = 0 + DO m1 = 1, SIZE(active_space_env%active_orbitals, 1) + i1 = active_space_env%active_orbitals(m1, 2) + DO m2 = m1, SIZE(active_space_env%active_orbitals, 1) + i2 = active_space_env%active_orbitals(m2, 2) + checksum = checksum + ABS(fmat(i1, i2)) + WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4 END DO - END IF - DEALLOCATE (fmat) - ! Print energy - esub = active_space_env%energy_inactive - i1 = 0; i2 = 0; i3 = 0; i4 = 0 - checksum = checksum + ABS(esub) - IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4 - END ASSOCIATE - END IF - ! - CALL cp_print_key_finished_output(iw, logger, as_input, "FCIDUMP") - - !>> - iw = cp_logger_get_default_io_unit(logger) - IF (iw > 0) WRITE (iw, '(T4,A,T66,F12.8)') "FCIDUMP| Checksum:", eri_checksum%checksum + checksum - !<< - + END DO + END IF + DEALLOCATE (fmat) + ! Print energy + esub = active_space_env%energy_inactive + i1 = 0; i2 = 0; i3 = 0; i4 = 0 + checksum = checksum + ABS(esub) + IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4 + END ASSOCIATE END IF + ! + CALL cp_print_key_finished_output(iw, logger, as_input, "FCIDUMP") + + !>> + iw = cp_logger_get_default_io_unit(logger) + IF (iw > 0) WRITE (iw, '(T4,A,T66,F12.8)') "FCIDUMP| Checksum:", eri_checksum%checksum + checksum + !<< END SUBROUTINE fcidump @@ -1987,9 +1974,7 @@ CONTAINS nspins = active_space_env%nspins IF (nspins == 1) THEN - CALL get_mo_set(active_space_env%mos_active(1), nmo=norb) - NULLIFY (mo_coef) - CALL get_mo_set(active_space_env%mos_active(1), mo_coeff=mo_coef) + CALL get_mo_set(active_space_env%mos_active(1), nmo=norb, mo_coeff=mo_coef) ! ! Loop over ERI, calculate subspace HF energy and Fock matrix ! @@ -1997,21 +1982,30 @@ CONTAINS ALLOCATE (ks_mat(norb, norb), ks_ref(norb, norb), p_mat(norb, norb)) ks_ref = 0.0_dp - CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_mat) + ! ks_mat contains the KS/Fock matrix (of full density) projected onto the AS MO subspace (f_ref in eq. 19) CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_mat) - ! - ! + CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_mat) + + ! compute ks_ref = V_H[rho^A] + V_HFX[rho^A] eri => active_space_env%eri%eri(1)%csr_mat CALL build_subspace_fock_matrix(active_space_env%active_orbitals, eri, p_mat, ks_ref) - ! - ! calculate energy - eeri = 0.0_dp + + ! compute eeri = E_H[rho^A] + E_HFX[rho^A] as + ! eeri = 1/2 * (SUM_pq (V_H[rho^A] + V_HFX[rho^A])_pq * D^A_pq) eeri = 0.5_dp*SUM(ks_ref*p_mat) + + ! now calculate the inactive energy acoording to eq. 19, that is + ! esub = E^I = E_ref - f_ref .* D^A + E_H[rho^A] + E_HFX[rho^A] + ! where f^ref = ks_mat, which is the KS/Fock matrix in MO basis, transformed previously + ! and is equal to ks_mat = h^0 + V_core + V_H[rho] + V_HFX[rho] esub = eref - SUM(ks_mat(1:norb, 1:norb)*p_mat(1:norb, 1:norb)) + eeri + + ! reuse ks_mat to store f^I = f^ref - (V_H[rho^A] + V_HFX[rho^A]) according to eq. 20 ks_mat(1:norb, 1:norb) = ks_mat(1:norb, 1:norb) - ks_ref(1:norb, 1:norb) - ! + ! this is now the embedding potential for the AS calculation! + active_space_env%energy_inactive = esub - ! + CALL cp_fm_release(active_space_env%fock_sub) ALLOCATE (active_space_env%fock_sub(nspins)) DO is = 1, nspins @@ -2343,172 +2337,51 @@ CONTAINS END SUBROUTINE create_pro_basis ! ************************************************************************************************** -!> \brief ... -!> \param qs_env ... -!> \param active_space_env ... -!> \param as_input ... +!> \brief Update the density matrix in AO basis with the active density contribution +!> \param active_space_env the active space environment +!> \param rho_ao the density matrix in AO basis ! ************************************************************************************************** - SUBROUTINE update_active_space(qs_env, active_space_env, as_input) - TYPE(qs_environment_type), POINTER :: qs_env + SUBROUTINE update_density_ao(active_space_env, rho_ao) TYPE(active_space_type), POINTER :: active_space_env - TYPE(section_vals_type), POINTER :: as_input - - CHARACTER(len=*), PARAMETER :: routineN = 'update_active_space' - - INTEGER :: handle, ispin, nao, nmo, nspins - LOGICAL :: should_stop - TYPE(cp_fm_type) :: lamat, vec - TYPE(cp_fm_type), DIMENSION(:), POINTER :: p_ref - TYPE(cp_fm_type), POINTER :: fm_active, pmat TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao - TYPE(dbcsr_type), POINTER :: pinact - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos - TYPE(mp_para_env_type), POINTER :: para_env - TYPE(qs_energy_type), POINTER :: energy - TYPE(qs_rho_type), POINTER :: rho - CALL timeset(routineN, handle) - - ! density matrix - CALL get_qs_env(qs_env, rho=rho) - CALL qs_rho_get(rho, rho_ao=rho_ao) + INTEGER :: ispin, nao, nmo, nspins + TYPE(cp_fm_type) :: R, U + TYPE(cp_fm_type), POINTER :: C_active, p_active_mo + TYPE(dbcsr_type), POINTER :: p_inactive_ao + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active + ! Transform the AS density matrix P_MO to the atomic orbital basis, + ! this is simply C * P_MO * C^T nspins = active_space_env%nspins - p_ref => active_space_env%p_active - mos => active_space_env%mos_active + mos_active => active_space_env%mos_active DO ispin = 1, nspins - pinact => active_space_env%pmat_inactive(ispin)%matrix - CALL dbcsr_copy(rho_ao(ispin)%matrix, pinact) - ! create active density matrix in AO basis - pmat => p_ref(ispin) + ! size of p_inactive_ao is (nao x nao) + p_inactive_ao => active_space_env%pmat_inactive(ispin)%matrix + + ! copy p_inactive_ao to rho_ao + CALL dbcsr_copy(rho_ao(ispin)%matrix, p_inactive_ao) + + ! size of p_active_mo is (nmo x nmo) + p_active_mo => active_space_env%p_active(ispin) + + ! calculate R = p_mo + CALL cp_fm_create(R, p_active_mo%matrix_struct) + CALL cp_fm_to_fm(p_active_mo, R) + + ! calculate U = C * p_mo + CALL get_mo_set(mos_active(ispin), mo_coeff=C_active, nao=nao, nmo=nmo) + CALL cp_fm_create(U, C_active%matrix_struct) + CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, C_active, R, 0.0_dp, U) - CALL cp_fm_create(lamat, pmat%matrix_struct) - CALL cp_fm_to_fm(pmat, lamat) - ! calclulate SQRT(P): P = U*U - CALL cp_fm_sqrt(lamat) - ! R = C * U - CALL get_mo_set(mos(ispin), mo_coeff=fm_active, nao=nao, nmo=nmo) - CALL cp_fm_create(vec, fm_active%matrix_struct) - CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, fm_active, lamat, 0.0_dp, vec) CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(ispin)%matrix, & - matrix_v=vec, ncol=nmo, alpha=1.0_dp) - CALL cp_fm_release(vec) - CALL cp_fm_release(lamat) + matrix_v=U, matrix_g=C_active, ncol=nmo, alpha=1.0_dp) + + CALL cp_fm_release(R) + CALL cp_fm_release(U) END DO - ! let's also calculate Vxc - qs_env%requires_matrix_vxc = .TRUE. - - ! new density and Kohn-Sham matrix - CALL qs_rho_update_rho(rho, qs_env) - CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) - CALL qs_ks_update_qs_env(qs_env) - - ! reference energy - CALL get_qs_env(qs_env, energy=energy) - active_space_env%energy_ref = energy%total - - ! update operators - mos => active_space_env%mos_active - CALL calculate_operators(mos, qs_env, active_space_env) - - ! Read active space density matix if requested - IF (active_space_env%read_p_act) THEN - ! inactive subspace energy and active Fock operator - CALL subspace_fock_matrix(active_space_env) - CALL fcidump(active_space_env, as_input) - CALL external_control(should_stop, "AS", target_time=qs_env%target_time, start_time=m_walltime(), force_check=.TRUE.) - - DO WHILE (.NOT. should_stop) - - CALL get_qs_env(qs_env, para_env=para_env) - CALL read_active_density(active_space_env, para_env) - - ! Debug: remove or refactor - nspins = active_space_env%nspins - p_ref => active_space_env%p_active - mos => active_space_env%mos_active - DO ispin = 1, nspins - pinact => active_space_env%pmat_inactive(ispin)%matrix - CALL dbcsr_copy(rho_ao(ispin)%matrix, pinact) - ! create active density matrix in AO basis - pmat => p_ref(ispin) - - CALL cp_fm_create(lamat, pmat%matrix_struct) - CALL cp_fm_to_fm(pmat, lamat) - ! calclulate SQRT(P): P = U*U - CALL cp_fm_sqrt(lamat) - ! R = C * U - CALL get_mo_set(mos(ispin), mo_coeff=fm_active, nao=nao, nmo=nmo) - CALL cp_fm_create(vec, fm_active%matrix_struct) - CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, fm_active, lamat, 0.0_dp, vec) - CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(ispin)%matrix, & - matrix_v=vec, ncol=nmo, alpha=1.0_dp) - CALL cp_fm_release(vec) - CALL cp_fm_release(lamat) - END DO - - ! Compute total energy with update density - ! let's also calculate Vxc - qs_env%requires_matrix_vxc = .TRUE. - - !CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.) - CALL qs_rho_update_rho(rho, qs_env) - CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) - CALL qs_ks_update_qs_env(qs_env) - - ! Updated reference energy - CALL get_qs_env(qs_env, energy=energy) - - ! inactive subspace energy and active Fock operator - CALL subspace_fock_matrix(active_space_env) - CALL fcidump(active_space_env, as_input) - CALL external_control(should_stop, "AS", target_time=qs_env%target_time, start_time=m_walltime(), force_check=.TRUE.) - active_space_env%energy_ref = energy%total - END DO - ELSE - ! inactive subspace energy and active Fock operator - CALL subspace_fock_matrix(active_space_env) - CALL fcidump(active_space_env, as_input) - END IF - - CALL timestop(handle) - - END SUBROUTINE update_active_space - -! ************************************************************************************************** -!> \brief Calculate the SQRT of a FM matrix -!> \param amat ... -! ************************************************************************************************** - SUBROUTINE cp_fm_sqrt(amat) - TYPE(cp_fm_type), INTENT(IN) :: amat - - INTEGER :: i, n - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigv - TYPE(cp_fm_type) :: bmat, vmat - - ! eigenvalues - CALL cp_fm_get_info(amat, nrow_global=n) - ALLOCATE (eigv(n)) - ! scratch arrays - CALL cp_fm_create(vmat, amat%matrix_struct) - CALL cp_fm_create(bmat, amat%matrix_struct) - ! diag - CALL cp_fm_syevd(amat, vmat, eigv) - ! SQRT - DO i = 1, n - eigv(i) = MAX(eigv(i), 0.0_dp) - eigv(i) = SQRT(eigv(i)) - END DO - CALL cp_fm_to_fm(vmat, bmat) - CALL cp_fm_column_scale(bmat, eigv) - CALL parallel_gemm("N", "T", n, n, n, 1.0_dp, bmat, vmat, 0.0_dp, amat) - ! clean up - DEALLOCATE (eigv) - CALL cp_fm_release(vmat) - CALL cp_fm_release(bmat) - - END SUBROUTINE cp_fm_sqrt + END SUBROUTINE update_density_ao ! ************************************************************************************************** !> \brief Print each value on the master node @@ -2559,135 +2432,485 @@ CONTAINS END FUNCTION eri_fcidump_checksum_func ! ************************************************************************************************** -!> \brief Reads active space density from external file -!> \param active_space_env ... -!> \param para_env ... +!> \brief Update active space density matrix from a fortran array +!> \param p_act_mo density matrix in active space MO basis +!> \param active_space_env active space environment +!> \param ispin spin index !> \author Vladimir Rybkin ! ************************************************************************************************** - SUBROUTINE read_active_density(active_space_env, para_env) + SUBROUTINE update_active_density(p_act_mo, active_space_env, ispin) + REAL(KIND=dp), DIMENSION(:) :: p_act_mo TYPE(active_space_type), POINTER :: active_space_env - TYPE(mp_para_env_type), POINTER :: para_env + INTEGER :: ispin - CHARACTER(LEN=default_path_length) :: p_act_filename - INTEGER :: i1, i2, ispin, l_global, LLL, m1, m2, & - m_global, MMM, nactive, ncol_local, & - nmo, nrow_local, nspins, p_active_unit - INTEGER, ALLOCATABLE, DIMENSION(:) :: nmos - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices - LOGICAL :: exist - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: p_act, p_act_beta, p_act_beta_read, & - p_act_read + INTEGER :: i1, i2, m1, m2, nmo_active + REAL(KIND=dp) :: mval + TYPE(cp_fm_type), POINTER :: p_active - ! Find out the spin states - nspins = SIZE(active_space_env%p_active) - ALLOCATE (nmos(nspins)) + p_active => active_space_env%p_active(ispin) + nmo_active = active_space_env%nmo_active - ! Find active DM dimensions - DO ispin = 1, nspins - CALL cp_fm_get_info(matrix=active_space_env%p_active(ispin), ncol_global=nmo) - nmos(ispin) = nmo + DO i1 = 1, nmo_active + m1 = active_space_env%active_orbitals(i1, ispin) + DO i2 = 1, nmo_active + m2 = active_space_env%active_orbitals(i2, ispin) + mval = p_act_mo(i2 + (i1 - 1)*nmo_active) + CALL cp_fm_set_element(p_active, m1, m2, mval) + END DO END DO - ! Allocate active DM - nmo = nmos(1) - ALLOCATE (p_act(nmo, nmo)) - p_act = 0.0_dp - IF (nspins .EQ. 2) THEN - nmo = nmos(2) - ALLOCATE (p_act_beta(nmo, nmo)) - p_act_beta = 0.0_dp + END SUBROUTINE update_active_density + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param active_space_env ... +!> \param as_input ... +! ************************************************************************************************** + SUBROUTINE rsdft_embedding(qs_env, active_space_env, as_input) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(active_space_type), POINTER :: active_space_env + TYPE(section_vals_type), POINTER :: as_input + + CHARACTER(len=*), PARAMETER :: routineN = 'rsdft_embedding' + +#ifdef __NO_SOCKETS + INTEGER :: handle + CALL timeset(routineN, handle) + CPABORT("CP2K was compiled with the __NO_SOCKETS option!") + MARK_USED(qs_env) + MARK_USED(active_space_env) + MARK_USED(as_input) +#else + + INTEGER :: client_fd, handle, iter, iw, max_iter, & + socket_fd + LOGICAL :: converged, do_scf_embedding, ionode + REAL(KIND=dp) :: alpha, delta_E, energy_corr, energy_new, & + energy_old, energy_scf, eps_iter, t1, & + t2 + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_rho_type), POINTER :: rho + + CALL timeset(routineN, handle) + + t1 = m_walltime() + + logger => cp_get_default_logger() + iw = cp_logger_get_default_io_unit(logger) + + CALL get_qs_env(qs_env, para_env=para_env) + ionode = para_env%is_source() + + ! get info from the input + CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding) + active_space_env%do_scf_embedding = do_scf_embedding + CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter) + CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter) + alpha = 0.0 + + ! create the socket and wait for the client to connect + CALL initialize_socket(socket_fd, client_fd, as_input, ionode) + CALL para_env%sync() + + ! send two-electron integrals to the client + CALL send_eri_to_client(client_fd, active_space_env, para_env) + + ! get pointer to density in ao basis + CALL get_qs_env(qs_env, rho=rho, energy=energy) + CALL qs_rho_get(rho, rho_ao=rho_ao) + + IF ((iw > 0)) THEN + WRITE (UNIT=iw, FMT="(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") & + "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", REPEAT("-", 78) + END IF + ! CALL cp_add_iter_level(logger%iter_info, "QS_SCF") + + iter = 0 + converged = .FALSE. + ! store the scf energy + energy_scf = active_space_env%energy_ref + energy_new = energy_scf + mos_active => active_space_env%mos_active + ! CALL set_qs_env(qs_env, active_space=active_space_env) + + ! start the self-consistent embedding loop + DO WHILE (iter < max_iter) + iter = iter + 1 + + ! send V_emb and E_ina to the active space solver and update + ! the active space environment with the new active energy and density + CALL send_fock_to_client(client_fd, active_space_env, para_env) + + ! update energies + energy_old = energy_new + energy_new = active_space_env%energy_total + energy_corr = energy_new - energy_scf + delta_E = energy_new - energy_old + + ! get timer + t2 = m_walltime() + ! print out progress + IF ((iw > 0)) THEN + WRITE (UNIT=iw, & + FMT="(T3,I4,T11,A,T21,F4.2,T28,F18.10,T49,F18.10,T70,ES11.2)") & + iter, 'P_Mix', t2 - t1, energy_corr, energy_new, delta_E + END IF + + ! update total density in AO basis with the AS contribution + CALL update_density_ao(active_space_env, rho_ao) ! rho_ao is updated + + ! calculate F_ks in AO basis (which contains Vxc) with the new density + qs_env%requires_matrix_vxc = .TRUE. + CALL qs_rho_update_rho(rho, qs_env) ! updates rho_r and rho_g using rho_ao + CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) ! set flags about the change + CALL qs_ks_update_qs_env(qs_env) ! this call actually calculates F_ks + + ! update the reference energy + active_space_env%energy_ref = energy%total + + ! transform KS/Fock, Vxc and Hcore from AO to MO basis + CALL calculate_operators(mos_active, qs_env, active_space_env) + + ! calculate the new inactive energy and embedding potential + CALL subspace_fock_matrix(active_space_env) + + ! check if it is a one-shot correction + IF (.NOT. active_space_env%do_scf_embedding) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") & + "*** one-shot embedding correction finished ***" + END IF + converged = .TRUE. + EXIT + ! check for convergence + ELSEIF (ABS(delta_E) <= eps_iter) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") & + "*** rsDFT embedding run converged in ", iter, " iteration(s) ***" + END IF + converged = .TRUE. + EXIT + END IF + + t1 = m_walltime() + END DO + + IF (.NOT. converged) THEN + IF (iw > 0) THEN + WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") & + "*** rsDFT embedding did not converged after ", iter, " iteration(s) ***" + END IF END IF - ! File name - p_act_filename = active_space_env%p_act_filename + ! update qs total energy to the final embedding energy + energy%total = active_space_env%energy_total - ! Read matrices into fm objects - IF (para_env%is_source()) THEN - ! Check, whether the file exists - INQUIRE (FILE=p_act_filename, exist=exist) - IF (.NOT. exist) & - CPABORT("Active density file not found. ") + CALL finalize_socket(socket_fd, client_fd, as_input, ionode) + CALL para_env%sync() +#endif - ! the actual active size may be different - nactive = SIZE(active_space_env%active_orbitals, 1) - ! Allocate active DM to read - ALLOCATE (p_act_read(nactive, nactive)) - p_act_read = 0.0_dp - IF (nspins .EQ. 2) THEN - ALLOCATE (p_act_beta_read(nactive, nactive)) - p_act_beta_read = 0.0_dp - END IF + CALL timestop(handle) - p_active_unit = -1 - CALL open_file(file_name=p_act_filename, & - file_action="READ", & - file_form="UNFORMATTED", & - file_status="UNKNOWN", & - file_access='STREAM', & - unit_number=p_active_unit) - ! Now read - READ (p_active_unit) p_act_read - DO i1 = 1, SIZE(active_space_env%active_orbitals, 1) - m1 = active_space_env%active_orbitals(i1, 1) - DO i2 = 1, SIZE(active_space_env%active_orbitals, 1) - m2 = active_space_env%active_orbitals(i2, 1) - p_act(m1, m2) = p_act_read(i1, i2) - END DO - END DO - DEALLOCATE (p_act_read) + END SUBROUTINE rsdft_embedding - IF (nspins .EQ. 2) THEN - READ (p_active_unit) p_act_beta_read - DO i1 = 1, SIZE(active_space_env%active_orbitals, 1) - m1 = active_space_env%active_orbitals(i1, 2) - DO i2 = 1, SIZE(active_space_env%active_orbitals, 1) - m2 = active_space_env%active_orbitals(i2, 2) - p_act_beta(m1, m2) = p_act_beta_read(i1, i2) - END DO - END DO - DEALLOCATE (p_act_beta_read) - END IF +#ifndef __NO_SOCKETS +! ************************************************************************************************** +!> \brief Creates the socket, spawns the client and connects to it +!> \param socket_fd the socket file descriptor +!> \param client_fd the client file descriptor +!> \param as_input active space inpute section +!> \param ionode logical flag indicating if the process is the master +! ************************************************************************************************** + SUBROUTINE initialize_socket(socket_fd, client_fd, as_input, ionode) + INTEGER, INTENT(OUT) :: socket_fd, client_fd + TYPE(section_vals_type), INTENT(IN), POINTER :: as_input + LOGICAL, INTENT(IN) :: ionode - ! Close restart file - IF (para_env%is_source()) THEN - CALL close_file(unit_number=p_active_unit) - END IF + CHARACTER(len=*), PARAMETER :: routineN = 'initialize_socket' + INTEGER, PARAMETER :: backlog = 10 + CHARACTER(len=default_path_length) :: hostname + INTEGER :: handle, iw, port, protocol + LOGICAL :: unix + TYPE(cp_logger_type), POINTER :: logger + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + iw = cp_logger_get_default_io_unit(logger) + + socket_fd = 0 + protocol = 0 + + ! protocol == 0 for UNIX, protocol > 0 for INET + CALL section_vals_val_get(as_input, "SOCKET%UNIX", l_val=unix) + IF (.NOT. unix) protocol = 1 + CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname) + CALL section_vals_val_get(as_input, "SOCKET%PORT", i_val=port) + + IF (ionode) THEN + CALL open_bind_socket(socket_fd, protocol, port, TRIM(hostname)//C_NULL_CHAR) + WRITE (iw, '(/,T3,A,A)') "@SERVER: Created socket with address ", TRIM(hostname) + CALL listen_socket(socket_fd, backlog) + + ! wait until a connetion request arrives + WRITE (iw, '(T3,A)') "@SERVER: Waiting for requests..." + CALL accept_socket(socket_fd, client_fd) + WRITE (iw, '(T3,A,I2)') "@SERVER: Accepted socket with fd ", client_fd END IF - ! Copy to full matrix structure - ! Broadcast the DM on all processes - CALL para_env%bcast(p_act, para_env%source) - IF (nspins .EQ. 2) CALL para_env%bcast(p_act_beta, para_env%source) + CALL timestop(handle) - DO ispin = 1, nspins - ! Copy to fm_type structure - CALL cp_fm_get_info(matrix=active_space_env%p_active(ispin), & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - col_indices=col_indices, & - row_indices=row_indices) + END SUBROUTINE initialize_socket - DO LLL = 1, nrow_local - l_global = row_indices(LLL) - DO MMM = 1, ncol_local - m_global = col_indices(MMM) - IF (ispin .EQ. 1) THEN - active_space_env%p_active(ispin)%local_data(LLL, MMM) = & - p_act(l_global, m_global) - ELSE - active_space_env%p_active(ispin)%local_data(LLL, MMM) = & - p_act_beta(l_global, m_global) +! ************************************************************************************************** +!> \brief Closes the connection to the socket and deletes the file +!> \param socket_fd the socket file descriptor +!> \param client_fd the client file descriptor +!> \param as_input active space inpute section +!> \param ionode logical flag indicating if the process is the master +! ************************************************************************************************** + SUBROUTINE finalize_socket(socket_fd, client_fd, as_input, ionode) + INTEGER, INTENT(IN) :: socket_fd, client_fd + TYPE(section_vals_type), INTENT(IN), POINTER :: as_input + LOGICAL, INTENT(IN) :: ionode + + CHARACTER(len=*), PARAMETER :: routineN = 'finalize_socket' + INTEGER, PARAMETER :: header_len = 12 + + CHARACTER(len=default_path_length) :: hostname + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname) + + IF (ionode) THEN + ! signal the client to quit + CALL writebuffer(client_fd, "QUIT ", header_len) + ! close the connection + CALL close_socket(client_fd) + CALL close_socket(socket_fd) + + ! delete the socket file + IF (file_exists(TRIM(hostname))) THEN + CALL remove_socket_file(TRIM(hostname)//C_NULL_CHAR) + END IF + END IF + + CALL timestop(handle) + + END SUBROUTINE finalize_socket + +! ************************************************************************************************** +!> \brief Sends the two-electron integrals to the client vie the socket +!> \param client_fd the client file descriptor +!> \param active_space_env active space environment +!> \param para_env parallel environment +! ************************************************************************************************** + SUBROUTINE send_eri_to_client(client_fd, active_space_env, para_env) + INTEGER, INTENT(IN) :: client_fd + TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env + TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env + + CHARACTER(len=*), PARAMETER :: routineN = 'send_eri_to_client' + INTEGER, PARAMETER :: header_len = 12 + + CHARACTER(len=default_string_length) :: header + INTEGER :: handle, iw + LOGICAL :: ionode + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb + TYPE(cp_logger_type), POINTER :: logger + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + iw = cp_logger_get_default_io_unit(logger) + ionode = para_env%is_source() + + ! TODO: do we really need to allocate the arrays on every process? + ALLOCATE (eri_aa(active_space_env%nmo_active**4)) + CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1) + IF (active_space_env%nspins == 2) THEN + ALLOCATE (eri_ab(active_space_env%nmo_active**4)) + CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2) + ALLOCATE (eri_bb(active_space_env%nmo_active**4)) + CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2) + END IF + + ! ask the status of the client + IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len) + DO + header = "" + CALL para_env%sync() + IF (ionode) THEN + ! IF (iw > 0) WRITE(iw, *) "@SERVER: Waiting for messages..." + CALL readbuffer(client_fd, header, header_len) + END IF + CALL para_env%bcast(header, para_env%source) + + ! IF (iw > 0) WRITE(iw, *) "@SERVER: Message from client: ", TRIM(header) + + IF (TRIM(header) == "READY") THEN + ! if the client is ready, send the data + CALL para_env%sync() + IF (ionode) THEN + CALL writebuffer(client_fd, "TWOBODY ", header_len) + CALL writebuffer(client_fd, active_space_env%nspins) + CALL writebuffer(client_fd, active_space_env%nmo_active) + CALL writebuffer(client_fd, active_space_env%nelec_active) + CALL writebuffer(client_fd, active_space_env%multiplicity) + ! send the alpha component + CALL writebuffer(client_fd, eri_aa, SIZE(eri_aa)) + ! send the beta part for unrestricted calculations + IF (active_space_env%nspins == 2) THEN + CALL writebuffer(client_fd, eri_ab, SIZE(eri_ab)) + CALL writebuffer(client_fd, eri_bb, SIZE(eri_bb)) END IF - END DO ! MMM - END DO ! LLL - END DO ! ispin + END IF + ELSE IF (TRIM(header) == "RECEIVED") THEN + EXIT + END IF + END DO - ! Release memory - DEALLOCATE (nmos) - DEALLOCATE (p_act) - IF (nspins .EQ. 2) DEALLOCATE (p_act_beta) + DEALLOCATE (eri_aa) + IF (active_space_env%nspins == 2) THEN + DEALLOCATE (eri_ab) + DEALLOCATE (eri_bb) + END IF - END SUBROUTINE read_active_density + CALL para_env%sync() + + CALL timestop(handle) + + END SUBROUTINE send_eri_to_client + +! ************************************************************************************************** +!> \brief Sends the one-electron embedding potential and the inactive energy to the client +!> \param client_fd the client file descriptor +!> \param active_space_env active space environment +!> \param para_env parallel environment +! ************************************************************************************************** + SUBROUTINE send_fock_to_client(client_fd, active_space_env, para_env) + INTEGER, INTENT(IN) :: client_fd + TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env + TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env + + CHARACTER(len=*), PARAMETER :: routineN = 'send_fock_to_client' + INTEGER, PARAMETER :: header_len = 12 + + CHARACTER(len=default_string_length) :: header + INTEGER :: handle, iw + LOGICAL :: debug = .FALSE., ionode + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fock_a, fock_b, p_act_mo_a, p_act_mo_b + TYPE(cp_logger_type), POINTER :: logger + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + iw = cp_logger_get_default_io_unit(logger) + ionode = para_env%is_source() + + ALLOCATE (p_act_mo_a(active_space_env%nmo_active**2)) + ALLOCATE (fock_a(active_space_env%nmo_active**2)) + IF (active_space_env%nspins == 2) THEN + ALLOCATE (p_act_mo_b(active_space_env%nmo_active**2)) + ALLOCATE (fock_b(active_space_env%nmo_active**2)) + END IF + + ! get the fock matrix into Fortran arrays + ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 1)) + CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices) + END ASSOCIATE + + IF (active_space_env%nspins == 2) THEN + ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 2)) + CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b, act_indices, act_indices) + END ASSOCIATE + END IF + + ! ask the status of the client + IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len) + DO + header = "" + + CALL para_env%sync() + IF (ionode) THEN + IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Waiting for messages..." + CALL readbuffer(client_fd, header, header_len) + END IF + CALL para_env%bcast(header, para_env%source) + + IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Message from client: ", TRIM(header) + + IF (TRIM(header) == "READY") THEN + ! if the client is ready, send the data + CALL para_env%sync() + IF (ionode) THEN + CALL writebuffer(client_fd, "ONEBODY ", header_len) + CALL writebuffer(client_fd, active_space_env%energy_inactive) + ! send the alpha component + CALL writebuffer(client_fd, fock_a, SIZE(fock_a)) + ! send the beta part for unrestricted calculations + IF (active_space_env%nspins == 2) THEN + CALL writebuffer(client_fd, fock_b, SIZE(fock_b)) + END IF + END IF + + ELSE IF (TRIM(header) == "HAVEDATA") THEN + ! qiskit has data to transfer, let them know we want it and wait for it + CALL para_env%sync() + IF (ionode) THEN + IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Qiskit has data to transfer" + CALL writebuffer(client_fd, "GETDENSITY ", header_len) + + ! read the active energy and density + CALL readbuffer(client_fd, active_space_env%energy_active) + CALL readbuffer(client_fd, p_act_mo_a, SIZE(p_act_mo_a)) + IF (active_space_env%nspins == 2) THEN + CALL readbuffer(client_fd, p_act_mo_b, SIZE(p_act_mo_b)) + END IF + END IF + + ! broadcast the data to all processors + CALL para_env%bcast(active_space_env%energy_active, para_env%source) + CALL para_env%bcast(p_act_mo_a, para_env%source) + IF (active_space_env%nspins == 2) THEN + CALL para_env%bcast(p_act_mo_b, para_env%source) + END IF + + ! update total and reference energies in active space enviornment + active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active + + ! update the active density matrix in the active space environment + CALL update_active_density(p_act_mo_a, active_space_env, 1) + IF (active_space_env%nspins == 2) THEN + CALL update_active_density(p_act_mo_b, active_space_env, 2) + END IF + + ! the non-iterative part is done, we can continue + EXIT + END IF + + END DO + + DEALLOCATE (p_act_mo_a) + DEALLOCATE (fock_a) + IF (active_space_env%nspins == 2) THEN + DEALLOCATE (p_act_mo_b) + DEALLOCATE (fock_b) + END IF + + CALL para_env%sync() + + CALL timestop(handle) + + END SUBROUTINE send_fock_to_client +#endif END MODULE qs_active_space_methods diff --git a/src/qs_active_space_types.F b/src/qs_active_space_types.F index 5dad65f69a..fa1fed00de 100644 --- a/src/qs_active_space_types.F +++ b/src/qs_active_space_types.F @@ -73,31 +73,34 @@ MODULE qs_active_space_types END TYPE eri_type_eri_element_func TYPE active_space_type - INTEGER :: nelec_active = 0 - INTEGER, POINTER, DIMENSION(:, :) :: active_orbitals => NULL() - INTEGER, POINTER, DIMENSION(:, :) :: inactive_orbitals => NULL() - INTEGER :: nelec_inactive = 0 - INTEGER, DIMENSION(2) :: nelec_inactive_spinwise = 0 - INTEGER, DIMENSION(2) :: nelec_total = 0 - INTEGER :: multiplicity = 0 - INTEGER :: nspins = 0 - LOGICAL :: molecule = .FALSE. - INTEGER :: model = 0 - REAL(KIND=dp) :: energy_total = 0.0_dp - REAL(KIND=dp) :: energy_ref = 0.0_dp - REAL(KIND=dp) :: energy_inactive = 0.0_dp - REAL(KIND=dp) :: energy_active = 0.0_dp - LOGICAL :: read_p_act = .FALSE. - CHARACTER(LEN=default_path_length) :: p_act_filename = "" - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active => NULL() - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_inactive => NULL() - TYPE(eri_type) :: eri = eri_type() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: p_active => NULL() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: ks_sub => NULL() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: vxc_sub => NULL() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: h_sub => NULL() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: fock_sub => NULL() - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat_inactive => NULL() + INTEGER :: nelec_active = 0 + INTEGER :: nelec_inactive = 0 + INTEGER :: nelec_total = 0 + INTEGER, POINTER, DIMENSION(:, :) :: active_orbitals => NULL() + INTEGER, POINTER, DIMENSION(:, :) :: inactive_orbitals => NULL() + INTEGER :: nmo_active = 0 + INTEGER :: nmo_inactive = 0 + INTEGER :: multiplicity = 0 + INTEGER :: nspins = 0 + LOGICAL :: molecule = .FALSE. + INTEGER :: model = 0 + REAL(KIND=dp) :: energy_total = 0.0_dp + REAL(KIND=dp) :: energy_ref = 0.0_dp + REAL(KIND=dp) :: energy_inactive = 0.0_dp + REAL(KIND=dp) :: energy_active = 0.0_dp + LOGICAL :: do_scf_embedding = .FALSE. + LOGICAL :: qcschema = .FALSE. + LOGICAL :: fcidump = .FALSE. + CHARACTER(LEN=default_path_length) :: qcschema_filename = '' + TYPE(eri_type) :: eri = eri_type() + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active => NULL() + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_inactive => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: p_active => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: ks_sub => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: vxc_sub => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: h_sub => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: fock_sub => NULL() + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat_inactive => NULL() END TYPE active_space_type ABSTRACT INTERFACE diff --git a/src/qs_active_space_utils.F b/src/qs_active_space_utils.F new file mode 100644 index 0000000000..2e3a78bcd3 --- /dev/null +++ b/src/qs_active_space_utils.F @@ -0,0 +1,158 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2023 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Contains utility routines for the active space module +!> \par History +!> 04.2023 created [SB] +!> \author SB +! ************************************************************************************************** +MODULE qs_active_space_utils + + USE cp_fm_types, ONLY: cp_fm_get_element,& + cp_fm_get_info,& + cp_fm_type + USE dbcsr_api, ONLY: dbcsr_csr_type + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_comm_type + USE qs_active_space_types, ONLY: csr_idx_from_combined,& + csr_idx_to_combined,& + eri_type,& + get_irange_csr +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_utils' + + PUBLIC :: subspace_matrix_to_array, eri_to_array + +CONTAINS + +! ************************************************************************************************** +!> \brief Copy a (square portion) of a `cp_fm_type` matrix to a standard 1D Fortran array +!> \param source_matrix the matrix from where the data is taken +!> \param target_array the array were the data is copied to +!> \param row_index a list containing the row subspace indices +!> \param col_index a list containing the column subspace indices +! ************************************************************************************************** + SUBROUTINE subspace_matrix_to_array(source_matrix, target_array, row_index, col_index) + TYPE(cp_fm_type), INTENT(IN) :: source_matrix + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: target_array + INTEGER, DIMENSION(:), INTENT(IN) :: row_index, col_index + + INTEGER :: i, i_sub, j, j_sub, max_col, max_row, & + ncols, nrows + REAL(KIND=dp) :: mval + + CALL cp_fm_get_info(source_matrix, nrow_global=max_row, ncol_global=max_col) + nrows = SIZE(row_index) + ncols = SIZE(col_index) + + CPASSERT(MAXVAL(row_index) <= max_row) + CPASSERT(MAXVAL(col_index) <= max_col) + CPASSERT(MINVAL(row_index) > 0) + CPASSERT(MINVAL(col_index) > 0) + CPASSERT(nrows <= max_row) + CPASSERT(ncols <= max_col) + + CPASSERT(SIZE(target_array) == nrows*ncols) + + DO j = 1, ncols + j_sub = col_index(j) + DO i = 1, nrows + i_sub = row_index(i) + CALL cp_fm_get_element(source_matrix, i_sub, j_sub, mval) + target_array(i + (j - 1)*nrows) = mval + END DO + END DO + END SUBROUTINE subspace_matrix_to_array + +! ************************************************************************************************** +!> \brief Copy the eri tensor for spins isp1 and isp2 to a standard 1D Fortran array +!> \param eri_env the eri environment +!> \param array the 1D Fortran array where the eri are copied to +!> \param active_orbitals a list containing the active orbitals indices +!> \param spin1 the spin of the bra +!> \param spin2 the spin of the ket +! ************************************************************************************************** + SUBROUTINE eri_to_array(eri_env, array, active_orbitals, spin1, spin2) + TYPE(eri_type), INTENT(IN) :: eri_env + REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: array + INTEGER, DIMENSION(:, :), INTENT(IN) :: active_orbitals + INTEGER, INTENT(IN) :: spin1, spin2 + + INTEGER :: i, i1, i12, i12l, i2, i3, i34, i34l, i4, & + ijkl, ijlk, irptr, j, jikl, jilk, k, & + klij, klji, l, lkij, lkji, nindex, & + nmo_active, nmo_max + INTEGER, DIMENSION(2) :: irange + REAL(KIND=dp) :: erival + TYPE(dbcsr_csr_type), POINTER :: eri + TYPE(mp_comm_type) :: mp_group + + nmo_active = SIZE(active_orbitals, 1) + nmo_max = eri_env%norb + nindex = (nmo_max*(nmo_max + 1))/2 + IF (spin1 == 1 .AND. spin2 == 1) THEN + eri => eri_env%eri(1)%csr_mat + ELSE IF ((spin1 == 1 .AND. spin2 == 2) .OR. (spin1 == 2 .AND. spin2 == 1)) THEN + eri => eri_env%eri(2)%csr_mat + ELSE + eri => eri_env%eri(3)%csr_mat + END IF + + CALL mp_group%set_handle(eri%mp_group%get_handle()) + irange = get_irange_csr(nindex, mp_group) + + array = 0.0_dp + + DO i = 1, nmo_active + i1 = active_orbitals(i, spin1) + DO j = i, nmo_active + i2 = active_orbitals(j, spin1) + i12 = csr_idx_to_combined(i1, i2, nmo_max) + IF (i12 >= irange(1) .AND. i12 <= irange(2)) THEN + i12l = i12 - irange(1) + 1 + irptr = eri%rowptr_local(i12l) - 1 + DO i34l = 1, eri%nzerow_local(i12l) + i34 = eri%colind_local(irptr + i34l) + CALL csr_idx_from_combined(i34, nmo_max, i3, i4) + k = FINDLOC(active_orbitals(:, spin2), i3, dim=1) + l = FINDLOC(active_orbitals(:, spin2), i4, dim=1) + erival = eri%nzval_local%r_dp(irptr + i34l) + + ! 8-fold permutational symmetry + ijkl = i + (j - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3 + jikl = j + (i - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3 + ijlk = i + (j - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3 + jilk = j + (i - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3 + array(ijkl) = erival + array(jikl) = erival + array(ijlk) = erival + array(jilk) = erival + IF (spin1 == spin2) THEN + klij = k + (l - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3 + lkij = l + (k - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3 + klji = k + (l - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3 + lkji = l + (k - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3 + array(klij) = erival + array(lkij) = erival + array(klji) = erival + array(lkji) = erival + END IF + END DO + END IF + END DO + END DO + CALL mp_group%sum(array) + + END SUBROUTINE eri_to_array + +END MODULE qs_active_space_utils diff --git a/src/qs_energy.F b/src/qs_energy.F index 0b207cc0fc..2861dd2b8f 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -24,6 +24,7 @@ MODULE qs_energy section_vals_val_get USE lri_environment_methods, ONLY: lri_print_stat USE mp2, ONLY: mp2_main + USE qs_active_space_methods, ONLY: active_space_main USE qs_energy_init, ONLY: qs_energies_init USE qs_energy_types, ONLY: qs_energy_type USE qs_energy_utils, ONLY: qs_energies_properties @@ -127,6 +128,9 @@ CONTAINS END IF + ! Do active space calculation + CALL active_space_main(qs_env) + ! Check for energy correction IF (qs_env%energy_correction) THEN CALL energy_correction(qs_env, ec_init=.TRUE., calculate_forces=.FALSE.) diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index c13516ddb2..69f450eb78 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -126,7 +126,6 @@ MODULE qs_scf_post_gpw REALSPACE,& RECIPROCALSPACE,& pw_type - USE qs_active_space_methods, ONLY: active_space_main USE qs_chargemol, ONLY: write_wfx USE qs_collocate_density, ONLY: calculate_wavefunction USE qs_commutators, ONLY: build_com_hr_matrix @@ -822,9 +821,6 @@ CONTAINS ! Use Wannier90 interface CALL wannier90_interface(input, logger, qs_env) - ! Do active space calculation - CALL active_space_main(input, logger, qs_env) - IF (my_do_mp2) THEN ! Get everything back DO ispin = 1, dft_control%nspins diff --git a/src/sockets.c b/src/sockets.c index 26f425a271..00e7e28555 100644 --- a/src/sockets.c +++ b/src/sockets.c @@ -36,7 +36,7 @@ * support sockets natively. * \author Joshua More and Michele Ceriotti ******************************************************************************/ -#ifndef __NO_IPI_DRIVER +#ifndef __NO_SOCKETS #define _POSIX_C_SOURCE 200809L @@ -46,6 +46,7 @@ #include #include #include +#include #include #include #include @@ -53,7 +54,7 @@ #include /******************************************************************************* - * \brief Opens a socket. + * \brief Opens and connects 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. @@ -64,7 +65,7 @@ * \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) { +void open_connect_socket(int *psockfd, int *inet, int *port, char *host) { int sockfd, ai_err; if (*inet > 0) { // creates an internet socket @@ -104,8 +105,8 @@ void open_socket(int *psockfd, int *inet, int *port, char *host) { // 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); + strcpy(serv_addr.sun_path, "/tmp/qiskit_"); + strcpy(serv_addr.sun_path + 12, host); // creates the socket sockfd = socket(AF_UNIX, SOCK_STREAM, 0); @@ -121,6 +122,76 @@ void open_socket(int *psockfd, int *inet, int *port, char *host) { *psockfd = sockfd; } +/******************************************************************************* + * \brief Opens and binds 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_bind_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); + } + + // binds + if (bind(sockfd, res->ai_addr, res->ai_addrlen) < 0) { + perror("Error binding 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, host); + + // creates the socket + sockfd = socket(AF_UNIX, SOCK_STREAM, 0); + + remove(serv_addr.sun_path); + + // binds + if (bind(sockfd, (struct sockaddr *)&serv_addr, sizeof(serv_addr)) < 0) { + perror( + "Error binding 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. @@ -163,6 +234,44 @@ void readbuffer(int *psockfd, char *data, int *plen) { } } +/******************************************************************************* + * \brief Listens to a socket. + * \param psockfd The id of the socket to listen. + * \param n An integer that determines the number of requests that will + * be queued before further requests are refused. + ******************************************************************************/ +void listen_socket(int *psockfd, int *backlog) { + + if (listen(*psockfd, *backlog) < 0) { + perror("Error listening socket"); + exit(-1); + }; +} + +/******************************************************************************* + * \brief Listens to a socket. + * \param psockfd The id of the socket to listen. + * \param pclientfd The id of the accepted socket. + ******************************************************************************/ +void accept_socket(int *psockfd, int *pclientfd) { + + int client_fd = accept(*psockfd, NULL, NULL); + + *pclientfd = client_fd; +} + +/******************************************************************************* + * \brief Closes a socket. + * \param psockfd The id of the socket to close. + ******************************************************************************/ +void close_socket(int *psockfd) { close(*psockfd); } + +/******************************************************************************* + * \brief Removes a socket file. + * \param hostname The name of the socket file to remove. + ******************************************************************************/ +void remove_socket_file(char *host) { remove(host); } + /******************************************************************************* * \brief Mini-wrapper to nanosleep * \param dsec number of seconds to wait (float values accepted) diff --git a/src/sockets_interface.F b/src/sockets_interface.F new file mode 100644 index 0000000000..95e54b2742 --- /dev/null +++ b/src/sockets_interface.F @@ -0,0 +1,304 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2023 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Implements UNIX and INET sockets +!> \par History +!> 08.2023 moved here and expanded for AS module by S. Battaglia +!> 03.2012 created by MC in ipi_driver.F +!> \author M. Ceriotti +! ************************************************************************************************** + +MODULE sockets_interface + + USE ISO_C_BINDING, ONLY: C_CHAR,& + C_DOUBLE,& + C_INT,& + C_LOC,& + C_PTR + USE kinds, ONLY: dp +#include "./base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sockets_interface' + +#ifndef __NO_SOCKETS + PUBLIC :: writebuffer, readbuffer, open_connect_socket, & + uwait, open_bind_socket, listen_socket, & + accept_socket, close_socket, remove_socket_file + + INTERFACE writebuffer + MODULE PROCEDURE writebuffer_s, & + writebuffer_d, writebuffer_dv, & + writebuffer_i + + END INTERFACE + + INTERFACE readbuffer + MODULE PROCEDURE readbuffer_s, & + readbuffer_dv, readbuffer_d, & + readbuffer_i + + END INTERFACE + + INTERFACE + SUBROUTINE uwait(sec) BIND(C, NAME="uwait") + USE ISO_C_BINDING, ONLY: C_DOUBLE + REAL(C_DOUBLE) :: sec + + END SUBROUTINE + END INTERFACE + + INTERFACE + SUBROUTINE open_connect_socket(psockfd, inet, port, host) BIND(C) + IMPORT + INTEGER(KIND=C_INT) :: psockfd, inet, port + CHARACTER(KIND=C_CHAR), DIMENSION(*) :: host + + END SUBROUTINE open_connect_socket + + SUBROUTINE open_bind_socket(psockfd, inet, port, host) BIND(C) + IMPORT + INTEGER(KIND=C_INT) :: psockfd, inet, port + CHARACTER(KIND=C_CHAR), DIMENSION(*) :: host + + END SUBROUTINE open_bind_socket + + SUBROUTINE listen_socket(psockfd, backlog) BIND(C) + IMPORT + INTEGER(KIND=C_INT) :: psockfd, backlog + + END SUBROUTINE listen_socket + + SUBROUTINE accept_socket(psockfd, pclientfd) BIND(C) + IMPORT + INTEGER(KIND=C_INT) :: psockfd, pclientfd + + END SUBROUTINE accept_socket + + SUBROUTINE close_socket(psockfd) BIND(C) + IMPORT + INTEGER(KIND=C_INT) :: psockfd + + END SUBROUTINE close_socket + + SUBROUTINE remove_socket_file(host) BIND(C) + IMPORT + CHARACTER(KIND=C_CHAR), DIMENSION(*) :: host + + END SUBROUTINE remove_socket_file + + SUBROUTINE writebuffer_csocket(psockfd, pdata, plen) BIND(C, name="writebuffer") + IMPORT + INTEGER(KIND=C_INT) :: psockfd + TYPE(C_PTR), VALUE :: pdata + INTEGER(KIND=C_INT) :: plen + + END SUBROUTINE writebuffer_csocket + + SUBROUTINE readbuffer_csocket(psockfd, pdata, plen) BIND(C, name="readbuffer") + IMPORT + INTEGER(KIND=C_INT) :: psockfd + TYPE(C_PTR), VALUE :: pdata + INTEGER(KIND=C_INT) :: plen + + END SUBROUTINE readbuffer_csocket + END INTERFACE +#endif + +CONTAINS + +#ifndef __NO_SOCKETS +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +! ************************************************************************************************** + SUBROUTINE writebuffer_d(psockfd, fdata) + INTEGER, INTENT(IN) :: psockfd + REAL(KIND=dp), INTENT(IN) :: fdata + + CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_d' + + INTEGER :: handle + REAL(KIND=C_DOUBLE), TARGET :: cdata + + CALL timeset(routineN, handle) + + cdata = fdata + CALL writebuffer_csocket(psockfd, c_loc(cdata), 8) + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +! ************************************************************************************************** + SUBROUTINE writebuffer_i(psockfd, fdata) + INTEGER, INTENT(IN) :: psockfd, fdata + + CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_i' + + INTEGER :: handle + INTEGER(KIND=C_INT), TARGET :: cdata + + CALL timeset(routineN, handle) + + cdata = fdata + CALL writebuffer_csocket(psockfd, c_loc(cdata), 4) + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fstring ... +!> \param plen ... +! ************************************************************************************************** + SUBROUTINE writebuffer_s(psockfd, fstring, plen) + INTEGER, INTENT(IN) :: psockfd + CHARACTER(LEN=*), INTENT(IN) :: fstring + INTEGER, INTENT(IN) :: plen + + CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_s' + + INTEGER :: handle, i + CHARACTER(LEN=1, KIND=C_CHAR), TARGET :: cstring(plen) + + CALL timeset(routineN, handle) + + DO i = 1, plen + cstring(i) = fstring(i:i) + END DO + CALL writebuffer_csocket(psockfd, c_loc(cstring(1)), plen) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +!> \param plen ... +! ************************************************************************************************** + SUBROUTINE writebuffer_dv(psockfd, fdata, plen) + INTEGER, INTENT(IN) :: psockfd, plen + REAL(KIND=dp), INTENT(IN), TARGET :: fdata(plen) + + CHARACTER(len=*), PARAMETER :: routineN = 'writebuffer_dv' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL writebuffer_csocket(psockfd, c_loc(fdata(1)), 8*plen) + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +! ************************************************************************************************** + SUBROUTINE readbuffer_d(psockfd, fdata) + INTEGER, INTENT(IN) :: psockfd + REAL(KIND=dp), INTENT(OUT) :: fdata + + CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_d' + + INTEGER :: handle + REAL(KIND=C_DOUBLE), TARGET :: cdata + + CALL timeset(routineN, handle) + + CALL readbuffer_csocket(psockfd, c_loc(cdata), 8) + fdata = cdata + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +! ************************************************************************************************** + SUBROUTINE readbuffer_i(psockfd, fdata) + INTEGER, INTENT(IN) :: psockfd + INTEGER, INTENT(OUT) :: fdata + + CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_i' + + INTEGER :: handle + INTEGER(KIND=C_INT), TARGET :: cdata + + CALL timeset(routineN, handle) + + CALL readbuffer_csocket(psockfd, c_loc(cdata), 4) + fdata = cdata + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fstring ... +!> \param plen ... +! ************************************************************************************************** + SUBROUTINE readbuffer_s(psockfd, fstring, plen) + INTEGER, INTENT(IN) :: psockfd + CHARACTER(LEN=*), INTENT(OUT) :: fstring + INTEGER, INTENT(IN) :: plen + + CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_s' + + INTEGER :: handle, i + CHARACTER(LEN=1, KIND=C_CHAR), TARGET :: cstring(plen) + + CALL timeset(routineN, handle) + + CALL readbuffer_csocket(psockfd, c_loc(cstring(1)), plen) + fstring = "" + DO i = 1, plen + fstring(i:i) = cstring(i) + END DO + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param psockfd ... +!> \param fdata ... +!> \param plen ... +! ************************************************************************************************** + SUBROUTINE readbuffer_dv(psockfd, fdata, plen) + INTEGER, INTENT(IN) :: psockfd, plen + REAL(KIND=dp), INTENT(OUT), TARGET :: fdata(plen) + + CHARACTER(len=*), PARAMETER :: routineN = 'readbuffer_dv' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL readbuffer_csocket(psockfd, c_loc(fdata(1)), 8*plen) + + CALL timestop(handle) + + END SUBROUTINE +#endif + +END MODULE sockets_interface diff --git a/tests/QS/regtest-as-1/TEST_FILES b/tests/QS/regtest-as-1/TEST_FILES index 5c1b0aab1b..71736a8d62 100644 --- a/tests/QS/regtest-as-1/TEST_FILES +++ b/tests/QS/regtest-as-1/TEST_FILES @@ -5,29 +5,29 @@ # for details see cp2k/tools/do_regtest # h2_gapw_2-2.inp 1 1e-12 -1.12703481947359 -h2_gapw_2-2.inp 92 1e-8 4.08762828 +h2_gapw_2-2.inp 92 1e-8 4.08763868 h2_gapw_2-3.inp 1 1e-12 -1.12703481947359 -h2_gapw_2-3.inp 92 1e-8 6.34130225 +h2_gapw_2-3.inp 92 1e-8 6.34131299 h2_gapw_2-4.inp 1 1e-12 -1.12703481947359 -h2_gapw_2-4.inp 92 1e-8 10.57997331 +h2_gapw_2-4.inp 92 1e-8 10.57998079 h2_gapw_pp_2-2.inp 1 1e-12 -1.12609343153655 -h2_gapw_pp_2-2.inp 92 1e-8 4.08439624 +h2_gapw_pp_2-2.inp 92 1e-8 4.08440577 h2_gapw_pp_2-3.inp 1 1e-12 -1.12609343153655 -h2_gapw_pp_2-3.inp 92 1e-8 6.33490647 +h2_gapw_pp_2-3.inp 92 1e-8 6.33491633 h2_gapw_pp_2-4.inp 1 1e-12 -1.12609343153655 -h2_gapw_pp_2-4.inp 92 1e-8 10.58497633 +h2_gapw_pp_2-4.inp 92 1e-8 10.58498318 h2_gpw_pp_2-2.inp 1 1e-12 -1.12622646044780 -h2_gpw_pp_2-2.inp 92 1e-8 4.08480786 +h2_gpw_pp_2-2.inp 92 1e-8 4.08481739 h2_gpw_pp_2-3.inp 1 1e-12 -1.12622646044780 -h2_gpw_pp_2-3.inp 92 1e-8 6.33542961 +h2_gpw_pp_2-3.inp 92 1e-8 6.33543946 h2_gpw_pp_2-4.inp 1 1e-12 -1.12622646044780 -h2_gpw_pp_2-4.inp 92 1e-8 10.58532520 +h2_gpw_pp_2-4.inp 92 1e-8 10.58533205 h2o_gapw_2-2.inp 1 1e-12 -76.011970410506123 -h2o_gapw_2-2.inp 92 1e-8 77.48403223 +h2o_gapw_2-2.inp 92 1e-8 77.48403318 ch2_gapw_2-3.inp 1 1e-12 -38.91914381184247 ch2_gapw_2-3.inp 92 1e-8 50.38752368 ch2_gapw_pp_2-3.inp 1 1e-12 -6.51646227467145 -ch2_gapw_pp_2-3.inp 92 1e-8 17.95445556 +ch2_gapw_pp_2-3.inp 92 1e-8 17.95445524 ch2_gpw_pp_2-3.inp 1 1e-12 -6.51682932586237 -ch2_gpw_pp_2-3.inp 92 1e-8 17.95513322 +ch2_gpw_pp_2-3.inp 92 1e-8 17.95513292 #EOF diff --git a/tests/QS/regtest-as-1/ch2_gapw_2-3.inp b/tests/QS/regtest-as-1/ch2_gapw_2-3.inp index 8ae849f0ad..dec71d907b 100644 --- a/tests/QS/regtest-as-1/ch2_gapw_2-3.inp +++ b/tests/QS/regtest-as-1/ch2_gapw_2-3.inp @@ -33,23 +33,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 3 - INACTIVE_ELECTRONS 3 3 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/ch2_gapw_pp_2-3.inp b/tests/QS/regtest-as-1/ch2_gapw_pp_2-3.inp index 9a84505bf9..0f5b854b11 100644 --- a/tests/QS/regtest-as-1/ch2_gapw_pp_2-3.inp +++ b/tests/QS/regtest-as-1/ch2_gapw_pp_2-3.inp @@ -33,23 +33,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 3 - INACTIVE_ELECTRONS 2 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/ch2_gpw_pp_2-3.inp b/tests/QS/regtest-as-1/ch2_gpw_pp_2-3.inp index 383fc81e5b..c7b9099942 100644 --- a/tests/QS/regtest-as-1/ch2_gpw_pp_2-3.inp +++ b/tests/QS/regtest-as-1/ch2_gpw_pp_2-3.inp @@ -33,23 +33,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 3 - INACTIVE_ELECTRONS 2 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_2-2.inp b/tests/QS/regtest-as-1/h2_gapw_2-2.inp index 878d933473..ff5eab102b 100644 --- a/tests/QS/regtest-as-1/h2_gapw_2-2.inp +++ b/tests/QS/regtest-as-1/h2_gapw_2-2.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_2-3.inp b/tests/QS/regtest-as-1/h2_gapw_2-3.inp index 156bf7cc8e..03436d723a 100644 --- a/tests/QS/regtest-as-1/h2_gapw_2-3.inp +++ b/tests/QS/regtest-as-1/h2_gapw_2-3.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_2-4.inp b/tests/QS/regtest-as-1/h2_gapw_2-4.inp index dfeaea5774..a78978d600 100644 --- a/tests/QS/regtest-as-1/h2_gapw_2-4.inp +++ b/tests/QS/regtest-as-1/h2_gapw_2-4.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 4 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 4 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_pp_2-2.inp b/tests/QS/regtest-as-1/h2_gapw_pp_2-2.inp index 4b23fe89c6..5b356d7889 100644 --- a/tests/QS/regtest-as-1/h2_gapw_pp_2-2.inp +++ b/tests/QS/regtest-as-1/h2_gapw_pp_2-2.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_pp_2-3.inp b/tests/QS/regtest-as-1/h2_gapw_pp_2-3.inp index 26f3821f4b..5481730335 100644 --- a/tests/QS/regtest-as-1/h2_gapw_pp_2-3.inp +++ b/tests/QS/regtest-as-1/h2_gapw_pp_2-3.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gapw_pp_2-4.inp b/tests/QS/regtest-as-1/h2_gapw_pp_2-4.inp index c846e33ec4..35180f0a76 100644 --- a/tests/QS/regtest-as-1/h2_gapw_pp_2-4.inp +++ b/tests/QS/regtest-as-1/h2_gapw_pp_2-4.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 4 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 4 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gpw_pp_2-2.inp b/tests/QS/regtest-as-1/h2_gpw_pp_2-2.inp index b205145eab..b899b1e01c 100644 --- a/tests/QS/regtest-as-1/h2_gpw_pp_2-2.inp +++ b/tests/QS/regtest-as-1/h2_gpw_pp_2-2.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gpw_pp_2-3.inp b/tests/QS/regtest-as-1/h2_gpw_pp_2-3.inp index 2230ec9214..7637b951d5 100644 --- a/tests/QS/regtest-as-1/h2_gpw_pp_2-3.inp +++ b/tests/QS/regtest-as-1/h2_gpw_pp_2-3.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 3 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 3 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2_gpw_pp_2-4.inp b/tests/QS/regtest-as-1/h2_gpw_pp_2-4.inp index 9f599f5f6a..f33b1e348e 100644 --- a/tests/QS/regtest-as-1/h2_gpw_pp_2-4.inp +++ b/tests/QS/regtest-as-1/h2_gpw_pp_2-4.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 4 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 4 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-1/h2o_gapw_2-2.inp b/tests/QS/regtest-as-1/h2o_gapw_2-2.inp index ab59b1ad2e..51cd17839c 100644 --- a/tests/QS/regtest-as-1/h2o_gapw_2-2.inp +++ b/tests/QS/regtest-as-1/h2o_gapw_2-2.inp @@ -31,22 +31,22 @@ KINETIC_ENERGY TRUE POTENTIAL_ENERGY TRUE &END AO_MATRICES - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - METHOD FULL_GPW - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 500 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + METHOD FULL_GPW + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 500 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as-dft/h2_2-2.inp b/tests/QS/regtest-as-dft/h2_2-2.inp index 31151ce7fd..6a9ce26597 100644 --- a/tests/QS/regtest-as-dft/h2_2-2.inp +++ b/tests/QS/regtest-as-dft/h2_2-2.inp @@ -79,26 +79,24 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_iterative_density.inp b/tests/QS/regtest-as-dft/h2_iterative_density.inp index 0b48eac853..3aa702b8d6 100644 --- a/tests/QS/regtest-as-dft/h2_iterative_density.inp +++ b/tests/QS/regtest-as-dft/h2_iterative_density.inp @@ -79,28 +79,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - READ_P_ACTIVE .TRUE. - P_ACTIVE_FILE_NAME h2_iterative_density.bin - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + !READ_P_ACTIVE .TRUE. + !P_ACTIVE_FILE_NAME h2_iterative_density.bin + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_1+2.inp b/tests/QS/regtest-as-dft/h2_manual_1+2.inp index 15f3b8d91f..de88dbf59c 100644 --- a/tests/QS/regtest-as-dft/h2_manual_1+2.inp +++ b/tests/QS/regtest-as-dft/h2_manual_1+2.inp @@ -79,28 +79,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 2 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 2 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_1+2_1+3.inp b/tests/QS/regtest-as-dft/h2_manual_1+2_1+3.inp index 629c92c55e..78b4de764f 100644 --- a/tests/QS/regtest-as-dft/h2_manual_1+2_1+3.inp +++ b/tests/QS/regtest-as-dft/h2_manual_1+2_1+3.inp @@ -81,29 +81,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 2 - INACTIVE_ELECTRONS 0 0 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 2 1 3 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 2 1 3 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_1+3.inp b/tests/QS/regtest-as-dft/h2_manual_1+3.inp index 6d5bec9521..8d80593387 100644 --- a/tests/QS/regtest-as-dft/h2_manual_1+3.inp +++ b/tests/QS/regtest-as-dft/h2_manual_1+3.inp @@ -79,28 +79,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 3 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 3 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_1+3_1+2.inp b/tests/QS/regtest-as-dft/h2_manual_1+3_1+2.inp index a73ef407ae..3fc8364cde 100644 --- a/tests/QS/regtest-as-dft/h2_manual_1+3_1+2.inp +++ b/tests/QS/regtest-as-dft/h2_manual_1+3_1+2.inp @@ -81,29 +81,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 2 - INACTIVE_ELECTRONS 0 0 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 3 1 2 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 3 1 2 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_1+4.inp b/tests/QS/regtest-as-dft/h2_manual_1+4.inp index 50d9fb8ddf..28101e7631 100644 --- a/tests/QS/regtest-as-dft/h2_manual_1+4.inp +++ b/tests/QS/regtest-as-dft/h2_manual_1+4.inp @@ -79,28 +79,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 4 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 4 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as-dft/h2_manual_UKS_1+2.inp b/tests/QS/regtest-as-dft/h2_manual_UKS_1+2.inp index 0a8b4a0628..73255279fb 100644 --- a/tests/QS/regtest-as-dft/h2_manual_UKS_1+2.inp +++ b/tests/QS/regtest-as-dft/h2_manual_UKS_1+2.inp @@ -81,29 +81,26 @@ &END PRINT &END SCF - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 2 - ACTIVE_ORBITALS 2 2 - INACTIVE_ELECTRONS 0 0 - ISOLATED_SYSTEM TRUE - ORBITAL_SELECTION MANUAL - ACTIVE_ORBITAL_INDICES 1 2 - &ERI - EPS_INTEGRAL 1E-10 - METHOD FULL_GPW - OPERATOR - OPERATOR_PARAMETER 1.0 - PERIODICITY 0 0 0 - &END ERI - &ERI_GPW - CUTOFF 250 - &END ERI_GPW - &FCIDUMP - FILENAME __STD_OUT__ - &END FCIDUMP - &END ACTIVE_SPACE - &END PRINT + &ACTIVE_SPACE + ACTIVE_ELECTRONS 2 + ACTIVE_ORBITALS 2 + ISOLATED_SYSTEM TRUE + ORBITAL_SELECTION MANUAL + ACTIVE_ORBITAL_INDICES 1 2 1 2 + &ERI + EPS_INTEGRAL 1E-10 + METHOD FULL_GPW + OPERATOR + OPERATOR_PARAMETER 1.0 + PERIODICITY 0 0 0 + &END ERI + &ERI_GPW + CUTOFF 250 + &END ERI_GPW + &FCIDUMP + FILENAME __STD_OUT__ + &END FCIDUMP + &END ACTIVE_SPACE &END DFT &END FORCE_EVAL diff --git a/tests/QS/regtest-as/be.inp b/tests/QS/regtest-as/be.inp index b9942a5551..33ce19fef6 100644 --- a/tests/QS/regtest-as/be.inp +++ b/tests/QS/regtest-as/be.inp @@ -29,26 +29,24 @@ &END &END XC &POISSON - PERIODIC NONE - POISSON_SOLVER ANALYTIC - &END - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 4 - ACTIVE_ORBITALS 3 - ORBITAL_SELECTION CANONICAL - &ERI - METHOD FULL_GPW - OPERATOR <1/R> - &END - &ERI_GPW - STORE_WFN T - PRINT_LEVEL LOW - &END - &FCIDUMP - &END - &END ACTIVE_SPACE + PERIODIC NONE + POISSON_SOLVER ANALYTIC &END + &ACTIVE_SPACE + ACTIVE_ELECTRONS 4 + ACTIVE_ORBITALS 3 + ORBITAL_SELECTION CANONICAL + &ERI + METHOD FULL_GPW + OPERATOR <1/R> + &END + &ERI_GPW + STORE_WFN T + PRINT_LEVEL LOW + &END + &FCIDUMP + &END + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tests/QS/regtest-as/h2o.inp b/tests/QS/regtest-as/h2o.inp index eea0614e34..e7043c4a1a 100644 --- a/tests/QS/regtest-as/h2o.inp +++ b/tests/QS/regtest-as/h2o.inp @@ -16,29 +16,27 @@ &XC_FUNCTIONAL Pade &END XC_FUNCTIONAL &END XC - &PRINT - &ACTIVE_SPACE - ACTIVE_ELECTRONS 4 - ACTIVE_ORBITALS 6 - ORBITAL_SELECTION CANONICAL - &ERI - METHOD FULL_GPW - OPERATOR <1/R> - &END - &ERI_GPW - STORE_WFN T - PRINT_LEVEL MEDIUM - &END - &PRINT_ORBITAL_CUBES - ALIST 3 4 1 - BLIST - STRIDE 3 - FILENAME Active_orbital - &END - &FCIDUMP - &END - &END ACTIVE_SPACE - &END + &ACTIVE_SPACE + ACTIVE_ELECTRONS 4 + ACTIVE_ORBITALS 6 + ORBITAL_SELECTION CANONICAL + &ERI + METHOD FULL_GPW + OPERATOR <1/R> + &END + &ERI_GPW + STORE_WFN T + PRINT_LEVEL MEDIUM + &END + &PRINT_ORBITAL_CUBES + ALIST 3 4 1 + BLIST + STRIDE 3 + FILENAME Active_orbital + &END + &FCIDUMP + &END + &END ACTIVE_SPACE &END DFT &SUBSYS &CELL diff --git a/tools/toolchain/scripts/stage7/install_hdf5.sh b/tools/toolchain/scripts/stage7/install_hdf5.sh index 9907c97fdd..9abdaff853 100755 --- a/tools/toolchain/scripts/stage7/install_hdf5.sh +++ b/tools/toolchain/scripts/stage7/install_hdf5.sh @@ -40,7 +40,7 @@ case "$with_hdf5" in ./configure \ --prefix="${pkg_install_dir}" \ --libdir="${pkg_install_dir}/lib" \ - --disable-shared \ + --enable-fortran \ > configure.log 2>&1 || tail -n ${LOG_LINES} configure.log make -j $(get_nprocs) > make.log 2>&1 || tail -n ${LOG_LINES} make.log make -j $(get_nprocs) install > install.log 2>&1 || tail -n ${LOG_LINES} install.log @@ -63,7 +63,7 @@ case "$with_hdf5" in esac if [ "$with_hdf5" != "__DONTUSE__" ]; then - HDF5_LIBS="-lhdf5 -lhdf5_hl -lz" + HDF5_LIBS="-lhdf5 -lhdf5_hl -lhdf5_fortran -lz" if [ "$with_hdf5" != "__SYSTEM__" ]; then cat << EOF > "${BUILDDIR}/setup_hdf5" prepend_path LD_LIBRARY_PATH "$pkg_install_dir/lib" @@ -77,17 +77,10 @@ EOF cat << EOF >> "${BUILDDIR}/setup_hdf5" export HDF5_CFLAGS="${HDF5_CFLAGS}" export HDF5_LDFLAGS="${HDF5_LDFLAGS}" -export CP_DFLAGS="\${CP_DFLAGS} IF_MPI(-D__HDF5|)" +export CP_DFLAGS="\${CP_DFLAGS} -D__HDF5" export CP_CFLAGS="\${CP_CFLAGS} ${HDF5_CFLAGS}" export CP_LDFLAGS="\${CP_LDFLAGS} ${HDF5_LDFLAGS}" -#################################################### -# -# include hdf5 only if sirius is activated and build -# depends them on mpi -# -#################################################### - -export CP_LIBS="IF_MPI(${HDF5_LIBS}|) \${CP_LIBS}" +export CP_LIBS="${HDF5_LIBS} \${CP_LIBS}" export HDF5_ROOT="$pkg_install_dir" export HDF5_LIBRARIES="$HDF5_LIBS" export HDF5_HL_LIBRARIES="$HDF5_LIBS"