diff --git a/src/ec_external.F b/src/ec_external.F index c6ae525e1f..08873fc1d3 100644 --- a/src/ec_external.F +++ b/src/ec_external.F @@ -58,12 +58,14 @@ MODULE ec_external USE qs_kind_types, ONLY: qs_kind_type USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_types, ONLY: qs_ks_env_type - USE qs_mo_types, ONLY: get_mo_set,& + USE qs_mo_types, ONLY: deallocate_mo_set,& + get_mo_set,& mo_set_type USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type USE qs_overlap, ONLY: build_overlap_matrix USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type + USE trexio_utils, ONLY: read_trexio USE virial_methods, ONLY: one_third_sum_diag USE virial_types, ONLY: virial_type #include "./base/base_uses.f90" @@ -120,23 +122,24 @@ CONTAINS CALL get_qs_env(qs_env, mos=mos) ALLOCATE (cpmos(nspins), mo_ref(nspins), mo_occ(nspins)) - DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) - NULLIFY (fm_struct) - CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & - template_fmstruct=mo_coeff%matrix_struct) - CALL cp_fm_create(cpmos(ispin), fm_struct) - CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) - CALL cp_fm_create(mo_ref(ispin), fm_struct) - CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp) - CALL cp_fm_create(mo_occ(ispin), fm_struct) - CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) - CALL cp_fm_struct_release(fm_struct) - END DO CALL cp_fm_release(ec_env%mo_occ) CALL cp_fm_release(ec_env%cpmos) IF (ec_env%debug_external) THEN + ! + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & + template_fmstruct=mo_coeff%matrix_struct) + CALL cp_fm_create(cpmos(ispin), fm_struct) + CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) + CALL cp_fm_create(mo_ref(ispin), fm_struct) + CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp) + CALL cp_fm_create(mo_occ(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) + CALL cp_fm_struct_release(fm_struct) + END DO ! ec_env%mo_occ => mo_ref CALL ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr) @@ -153,8 +156,22 @@ CONTAINS END IF ec_env%cpmos => cpmos ELSE + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & + template_fmstruct=mo_coeff%matrix_struct) + CALL cp_fm_create(cpmos(ispin), fm_struct) + CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) + CALL cp_fm_create(mo_occ(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) + CALL cp_fm_create(mo_ref(ispin), fm_struct) + CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp) + CALL cp_fm_struct_release(fm_struct) + END DO + ! get external information - CALL ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr) + CALL ec_ext_interface(qs_env, ec_env%exresp_fn, mo_occ, mo_ref, cpmos, calculate_forces, unit_nr) ec_env%mo_occ => mo_ref ec_env%cpmos => cpmos END IF @@ -186,47 +203,72 @@ CONTAINS ! ************************************************************************************************** !> \brief ... !> \param qs_env ... -!> \param ec_env ... +!> \param trexio_fn ... +!> \param mo_occ ... !> \param mo_ref ... !> \param cpmos ... !> \param calculate_forces ... !> \param unit_nr ... ! ************************************************************************************************** - SUBROUTINE ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr) + SUBROUTINE ec_ext_interface(qs_env, trexio_fn, mo_occ, mo_ref, cpmos, calculate_forces, unit_nr) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(energy_correction_type), POINTER :: ec_env - TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_ref, cpmos + CHARACTER(LEN=*) :: trexio_fn + TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_occ, mo_ref, cpmos LOGICAL, INTENT(IN) :: calculate_forces INTEGER, INTENT(IN) :: unit_nr CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_interface' - INTEGER :: handle, ispin, nao, nocc(2), nspins + INTEGER :: handle, ispin, nao, nmos, nocc(2), nspins + REAL(KIND=dp) :: focc + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: energy_derivative TYPE(dft_control_type), POINTER :: dft_control - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_trexio TYPE(mp_para_env_type), POINTER :: para_env CALL timeset(routineN, handle) - MARK_USED(qs_env) - MARK_USED(ec_env) - MARK_USED(mo_ref) - MARK_USED(cpmos) - MARK_USED(calculate_forces) - MARK_USED(unit_nr) - CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env) nspins = dft_control%nspins - CALL get_qs_env(qs_env, mos=mos) nocc = 0 DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), nao=nao, homo=nocc(ispin)) + CALL cp_fm_get_info(mo_occ(ispin), nrow_global=nao, ncol_global=nocc(ispin)) END DO IF (unit_nr > 0) THEN - WRITE (unit_nr, '(T2,A)') "Read EXTERNAL Response from file: "//TRIM(ec_env%exresp_fn) + WRITE (unit_nr, '(T2,A)') " Read EXTERNAL Response from file: "//TRIM(trexio_fn) END IF + ALLOCATE (mos_trexio(nspins)) + IF (calculate_forces) THEN + NULLIFY (energy_derivative) + CALL dbcsr_allocate_matrix_set(energy_derivative, nspins) + ! + CALL read_trexio(qs_env, trexio_filename=trexio_fn, & + mo_set_trexio=mos_trexio, & + energy_derivative=energy_derivative) + ! + focc = 2.0_dp + IF (nspins == 1) focc = 4.0_dp + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos_trexio(ispin), mo_coeff=mo_coeff, homo=nmos) + CALL cp_dbcsr_sm_fm_multiply(energy_derivative(ispin)%matrix, mo_occ(ispin), & + cpmos(ispin), ncol=nmos, alpha=focc, beta=0.0_dp) + END DO + ! + CALL dbcsr_deallocate_matrix_set(energy_derivative) + ELSE + CALL read_trexio(qs_env, trexio_filename=trexio_fn, & + mo_set_trexio=mos_trexio) + END IF + ! + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos_trexio(ispin), mo_coeff=mo_coeff, homo=nmos) + CALL cp_fm_to_fm(mo_coeff, mo_ref(ispin), nmos) + CALL deallocate_mo_set(mos_trexio(ispin)) + END DO + DEALLOCATE (mos_trexio) CALL timestop(handle) diff --git a/src/energy_corrections.F b/src/energy_corrections.F index e5b9cbcaa8..8fb07e827c 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -608,7 +608,7 @@ CONTAINS WRITE (funit, *) WRITE (funit, *) " CELL // RESPONSE PRESSURE // ERRORS " DO ia = 1, 3 - WRITE (funit, "(3F15.8,5x,3G15.8,5x,3G15.8)") cell%hmat(ia, 1:3), & + WRITE (funit, "(3F15.8,5x,3F15.8,5x,3F15.8)") cell%hmat(ia, 1:3), & ec_env%rpv(ia, 1:3), ec_env%rpverror(ia, 1:3) END DO diff --git a/src/trexio_utils.F b/src/trexio_utils.F index c1c97b0593..98876e42ed 100644 --- a/src/trexio_utils.F +++ b/src/trexio_utils.F @@ -1018,7 +1018,7 @@ CONTAINS INTEGER :: ao_num, mo_num, nmo, nspins, ispin, nsgf, & save_cartesian, i, j, k, l, m, imo, ishell, & nshell, shell_num, nucleus_num, natoms, ikind, & - iatom, nrow_global, nelectron, nrows, ncols, & + iatom, nelectron, nrows, ncols, & row, col, row_size, col_size, & row_offset, col_offset, myprint INTEGER, DIMENSION(2) :: nmo_spin, electron_num @@ -1031,8 +1031,6 @@ CONTAINS TYPE(cp_logger_type), POINTER :: logger => Null() TYPE(cp_fm_type), POINTER :: mo_coeff_ref, mo_coeff_target - TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp - TYPE(cp_blacs_env_type), POINTER :: context TYPE(mp_para_env_type), POINTER :: para_env => Null() TYPE(dft_control_type), POINTER :: dft_control => Null() TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s => Null() @@ -1270,11 +1268,7 @@ CONTAINS CALL allocate_mo_set(mo_set_trexio(ispin), nsgf, nmo, nelectron, 0.0_dp, maxocc, 0.0_dp) CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff_ref) - CALL cp_fm_get_info(mo_coeff_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=nmo) - CALL init_mo_set(mo_set_trexio(ispin), fm_struct=fm_struct_tmp, name="TREXIO MOs") - CALL cp_fm_struct_release(fm_struct_tmp) + CALL init_mo_set(mo_set_trexio(ispin), fm_ref=mo_coeff_ref, name="TREXIO MOs") CALL get_mo_set(mo_set_trexio(ispin), mo_coeff=mo_coeff_target) DO j = 1, nmo @@ -1371,12 +1365,12 @@ CONTAINS DEALLOCATE (temp) - CALL get_qs_env(qs_env, matrix_ks=matrix_s) + CALL get_qs_env(qs_env, matrix_s=matrix_s) DO ispin = 1, nspins ALLOCATE (energy_derivative(ispin)%matrix) ! we use the overlap matrix as a template, copying it but removing the sparsity - CALL dbcsr_copy(energy_derivative(ispin)%matrix, matrix_s(ispin)%matrix, & + CALL dbcsr_copy(energy_derivative(ispin)%matrix, matrix_s(1)%matrix, & name='Energy Derivative', keep_sparsity=.FALSE.) CALL dbcsr_set(energy_derivative(ispin)%matrix, 0.0_dp) @@ -1426,6 +1420,7 @@ CONTAINS MARK_USED(mo_set_trexio) MARK_USED(energy_derivative) CPWARN('TREXIO support has not been enabled in this build.') + CPABORT('TREXIO Not Available') #endif CALL timestop(handle)