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 a QM/MM calculation 8!> \par History 9!> 5.2004 created [fawzi] 10!> \author Fawzi Mohamed 11! ************************************************************************************************** 12MODULE qmmm_init 13 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind,& 16 set_atomic_kind 17 USE cell_methods, ONLY: read_cell 18 USE cell_types, ONLY: cell_type,& 19 use_perd_xyz 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_eri_mme_interface, ONLY: cp_eri_mme_init_read_input,& 22 cp_eri_mme_set_params 23 USE cp_log_handling, ONLY: cp_get_default_logger,& 24 cp_logger_type 25 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 26 cp_print_key_unit_nr 27 USE cp_para_types, ONLY: cp_para_env_type 28 USE cp_subsys_types, ONLY: cp_subsys_get,& 29 cp_subsys_type 30 USE cp_units, ONLY: cp_unit_from_cp2k,& 31 cp_unit_to_cp2k 32 USE external_potential_types, ONLY: fist_potential_type,& 33 get_potential,& 34 set_potential 35 USE force_field_types, ONLY: input_info_type 36 USE force_fields_input, ONLY: read_gd_section,& 37 read_gp_section,& 38 read_lj_section,& 39 read_wl_section 40 USE input_constants, ONLY: & 41 RADIUS_QMMM_DEFAULT, calc_always, calc_once, do_eri_mme, do_qmmm_gauss, & 42 do_qmmm_image_calcmatrix, do_qmmm_image_iter, do_qmmm_link_gho, do_qmmm_link_imomm, & 43 do_qmmm_link_pseudo, do_qmmm_pcharge, do_qmmm_swave 44 USE input_section_types, ONLY: section_vals_get,& 45 section_vals_get_subs_vals,& 46 section_vals_type,& 47 section_vals_val_get 48 USE kinds, ONLY: default_string_length,& 49 dp 50 USE molecule_kind_types, ONLY: molecule_kind_type 51 USE pair_potential_types, ONLY: pair_potential_reallocate 52 USE particle_list_types, ONLY: particle_list_type 53 USE particle_types, ONLY: particle_type 54 USE pw_env_types, ONLY: pw_env_type 55 USE qmmm_elpot, ONLY: qmmm_potential_init 56 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 57 USE qmmm_gaussian_init, ONLY: qmmm_gaussian_initialize 58 USE qmmm_per_elpot, ONLY: qmmm_ewald_potential_init,& 59 qmmm_per_potential_init 60 USE qmmm_types_low, ONLY: add_set_type,& 61 add_shell_type,& 62 create_add_set_type,& 63 create_add_shell_type,& 64 qmmm_env_mm_type,& 65 qmmm_env_qm_type,& 66 qmmm_links_type 67 USE qs_environment_types, ONLY: get_qs_env,& 68 qs_environment_type 69 USE shell_potential_types, ONLY: get_shell,& 70 shell_kind_type 71#include "./base/base_uses.f90" 72 73 IMPLICIT NONE 74 PRIVATE 75 76 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_init' 78 79 PUBLIC :: assign_mm_charges_and_radius, & 80 print_qmmm_charges, & 81 print_qmmm_links, & 82 print_image_charge_info, & 83 qmmm_init_gaussian_type, & 84 qmmm_init_potential, & 85 qmmm_init_periodic_potential, & 86 setup_qmmm_vars_qm, & 87 setup_qmmm_vars_mm, & 88 setup_qmmm_links, & 89 move_or_add_atoms, & 90 setup_origin_mm_cell 91 92CONTAINS 93 94! ************************************************************************************************** 95!> \brief Assigns charges and radius to evaluate the MM electrostatic potential 96!> \param subsys the subsys containing the MM charges 97!> \param charges ... 98!> \param mm_atom_chrg ... 99!> \param mm_el_pot_radius ... 100!> \param mm_el_pot_radius_corr ... 101!> \param mm_atom_index ... 102!> \param mm_link_atoms ... 103!> \param mm_link_scale_factor ... 104!> \param added_shells ... 105!> \param shell_model ... 106!> \par History 107!> 06.2004 created [tlaino] 108!> \author Teodoro Laino 109! ************************************************************************************************** 110 SUBROUTINE assign_mm_charges_and_radius(subsys, charges, mm_atom_chrg, mm_el_pot_radius, & 111 mm_el_pot_radius_corr, mm_atom_index, mm_link_atoms, & 112 mm_link_scale_factor, added_shells, shell_model) 113 TYPE(cp_subsys_type), POINTER :: subsys 114 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 115 REAL(dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, & 116 mm_el_pot_radius_corr 117 INTEGER, DIMENSION(:), POINTER :: mm_atom_index, mm_link_atoms 118 REAL(dp), DIMENSION(:), POINTER :: mm_link_scale_factor 119 TYPE(add_shell_type), OPTIONAL, POINTER :: added_shells 120 LOGICAL :: shell_model 121 122 CHARACTER(LEN=*), PARAMETER :: routineN = 'assign_mm_charges_and_radius', & 123 routineP = moduleN//':'//routineN 124 125 INTEGER :: I, ilink, IndMM, IndShell, ishell 126 LOGICAL :: is_shell 127 REAL(dp) :: qcore, qi, qshell, rc, ri 128 TYPE(atomic_kind_type), POINTER :: my_kind 129 TYPE(fist_potential_type), POINTER :: my_potential 130 TYPE(particle_list_type), POINTER :: core_particles, particles, & 131 shell_particles 132 TYPE(particle_type), DIMENSION(:), POINTER :: core_set, particle_set, shell_set 133 TYPE(shell_kind_type), POINTER :: shell_kind 134 135 NULLIFY (particle_set, my_kind, added_shells) 136 CALL cp_subsys_get(subsys=subsys, particles=particles, core_particles=core_particles, & 137 shell_particles=shell_particles) 138 particle_set => particles%els 139 140 IF (ALL(particle_set(:)%shell_index .EQ. 0)) THEN 141 shell_model = .FALSE. 142 CALL create_add_shell_type(added_shells, ndim=0) 143 ELSE 144 shell_model = .TRUE. 145 END IF 146 147 IF (shell_model) THEN 148 shell_set => shell_particles%els 149 core_set => core_particles%els 150 ishell = SIZE(shell_set) 151 CALL create_add_shell_type(added_shells, ndim=ishell) 152 added_shells%added_particles => shell_set 153 added_shells%added_cores => core_set 154 END IF 155 156 DO I = 1, SIZE(mm_atom_index) 157 IndMM = mm_atom_index(I) 158 my_kind => particle_set(IndMM)%atomic_kind 159 CALL get_atomic_kind(atomic_kind=my_kind, fist_potential=my_potential, & 160 shell_active=is_shell, shell=shell_kind) 161 CALL get_potential(potential=my_potential, & 162 qeff=qi, & 163 qmmm_radius=ri, & 164 qmmm_corr_radius=rc) 165 IF (ASSOCIATED(charges)) qi = charges(IndMM) 166 mm_atom_chrg(I) = qi 167 mm_el_pot_radius(I) = ri 168 mm_el_pot_radius_corr(I) = rc 169 IF (is_shell) THEN 170 IndShell = particle_set(IndMM)%shell_index 171 IF (ASSOCIATED(shell_kind)) THEN 172 CALL get_shell(shell=shell_kind, & 173 charge_core=qcore, & 174 charge_shell=qshell) 175 mm_atom_chrg(I) = qcore 176 END IF 177 added_shells%mm_core_index(IndShell) = IndMM 178 added_shells%mm_core_chrg(IndShell) = qshell 179 added_shells%mm_el_pot_radius(Indshell) = ri*1.0_dp 180 added_shells%mm_el_pot_radius_corr(Indshell) = rc*1.0_dp 181 END IF 182 END DO 183 184 IF (ASSOCIATED(mm_link_atoms)) THEN 185 DO ilink = 1, SIZE(mm_link_atoms) 186 DO i = 1, SIZE(mm_atom_index) 187 IF (mm_atom_index(i) == mm_link_atoms(ilink)) EXIT 188 END DO 189 IndMM = mm_atom_index(I) 190 mm_atom_chrg(i) = mm_atom_chrg(i)*mm_link_scale_factor(ilink) 191 END DO 192 END IF 193 194 END SUBROUTINE assign_mm_charges_and_radius 195 196! ************************************************************************************************** 197!> \brief Print info on charges generating the qmmm potential.. 198!> \param mm_atom_index ... 199!> \param mm_atom_chrg ... 200!> \param mm_el_pot_radius ... 201!> \param mm_el_pot_radius_corr ... 202!> \param added_charges ... 203!> \param added_shells ... 204!> \param qmmm_section ... 205!> \param nocompatibility ... 206!> \param shell_model ... 207!> \par History 208!> 01.2005 created [tlaino] 209!> \author Teodoro Laino 210! ************************************************************************************************** 211 SUBROUTINE print_qmmm_charges(mm_atom_index, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, & 212 added_charges, added_shells, qmmm_section, nocompatibility, shell_model) 213 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 214 REAL(dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, & 215 mm_el_pot_radius_corr 216 TYPE(add_set_type), POINTER :: added_charges 217 TYPE(add_shell_type), POINTER :: added_shells 218 TYPE(section_vals_type), POINTER :: qmmm_section 219 LOGICAL, INTENT(IN) :: nocompatibility, shell_model 220 221 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_qmmm_charges', & 222 routineP = moduleN//':'//routineN 223 224 INTEGER :: I, ind1, ind2, IndMM, iw 225 REAL(KIND=dp) :: qi, qtot, rc, ri 226 TYPE(cp_logger_type), POINTER :: logger 227 228 qtot = 0.0_dp 229 logger => cp_get_default_logger() 230 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%QMMM_CHARGES", & 231 extension=".log") 232 IF (iw > 0) THEN 233 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79) 234 WRITE (iw, FMT='(/5X,A)') "MM POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL" 235 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79) 236 DO I = 1, SIZE(mm_atom_index) 237 IndMM = mm_atom_index(I) 238 qi = mm_atom_chrg(I) 239 qtot = qtot + qi 240 ri = mm_el_pot_radius(I) 241 rc = mm_el_pot_radius_corr(I) 242 IF (nocompatibility) THEN 243 WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6)') ' MM ATOM:', IndMM, ' RADIUS:', ri, & 244 ' CHARGE:', qi 245 ELSE 246 WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6,/,T56,A12,T69,F12.6)') & 247 ' MM ATOM:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, 'CORR. RADIUS', rc 248 END IF 249 END DO 250 IF (added_charges%num_mm_atoms /= 0) THEN 251 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79) 252 WRITE (iw, '(/5X,A)') "ADDED POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL" 253 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79) 254 DO I = 1, SIZE(added_charges%mm_atom_index) 255 IndMM = added_charges%mm_atom_index(I) 256 qi = added_charges%mm_atom_chrg(I) 257 qtot = qtot + qi 258 ri = added_charges%mm_el_pot_radius(I) 259 ind1 = added_charges%add_env(I)%Index1 260 ind2 = added_charges%add_env(I)%Index2 261 IF (nocompatibility) THEN 262 WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5)') 'MM POINT:', IndMM, ' RADIUS:', ri, & 263 ' CHARGE:', qi, ind1, ind2 264 ELSE 265 WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5,/,T56,A12,T69,F12.6)') & 266 'MM POINT:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, ind1, ind2, 'CORR. RADIUS', rc 267 END IF 268 END DO 269 270 END IF 271 272 IF (shell_model) THEN 273 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73) 274 WRITE (iw, '(/5X,A)') "ADDED SHELL CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL" 275 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73) 276 277 DO I = 1, SIZE(added_shells%mm_core_index) 278 IndMM = added_shells%mm_core_index(I) 279 qi = added_shells%mm_core_chrg(I) 280 qtot = qtot + qi 281 ri = added_shells%mm_el_pot_radius(I) 282 IF (nocompatibility) THEN 283 WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,3F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, & 284 ' CHARGE:', qi, added_shells%added_particles(I)%r 285 ELSE 286 WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,A,F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, & 287 ' CHARGE:', qi, ' CORR. RADIUS', rc 288 END IF 289 290 END DO 291 292 END IF 293 294 WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79) 295 WRITE (iw, '(/,T50,A,T69,F12.6)') ' TOTAL CHARGE:', qtot 296 WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79) 297 END IF 298 CALL cp_print_key_finished_output(iw, logger, qmmm_section, & 299 "PRINT%QMMM_CHARGES") 300 END SUBROUTINE print_qmmm_charges 301 302! ************************************************************************************************** 303!> \brief Print info on qm/mm links 304!> \param qmmm_section ... 305!> \param qmmm_links ... 306!> \par History 307!> 01.2005 created [tlaino] 308!> \author Teodoro Laino 309! ************************************************************************************************** 310 SUBROUTINE print_qmmm_links(qmmm_section, qmmm_links) 311 TYPE(section_vals_type), POINTER :: qmmm_section 312 TYPE(qmmm_links_type), POINTER :: qmmm_links 313 314 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_qmmm_links', & 315 routineP = moduleN//':'//routineN 316 317 INTEGER :: i, iw, mm_index, qm_index 318 REAL(KIND=dp) :: alpha 319 TYPE(cp_logger_type), POINTER :: logger 320 321 logger => cp_get_default_logger() 322 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%qmmm_link_info", extension=".log") 323 IF (iw > 0) THEN 324 IF (ASSOCIATED(qmmm_links)) THEN 325 WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73) 326 WRITE (iw, FMT="(/,T31,A)") " QM/MM LINKS " 327 WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73) 328 IF (ASSOCIATED(qmmm_links%imomm)) THEN 329 WRITE (iw, FMT="(/,T31,A)") " IMOMM TYPE LINK " 330 DO I = 1, SIZE(qmmm_links%imomm) 331 qm_index = qmmm_links%imomm(I)%link%qm_index 332 mm_index = qmmm_links%imomm(I)%link%mm_index 333 alpha = qmmm_links%imomm(I)%link%alpha 334 WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8,T62,A6,F12.6)") "TYPE: IMOMM", & 335 "QM INDEX:", qm_index, "MM INDEX:", mm_index, "ALPHA:", alpha 336 END DO 337 END IF 338 IF (ASSOCIATED(qmmm_links%pseudo)) THEN 339 WRITE (iw, FMT="(/,T31,A)") " PSEUDO TYPE LINK " 340 DO I = 1, SIZE(qmmm_links%pseudo) 341 qm_index = qmmm_links%pseudo(I)%link%qm_index 342 mm_index = qmmm_links%pseudo(I)%link%mm_index 343 WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8)") "TYPE: PSEUDO", & 344 "QM INDEX:", qm_index, "MM INDEX:", mm_index 345 END DO 346 END IF 347 WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 73) 348 ELSE 349 WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73) 350 WRITE (iw, FMT="(/,T26,A)") " NO QM/MM LINKS DETECTED" 351 WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73) 352 END IF 353 END IF 354 CALL cp_print_key_finished_output(iw, logger, qmmm_section, & 355 "PRINT%qmmm_link_info") 356 END SUBROUTINE print_qmmm_links 357 358! ************************************************************************************************** 359!> \brief ... 360!> \param qmmm_env_qm ... 361!> \param para_env ... 362!> \param mm_atom_chrg ... 363!> \param qs_env ... 364!> \param added_charges ... 365!> \param added_shells ... 366!> \param print_section ... 367!> \param qmmm_section ... 368!> \par History 369!> 1.2005 created [tlaino] 370!> \author Teodoro Laino 371! ************************************************************************************************** 372 SUBROUTINE qmmm_init_gaussian_type(qmmm_env_qm, para_env, & 373 mm_atom_chrg, qs_env, added_charges, added_shells, & 374 print_section, qmmm_section) 375 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm 376 TYPE(cp_para_env_type), POINTER :: para_env 377 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg 378 TYPE(qs_environment_type), POINTER :: qs_env 379 TYPE(add_set_type), POINTER :: added_charges 380 TYPE(add_shell_type), POINTER :: added_shells 381 TYPE(section_vals_type), POINTER :: print_section, qmmm_section 382 383 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_init_gaussian_type', & 384 routineP = moduleN//':'//routineN 385 386 INTEGER :: i 387 REAL(KIND=dp) :: maxchrg 388 REAL(KIND=dp), DIMENSION(:), POINTER :: maxradius, maxradius2 389 TYPE(pw_env_type), POINTER :: pw_env 390 391 NULLIFY (maxradius, maxradius2, pw_env) 392 393 maxchrg = MAXVAL(ABS(mm_atom_chrg(:))) 394 CALL get_qs_env(qs_env, pw_env=pw_env) 395 IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:)))) 396 CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=qmmm_env_qm%pgfs, & 397 para_env=para_env, & 398 pw_env=pw_env, & 399 mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, & 400 mm_el_pot_radius_corr=qmmm_env_qm%mm_el_pot_radius_corr, & 401 qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 402 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 403 maxradius=maxradius, & 404 maxchrg=maxchrg, & 405 compatibility=qmmm_env_qm%compatibility, & 406 print_section=print_section, & 407 qmmm_section=qmmm_section) 408 409 IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN 410 CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_charges%pgfs, & 411 para_env=para_env, & 412 pw_env=pw_env, & 413 mm_el_pot_radius=added_charges%mm_el_pot_radius, & 414 mm_el_pot_radius_corr=added_charges%mm_el_pot_radius_corr, & 415 qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 416 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 417 maxradius=maxradius2, & 418 maxchrg=maxchrg, & 419 compatibility=qmmm_env_qm%compatibility, & 420 print_section=print_section, & 421 qmmm_section=qmmm_section) 422 423 SELECT CASE (qmmm_env_qm%qmmm_coupl_type) 424 CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge) 425 DO i = 1, SIZE(maxradius) 426 maxradius(i) = MAX(maxradius(i), maxradius2(i)) 427 END DO 428 END SELECT 429 430 IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2) 431 END IF 432 433 IF (qmmm_env_qm%added_shells%num_mm_atoms .GT. 0) THEN 434 435 maxchrg = MAXVAL(ABS(added_shells%mm_core_chrg(:))) 436 CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_shells%pgfs, & 437 para_env=para_env, & 438 pw_env=pw_env, & 439 mm_el_pot_radius=added_shells%mm_el_pot_radius, & 440 mm_el_pot_radius_corr=added_shells%mm_el_pot_radius_corr, & 441 qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 442 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 443 maxradius=maxradius2, & 444 maxchrg=maxchrg, & 445 compatibility=qmmm_env_qm%compatibility, & 446 print_section=print_section, & 447 qmmm_section=qmmm_section) 448 449 SELECT CASE (qmmm_env_qm%qmmm_coupl_type) 450 CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge) 451 DO i = 1, SIZE(maxradius) 452 maxradius(i) = MAX(maxradius(i), maxradius2(i)) 453 END DO 454 END SELECT 455 456 IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2) 457 458 END IF 459 460 qmmm_env_qm%maxradius => maxradius 461 462 END SUBROUTINE qmmm_init_gaussian_type 463 464! ************************************************************************************************** 465!> \brief ... 466!> \param qmmm_env_qm ... 467!> \param mm_cell ... 468!> \param added_charges ... 469!> \param added_shells ... 470!> \param print_section ... 471!> \par History 472!> 1.2005 created [tlaino] 473!> \author Teodoro Laino 474! ************************************************************************************************** 475 SUBROUTINE qmmm_init_potential(qmmm_env_qm, mm_cell, & 476 added_charges, added_shells, print_section) 477 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm 478 TYPE(cell_type), POINTER :: mm_cell 479 TYPE(add_set_type), POINTER :: added_charges 480 TYPE(add_shell_type), POINTER :: added_shells 481 TYPE(section_vals_type), POINTER :: print_section 482 483 CHARACTER(LEN=*), PARAMETER :: routineN = 'qmmm_init_potential', & 484 routineP = moduleN//':'//routineN 485 486 CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 487 mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, & 488 potentials=qmmm_env_qm%potentials, & 489 pgfs=qmmm_env_qm%pgfs, & 490 mm_cell=mm_cell, & 491 compatibility=qmmm_env_qm%compatibility, & 492 print_section=print_section) 493 494 IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN 495 496 CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 497 mm_el_pot_radius=added_charges%mm_el_pot_radius, & 498 potentials=added_charges%potentials, & 499 pgfs=added_charges%pgfs, & 500 mm_cell=mm_cell, & 501 compatibility=qmmm_env_qm%compatibility, & 502 print_section=print_section) 503 END IF 504 505 IF (qmmm_env_qm%added_shells%num_mm_atoms .GT. 0) THEN 506 507 CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 508 mm_el_pot_radius=added_shells%mm_el_pot_radius, & 509 potentials=added_shells%potentials, & 510 pgfs=added_shells%pgfs, & 511 mm_cell=mm_cell, & 512 compatibility=qmmm_env_qm%compatibility, & 513 print_section=print_section) 514 END IF 515 516 END SUBROUTINE qmmm_init_potential 517 518! ************************************************************************************************** 519!> \brief ... 520!> \param qmmm_env_qm ... 521!> \param qm_cell_small ... 522!> \param mm_cell ... 523!> \param para_env ... 524!> \param qs_env ... 525!> \param added_charges ... 526!> \param added_shells ... 527!> \param qmmm_periodic ... 528!> \param print_section ... 529!> \param mm_atom_chrg ... 530!> \par History 531!> 7.2005 created [tlaino] 532!> \author Teodoro Laino 533! ************************************************************************************************** 534 SUBROUTINE qmmm_init_periodic_potential(qmmm_env_qm, qm_cell_small, mm_cell, para_env, qs_env, & 535 added_charges, added_shells, qmmm_periodic, print_section, mm_atom_chrg) 536 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm 537 TYPE(cell_type), POINTER :: qm_cell_small, mm_cell 538 TYPE(cp_para_env_type), POINTER :: para_env 539 TYPE(qs_environment_type), POINTER :: qs_env 540 TYPE(add_set_type), POINTER :: added_charges 541 TYPE(add_shell_type), POINTER :: added_shells 542 TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section 543 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg 544 545 CHARACTER(LEN=*), PARAMETER :: routineN = 'qmmm_init_periodic_potential', & 546 routineP = moduleN//':'//routineN 547 548 REAL(KIND=dp) :: maxchrg 549 TYPE(dft_control_type), POINTER :: dft_control 550 551 IF (qmmm_env_qm%periodic) THEN 552 553 NULLIFY (dft_control) 554 CALL get_qs_env(qs_env, dft_control=dft_control) 555 556 IF (dft_control%qs_control%semi_empirical) THEN 557 CPABORT("QM/MM periodic calculations not implemented for semi empirical methods") 558 ELSE IF (dft_control%qs_control%dftb) THEN 559 CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, & 560 qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, & 561 para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section) 562 ELSE IF (dft_control%qs_control%xtb) THEN 563 CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, & 564 qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, & 565 para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section) 566 ELSE 567 568 ! setup for GPW/GPAW 569 maxchrg = MAXVAL(ABS(mm_atom_chrg(:))) 570 IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:)))) 571 572 CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 573 per_potentials=qmmm_env_qm%per_potentials, & 574 potentials=qmmm_env_qm%potentials, & 575 pgfs=qmmm_env_qm%pgfs, & 576 qm_cell_small=qm_cell_small, & 577 mm_cell=mm_cell, & 578 para_env=para_env, & 579 compatibility=qmmm_env_qm%compatibility, & 580 qmmm_periodic=qmmm_periodic, & 581 print_section=print_section, & 582 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 583 maxchrg=maxchrg, & 584 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, & 585 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local) 586 587 IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN 588 589 CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 590 per_potentials=added_charges%per_potentials, & 591 potentials=added_charges%potentials, & 592 pgfs=added_charges%pgfs, & 593 qm_cell_small=qm_cell_small, & 594 mm_cell=mm_cell, & 595 para_env=para_env, & 596 compatibility=qmmm_env_qm%compatibility, & 597 qmmm_periodic=qmmm_periodic, & 598 print_section=print_section, & 599 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 600 maxchrg=maxchrg, & 601 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, & 602 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local) 603 END IF 604 605 IF (qmmm_env_qm%added_shells%num_mm_atoms .GT. 0) THEN 606 607 CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, & 608 per_potentials=added_shells%per_potentials, & 609 potentials=added_shells%potentials, & 610 pgfs=added_shells%pgfs, & 611 qm_cell_small=qm_cell_small, & 612 mm_cell=mm_cell, & 613 para_env=para_env, & 614 compatibility=qmmm_env_qm%compatibility, & 615 qmmm_periodic=qmmm_periodic, & 616 print_section=print_section, & 617 eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, & 618 maxchrg=maxchrg, & 619 ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, & 620 ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local) 621 END IF 622 623 END IF 624 625 END IF 626 627 END SUBROUTINE qmmm_init_periodic_potential 628 629! ************************************************************************************************** 630!> \brief ... 631!> \param qmmm_section ... 632!> \param qmmm_env ... 633!> \param subsys_mm ... 634!> \param qm_atom_type ... 635!> \param qm_atom_index ... 636!> \param mm_atom_index ... 637!> \param qm_cell_small ... 638!> \param qmmm_coupl_type ... 639!> \param eps_mm_rspace ... 640!> \param qmmm_link ... 641!> \param para_env ... 642!> \par History 643!> 11.2004 created [tlaino] 644!> \author Teodoro Laino 645! ************************************************************************************************** 646 SUBROUTINE setup_qmmm_vars_qm(qmmm_section, qmmm_env, subsys_mm, qm_atom_type, & 647 qm_atom_index, mm_atom_index, qm_cell_small, qmmm_coupl_type, eps_mm_rspace, & 648 qmmm_link, para_env) 649 TYPE(section_vals_type), POINTER :: qmmm_section 650 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 651 TYPE(cp_subsys_type), POINTER :: subsys_mm 652 CHARACTER(len=default_string_length), & 653 DIMENSION(:), POINTER :: qm_atom_type 654 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_atom_index 655 TYPE(cell_type), POINTER :: qm_cell_small 656 INTEGER, INTENT(OUT) :: qmmm_coupl_type 657 REAL(KIND=dp), INTENT(OUT) :: eps_mm_rspace 658 LOGICAL, INTENT(OUT) :: qmmm_link 659 TYPE(cp_para_env_type), POINTER :: para_env 660 661 CHARACTER(len=*), PARAMETER :: routineN = 'setup_qmmm_vars_qm', & 662 routineP = moduleN//':'//routineN 663 664 CHARACTER(len=default_string_length) :: atmname, mm_atom_kind 665 INTEGER :: i, icount, ikind, ikindr, my_type, & 666 n_rep_val, nkind, size_mm_system 667 INTEGER, DIMENSION(:), POINTER :: mm_link_atoms 668 LOGICAL :: explicit, is_mm, is_qm 669 REAL(KIND=dp) :: tmp_radius, tmp_radius_c 670 REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_sph_cut 671 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 672 TYPE(atomic_kind_type), POINTER :: atomic_kind 673 TYPE(fist_potential_type), POINTER :: fist_potential 674 TYPE(section_vals_type), POINTER :: cell_section, eri_mme_section, & 675 image_charge_section, mm_kinds 676 677 NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut) 678 NULLIFY (image_charge_section) 679 qmmm_link = .FALSE. 680 681 CALL section_vals_get(qmmm_section, explicit=explicit) 682 IF (explicit) THEN 683 ! 684 ! QM_CELL 685 ! 686 cell_section => section_vals_get_subs_vals(qmmm_section, "CELL") 687 CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, & 688 check_for_ref=.FALSE., para_env=para_env) 689 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type) 690 CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace) 691 CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut) 692 CPASSERT(SIZE(tmp_sph_cut) == 2) 693 qmmm_env%spherical_cutoff = tmp_sph_cut 694 IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN 695 qmmm_env%spherical_cutoff(2) = 0.0_dp 696 ELSE 697 IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = EPSILON(0.0_dp) 698 tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2) 699 IF (tmp_radius <= 0.0_dp) & 700 CALL cp_abort(__LOCATION__, & 701 "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// & 702 "the Spherical Cutoff in order to satisfy the previous condition!") 703 END IF 704 ! 705 ! Initialization of arrays and core_charge_radius... 706 ! 707 tmp_radius = 0.0_dp 708 CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds) 709 DO Ikind = 1, SIZE(atomic_kinds%els) 710 atomic_kind => atomic_kinds%els(Ikind) 711 CALL get_atomic_kind(atomic_kind=atomic_kind, & 712 fist_potential=fist_potential) 713 CALL set_potential(potential=fist_potential, & 714 qmmm_radius=tmp_radius, & 715 qmmm_corr_radius=tmp_radius) 716 CALL set_atomic_kind(atomic_kind=atomic_kind, & 717 fist_potential=fist_potential) 718 END DO 719 CALL setup_qm_atom_list(qmmm_section=qmmm_section, & 720 qm_atom_index=qm_atom_index, & 721 qm_atom_type=qm_atom_type, & 722 mm_link_atoms=mm_link_atoms, & 723 qmmm_link=qmmm_link) 724 ! 725 ! MM_KINDS 726 ! 727 mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND") 728 CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind) 729 ! 730 ! Default 731 ! 732 tmp_radius = cp_unit_to_cp2k(RADIUS_QMMM_DEFAULT, "angstrom") 733 Set_Radius_Pot_0: DO IkindR = 1, SIZE(atomic_kinds%els) 734 atomic_kind => atomic_kinds%els(IkindR) 735 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname) 736 CALL get_atomic_kind(atomic_kind=atomic_kind, & 737 fist_potential=fist_potential) 738 CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, & 739 qmmm_corr_radius=tmp_radius) 740 CALL set_atomic_kind(atomic_kind=atomic_kind, & 741 fist_potential=fist_potential) 742 END DO Set_Radius_Pot_0 743 ! 744 ! If present overwrite the kind specified in input file... 745 ! 746 IF (explicit) THEN 747 DO ikind = 1, nkind 748 CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, & 749 c_val=mm_atom_kind) 750 CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius) 751 tmp_radius_c = tmp_radius 752 CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val) 753 IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, & 754 r_val=tmp_radius_c) 755 Set_Radius_Pot_1: DO IkindR = 1, SIZE(atomic_kinds%els) 756 atomic_kind => atomic_kinds%els(IkindR) 757 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname) 758 is_qm = qmmm_ff_precond_only_qm(atmname) 759 IF (TRIM(mm_atom_kind) == atmname) THEN 760 CALL get_atomic_kind(atomic_kind=atomic_kind, & 761 fist_potential=fist_potential) 762 CALL set_potential(potential=fist_potential, & 763 qmmm_radius=tmp_radius, & 764 qmmm_corr_radius=tmp_radius_c) 765 CALL set_atomic_kind(atomic_kind=atomic_kind, & 766 fist_potential=fist_potential) 767 END IF 768 END DO Set_Radius_Pot_1 769 END DO 770 END IF 771 772 !Image charge section 773 774 image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE") 775 CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge) 776 777 ELSE 778 CPABORT("QMMM section not present in input file!") 779 ENDIF 780 ! 781 ! Build MM atoms list 782 ! 783 size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index) 784 IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms) 785 ALLOCATE (mm_atom_index(size_mm_system)) 786 icount = 0 787 788 DO i = 1, SIZE(subsys_mm%particles%els) 789 is_mm = .TRUE. 790 IF (ANY(qm_atom_index == i)) THEN 791 is_mm = .FALSE. 792 END IF 793 IF (ASSOCIATED(mm_link_atoms)) THEN 794 IF (ANY(mm_link_atoms == i) .AND. qmmm_link) is_mm = .TRUE. 795 END IF 796 IF (is_mm) THEN 797 icount = icount + 1 798 IF (icount <= size_mm_system) mm_atom_index(icount) = i 799 END IF 800 END DO 801 CPASSERT(icount == size_mm_system) 802 IF (ASSOCIATED(mm_link_atoms)) THEN 803 DEALLOCATE (mm_link_atoms) 804 END IF 805 806 ! Build image charge atom list + set up variables 807 ! 808 IF (qmmm_env%image_charge) THEN 809 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", & 810 explicit=explicit) 811 IF (explicit) qmmm_env%image_charge_pot%all_mm = .FALSE. 812 813 IF (qmmm_env%image_charge_pot%all_mm) THEN 814 qmmm_env%image_charge_pot%image_mm_list => mm_atom_index 815 ELSE 816 CALL setup_image_atom_list(image_charge_section, qmmm_env, & 817 qm_atom_index, subsys_mm) 818 END IF 819 820 qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els 821 822 CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", & 823 r_val=qmmm_env%image_charge_pot%V0) 824 CALL section_vals_val_get(image_charge_section, "WIDTH", & 825 r_val=qmmm_env%image_charge_pot%eta) 826 CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", & 827 i_val=my_type) 828 SELECT CASE (my_type) 829 CASE (do_qmmm_image_calcmatrix) 830 qmmm_env%image_charge_pot%coeff_iterative = .FALSE. 831 CASE (do_qmmm_image_iter) 832 qmmm_env%image_charge_pot%coeff_iterative = .TRUE. 833 END SELECT 834 835 CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", & 836 l_val=qmmm_env%image_charge_pot%image_restart) 837 838 CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", & 839 i_val=qmmm_env%image_charge_pot%image_matrix_method) 840 841 IF (qmmm_env%image_charge_pot%image_matrix_method .EQ. do_eri_mme) THEN 842 eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME") 843 CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param) 844 CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, & 845 hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, & 846 zet_min=qmmm_env%image_charge_pot%eta, & 847 zet_max=qmmm_env%image_charge_pot%eta, & 848 l_max_zet=0, & 849 l_max=0, & 850 para_env=para_env) 851 852 ENDIF 853 END IF 854 855 END SUBROUTINE setup_qmmm_vars_qm 856 857! ************************************************************************************************** 858!> \brief ... 859!> \param qmmm_section ... 860!> \param qmmm_env ... 861!> \param qm_atom_index ... 862!> \param mm_link_atoms ... 863!> \param mm_link_scale_factor ... 864!> \param fist_scale_charge_link ... 865!> \param qmmm_coupl_type ... 866!> \param qmmm_link ... 867!> \par History 868!> 12.2004 created [tlaino] 869!> \author Teodoro Laino 870! ************************************************************************************************** 871 SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, & 872 mm_link_atoms, mm_link_scale_factor, & 873 fist_scale_charge_link, qmmm_coupl_type, & 874 qmmm_link) 875 TYPE(section_vals_type), POINTER :: qmmm_section 876 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 877 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, mm_link_atoms 878 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_link_scale_factor, & 879 fist_scale_charge_link 880 INTEGER, INTENT(OUT) :: qmmm_coupl_type 881 LOGICAL, INTENT(OUT) :: qmmm_link 882 883 CHARACTER(len=*), PARAMETER :: routineN = 'setup_qmmm_vars_mm', & 884 routineP = moduleN//':'//routineN 885 886 LOGICAL :: explicit 887 TYPE(section_vals_type), POINTER :: qmmm_ff_section 888 889 NULLIFY (qmmm_ff_section) 890 qmmm_link = .FALSE. 891 CALL section_vals_get(qmmm_section, explicit=explicit) 892 IF (explicit) THEN 893 CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type) 894 CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, & 895 mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, & 896 fist_scale_charge_link=fist_scale_charge_link) 897 ! 898 ! Do we want to use a different FF for the non-bonded QM/MM interactions? 899 ! 900 qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD") 901 CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff) 902 IF (qmmm_env%use_qmmm_ff) THEN 903 CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", & 904 l_val=qmmm_env%multiple_potential) 905 CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info) 906 END IF 907 END IF 908 END SUBROUTINE setup_qmmm_vars_mm 909 910! ************************************************************************************************** 911!> \brief reads information regarding the forcefield specific for the QM/MM 912!> interactions 913!> \param qmmm_ff_section ... 914!> \param inp_info ... 915!> \par History 916!> 12.2004 created [tlaino] 917!> \author Teodoro Laino 918! ************************************************************************************************** 919 SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info) 920 TYPE(section_vals_type), POINTER :: qmmm_ff_section 921 TYPE(input_info_type), POINTER :: inp_info 922 923 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_qmmm_ff_section', & 924 routineP = moduleN//':'//routineN 925 926 INTEGER :: n_gd, n_gp, n_lj, n_wl, np 927 TYPE(section_vals_type), POINTER :: gd_section, gp_section, lj_section, & 928 wl_section 929 930! 931! NONBONDED 932! 933 934 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES") 935 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS") 936 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN") 937 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT") 938 CALL section_vals_get(lj_section, n_repetition=n_lj) 939 np = n_lj 940 IF (n_lj /= 0) THEN 941 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.TRUE.) 942 CALL read_lj_section(inp_info%nonbonded, lj_section, start=0) 943 END IF 944 CALL section_vals_get(wl_section, n_repetition=n_wl) 945 np = n_lj + n_wl 946 IF (n_wl /= 0) THEN 947 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.TRUE.) 948 CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj) 949 END IF 950 CALL section_vals_get(gd_section, n_repetition=n_gd) 951 np = n_lj + n_wl + n_gd 952 IF (n_gd /= 0) THEN 953 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.TRUE.) 954 CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl) 955 END IF 956 CALL section_vals_get(gp_section, n_repetition=n_gp) 957 np = n_lj + n_wl + n_gd + n_gp 958 IF (n_gp /= 0) THEN 959 CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.TRUE.) 960 CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd) 961 END IF 962 ! 963 ! NONBONDED14 964 ! 965 lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES") 966 wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS") 967 gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN") 968 gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT") 969 CALL section_vals_get(lj_section, n_repetition=n_lj) 970 np = n_lj 971 IF (n_lj /= 0) THEN 972 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.TRUE.) 973 CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0) 974 END IF 975 CALL section_vals_get(wl_section, n_repetition=n_wl) 976 np = n_lj + n_wl 977 IF (n_wl /= 0) THEN 978 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.TRUE.) 979 CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj) 980 END IF 981 CALL section_vals_get(gd_section, n_repetition=n_gd) 982 np = n_lj + n_wl + n_gd 983 IF (n_gd /= 0) THEN 984 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.TRUE.) 985 CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl) 986 END IF 987 CALL section_vals_get(gp_section, n_repetition=n_gp) 988 np = n_lj + n_wl + n_gd + n_gp 989 IF (n_gp /= 0) THEN 990 CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.TRUE.) 991 CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd) 992 END IF 993 994 END SUBROUTINE read_qmmm_ff_section 995 996! ************************************************************************************************** 997!> \brief ... 998!> \param qmmm_section ... 999!> \param qm_atom_index ... 1000!> \param qm_atom_type ... 1001!> \param mm_link_atoms ... 1002!> \param mm_link_scale_factor ... 1003!> \param qmmm_link ... 1004!> \param fist_scale_charge_link ... 1005!> \par History 1006!> 12.2004 created [tlaino] 1007!> \author Teodoro Laino 1008! ************************************************************************************************** 1009 SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, & 1010 mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link) 1011 TYPE(section_vals_type), POINTER :: qmmm_section 1012 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: qm_atom_index 1013 CHARACTER(len=default_string_length), & 1014 DIMENSION(:), OPTIONAL, POINTER :: qm_atom_type 1015 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mm_link_atoms 1016 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: mm_link_scale_factor 1017 LOGICAL, INTENT(OUT), OPTIONAL :: qmmm_link 1018 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: fist_scale_charge_link 1019 1020 CHARACTER(len=*), PARAMETER :: routineN = 'setup_qm_atom_list', & 1021 routineP = moduleN//':'//routineN 1022 1023 CHARACTER(len=default_string_length) :: qm_atom_kind, qm_link_element 1024 INTEGER :: ikind, k, link_involv_mm, link_type, & 1025 mm_index, n_var, nkind, nlinks, & 1026 num_qm_atom_tot 1027 INTEGER, DIMENSION(:), POINTER :: mm_indexes 1028 LOGICAL :: explicit 1029 REAL(KIND=dp) :: scale_f 1030 TYPE(section_vals_type), POINTER :: qm_kinds, qmmm_links 1031 1032 num_qm_atom_tot = 0 1033 link_involv_mm = 0 1034 nlinks = 0 1035 ! 1036 ! QM_KINDS 1037 ! 1038 qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND") 1039 CALL section_vals_get(qm_kinds, n_repetition=nkind) 1040 DO ikind = 1, nkind 1041 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var) 1042 DO k = 1, n_var 1043 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, & 1044 i_vals=mm_indexes) 1045 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes) 1046 END DO 1047 END DO 1048 ! 1049 ! QM/MM LINKS 1050 ! 1051 qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK") 1052 CALL section_vals_get(qmmm_links, explicit=explicit) 1053 IF (explicit) THEN 1054 qmmm_link = .TRUE. 1055 CALL section_vals_get(qmmm_links, n_repetition=nlinks) 1056 ! Take care of the various link types 1057 DO ikind = 1, nlinks 1058 CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, & 1059 i_val=link_type) 1060 SELECT CASE (link_type) 1061 CASE (do_qmmm_link_imomm) 1062 num_qm_atom_tot = num_qm_atom_tot + 1 1063 link_involv_mm = link_involv_mm + 1 1064 CASE (do_qmmm_link_pseudo) 1065 num_qm_atom_tot = num_qm_atom_tot + 1 1066 CASE (do_qmmm_link_gho) 1067 ! do nothing for the moment 1068 CASE DEFAULT 1069 CPABORT("") 1070 END SELECT 1071 END DO 1072 END IF 1073 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) & 1074 ALLOCATE (mm_link_scale_factor(link_involv_mm)) 1075 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) & 1076 ALLOCATE (fist_scale_charge_link(link_involv_mm)) 1077 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) & 1078 ALLOCATE (mm_link_atoms(link_involv_mm)) 1079 IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot)) 1080 IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot)) 1081 IF (PRESENT(qm_atom_index)) qm_atom_index = 0 1082 IF (PRESENT(qm_atom_type)) qm_atom_type = " " 1083 num_qm_atom_tot = 1 1084 DO ikind = 1, nkind 1085 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var) 1086 DO k = 1, n_var 1087 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, & 1088 i_vals=mm_indexes) 1089 IF (PRESENT(qm_atom_index)) THEN 1090 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:) 1091 END IF 1092 IF (PRESENT(qm_atom_type)) THEN 1093 CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, & 1094 c_val=qm_atom_kind) 1095 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind 1096 END IF 1097 num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes) 1098 END DO 1099 END DO 1100 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp 1101 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp 1102 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0 1103 IF (explicit) THEN 1104 DO ikind = 1, nlinks 1105 IF (PRESENT(qm_atom_type)) THEN 1106 CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element) 1107 qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = TRIM(qm_link_element)//"_LINK" 1108 END IF 1109 IF (PRESENT(qm_atom_index)) THEN 1110 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index) 1111 CPASSERT(ALL(qm_atom_index /= mm_index)) 1112 qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index 1113 num_qm_atom_tot = num_qm_atom_tot + 1 1114 END IF 1115 IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN 1116 CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index) 1117 mm_link_atoms(ikind) = mm_index 1118 END IF 1119 IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN 1120 CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f) 1121 mm_link_scale_factor(ikind) = scale_f 1122 END IF 1123 IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN 1124 CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f) 1125 fist_scale_charge_link(ikind) = scale_f 1126 END IF 1127 END DO 1128 END IF 1129 CPASSERT(num_qm_atom_tot - 1 == SIZE(qm_atom_index)) 1130 1131 END SUBROUTINE setup_qm_atom_list 1132 1133! ************************************************************************************************** 1134!> \brief this routine sets up all variables to treat qmmm links 1135!> \param qmmm_section ... 1136!> \param qmmm_links ... 1137!> \param mm_el_pot_radius ... 1138!> \param mm_el_pot_radius_corr ... 1139!> \param mm_atom_index ... 1140!> \param iw ... 1141!> \par History 1142!> 12.2004 created [tlaino] 1143!> \author Teodoro Laino 1144! ************************************************************************************************** 1145 SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, & 1146 mm_atom_index, iw) 1147 TYPE(section_vals_type), POINTER :: qmmm_section 1148 TYPE(qmmm_links_type), POINTER :: qmmm_links 1149 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius, mm_el_pot_radius_corr 1150 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1151 INTEGER, INTENT(IN) :: iw 1152 1153 CHARACTER(len=*), PARAMETER :: routineN = 'setup_qmmm_links', & 1154 routineP = moduleN//':'//routineN 1155 1156 INTEGER :: ikind, link_type, mm_index, n_gho, & 1157 n_imomm, n_pseudo, n_rep_val, n_tot, & 1158 nlinks, qm_index 1159 INTEGER, DIMENSION(:), POINTER :: wrk_tmp 1160 REAL(KIND=dp) :: alpha, my_radius 1161 TYPE(section_vals_type), POINTER :: qmmm_link_section 1162 1163 NULLIFY (wrk_tmp) 1164 n_imomm = 0 1165 n_gho = 0 1166 n_pseudo = 0 1167 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK") 1168 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks) 1169 CPASSERT(nlinks /= 0) 1170 DO ikind = 1, nlinks 1171 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type) 1172 IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1 1173 IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1 1174 IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1 1175 END DO 1176 n_tot = n_imomm + n_gho + n_pseudo 1177 CPASSERT(n_tot /= 0) 1178 ALLOCATE (qmmm_links) 1179 NULLIFY (qmmm_links%imomm, & 1180 qmmm_links%pseudo) 1181 ! IMOMM 1182 IF (n_imomm /= 0) THEN 1183 ALLOCATE (qmmm_links%imomm(n_imomm)) 1184 ALLOCATE (wrk_tmp(n_imomm)) 1185 DO ikind = 1, n_imomm 1186 NULLIFY (qmmm_links%imomm(ikind)%link) 1187 ALLOCATE (qmmm_links%imomm(ikind)%link) 1188 END DO 1189 n_imomm = 0 1190 DO ikind = 1, nlinks 1191 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type) 1192 IF (link_type == do_qmmm_link_imomm) THEN 1193 n_imomm = n_imomm + 1 1194 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index) 1195 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index) 1196 CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha) 1197 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val) 1198 qmmm_links%imomm(n_imomm)%link%qm_index = qm_index 1199 qmmm_links%imomm(n_imomm)%link%mm_index = mm_index 1200 qmmm_links%imomm(n_imomm)%link%alpha = alpha 1201 wrk_tmp(n_imomm) = mm_index 1202 IF (n_rep_val == 1) THEN 1203 CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius) 1204 WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius 1205 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius 1206 END IF 1207 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val) 1208 IF (n_rep_val == 1) THEN 1209 CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius) 1210 WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius 1211 END IF 1212 END IF 1213 END DO 1214 ! 1215 ! Checking the link structure 1216 ! 1217 DO ikind = 1, SIZE(wrk_tmp) 1218 IF (COUNT(wrk_tmp == wrk_tmp(ikind)) > 1) THEN 1219 WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom." 1220 WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections." 1221 CPABORT("") 1222 END IF 1223 END DO 1224 DEALLOCATE (wrk_tmp) 1225 END IF 1226 ! PSEUDO 1227 IF (n_pseudo /= 0) THEN 1228 ALLOCATE (qmmm_links%pseudo(n_pseudo)) 1229 DO ikind = 1, n_pseudo 1230 NULLIFY (qmmm_links%pseudo(ikind)%link) 1231 ALLOCATE (qmmm_links%pseudo(ikind)%link) 1232 END DO 1233 n_pseudo = 0 1234 DO ikind = 1, nlinks 1235 CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type) 1236 IF (link_type == do_qmmm_link_pseudo) THEN 1237 n_pseudo = n_pseudo + 1 1238 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index) 1239 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index) 1240 qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index 1241 qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index 1242 END IF 1243 END DO 1244 END IF 1245 ! GHO 1246 IF (n_gho /= 0) THEN 1247 ! not yet implemented 1248 ! still to define : type, implementation into QS 1249 CPABORT("") 1250 END IF 1251 END SUBROUTINE setup_qmmm_links 1252 1253! ************************************************************************************************** 1254!> \brief this routine sets up all variables to treat qmmm links 1255!> \param qmmm_section ... 1256!> \param move_mm_charges ... 1257!> \param add_mm_charges ... 1258!> \param mm_atom_chrg ... 1259!> \param mm_el_pot_radius ... 1260!> \param mm_el_pot_radius_corr ... 1261!> \param added_charges ... 1262!> \param mm_atom_index ... 1263!> \par History 1264!> 12.2004 created [tlaino] 1265!> \author Teodoro Laino 1266! ************************************************************************************************** 1267 SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, & 1268 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, & 1269 added_charges, mm_atom_index) 1270 TYPE(section_vals_type), POINTER :: qmmm_section 1271 LOGICAL, INTENT(OUT) :: move_mm_charges, add_mm_charges 1272 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, & 1273 mm_el_pot_radius_corr 1274 TYPE(add_set_type), POINTER :: added_charges 1275 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1276 1277 CHARACTER(len=*), PARAMETER :: routineN = 'move_or_add_atoms', & 1278 routineP = moduleN//':'//routineN 1279 1280 INTEGER :: i_add, icount, ikind, ind1, Index1, & 1281 Index2, n_add_tot, n_adds, n_move_tot, & 1282 n_moves, n_rep_val, nlinks 1283 LOGICAL :: explicit 1284 REAL(KIND=dp) :: alpha, c_radius, charge, radius 1285 TYPE(section_vals_type), POINTER :: add_section, move_section, & 1286 qmmm_link_section 1287 1288 explicit = .FALSE. 1289 move_mm_charges = .FALSE. 1290 add_mm_charges = .FALSE. 1291 NULLIFY (qmmm_link_section, move_section, add_section) 1292 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK") 1293 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks) 1294 CPASSERT(nlinks /= 0) 1295 icount = 0 1296 n_move_tot = 0 1297 n_add_tot = 0 1298 DO ikind = 1, nlinks 1299 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", & 1300 i_rep_section=ikind) 1301 CALL section_vals_get(move_section, n_repetition=n_moves) 1302 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", & 1303 i_rep_section=ikind) 1304 CALL section_vals_get(add_section, n_repetition=n_adds) 1305 n_move_tot = n_move_tot + n_moves 1306 n_add_tot = n_add_tot + n_adds 1307 END DO 1308 icount = n_move_tot + n_add_tot 1309 IF (n_add_tot /= 0) add_mm_charges = .TRUE. 1310 IF (n_move_tot /= 0) move_mm_charges = .TRUE. 1311 ! 1312 ! create add_set_type 1313 ! 1314 CALL create_add_set_type(added_charges, ndim=icount) 1315 ! 1316 ! Fill in structures 1317 ! 1318 icount = 0 1319 DO ikind = 1, nlinks 1320 move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", & 1321 i_rep_section=ikind) 1322 CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves) 1323 ! 1324 ! Moving charge atoms 1325 ! 1326 IF (explicit) THEN 1327 DO i_add = 1, n_moves 1328 icount = icount + 1 1329 CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add) 1330 CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add) 1331 CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add) 1332 CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add) 1333 CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add) 1334 c_radius = radius 1335 IF (n_rep_val == 1) & 1336 CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add) 1337 1338 CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, & 1339 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, & 1340 mm_el_pot_radius_corr=mm_el_pot_radius_corr, & 1341 mm_atom_index=mm_atom_index, move=n_moves, Ind1=ind1) 1342 END DO 1343 mm_atom_chrg(ind1) = 0.0_dp 1344 END IF 1345 1346 add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", & 1347 i_rep_section=ikind) 1348 CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds) 1349 ! 1350 ! Adding charge atoms 1351 ! 1352 IF (explicit) THEN 1353 DO i_add = 1, n_adds 1354 icount = icount + 1 1355 CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add) 1356 CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add) 1357 CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add) 1358 CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add) 1359 CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add) 1360 CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add) 1361 c_radius = radius 1362 IF (n_rep_val == 1) & 1363 CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add) 1364 1365 CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, & 1366 mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, & 1367 mm_el_pot_radius_corr=mm_el_pot_radius_corr, & 1368 mm_atom_index=mm_atom_index) 1369 END DO 1370 END IF 1371 END DO 1372 1373 END SUBROUTINE move_or_add_atoms 1374 1375! ************************************************************************************************** 1376!> \brief this routine sets up all variables of the add_set_type type 1377!> \param added_charges ... 1378!> \param icount ... 1379!> \param Index1 ... 1380!> \param Index2 ... 1381!> \param alpha ... 1382!> \param radius ... 1383!> \param c_radius ... 1384!> \param charge ... 1385!> \param mm_atom_chrg ... 1386!> \param mm_el_pot_radius ... 1387!> \param mm_el_pot_radius_corr ... 1388!> \param mm_atom_index ... 1389!> \param move ... 1390!> \param ind1 ... 1391!> \par History 1392!> 12.2004 created [tlaino] 1393!> \author Teodoro Laino 1394! ************************************************************************************************** 1395 SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, & 1396 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1) 1397 TYPE(add_set_type), POINTER :: added_charges 1398 INTEGER, INTENT(IN) :: icount, Index1, Index2 1399 REAL(KIND=dp), INTENT(IN) :: alpha, radius, c_radius 1400 REAL(KIND=dp), INTENT(IN), OPTIONAL :: charge 1401 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, & 1402 mm_el_pot_radius_corr 1403 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1404 INTEGER, INTENT(in), OPTIONAL :: move 1405 INTEGER, INTENT(OUT), OPTIONAL :: ind1 1406 1407 CHARACTER(len=*), PARAMETER :: routineN = 'set_add_set_type', & 1408 routineP = moduleN//':'//routineN 1409 1410 INTEGER :: i, my_move 1411 REAL(KIND=dp) :: my_c_radius, my_charge, my_radius 1412 1413 my_move = 0 1414 my_radius = radius 1415 my_c_radius = c_radius 1416 IF (PRESENT(charge)) my_charge = charge 1417 IF (PRESENT(move)) my_move = move 1418 i = 1 1419 GetId: DO WHILE (i <= SIZE(mm_atom_index)) 1420 IF (Index1 == mm_atom_index(i)) EXIT GetId 1421 i = i + 1 1422 END DO GetId 1423 IF (PRESENT(ind1)) ind1 = i 1424 CPASSERT(i <= SIZE(mm_atom_index)) 1425 IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/REAL(my_move, KIND=dp) 1426 IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i) 1427 IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i) 1428 1429 added_charges%add_env(icount)%Index1 = Index1 1430 added_charges%add_env(icount)%Index2 = Index2 1431 added_charges%add_env(icount)%alpha = alpha 1432 added_charges%mm_atom_index(icount) = icount 1433 added_charges%mm_atom_chrg(icount) = my_charge 1434 added_charges%mm_el_pot_radius(icount) = my_radius 1435 added_charges%mm_el_pot_radius_corr(icount) = my_c_radius 1436 END SUBROUTINE set_add_set_type 1437 1438! ************************************************************************************************** 1439!> \brief this routine sets up the origin of the MM cell respect to the 1440!> origin of the QM cell. The origin of the QM cell is assumed to be 1441!> in (0.0,0.0,0.0)... 1442!> \param qmmm_section ... 1443!> \param qmmm_env ... 1444!> \param qm_cell_small ... 1445!> \param dr ... 1446!> \par History 1447!> 02.2005 created [tlaino] 1448!> \author Teodoro Laino 1449! ************************************************************************************************** 1450 SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, & 1451 dr) 1452 TYPE(section_vals_type), POINTER :: qmmm_section 1453 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 1454 TYPE(cell_type), POINTER :: qm_cell_small 1455 REAL(KIND=dp), DIMENSION(3), INTENT(in) :: dr 1456 1457 CHARACTER(len=*), PARAMETER :: routineN = 'setup_origin_mm_cell', & 1458 routineP = moduleN//':'//routineN 1459 1460 LOGICAL :: center_grid 1461 REAL(KIND=dp), DIMENSION(3) :: tmp 1462 REAL(KINd=dp), DIMENSION(:), POINTER :: vec 1463 1464! This is the vector that corrects position to apply properly the PBC 1465 1466 tmp(1) = qm_cell_small%hmat(1, 1) 1467 tmp(2) = qm_cell_small%hmat(2, 2) 1468 tmp(3) = qm_cell_small%hmat(3, 3) 1469 CPASSERT(ALL(tmp > 0)) 1470 qmmm_env%dOmmOqm = tmp/2.0_dp 1471 ! This is unit vector to translate the QM system in order to center it 1472 ! in QM cell 1473 CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid) 1474 IF (center_grid) THEN 1475 qmmm_env%utrasl = dr 1476 ELSE 1477 qmmm_env%utrasl = 1.0_dp 1478 ENDIF 1479 CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec) 1480 qmmm_env%transl_v = vec 1481 END SUBROUTINE setup_origin_mm_cell 1482 1483! ************************************************************************************************** 1484!> \brief this routine sets up list of MM atoms carrying an image charge 1485!> \param image_charge_section ... 1486!> \param qmmm_env ... 1487!> \param qm_atom_index ... 1488!> \param subsys_mm ... 1489!> \par History 1490!> 02.2012 created 1491!> \author Dorothea Golze 1492! ************************************************************************************************** 1493 SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, & 1494 qm_atom_index, subsys_mm) 1495 1496 TYPE(section_vals_type), POINTER :: image_charge_section 1497 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 1498 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 1499 TYPE(cp_subsys_type), POINTER :: subsys_mm 1500 1501 CHARACTER(len=*), PARAMETER :: routineN = 'setup_image_atom_list', & 1502 routineP = moduleN//':'//routineN 1503 1504 INTEGER :: atom_a, atom_b, i, j, k, max_index, & 1505 n_var, num_const_atom, & 1506 num_image_mm_atom 1507 INTEGER, DIMENSION(:), POINTER :: mm_indexes 1508 LOGICAL :: fix_xyz, imageind_in_range 1509 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind 1510 1511 NULLIFY (mm_indexes, molecule_kind) 1512 imageind_in_range = .FALSE. 1513 num_image_mm_atom = 0 1514 max_index = 0 1515 1516 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", & 1517 n_rep_val=n_var) 1518 DO i = 1, n_var 1519 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", & 1520 i_rep_val=i, i_vals=mm_indexes) 1521 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes) 1522 END DO 1523 1524 ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom)) 1525 1526 qmmm_env%image_charge_pot%image_mm_list = 0 1527 num_image_mm_atom = 1 1528 1529 DO i = 1, n_var 1530 CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", & 1531 i_rep_val=i, i_vals=mm_indexes) 1532 qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom & 1533 + SIZE(mm_indexes) - 1) = mm_indexes(:) 1534 num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes) 1535 END DO 1536 1537 ! checking, if in range, if list contains QM atoms or any atoms doubled 1538 num_image_mm_atom = num_image_mm_atom - 1 1539 1540 max_index = SIZE(subsys_mm%particles%els) 1541 1542 CPASSERT(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0) 1543 imageind_in_range = (MAXVAL(qmmm_env%image_charge_pot%image_mm_list) <= max_index) & 1544 .AND. (MINVAL(qmmm_env%image_charge_pot%image_mm_list) > 0) 1545 CPASSERT(imageind_in_range) 1546 1547 DO i = 1, num_image_mm_atom 1548 atom_a = qmmm_env%image_charge_pot%image_mm_list(i) 1549 IF (ANY(qm_atom_index == atom_a)) THEN 1550 CPABORT("Image atom list must only contain MM atoms") 1551 ENDIF 1552 DO j = i + 1, num_image_mm_atom 1553 atom_b = qmmm_env%image_charge_pot%image_mm_list(j) 1554 IF (atom_a == atom_b) & 1555 CPABORT("There are atoms doubled in image list.") 1556 ENDDO 1557 ENDDO 1558 1559 ! check if molecules in list carry constraints 1560 num_const_atom = 0 1561 fix_xyz = .TRUE. 1562 IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN 1563 IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN 1564 molecule_kind => subsys_mm%molecule_kinds%els 1565 DO i = 1, SIZE(molecule_kind) 1566 IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT 1567 IF (.NOT. fix_xyz) EXIT 1568 DO j = 1, SIZE(molecule_kind(i)%fixd_list) 1569 IF (.NOT. fix_xyz) EXIT 1570 DO k = 1, num_image_mm_atom 1571 atom_a = qmmm_env%image_charge_pot%image_mm_list(k) 1572 IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN 1573 num_const_atom = num_const_atom + 1 1574 IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN 1575 fix_xyz = .FALSE. 1576 EXIT 1577 ENDIF 1578 ENDIF 1579 ENDDO 1580 ENDDO 1581 ENDDO 1582 ENDIF 1583 ENDIF 1584 1585 ! if all image atoms are constrained, calculate image matrix only 1586 ! once for the first MD or GEO_OPT step (for non-iterative case) 1587 IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN 1588 qmmm_env%image_charge_pot%state_image_matrix = calc_once 1589 ELSE 1590 qmmm_env%image_charge_pot%state_image_matrix = calc_always 1591 ENDIF 1592 1593 END SUBROUTINE setup_image_atom_list 1594 1595! ************************************************************************************************** 1596!> \brief Print info on image charges 1597!> \param qmmm_env ... 1598!> \param qmmm_section ... 1599!> \par History 1600!> 03.2012 created 1601!> \author Dorothea Golze 1602! ************************************************************************************************** 1603 SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section) 1604 1605 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 1606 TYPE(section_vals_type), POINTER :: qmmm_section 1607 1608 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_image_charge_info', & 1609 routineP = moduleN//':'//routineN 1610 1611 INTEGER :: iw 1612 REAL(KIND=dp) :: eta, eta_conv, V0, V0_conv 1613 TYPE(cp_logger_type), POINTER :: logger 1614 1615 logger => cp_get_default_logger() 1616 iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", & 1617 extension=".log") 1618 eta = qmmm_env%image_charge_pot%eta 1619 eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2) 1620 V0 = qmmm_env%image_charge_pot%V0 1621 V0_conv = cp_unit_from_cp2k(V0, "volt") 1622 1623 IF (iw > 0) THEN 1624 WRITE (iw, FMT="(T25,A)") "IMAGE CHARGE PARAMETERS" 1625 WRITE (iw, FMT="(T25,A)") REPEAT("-", 23) 1626 WRITE (iw, FMT="(/)") 1627 WRITE (iw, FMT="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:" 1628 WRITE (iw, FMT="(/)") 1629 1630 WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list 1631 WRITE (iw, FMT="(/)") 1632 WRITE (iw, "(T2,A52,T69,F12.8)") & 1633 "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv 1634 WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", V0_conv 1635 WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79) 1636 END IF 1637 CALL cp_print_key_finished_output(iw, logger, qmmm_section, & 1638 "PRINT%PROGRAM_RUN_INFO") 1639 1640 END SUBROUTINE print_image_charge_info 1641 1642END MODULE qmmm_init 1643