1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 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_create 13 USE bibliography, ONLY: Golze2013,& 14 Laino2005,& 15 Laino2006,& 16 cite_reference 17 USE cell_methods, ONLY: write_cell 18 USE cell_types, ONLY: cell_clone,& 19 cell_release,& 20 cell_type,& 21 get_cell 22 USE cp_log_handling, ONLY: cp_get_default_logger,& 23 cp_logger_type 24 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 25 cp_print_key_unit_nr 26 USE cp_para_types, ONLY: cp_para_env_type 27 USE cp_subsys_methods, ONLY: create_small_subsys 28 USE cp_subsys_types, ONLY: cp_subsys_release,& 29 cp_subsys_type 30 USE fist_environment, ONLY: fist_init 31 USE fist_environment_types, ONLY: fist_env_create,& 32 fist_env_get,& 33 fist_env_set,& 34 fist_environment_type 35 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type 36 USE global_types, ONLY: global_environment_type 37 USE header, ONLY: qmmm_header 38 USE input_constants, ONLY: & 39 do_fist, do_multipole_section_off, do_multipole_section_on, do_qmmm, & 40 do_qmmm_center_every_step, do_qmmm_center_never, do_qmmm_center_pbc_aware, & 41 do_qmmm_center_setup_only, do_qmmm_none, do_qs 42 USE input_section_types, ONLY: section_vals_get,& 43 section_vals_get_subs_vals,& 44 section_vals_type,& 45 section_vals_val_get,& 46 section_vals_val_set 47 USE kinds, ONLY: default_string_length,& 48 dp 49 USE pw_env_types, ONLY: pw_env_type 50 USE qmmm_init, ONLY: & 51 assign_mm_charges_and_radius, move_or_add_atoms, print_image_charge_info, & 52 print_qmmm_charges, print_qmmm_links, qmmm_init_gaussian_type, & 53 qmmm_init_periodic_potential, qmmm_init_potential, setup_origin_mm_cell, setup_qmmm_links, & 54 setup_qmmm_vars_mm, setup_qmmm_vars_qm 55 USE qmmm_links_methods, ONLY: qmmm_link_Imomm_coord 56 USE qmmm_pw_grid, ONLY: qmmm_pw_grid_init 57 USE qmmm_types, ONLY: qmmm_env_type 58 USE qmmm_types_low, ONLY: & 59 add_set_release, add_set_type, add_shell_type, qmmm_env_mm_create, qmmm_env_mm_release, & 60 qmmm_env_mm_type, qmmm_env_qm_create, qmmm_env_qm_type, qmmm_links_type 61 USE qs_environment, ONLY: qs_init 62 USE qs_environment_types, ONLY: get_qs_env,& 63 qs_env_create,& 64 qs_environment_type,& 65 set_qs_env 66#include "./base/base_uses.f90" 67 68 IMPLICIT NONE 69 PRIVATE 70 71 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_create' 73 74 PUBLIC :: qmmm_env_create 75 76CONTAINS 77 78! ************************************************************************************************** 79!> \brief ... 80!> \param qmmm_env ... 81!> \param root_section ... 82!> \param para_env ... 83!> \param globenv ... 84!> \param force_env_section ... 85!> \param qmmm_section ... 86!> \param subsys_section ... 87!> \param use_motion_section ... 88!> \param prev_subsys ... 89!> \param ignore_outside_box ... 90!> \par History 91!> 05.2004 created [fawzi] 92!> \author Fawzi Mohamed 93! ************************************************************************************************** 94 SUBROUTINE qmmm_env_create(qmmm_env, root_section, para_env, globenv, & 95 force_env_section, qmmm_section, subsys_section, use_motion_section, prev_subsys, & 96 ignore_outside_box) 97 TYPE(qmmm_env_type), POINTER :: qmmm_env 98 TYPE(section_vals_type), POINTER :: root_section 99 TYPE(cp_para_env_type), POINTER :: para_env 100 TYPE(global_environment_type), POINTER :: globenv 101 TYPE(section_vals_type), POINTER :: force_env_section, qmmm_section, & 102 subsys_section 103 LOGICAL, INTENT(IN) :: use_motion_section 104 TYPE(cp_subsys_type), OPTIONAL, POINTER :: prev_subsys 105 LOGICAL, INTENT(in), OPTIONAL :: ignore_outside_box 106 107 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_env_create', & 108 routineP = moduleN//':'//routineN 109 110 CHARACTER(len=default_string_length), & 111 DIMENSION(:), POINTER :: qm_atom_type 112 INTEGER :: center_i, delta_charge, handle, iw, iw2, & 113 orig_charge, qmmm_coupl_type, & 114 use_multipole 115 INTEGER, DIMENSION(:), POINTER :: mm_atom_index, mm_link_atoms, & 116 qm_atom_index 117 LOGICAL :: add_mm_charges, explicit, & 118 move_mm_charges, nocompatibility, & 119 qmmm_link, qmmm_link_Imomm, shell_model 120 REAL(dp), DIMENSION(:), POINTER :: mm_atom_chrg, mm_el_pot_radius, & 121 mm_el_pot_radius_corr 122 REAL(KIND=dp) :: eps_mm_rspace 123 REAL(KIND=dp), DIMENSION(3) :: abc_mm, abc_qm 124 REAL(KIND=dp), DIMENSION(:), POINTER :: fist_scale_charge_link, & 125 mm_link_scale_factor 126 TYPE(add_set_type), POINTER :: added_charges 127 TYPE(add_shell_type), POINTER :: added_shells 128 TYPE(cell_type), POINTER :: mm_cell, qm_cell_small, super_cell 129 TYPE(cp_logger_type), POINTER :: logger 130 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 131 TYPE(fist_environment_type), POINTER :: fist_env 132 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 133 TYPE(pw_env_type), POINTER :: pw_env 134 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env_mm 135 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm 136 TYPE(qmmm_links_type), POINTER :: qmmm_links 137 TYPE(qs_environment_type), POINTER :: qs_env 138 TYPE(section_vals_type), POINTER :: multipole_section, print_gen, & 139 print_section, qmmm_periodic 140 141 CALL timeset(routineN, handle) 142 143 NULLIFY (qm_atom_index, mm_atom_index, qm_atom_type) 144 NULLIFY (qmmm_env_qm, subsys_mm, subsys_qm, mm_cell, qm_cell_small) 145 NULLIFY (mm_atom_chrg, mm_el_pot_radius, qmmm_env_mm, fist_env) 146 NULLIFY (qmmm_env) 147 NULLIFY (mm_link_atoms, mm_link_scale_factor, qmmm_links, added_charges, added_shells) 148 NULLIFY (fist_scale_charge_link, print_section, fist_nonbond_env) 149 NULLIFY (print_gen, logger, mm_el_pot_radius_corr, super_cell, pw_env) 150 151 logger => cp_get_default_logger() 152 153 ! citations 154 CALL cite_reference(Laino2005) 155 156 ! Input section... 157 IF (.NOT. ASSOCIATED(subsys_section)) THEN 158 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS") 159 END IF 160 qmmm_periodic => section_vals_get_subs_vals(qmmm_section, "PERIODIC") 161 multipole_section => section_vals_get_subs_vals(qmmm_section, "PERIODIC%MULTIPOLE") 162 print_section => section_vals_get_subs_vals(qmmm_section, "PRINT") 163 print_gen => section_vals_get_subs_vals(print_section, "PROGRAM_RUN_INFO") 164 iw = cp_print_key_unit_nr(logger, print_gen, "", extension=".log") 165 166 ! Create QM/MM Environments.. 167 CALL qmmm_env_qm_create(qmmm_env_qm) 168 CALL qmmm_env_mm_create(qmmm_env_mm) 169 170 ! Set up QM/MM Options 171 CALL setup_qmmm_vars_mm(qmmm_section, & 172 qmmm_env_mm, & 173 qm_atom_index, & 174 mm_link_atoms, & 175 mm_link_scale_factor, & 176 fist_scale_charge_link, & 177 qmmm_coupl_type, & 178 qmmm_link) 179 180 qmmm_env_mm%qm_atom_index => qm_atom_index 181 qmmm_env_mm%mm_link_atoms => mm_link_atoms 182 qmmm_env_mm%mm_link_scale_factor => mm_link_scale_factor 183 qmmm_env_mm%fist_scale_charge_link => fist_scale_charge_link 184 qmmm_env_mm%qmmm_coupl_type = qmmm_coupl_type 185 qmmm_env_mm%qmmm_link = qmmm_link 186 ! Center the qm subsys into the qm box 187 CALL section_vals_val_get(qmmm_section, "CENTER", i_val=center_i) 188 IF (center_i == do_qmmm_center_never) THEN 189 qmmm_env_qm%center_qm_subsys = .FALSE. 190 qmmm_env_qm%center_qm_subsys0 = .FALSE. 191 ELSE IF (center_i == do_qmmm_center_setup_only) THEN 192 qmmm_env_qm%center_qm_subsys = .FALSE. 193 qmmm_env_qm%center_qm_subsys0 = .TRUE. 194 ELSE IF (center_i == do_qmmm_center_every_step) THEN 195 qmmm_env_qm%center_qm_subsys = .TRUE. 196 qmmm_env_qm%center_qm_subsys0 = .TRUE. 197 ELSE 198 CPABORT("Unknown type of CENTER! ") 199 ENDIF 200 201 CALL section_vals_val_get(qmmm_section, "CENTER_TYPE", i_val=center_i) 202 qmmm_env_qm%center_qm_subsys_pbc_aware = (center_i == do_qmmm_center_pbc_aware) 203 204 ! Compatibility with the QM/MM in CPMD code 205 CALL section_vals_val_get(qmmm_section, "NOCOMPATIBILITY", l_val=nocompatibility) 206 qmmm_env_qm%compatibility = .NOT. nocompatibility 207 208 ! Parallel scheme for the long range 209 CALL section_vals_val_get(qmmm_section, "PARALLEL_SCHEME", & 210 i_val=qmmm_env_qm%par_scheme) 211 212 ! Periodic boundary condition calculation 213 CALL section_vals_get(qmmm_periodic, explicit=explicit) 214 qmmm_env_qm%periodic = explicit 215 !multipole section is switched on by default; switched off only if explicitly stated 216 IF (qmmm_env_qm%periodic) qmmm_env_qm%multipole = .TRUE. 217 CALL section_vals_get(multipole_section, explicit=explicit) 218 CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", i_val=use_multipole) 219 IF (explicit .AND. use_multipole == do_multipole_section_off) qmmm_env_qm%multipole = .FALSE. 220 IF (explicit .AND. use_multipole == do_multipole_section_on) qmmm_env_qm%multipole = .TRUE. 221 IF (qmmm_env_qm%periodic .AND. qmmm_env_qm%multipole) CALL cite_reference(Laino2006) 222 IF (qmmm_coupl_type == do_qmmm_none) THEN 223 IF (qmmm_env_qm%periodic) & 224 CALL cp_warn(__LOCATION__, & 225 "QMMM periodic calculation with coupling NONE was requested! "// & 226 "Switching off the periodic keyword since periodic and non-periodic "// & 227 "calculation with coupling NONE represent the same method! ") 228 qmmm_env_qm%periodic = .FALSE. 229 END IF 230 231 ! First Initialize Fist... 232 CALL section_vals_val_set(force_env_section, "METHOD", i_val=do_fist) 233 CALL fist_env_create(fist_env, para_env=para_env) 234 CALL fist_env_set(fist_env, qmmm=.TRUE., qmmm_env=qmmm_env_mm) 235 CALL fist_init(fist_env, root_section, para_env, force_env_section, & 236 subsys_section, use_motion_section, prev_subsys=prev_subsys) 237 CALL fist_env_get(fist_env, subsys=subsys_mm, cell=mm_cell) 238 239 ! Set up QM/MM Options 240 CALL setup_qmmm_vars_qm(qmmm_section, & 241 qmmm_env_qm, & 242 subsys_mm, & 243 qm_atom_type, & 244 qm_atom_index, & 245 mm_atom_index, & 246 qm_cell_small, & 247 qmmm_coupl_type, & 248 eps_mm_rspace, & 249 qmmm_link, & 250 para_env) 251 252 qmmm_env_qm%qm_atom_index => qm_atom_index 253 qmmm_env_qm%mm_atom_index => mm_atom_index 254 qmmm_env_qm%eps_mm_rspace = eps_mm_rspace 255 qmmm_env_qm%qmmm_coupl_type = qmmm_coupl_type 256 qmmm_env_qm%qmmm_link = qmmm_link 257 qmmm_env_qm%num_qm_atoms = SIZE(qm_atom_index) 258 qmmm_env_qm%num_mm_atoms = SIZE(mm_atom_index) 259 IF (qmmm_env_qm%image_charge) THEN 260 qmmm_env_qm%num_image_mm_atoms = SIZE(qmmm_env_qm%image_charge_pot%image_mm_list) 261 CALL cite_reference(Golze2013) 262 END IF 263 264 ! Duplicate structure for link atoms 265 IF (qmmm_link) THEN 266 IF (ASSOCIATED(mm_link_atoms)) THEN 267 ALLOCATE (qmmm_env_qm%mm_link_atoms(SIZE(mm_link_atoms))) 268 qmmm_env_qm%mm_link_atoms = mm_link_atoms 269 END IF 270 END IF 271 IF (iw > 0) THEN 272 WRITE (iw, '(A,I26)') " Number of QM atoms: ", qmmm_env_qm%num_qm_atoms 273 WRITE (iw, '(A,I26)') " Number of MM atoms: ", qmmm_env_qm%num_mm_atoms 274 IF (qmmm_env_qm%image_charge) THEN 275 WRITE (iw, '(A,I8)') " Number of MM atoms with image charge: ", & 276 qmmm_env_qm%num_image_mm_atoms 277 ENDIF 278 WRITE (iw, '(A)') " QM cell ::" 279 CALL write_cell(qm_cell_small, subsys_section) 280 END IF 281 CALL get_cell(qm_cell_small, abc=abc_qm) 282 CALL get_cell(mm_cell, abc=abc_mm) 283 284 IF (qmmm_env_qm%image_charge) THEN 285 IF (ANY(ABS(abc_mm - abc_qm) > 1.0E-12)) & 286 CPABORT("QM and MM box need to have the same size when using image charges") 287 ENDIF 288 289 ! Assign charges and mm_el_pot_radius from fist_topology 290 CALL fist_env_get(fist_env, fist_nonbond_env=fist_nonbond_env) 291 ALLOCATE (mm_atom_chrg(SIZE(mm_atom_index))) 292 ALLOCATE (mm_el_pot_radius(SIZE(mm_atom_index))) 293 ALLOCATE (mm_el_pot_radius_corr(SIZE(mm_atom_index))) 294 mm_atom_chrg = 0.0_dp 295 mm_el_pot_radius = 0.0_dp 296 mm_el_pot_radius_corr = 0.0_dp 297 298 CALL assign_mm_charges_and_radius(subsys=subsys_mm, & 299 charges=fist_nonbond_env%charges, & 300 mm_atom_chrg=mm_atom_chrg, & 301 mm_el_pot_radius=mm_el_pot_radius, & 302 mm_el_pot_radius_corr=mm_el_pot_radius_corr, & 303 mm_atom_index=mm_atom_index, & 304 mm_link_atoms=mm_link_atoms, & 305 mm_link_scale_factor=mm_link_scale_factor, & 306 added_shells=added_shells, & 307 shell_model=shell_model) 308 309 qmmm_env_qm%mm_atom_chrg => mm_atom_chrg 310 qmmm_env_qm%mm_el_pot_radius => mm_el_pot_radius 311 qmmm_env_qm%mm_el_pot_radius_corr => mm_el_pot_radius_corr 312 qmmm_env_qm%added_shells => added_shells 313 314 qmmm_link_Imomm = .FALSE. 315 IF (qmmm_link) THEN 316 CALL setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, & 317 mm_el_pot_radius_corr, mm_atom_index, iw) 318 qmmm_env_qm%qmmm_links => qmmm_links 319 320 CALL print_qmmm_links(qmmm_section, qmmm_links) 321 322 CALL add_set_release(qmmm_env_qm%added_charges) 323 CALL move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, & 324 mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, & 325 added_charges, mm_atom_index) 326 qmmm_env_qm%move_mm_charges = move_mm_charges 327 qmmm_env_qm%add_mm_charges = add_mm_charges 328 qmmm_env_qm%added_charges => added_charges 329 IF (ASSOCIATED(qmmm_links%imomm)) qmmm_link_imomm = (SIZE(qmmm_links%imomm) /= 0) 330 END IF 331 332 CALL print_qmmm_charges(mm_atom_index, mm_atom_chrg, mm_el_pot_radius, & 333 mm_el_pot_radius_corr, qmmm_env_qm%added_charges, & 334 qmmm_env_qm%added_shells, qmmm_section, nocompatibility, shell_model) 335 IF (qmmm_env_qm%image_charge) THEN 336 CALL print_image_charge_info(qmmm_env_qm, qmmm_section) 337 ENDIF 338 339 CALL section_vals_val_get(qmmm_section, "DELTA_CHARGE", i_val=delta_charge) 340 CALL section_vals_val_get(force_env_section, "DFT%CHARGE", i_val=orig_charge) 341 CALL section_vals_val_set(force_env_section, "DFT%CHARGE", i_val=orig_charge + delta_charge) 342 343 CALL section_vals_val_set(force_env_section, "METHOD", i_val=do_qs) 344 CALL create_small_subsys(subsys_qm, & 345 big_subsys=subsys_mm, small_para_env=para_env, & 346 small_cell=qm_cell_small, sub_atom_index=qm_atom_index, & 347 sub_atom_kind_name=qm_atom_type, para_env=para_env, & 348 force_env_section=force_env_section, subsys_section=subsys_section, & 349 ignore_outside_box=ignore_outside_box) 350 IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, subsys_qm%particles%els, & 351 qm_atom_index) 352 CALL qs_env_create(qs_env, globenv) 353 CALL qs_init(qs_env, para_env, root_section, globenv=globenv, cp_subsys=subsys_qm, & 354 cell=qm_cell_small, qmmm=.TRUE., qmmm_env_qm=qmmm_env_qm, & 355 force_env_section=force_env_section, subsys_section=subsys_section, & 356 use_motion_section=use_motion_section) 357 CALL cp_subsys_release(subsys_qm) 358 359 IF (qmmm_env_qm%periodic) THEN 360 IF (.NOT. ASSOCIATED(super_cell)) THEN 361 ALLOCATE (super_cell) 362 END IF 363 CALL cell_clone(mm_cell, super_cell) 364 CALL set_qs_env(qs_env, super_cell=super_cell, qmmm_periodic=qmmm_env_qm%periodic) 365 CALL cell_release(super_cell) 366 END IF 367 CALL section_vals_val_set(force_env_section, "DFT%CHARGE", i_val=orig_charge) 368 CALL cp_print_key_finished_output(iw, logger, print_gen, "") 369 iw2 = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_BANNER", & 370 extension=".qmmmLog") 371 CALL qmmm_header(iw2) 372 CALL cp_print_key_finished_output(iw2, logger, qmmm_section, & 373 "PRINT%PROGRAM_BANNER") 374 ! 375 ! Initialize MM Potential fitted with Gaussian 376 ! 377 CALL qmmm_init_gaussian_type(qmmm_env_qm=qmmm_env_qm, & 378 para_env=para_env, & 379 qs_env=qs_env, & 380 mm_atom_chrg=mm_atom_chrg, & 381 added_charges=qmmm_env_qm%added_charges, & 382 added_shells=qmmm_env_qm%added_shells, & 383 print_section=print_section, & 384 qmmm_section=qmmm_section) 385 ! 386 ! Initialize the MM potential stored on vector 387 ! 388 CALL qmmm_init_potential(qmmm_env_qm=qmmm_env_qm, & 389 mm_cell=mm_cell, & 390 added_charges=qmmm_env_qm%added_charges, & 391 added_shells=qmmm_env_qm%added_shells, & 392 print_section=print_section) 393 ! 394 ! Initialize the qmmm_pw_grid 395 ! 396 CALL get_qs_env(qs_env, pw_env=pw_env) 397 CALL qmmm_pw_grid_init(qmmm_env=qmmm_env_qm, & 398 pw_env=pw_env) 399 ! 400 ! Initialize the MM periodic potential 401 ! 402 CALL qmmm_init_periodic_potential(qmmm_env_qm=qmmm_env_qm, & 403 qm_cell_small=qm_cell_small, & 404 mm_cell=mm_cell, & 405 para_env=para_env, & 406 qs_env=qs_env, & 407 added_charges=qmmm_env_qm%added_charges, & 408 added_shells=qmmm_env_qm%added_shells, & 409 qmmm_periodic=qmmm_periodic, & 410 print_section=print_section, & 411 mm_atom_chrg=mm_atom_chrg) 412 ! 413 ! Preparing for PBC... 414 ! 415 CALL setup_origin_mm_cell(qmmm_section, qmmm_env_qm, qm_cell_small, & 416 dr=pw_env%pw_pools(pw_env%auxbas_grid)%pool%pw_grid%dr) 417 418 CALL cell_release(qm_cell_small) 419 420 ! assemble the actual qmmm_env 421 ALLOCATE (qmmm_env) 422 qmmm_env%qs_env => qs_env 423 qmmm_env%fist_env => fist_env 424 qmmm_env%qm => qmmm_env_qm 425 426 ! The qmmm_env inherits our ref_cout for qmmm_env_qm, fist_env, qs_env 427 ! An exception is qmmm_env_mm, because it's buried in the fist_env 428 CALL qmmm_env_mm_release(qmmm_env_mm) 429 430 CALL section_vals_val_set(force_env_section, "METHOD", i_val=do_qmmm) 431 DEALLOCATE (qm_atom_type) 432 433 CALL timestop(handle) 434 435 END SUBROUTINE qmmm_env_create 436 437END MODULE qmmm_create 438