1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Initialize the use of the gaussians to treat the QMMM 8!> coupling potential 9!> \par History 10!> 6.2004 created [tlaino] 11!> \author Teodoro Laino 12! ************************************************************************************************** 13MODULE qmmm_gaussian_init 14 USE ao_util, ONLY: exp_radius 15 USE cp_log_handling, ONLY: cp_get_default_logger,& 16 cp_logger_type 17 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 18 cp_print_key_unit_nr 19 USE cp_para_types, ONLY: cp_para_env_type 20 USE gaussian_gridlevels, ONLY: gaussian_gridlevel,& 21 gridlevel_info_type 22 USE input_constants, ONLY: do_qmmm_gauss,& 23 do_qmmm_swave 24 USE input_section_types, ONLY: section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE memory_utilities, ONLY: reallocate 28 USE pw_env_types, ONLY: pw_env_get,& 29 pw_env_type 30 USE pw_pool_types, ONLY: pw_pool_p_type 31 USE qmmm_gaussian_data, ONLY: max_geep_lib_gauss,& 32 min_geep_lib_gauss 33 USE qmmm_gaussian_input, ONLY: read_mm_potential,& 34 set_mm_potential_erf,& 35 set_mm_potential_swave 36 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,& 37 qmmm_gaussian_type 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 PRIVATE 42 43 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gaussian_init' 45 PUBLIC :: qmmm_gaussian_initialize 46 47CONTAINS 48 49! ************************************************************************************************** 50!> \brief Initialize the Gaussian QMMM Environment 51!> \param qmmm_gaussian_fns ... 52!> \param para_env ... 53!> \param pw_env ... 54!> \param mm_el_pot_radius ... 55!> \param mm_el_pot_radius_corr ... 56!> \param qmmm_coupl_type ... 57!> \param eps_mm_rspace ... 58!> \param maxradius ... 59!> \param maxchrg ... 60!> \param compatibility ... 61!> \param print_section ... 62!> \param qmmm_section ... 63!> \par History 64!> 06.2004 created [tlaino] 65!> \author Teodoro Laino 66! ************************************************************************************************** 67 SUBROUTINE qmmm_gaussian_initialize(qmmm_gaussian_fns, para_env, pw_env, & 68 mm_el_pot_radius, mm_el_pot_radius_corr, & 69 qmmm_coupl_type, eps_mm_rspace, maxradius, maxchrg, compatibility, & 70 print_section, qmmm_section) 71 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: qmmm_gaussian_fns 72 TYPE(cp_para_env_type), POINTER :: para_env 73 TYPE(pw_env_type), POINTER :: pw_env 74 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr 75 INTEGER, INTENT(IN) :: qmmm_coupl_type 76 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 77 REAL(KIND=dp), DIMENSION(:), POINTER :: maxradius 78 REAL(KIND=dp), INTENT(IN) :: maxchrg 79 LOGICAL, INTENT(IN) :: compatibility 80 TYPE(section_vals_type), POINTER :: print_section, qmmm_section 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_gaussian_initialize', & 83 routineP = moduleN//':'//routineN 84 85 INTEGER :: i, ilevel, j, Ndim, num_geep_gauss, & 86 output_unit 87 LOGICAL :: Found, use_geep_lib 88 REAL(KIND=dp) :: alpha, mymaxradius, Prefactor 89 REAL(KIND=dp), DIMENSION(:), POINTER :: c_radius, radius 90 TYPE(cp_logger_type), POINTER :: logger 91 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 92 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pools 93 TYPE(qmmm_gaussian_type), POINTER :: mypgf 94 95! Statements 96 97 NULLIFY (mypgf, gridlevel_info, radius, c_radius, logger) 98 logger => cp_get_default_logger() 99 CALL section_vals_val_get(qmmm_section, "USE_GEEP_LIB", i_val=num_geep_gauss) 100 IF (num_geep_gauss == 0) THEN 101 use_geep_lib = .FALSE. 102 ELSE 103 use_geep_lib = .TRUE. 104 CPASSERT(num_geep_gauss >= min_geep_lib_gauss) 105 CPASSERT(num_geep_gauss <= max_geep_lib_gauss) 106 END IF 107 SELECT CASE (qmmm_coupl_type) 108 CASE (do_qmmm_gauss, do_qmmm_swave) 109 ! 110 ! Preprocessing... 111 ! 112 ALLOCATE (radius(1)) 113 ALLOCATE (c_radius(1)) 114 Ndim = SIZE(radius) 115 Loop_on_all_values: DO I = 1, SIZE(mm_el_pot_radius) 116 Found = .FALSE. 117 Loop_on_found_values: DO J = 1, SIZE(radius) - 1 118 IF (mm_el_pot_radius(i) .EQ. radius(j)) THEN 119 Found = .TRUE. 120 EXIT Loop_on_found_values 121 END IF 122 END DO Loop_on_found_values 123 IF (.NOT. Found) THEN 124 Ndim = SIZE(radius) 125 radius(Ndim) = mm_el_pot_radius(i) 126 c_radius(Ndim) = mm_el_pot_radius_corr(i) 127 Ndim = Ndim + 1 128 CALL REALLOCATE(radius, 1, Ndim) 129 CALL REALLOCATE(c_radius, 1, Ndim) 130 END IF 131 END DO Loop_on_all_values 132 ! 133 IF (Ndim - 1 > 0) THEN 134 CALL REALLOCATE(radius, 1, Ndim - 1) 135 CALL REALLOCATE(c_radius, 1, Ndim - 1) 136 ELSE IF (Ndim - 1 == 0) THEN 137 DEALLOCATE (radius) 138 DEALLOCATE (c_radius) 139 ELSE 140 CPABORT("") 141 END IF 142 ! 143 ALLOCATE (qmmm_gaussian_fns(Ndim - 1)) 144 DO I = 1, Ndim - 1 145 NULLIFY (qmmm_gaussian_fns(I)%pgf) 146 ALLOCATE (qmmm_gaussian_fns(I)%pgf) 147 NULLIFY (qmmm_gaussian_fns(I)%pgf%Ak) 148 NULLIFY (qmmm_gaussian_fns(I)%pgf%Gk) 149 NULLIFY (qmmm_gaussian_fns(I)%pgf%grid_level) 150 ! 151 ! Default Values 152 ! 153 qmmm_gaussian_fns(I)%pgf%Elp_Radius = radius(I) 154 qmmm_gaussian_fns(I)%pgf%Elp_Radius_corr = c_radius(I) 155 END DO 156 IF (ASSOCIATED(radius)) THEN 157 DEALLOCATE (radius) 158 END IF 159 IF (ASSOCIATED(c_radius)) THEN 160 DEALLOCATE (c_radius) 161 END IF 162 ! 163 IF (use_geep_lib) THEN 164 IF (qmmm_coupl_type == do_qmmm_gauss) THEN 165 CALL set_mm_potential_erf(qmmm_gaussian_fns, & 166 compatibility, num_geep_gauss) 167 ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN 168 CALL set_mm_potential_swave(qmmm_gaussian_fns, & 169 num_geep_gauss) 170 END IF 171 ELSE 172 CALL read_mm_potential(para_env, qmmm_gaussian_fns, & 173 (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)), qmmm_section) 174 END IF 175 ! 176 CALL pw_env_get(pw_env, pw_pools=pools, gridlevel_info=gridlevel_info) 177 ALLOCATE (maxradius(SIZE(pools))) 178 maxradius = 0.0_dp 179 DO J = 1, SIZE(qmmm_gaussian_fns) 180 mypgf => qmmm_gaussian_fns(J)%pgf 181 ALLOCATE (mypgf%grid_level(SIZE(mypgf%Ak))) 182 mypgf%grid_level = 0 183 mymaxradius = 0.0_dp 184 DO I = 1, mypgf%number_of_gaussians 185 IF (mypgf%Gk(I) /= 0.0_dp) THEN 186 alpha = 1.0_dp/mypgf%Gk(I) 187 alpha = alpha*alpha 188 ilevel = gaussian_gridlevel(gridlevel_info, alpha) 189 Prefactor = mypgf%Ak(I)*maxchrg 190 mymaxradius = exp_radius(0, alpha, eps_mm_rspace, Prefactor) 191 maxradius(ilevel) = MAX(maxradius(ilevel), mymaxradius) 192 mypgf%grid_level(i) = ilevel 193 END IF 194 END DO 195 END DO 196 ! 197 ! End of gaussian initialization... 198 CASE DEFAULT 199 output_unit = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", & 200 extension=".qmmmLog") 201 IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Gaussian Data Not Initialized!" 202 CALL cp_print_key_finished_output(output_unit, logger, print_section, "PROGRAM_RUN_INFO") 203 END SELECT 204 END SUBROUTINE qmmm_gaussian_initialize 205 206END MODULE qmmm_gaussian_init 207