1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \author T.Laino 8! ************************************************************************************************** 9MODULE qs_ks_qmmm_methods 10 USE cp_control_types, ONLY: dft_control_type 11 USE cp_log_handling, ONLY: cp_get_default_logger,& 12 cp_logger_type 13 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 14 cp_print_key_unit_nr 15 USE cube_utils, ONLY: cube_info_type,& 16 init_cube_info 17 USE d3_poly, ONLY: init_d3_poly_module 18 USE input_constants, ONLY: do_qmmm_gauss,& 19 do_qmmm_swave 20 USE input_section_types, ONLY: section_vals_type 21 USE kinds, ONLY: dp 22 USE pw_env_types, ONLY: pw_env_get,& 23 pw_env_retain 24 USE pw_methods, ONLY: pw_integral_ab 25 USE pw_pool_types, ONLY: pw_pool_create_pw,& 26 pw_pool_p_type,& 27 pw_pool_type 28 USE pw_types, ONLY: REALDATA3D,& 29 REALSPACE,& 30 pw_p_type 31 USE qmmm_types_low, ONLY: qmmm_env_qm_type 32 USE qs_environment_types, ONLY: get_qs_env,& 33 qs_environment_type,& 34 set_qs_env 35 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type,& 36 qs_ks_qmmm_release 37#include "./base/base_uses.f90" 38 39 IMPLICIT NONE 40 41 PRIVATE 42 43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_qmmm_methods' 44 INTEGER, SAVE, PRIVATE :: last_ks_qmmm_nr = 0 45 46 PUBLIC :: ks_qmmm_env_rebuild, qmmm_calculate_energy, & 47 qmmm_modify_hartree_pot 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief Initialize the ks_qmmm_env 53!> \param qs_env ... 54!> \param qmmm_env ... 55!> \par History 56!> 05.2004 created [tlaino] 57!> \author Teodoro Laino 58! ************************************************************************************************** 59 SUBROUTINE ks_qmmm_env_rebuild(qs_env, qmmm_env) 60 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env 61 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 62 63 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env 64 65 NULLIFY (ks_qmmm_env) 66 CALL get_qs_env(qs_env=qs_env, & 67 ks_qmmm_env=ks_qmmm_env) 68 69 ! *** allocate the ks_qmmm env if not allocated yet!** 70 IF (.NOT. ASSOCIATED(ks_qmmm_env)) THEN 71 CALL qs_ks_qmmm_create(ks_qmmm_env=ks_qmmm_env, qs_env=qs_env, & 72 qmmm_env=qmmm_env) 73 CALL set_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env) 74 CALL qs_ks_qmmm_release(ks_qmmm_env=ks_qmmm_env) 75 END IF 76 END SUBROUTINE ks_qmmm_env_rebuild 77 78! ************************************************************************************************** 79!> \brief allocates and initializes the given ks_qmmm_env. 80!> \param ks_qmmm_env the ks_qmmm env to be initialized 81!> \param qs_env the qs environment 82!> \param qmmm_env ... 83!> \par History 84!> 05.2004 created [tlaino] 85!> \author Teodoro Laino 86! ************************************************************************************************** 87 SUBROUTINE qs_ks_qmmm_create(ks_qmmm_env, qs_env, qmmm_env) 88 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env 89 TYPE(qs_environment_type), POINTER :: qs_env 90 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 91 92 CHARACTER(len=*), PARAMETER :: routineN = 'qs_ks_qmmm_create', & 93 routineP = moduleN//':'//routineN 94 95 INTEGER :: handle, igrid 96 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 97 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pools 98 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 99 100 CALL timeset(routineN, handle) 101 102 CPASSERT(.NOT. ASSOCIATED(ks_qmmm_env)) 103 ALLOCATE (ks_qmmm_env) 104 105 NULLIFY (ks_qmmm_env%pw_env, & 106 ks_qmmm_env%cube_info) 107 NULLIFY (auxbas_pw_pool) 108 CALL get_qs_env(qs_env=qs_env, & 109 pw_env=ks_qmmm_env%pw_env) 110 CALL pw_env_get(ks_qmmm_env%pw_env, auxbas_pw_pool=auxbas_pw_pool) 111 CALL pw_env_retain(ks_qmmm_env%pw_env) 112 113 ks_qmmm_env%pc_ener = 0.0_dp 114 ks_qmmm_env%n_evals = 0 115 ks_qmmm_env%ref_count = 1 116 last_ks_qmmm_nr = last_ks_qmmm_nr + 1 117 ks_qmmm_env%id_nr = last_ks_qmmm_nr 118 119 CALL pw_pool_create_pw(auxbas_pw_pool, ks_qmmm_env%v_qmmm_rspace%pw, & 120 use_data=REALDATA3D, in_space=REALSPACE) 121 122 IF ((qmmm_env%qmmm_coupl_type == do_qmmm_gauss) .OR. (qmmm_env%qmmm_coupl_type == do_qmmm_swave)) THEN 123 CALL init_d3_poly_module() ! a fairly arbitrary but sufficient spot to do this 124 CALL pw_env_get(ks_qmmm_env%pw_env, pw_pools=pools) 125 ALLOCATE (cube_info(SIZE(pools))) 126 DO igrid = 1, SIZE(pools) 127 CALL init_cube_info(cube_info(igrid), & 128 pools(igrid)%pool%pw_grid%dr(:), & 129 pools(igrid)%pool%pw_grid%dh(:, :), & 130 pools(igrid)%pool%pw_grid%dh_inv(:, :), & 131 pools(igrid)%pool%pw_grid%orthorhombic, & 132 qmmm_env%maxRadius(igrid)) 133 END DO 134 ks_qmmm_env%cube_info => cube_info 135 END IF 136 NULLIFY (ks_qmmm_env%matrix_h) 137 ! 138 139 CALL timestop(handle) 140 141 END SUBROUTINE qs_ks_qmmm_create 142 143! ************************************************************************************************** 144!> \brief Computes the contribution to the total energy of the QM/MM 145!> electrostatic coupling 146!> \param qs_env ... 147!> \param rho ... 148!> \param v_qmmm ... 149!> \param qmmm_energy ... 150!> \par History 151!> 05.2004 created [tlaino] 152!> \author Teodoro Laino 153! ************************************************************************************************** 154 SUBROUTINE qmmm_calculate_energy(qs_env, rho, v_qmmm, qmmm_energy) 155 TYPE(qs_environment_type), POINTER :: qs_env 156 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho 157 TYPE(pw_p_type), INTENT(IN) :: v_qmmm 158 REAL(KIND=dp), INTENT(INOUT) :: qmmm_energy 159 160 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_calculate_energy', & 161 routineP = moduleN//':'//routineN 162 163 INTEGER :: handle, ispin, output_unit 164 TYPE(cp_logger_type), POINTER :: logger 165 TYPE(dft_control_type), POINTER :: dft_control 166 TYPE(pw_p_type), POINTER :: rho0_s_rs 167 TYPE(section_vals_type), POINTER :: input 168 169 CALL timeset(routineN, handle) 170 CPASSERT(ASSOCIATED(rho)) 171 NULLIFY (dft_control, input, logger) 172 logger => cp_get_default_logger() 173 174 CALL get_qs_env(qs_env=qs_env, & 175 dft_control=dft_control, & 176 input=input) 177 178 output_unit = cp_print_key_unit_nr(logger, input, "QMMM%PRINT%PROGRAM_RUN_INFO", & 179 extension=".qmmmLog") 180 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") & 181 "Adding QM/MM electrostatic potential to the Kohn-Sham potential." 182 183 qmmm_energy = 0.0_dp 184 DO ispin = 1, dft_control%nspins 185 qmmm_energy = qmmm_energy + pw_integral_ab(rho(ispin)%pw, v_qmmm%pw) 186 END DO 187 IF (dft_control%qs_control%gapw) THEN 188 CALL get_qs_env(qs_env=qs_env, & 189 rho0_s_rs=rho0_s_rs) 190 CPASSERT(ASSOCIATED(rho0_s_rs)) 191 qmmm_energy = qmmm_energy + pw_integral_ab(rho0_s_rs%pw, v_qmmm%pw) 192 END IF 193 194 CALL cp_print_key_finished_output(output_unit, logger, input, & 195 "QMMM%PRINT%PROGRAM_RUN_INFO") 196 197 CALL timestop(handle) 198 END SUBROUTINE qmmm_calculate_energy 199 200! ************************************************************************************************** 201!> \brief Modify the hartree potential in order to include the QM/MM correction 202!> \param v_hartree ... 203!> \param v_qmmm ... 204!> \param scale ... 205!> \par History 206!> 05.2004 created [tlaino] 207!> \author Teodoro Laino 208! ************************************************************************************************** 209 SUBROUTINE qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale) 210 TYPE(pw_p_type), INTENT(INOUT) :: v_hartree 211 TYPE(pw_p_type), INTENT(IN) :: v_qmmm 212 REAL(KIND=dp), INTENT(IN) :: scale 213 214 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_modify_hartree_pot', & 215 routineP = moduleN//':'//routineN 216 217 INTEGER :: handle 218 REAL(KIND=dp) :: fac 219 220 CALL timeset(routineN, handle) 221 222 fac = v_qmmm%pw%pw_grid%dvol*scale 223 v_hartree%pw%cr3d = v_hartree%pw%cr3d+fac*v_qmmm%pw%cr3d 224 225 CALL timestop(handle) 226 END SUBROUTINE qmmm_modify_hartree_pot 227 228END MODULE qs_ks_qmmm_methods 229