1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Working with the semi empirical parameter types. 8!> \author JGH (14.08.2004) 9! ************************************************************************************************** 10MODULE semi_empirical_utils 11 USE basis_set_types, ONLY: allocate_sto_basis_set,& 12 create_gto_from_sto_basis,& 13 gto_basis_set_type,& 14 set_sto_basis_set 15 USE cell_types, ONLY: cell_type,& 16 plane_distance 17 USE cp_control_types, ONLY: semi_empirical_control_type 18 USE input_constants, ONLY: & 19 do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, & 20 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1 21 USE input_section_types, ONLY: section_vals_type,& 22 section_vals_val_get 23 USE kinds, ONLY: dp 24 USE semi_empirical_mpole_methods, ONLY: semi_empirical_mpole_p_setup 25 USE semi_empirical_par_utils, ONLY: get_se_basis,& 26 setup_1c_2el_int 27 USE semi_empirical_parameters, ONLY: & 28 am1_default_parameter, mndo_default_parameter, pcharge_default_parameter, & 29 pdg_default_parameter, pm3_default_parameter, pm6_default_parameter, & 30 pm6fm_default_parameter, pnnl_default_parameter, rm1_default_parameter 31 USE semi_empirical_types, ONLY: se_taper_type,& 32 semi_empirical_type 33#include "./base/base_uses.f90" 34 35 IMPLICIT NONE 36 37 PRIVATE 38 39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_utils' 40 41 PUBLIC :: init_se_param, se_param_set_default, get_se_type, & 42 initialize_se_taper, finalize_se_taper, se_cutoff_compatible 43 44CONTAINS 45! ************************************************************************************************** 46!> \brief Reset cutoffs trying to be somehow a bit smarter 47!> \param se_control ... 48!> \param se_section ... 49!> \param cell ... 50!> \param output_unit ... 51!> \author Teodoro Laino [tlaino] - 03.2009 52! ************************************************************************************************** 53 SUBROUTINE se_cutoff_compatible(se_control, se_section, cell, output_unit) 54 TYPE(semi_empirical_control_type), POINTER :: se_control 55 TYPE(section_vals_type), POINTER :: se_section 56 TYPE(cell_type), POINTER :: cell 57 INTEGER, INTENT(IN) :: output_unit 58 59 CHARACTER(len=*), PARAMETER :: routineN = 'se_cutoff_compatible', & 60 routineP = moduleN//':'//routineN 61 62 LOGICAL :: explicit1, explicit2 63 REAL(KIND=dp) :: rc 64 65! Coulomb Cutoff Taper 66 67 CALL section_vals_val_get(se_section, "COULOMB%CUTOFF", explicit=explicit1) 68 CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", explicit=explicit2) 69 IF ((.NOT. explicit1) .AND. se_control%do_ewald_gks) THEN 70 rc = MAX(0.5*plane_distance(1, 0, 0, cell), & 71 0.5*plane_distance(0, 1, 0, cell), & 72 0.5*plane_distance(0, 0, 1, cell)) 73 IF (rc /= se_control%cutoff_cou) THEN 74 IF (output_unit > 0) THEN 75 WRITE (output_unit, *) 76 WRITE (output_unit, '(A,T37,A)') " SEMIEMPIRICAL|", & 77 " Coulomb Integral cutoff/taper was redefined" 78 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Old value [a.u.]", & 79 se_control%cutoff_cou 80 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc 81 WRITE (output_unit, *) 82 END IF 83 END IF 84 se_control%cutoff_cou = rc 85 IF (.NOT. explicit2) se_control%taper_cou = rc 86 ELSE IF ((.NOT. explicit1) .AND. (ALL(cell%perd == 0))) THEN 87 rc = MAX(plane_distance(1, 0, 0, cell), & 88 plane_distance(0, 1, 0, cell), & 89 plane_distance(0, 0, 1, cell)) 90 IF (rc /= se_control%cutoff_cou) THEN 91 IF (output_unit > 0) THEN 92 WRITE (output_unit, *) 93 WRITE (output_unit, '(A,T37,A)') " SEMIEMPIRICAL|", & 94 " Coulomb Integral cutoff/taper was redefined" 95 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Old value [a.u.]", & 96 se_control%cutoff_cou 97 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc 98 WRITE (output_unit, *) 99 END IF 100 END IF 101 se_control%cutoff_cou = rc 102 IF (.NOT. explicit2) se_control%taper_cou = rc 103 END IF 104 IF (output_unit > 0) THEN 105 WRITE (output_unit, *) 106 WRITE (output_unit, '(A,T44,A)') " SEMIEMPIRICAL|", & 107 " Coulomb Integral cutoff/taper values" 108 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Cutoff [a.u.]", & 109 se_control%cutoff_cou 110 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Taper [a.u.]", & 111 se_control%taper_cou 112 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Range [a.u.]", & 113 se_control%range_cou 114 WRITE (output_unit, *) 115 END IF 116 ! Exchange Cutoff Taper 117 CALL section_vals_val_get(se_section, "EXCHANGE%CUTOFF", explicit=explicit1) 118 CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", explicit=explicit2) 119 rc = se_control%cutoff_exc 120 IF (.NOT. explicit1) THEN 121 rc = MIN(rc, MAX(0.25_dp*plane_distance(1, 0, 0, cell), & 122 0.25_dp*plane_distance(0, 1, 0, cell), & 123 0.25_dp*plane_distance(0, 0, 1, cell))) 124 125 IF (rc /= se_control%cutoff_exc) THEN 126 IF (output_unit > 0) THEN 127 WRITE (output_unit, *) 128 WRITE (output_unit, '(A,T36,A)') " SEMIEMPIRICAL|", & 129 " Exchange Integral cutoff/taper was redefined" 130 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Default value [a.u.]", & 131 se_control%cutoff_exc 132 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc 133 WRITE (output_unit, *) 134 END IF 135 END IF 136 END IF 137 se_control%cutoff_exc = rc 138 IF (.NOT. explicit2) se_control%taper_exc = rc 139 140 IF (output_unit > 0) THEN 141 WRITE (output_unit, *) 142 WRITE (output_unit, '(A,T43,A)') " SEMIEMPIRICAL|", & 143 " Exchange Integral cutoff/taper values" 144 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Cutoff [a.u.]", & 145 se_control%cutoff_exc 146 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Taper [a.u.]", & 147 se_control%taper_exc 148 WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Range [a.u.]", & 149 se_control%range_exc 150 WRITE (output_unit, *) 151 END IF 152 153 END SUBROUTINE se_cutoff_compatible 154 155! ************************************************************************************************** 156!> \brief Initializes the semi-empirical taper for a chunk calculation 157!> \param se_taper ... 158!> \param coulomb ... 159!> \param exchange ... 160!> \param lr_corr ... 161!> \author Teodoro Laino [tlaino] - 03.2009 162! ************************************************************************************************** 163 SUBROUTINE initialize_se_taper(se_taper, coulomb, exchange, lr_corr) 164 TYPE(se_taper_type), POINTER :: se_taper 165 LOGICAL, INTENT(IN), OPTIONAL :: coulomb, exchange, lr_corr 166 167 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_se_taper', & 168 routineP = moduleN//':'//routineN 169 170 LOGICAL :: check, l_coulomb, l_exchange, l_lrc 171 172 check = .NOT. ASSOCIATED(se_taper%taper) 173 CPASSERT(check) 174 l_coulomb = .FALSE. 175 l_exchange = .FALSE. 176 l_lrc = .FALSE. 177 IF (PRESENT(coulomb)) l_coulomb = coulomb 178 IF (PRESENT(exchange)) l_exchange = exchange 179 IF (PRESENT(lr_corr)) l_lrc = lr_corr 180 IF (l_coulomb) THEN 181 check = (.NOT. l_exchange) .AND. (.NOT. l_lrc) 182 CPASSERT(check) 183 se_taper%taper => se_taper%taper_cou 184 END IF 185 IF (l_exchange) THEN 186 check = (.NOT. l_coulomb) .AND. (.NOT. l_lrc) 187 CPASSERT(check) 188 se_taper%taper => se_taper%taper_exc 189 END IF 190 IF (l_lrc) THEN 191 check = (.NOT. l_coulomb) .AND. (.NOT. l_exchange) 192 CPASSERT(check) 193 se_taper%taper => se_taper%taper_lrc 194 END IF 195 END SUBROUTINE initialize_se_taper 196 197! ************************************************************************************************** 198!> \brief Finalizes the semi-empirical taper for a chunk calculation 199!> \param se_taper ... 200!> \author Teodoro Laino [tlaino] - 03.2009 201! ************************************************************************************************** 202 SUBROUTINE finalize_se_taper(se_taper) 203 TYPE(se_taper_type), POINTER :: se_taper 204 205 CHARACTER(len=*), PARAMETER :: routineN = 'finalize_se_taper', & 206 routineP = moduleN//':'//routineN 207 208 LOGICAL :: check 209 210 check = ASSOCIATED(se_taper%taper) 211 CPASSERT(check) 212 NULLIFY (se_taper%taper) 213 END SUBROUTINE finalize_se_taper 214 215! ************************************************************************************************** 216!> \brief Initialize semi_empirical type 217!> \param sep ... 218!> \param orb_basis_set ... 219!> \param ngauss ... 220! ************************************************************************************************** 221 SUBROUTINE init_se_param(sep, orb_basis_set, ngauss) 222 TYPE(semi_empirical_type), POINTER :: sep 223 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 224 INTEGER, INTENT(IN) :: ngauss 225 226 CHARACTER(len=*), PARAMETER :: routineN = 'init_se_param', routineP = moduleN//':'//routineN 227 228 CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol 229 INTEGER :: l, nshell 230 INTEGER, DIMENSION(:), POINTER :: lq, nq 231 REAL(KIND=dp), DIMENSION(:), POINTER :: zet 232 233 IF (ASSOCIATED(sep)) THEN 234 CALL allocate_sto_basis_set(sep%basis) 235 nshell = 0 236 IF (sep%natorb == 1) nshell = 1 237 IF (sep%natorb == 4) nshell = 2 238 IF (sep%natorb == 9) nshell = 3 239 ALLOCATE (nq(0:3), lq(0:3), zet(0:3)) 240 241 ALLOCATE (symbol(0:3)) 242 243 symbol = "" 244 nq = 0 245 lq = 0 246 zet = 0._dp 247 DO l = 0, nshell - 1 248 nq(l) = get_se_basis(sep, l) 249 lq(l) = l 250 zet(l) = sep%sto_exponents(l) 251 IF (l == 0) WRITE (symbol(0), '(I1,A1)') nq(l), "S" 252 IF (l == 1) WRITE (symbol(1), '(I1,A1)') nq(l), "P" 253 IF (l == 2) WRITE (symbol(2), '(I1,A1)') nq(l), "D" 254 END DO 255 256 IF (nshell > 0) THEN 257 sep%ngauss = ngauss 258 CALL set_sto_basis_set(sep%basis, name=sep%name, nshell=nshell, symbol=symbol, & 259 nq=nq, lq=lq, zet=zet) 260 CALL create_gto_from_sto_basis(sep%basis, orb_basis_set, sep%ngauss) 261 END IF 262 263 DEALLOCATE (nq) 264 DEALLOCATE (lq) 265 DEALLOCATE (zet) 266 DEALLOCATE (symbol) 267 ELSE 268 CPABORT("The pointer sep is not associated") 269 END IF 270 271 END SUBROUTINE init_se_param 272 273! ************************************************************************************************** 274!> \brief Initialize parameter for a semi_empirival type 275!> \param sep ... 276!> \param z ... 277!> \param method ... 278! ************************************************************************************************** 279 SUBROUTINE se_param_set_default(sep, z, method) 280 281 TYPE(semi_empirical_type), POINTER :: sep 282 INTEGER, INTENT(IN) :: z, method 283 284 CHARACTER(len=*), PARAMETER :: routineN = 'se_param_set_default', & 285 routineP = moduleN//':'//routineN 286 287 IF (ASSOCIATED(sep)) THEN 288 IF (z < 0) THEN 289 CPABORT("Atomic number < 0") 290 END IF 291 SELECT CASE (method) 292 CASE (do_method_am1) 293 CALL am1_default_parameter(sep, z) 294 CASE (do_method_rm1) 295 CALL rm1_default_parameter(sep, z) 296 CASE (do_method_pm3) 297 CALL pm3_default_parameter(sep, z) 298 CASE (do_method_pm6) 299 CALL pm6_default_parameter(sep, z) 300 CASE (do_method_pm6fm) 301 CALL pm6fm_default_parameter(sep, z) 302 CASE (do_method_pdg) 303 CALL pdg_default_parameter(sep, z) 304 CASE (do_method_mndo) 305 CALL mndo_default_parameter(sep, z, do_method_mndo) 306 CASE (do_method_mndod) 307 CALL mndo_default_parameter(sep, z, do_method_mndod) 308 CASE (do_method_pnnl) 309 CALL pnnl_default_parameter(sep, z) 310 CASE (do_method_pchg) 311 CALL pcharge_default_parameter(sep, z) 312 CASE DEFAULT 313 CPABORT("Semiempirical method unknown") 314 END SELECT 315 ELSE 316 CPABORT("The pointer sep is not associated") 317 END IF 318 319 ! Check if the element has been defined.. 320 IF (.NOT. sep%defined) & 321 CALL cp_abort(__LOCATION__, & 322 "Semiempirical type ("//TRIM(sep%name)//") cannot be defined for "// & 323 "the requested parameterization.") 324 325 ! Fill 1 center - 2 electron integrals 326 CALL setup_1c_2el_int(sep) 327 328 ! Fill multipolar expansion of atomic orbitals charge distributions 329 CALL semi_empirical_mpole_p_setup(sep%w_mpole, sep, method) 330 331 ! Get the value of the size of CORE integral array 332 sep%core_size = 0 333 IF (sep%natorb > 0) sep%core_size = 1 334 IF (sep%natorb > 1) sep%core_size = 4 335 IF (sep%dorb) sep%core_size = 10 336 337 ! Get size of the all possible combinations of atomic orbitals 338 sep%atm_int_size = (sep%natorb + 1)*sep%natorb/2 339 340 END SUBROUTINE se_param_set_default 341 342! ************************************************************************************************** 343!> \brief Gives back the unique semi_empirical METHOD type 344!> \param se_method ... 345!> \return ... 346! ************************************************************************************************** 347 FUNCTION get_se_type(se_method) RESULT(se_type) 348 349 INTEGER, INTENT(IN) :: se_method 350 INTEGER :: se_type 351 352 CHARACTER(len=*), PARAMETER :: routineN = 'get_se_type', routineP = moduleN//':'//routineN 353 354 SELECT CASE (se_method) 355 CASE DEFAULT 356 se_type = se_method 357 CASE (do_method_am1, do_method_rm1) 358 se_type = do_method_am1 359 END SELECT 360 361 END FUNCTION get_se_type 362 363END MODULE semi_empirical_utils 364 365