1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> - Merged with the Quickstep MODULE method_specification (17.01.2002,MK) 9!> - USE statements cleaned, added 10!> (25.09.2002,MK) 11!> - Added more LSD structure (01.2003,Joost VandeVondele) 12!> - New molecule data types introduced (Sep. 2003,MK) 13!> - Cleaning; getting rid of pnode (02.10.2003,MK) 14!> - Sub-system setup added (08.10.2003,MK) 15!> \author MK (18.05.2000) 16! ************************************************************************************************** 17MODULE qs_environment 18 USE almo_scf_env_methods, ONLY: almo_scf_env_create 19 USE atom_kind_orbitals, ONLY: calculate_atomic_relkin 20 USE atomic_kind_types, ONLY: atomic_kind_type 21 USE auto_basis, ONLY: create_aux_fit_basis_set,& 22 create_lri_aux_basis_set,& 23 create_ri_aux_basis_set 24 USE basis_set_container_types, ONLY: add_basis_set_to_container 25 USE basis_set_types, ONLY: get_gto_basis_set,& 26 gto_basis_set_type 27 USE bibliography, ONLY: Iannuzzi2006,& 28 Iannuzzi2007,& 29 cite_reference 30 USE cell_types, ONLY: cell_type 31 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 32 cp_blacs_env_release,& 33 cp_blacs_env_retain,& 34 cp_blacs_env_type 35 USE cp_control_types, ONLY: dft_control_release,& 36 dft_control_type,& 37 dftb_control_type,& 38 gapw_control_type,& 39 qs_control_type,& 40 semi_empirical_control_type,& 41 xtb_control_type 42 USE cp_control_utils, ONLY: read_ddapc_section,& 43 read_dft_control,& 44 read_mgrid_section,& 45 read_qs_section,& 46 read_tddfpt2_control,& 47 read_tddfpt_control,& 48 write_dft_control,& 49 write_qs_control 50 USE cp_ddapc_types, ONLY: cp_ddapc_ewald_create 51 USE cp_log_handling, ONLY: cp_get_default_logger,& 52 cp_logger_get_default_io_unit,& 53 cp_logger_type 54 USE cp_output_handling, ONLY: cp_p_file,& 55 cp_print_key_finished_output,& 56 cp_print_key_should_output,& 57 cp_print_key_unit_nr 58 USE cp_para_env, ONLY: cp_para_env_retain 59 USE cp_para_types, ONLY: cp_para_env_type 60 USE cp_subsys_types, ONLY: cp_subsys_type 61 USE cp_symmetry, ONLY: write_symmetry 62 USE distribution_1d_types, ONLY: distribution_1d_release,& 63 distribution_1d_type 64 USE distribution_methods, ONLY: distribute_molecules_1d 65 USE dm_ls_scf_create, ONLY: ls_scf_create 66 USE et_coupling_types, ONLY: et_coupling_create 67 USE ewald_environment_types, ONLY: ewald_env_create,& 68 ewald_env_get,& 69 ewald_env_release,& 70 ewald_env_set,& 71 ewald_environment_type,& 72 read_ewald_section,& 73 read_ewald_section_tb 74 USE ewald_pw_methods, ONLY: ewald_pw_grid_update 75 USE ewald_pw_types, ONLY: ewald_pw_create,& 76 ewald_pw_release,& 77 ewald_pw_type 78 USE external_potential_types, ONLY: get_potential,& 79 init_potential,& 80 set_potential 81 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_create,& 82 fist_nonbond_env_type 83 USE gamma, ONLY: init_md_ftable 84 USE global_types, ONLY: global_environment_type 85 USE hartree_local_methods, ONLY: init_coulomb_local 86 USE header, ONLY: dftb_header,& 87 qs_header,& 88 se_header,& 89 xtb_header 90 USE hfx_types, ONLY: hfx_create 91 USE input_constants, ONLY: & 92 dispersion_d3, do_et_ddapc, do_method_am1, do_method_dftb, do_method_gapw, & 93 do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, & 94 do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, do_method_pm6fm, & 95 do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, do_qmmm_gauss, & 96 do_qmmm_swave, general_roks, kg_tnadd_embed_ri, rel_none, rel_trans_atom, & 97 vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, wfi_aspc_nr, wfi_linear_ps_method_nr, & 98 wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, xc_vdw_fun_none, & 99 xc_vdw_fun_nonloc, xc_vdw_fun_pairpot 100 USE input_section_types, ONLY: section_vals_get,& 101 section_vals_get_subs_vals,& 102 section_vals_type,& 103 section_vals_val_get 104 USE kg_environment, ONLY: kg_env_create 105 USE kinds, ONLY: dp 106 USE kpoint_methods, ONLY: kpoint_env_initialize,& 107 kpoint_initialize,& 108 kpoint_initialize_mos 109 USE kpoint_types, ONLY: kpoint_create,& 110 kpoint_type,& 111 read_kpoint_section,& 112 write_kpoint_info 113 USE lri_environment_init, ONLY: lri_env_basis,& 114 lri_env_init 115 USE lri_environment_types, ONLY: lri_environment_type 116 USE machine, ONLY: m_flush 117 USE mathconstants, ONLY: pi 118 USE molecule_kind_types, ONLY: molecule_kind_type,& 119 write_molecule_kind_set 120 USE molecule_types, ONLY: molecule_type 121 USE mp2_setup, ONLY: read_mp2_section 122 USE mp2_types, ONLY: mp2_env_create 123 USE multipole_types, ONLY: do_multipole_none 124 USE orbital_pointers, ONLY: init_orbital_pointers 125 USE orbital_transformation_matrices, ONLY: init_spherical_harmonics 126 USE particle_methods, ONLY: write_particle_distances,& 127 write_qs_particle_coordinates,& 128 write_structure_data 129 USE particle_types, ONLY: particle_type 130 USE qmmm_types_low, ONLY: qmmm_env_qm_type 131 USE qs_dftb_parameters, ONLY: qs_dftb_param_init 132 USE qs_dftb_types, ONLY: qs_dftb_pairpot_type 133 USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init 134 USE qs_dispersion_pairpot, ONLY: qs_dispersion_pairpot_init 135 USE qs_dispersion_types, ONLY: qs_dispersion_type 136 USE qs_dispersion_utils, ONLY: qs_dispersion_env_set,& 137 qs_write_dispersion 138 USE qs_energy_types, ONLY: allocate_qs_energy,& 139 qs_energy_type 140 USE qs_environment_methods, ONLY: qs_env_setup 141 USE qs_environment_types, ONLY: get_qs_env,& 142 qs_environment_type,& 143 set_qs_env 144 USE qs_force_types, ONLY: qs_force_type 145 USE qs_gcp_types, ONLY: qs_gcp_type 146 USE qs_gcp_utils, ONLY: qs_gcp_env_set,& 147 qs_gcp_init 148 USE qs_interactions, ONLY: init_interaction_radii,& 149 init_se_nlradius,& 150 write_core_charge_radii,& 151 write_paw_radii,& 152 write_pgf_orb_radii,& 153 write_ppl_radii,& 154 write_ppnl_radii 155 USE qs_kind_types, ONLY: & 156 check_qs_kind_set, get_qs_kind, get_qs_kind_set, init_gapw_basis_set, init_gapw_nlcc, & 157 init_qs_kind_set, qs_kind_type, set_qs_kind, write_gto_basis_sets, write_qs_kind_set 158 USE qs_ks_types, ONLY: qs_ks_env_create,& 159 qs_ks_env_type,& 160 qs_ks_release,& 161 set_ks_env 162 USE qs_mo_types, ONLY: allocate_mo_set,& 163 mo_set_p_type 164 USE qs_rho0_ggrid, ONLY: rho0_s_grid_create 165 USE qs_rho0_methods, ONLY: init_rho0 166 USE qs_rho0_types, ONLY: rho0_mpole_type 167 USE qs_rho_atom_methods, ONLY: init_rho_atom 168 USE qs_subsys_methods, ONLY: qs_subsys_create 169 USE qs_subsys_types, ONLY: qs_subsys_get,& 170 qs_subsys_release,& 171 qs_subsys_set,& 172 qs_subsys_type 173 USE qs_wf_history_methods, ONLY: wfi_create,& 174 wfi_create_for_kp 175 USE qs_wf_history_types, ONLY: qs_wf_history_type,& 176 wfi_release 177 USE rel_control_types, ONLY: rel_c_create,& 178 rel_c_read_parameters,& 179 rel_c_release,& 180 rel_control_type 181 USE scf_control_types, ONLY: scf_c_create,& 182 scf_c_read_parameters,& 183 scf_c_release,& 184 scf_c_write_parameters,& 185 scf_control_type 186 USE semi_empirical_expns3_methods, ONLY: semi_empirical_expns3_setup 187 USE semi_empirical_int_arrays, ONLY: init_se_intd_array 188 USE semi_empirical_mpole_methods, ONLY: nddo_mpole_setup 189 USE semi_empirical_mpole_types, ONLY: nddo_mpole_type 190 USE semi_empirical_store_int_types, ONLY: semi_empirical_si_create,& 191 semi_empirical_si_type 192 USE semi_empirical_types, ONLY: se_taper_create,& 193 se_taper_type 194 USE semi_empirical_utils, ONLY: se_cutoff_compatible 195 USE transport, ONLY: transport_env_create 196 USE xtb_parameters, ONLY: init_xtb_basis,& 197 xtb_parameters_init,& 198 xtb_parameters_read,& 199 xtb_parameters_set 200 USE xtb_types, ONLY: allocate_xtb_atom_param 201#include "./base/base_uses.f90" 202 203 IMPLICIT NONE 204 205 PRIVATE 206 207 ! *** Global parameters *** 208 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_environment' 209 210 ! *** Public subroutines *** 211 PUBLIC :: qs_init 212 213CONTAINS 214 215! ************************************************************************************************** 216!> \brief Read the input and the database files for the setup of the 217!> QUICKSTEP environment. 218!> \param qs_env ... 219!> \param para_env ... 220!> \param root_section ... 221!> \param globenv ... 222!> \param cp_subsys ... 223!> \param kpoint_env ... 224!> \param cell ... 225!> \param cell_ref ... 226!> \param qmmm ... 227!> \param qmmm_env_qm ... 228!> \param force_env_section ... 229!> \param subsys_section ... 230!> \param use_motion_section ... 231!> \author Creation (22.05.2000,MK) 232! ************************************************************************************************** 233 SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, cell, cell_ref, & 234 qmmm, qmmm_env_qm, force_env_section, subsys_section, use_motion_section) 235 236 TYPE(qs_environment_type), POINTER :: qs_env 237 TYPE(cp_para_env_type), POINTER :: para_env 238 TYPE(section_vals_type), OPTIONAL, POINTER :: root_section 239 TYPE(global_environment_type), OPTIONAL, POINTER :: globenv 240 TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys 241 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoint_env 242 TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref 243 LOGICAL, INTENT(IN), OPTIONAL :: qmmm 244 TYPE(qmmm_env_qm_type), OPTIONAL, POINTER :: qmmm_env_qm 245 TYPE(section_vals_type), POINTER :: force_env_section, subsys_section 246 LOGICAL, INTENT(IN) :: use_motion_section 247 248 CHARACTER(len=*), PARAMETER :: routineN = 'qs_init', routineP = moduleN//':'//routineN 249 250 INTEGER :: ikind, method_id, natom, nkind 251 LOGICAL :: do_admm_rpa, do_et, do_exx, do_hfx, & 252 do_kpoints, is_semi, mp2_present, & 253 my_qmmm, qmmm_decoupl, silent, & 254 use_ref_cell 255 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat 256 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 257 TYPE(cell_type), POINTER :: my_cell, my_cell_ref 258 TYPE(cp_blacs_env_type), POINTER :: blacs_env 259 TYPE(dft_control_type), POINTER :: dft_control 260 TYPE(kpoint_type), POINTER :: kpoints 261 TYPE(lri_environment_type), POINTER :: lri_env 262 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 263 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 264 TYPE(qs_ks_env_type), POINTER :: ks_env 265 TYPE(qs_subsys_type), POINTER :: subsys 266 TYPE(qs_wf_history_type), POINTER :: wf_history 267 TYPE(rel_control_type), POINTER :: rel_control 268 TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, & 269 hfx_section, kpoint_section, & 270 mp2_section, transport_section 271 272 NULLIFY (my_cell, my_cell_ref, atomic_kind_set, particle_set, & 273 qs_kind_set, kpoint_section, dft_section, & 274 subsys, ks_env, dft_control, blacs_env) 275 276 CALL set_qs_env(qs_env, input=force_env_section) 277 IF (.NOT. ASSOCIATED(subsys_section)) THEN 278 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS") 279 END IF 280 281 ! QMMM 282 my_qmmm = .FALSE. 283 IF (PRESENT(qmmm)) my_qmmm = qmmm 284 qmmm_decoupl = .FALSE. 285 IF (PRESENT(qmmm_env_qm)) THEN 286 IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. & 287 qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN 288 ! For GAUSS/SWAVE methods there could be a DDAPC decoupling requested 289 qmmm_decoupl = my_qmmm .AND. qmmm_env_qm%periodic .AND. qmmm_env_qm%multipole 290 END IF 291 qs_env%qmmm_env_qm => qmmm_env_qm 292 END IF 293 CALL set_qs_env(qs_env=qs_env, qmmm=my_qmmm) 294 295 ! Possibly initialize arrays for SE 296 CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=method_id) 297 SELECT CASE (method_id) 298 CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, & 299 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl) 300 CALL init_se_intd_array() 301 is_semi = .TRUE. 302 CASE (do_method_xtb, do_method_dftb) 303 is_semi = .TRUE. 304 CASE DEFAULT 305 is_semi = .FALSE. 306 END SELECT 307 308 CALL qs_subsys_create(subsys, para_env, & 309 force_env_section=force_env_section, & 310 subsys_section=subsys_section, & 311 use_motion_section=use_motion_section, & 312 root_section=root_section, & 313 cp_subsys=cp_subsys, cell=cell, cell_ref=cell_ref, & 314 elkind=is_semi) 315 316 CALL qs_ks_env_create(ks_env) 317 CALL set_ks_env(ks_env, subsys=subsys) 318 CALL set_qs_env(qs_env, ks_env=ks_env) 319 320 CALL qs_subsys_get(subsys, & 321 cell=my_cell, & 322 cell_ref=my_cell_ref, & 323 use_ref_cell=use_ref_cell, & 324 atomic_kind_set=atomic_kind_set, & 325 qs_kind_set=qs_kind_set, & 326 particle_set=particle_set) 327 328 CALL set_ks_env(ks_env, para_env=para_env) 329 IF (PRESENT(globenv)) THEN 330 CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, & 331 globenv%blacs_repeatable) 332 ELSE 333 CALL cp_blacs_env_create(blacs_env, para_env) 334 END IF 335 CALL set_ks_env(ks_env, blacs_env=blacs_env) 336 CALL cp_blacs_env_release(blacs_env) 337 338 ! *** Setup the grids for the G-space Interpolation if any 339 CALL cp_ddapc_ewald_create(qs_env%cp_ddapc_ewald, qmmm_decoupl, my_cell, & 340 force_env_section, subsys_section, para_env) 341 342 ! kpoints 343 IF (PRESENT(kpoint_env)) THEN 344 silent = .TRUE. 345 kpoints => kpoint_env 346 CALL set_qs_env(qs_env=qs_env, kpoints=kpoints) 347 CALL kpoint_initialize(kpoints, particle_set, my_cell) 348 ELSE 349 silent = .FALSE. 350 NULLIFY (kpoints) 351 CALL kpoint_create(kpoints) 352 CALL set_qs_env(qs_env=qs_env, kpoints=kpoints) 353 kpoint_section => section_vals_get_subs_vals(qs_env%input, "DFT%KPOINTS") 354 CALL read_kpoint_section(kpoints, kpoint_section, my_cell%hmat) 355 CALL kpoint_initialize(kpoints, particle_set, my_cell) 356 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") 357 CALL write_kpoint_info(kpoints, dft_section) 358 END IF 359 360 CALL qs_init_subsys(qs_env, para_env, subsys, my_cell, my_cell_ref, use_ref_cell, & 361 subsys_section, silent=silent) 362 363 CALL get_qs_env(qs_env, dft_control=dft_control) 364 IF (method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN 365 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env) 366 CALL lri_env_basis("LRI", qs_env, lri_env, qs_kind_set) 367 ELSE IF (method_id == do_method_rigpw) THEN 368 CALL cp_warn(__LOCATION__, "Experimental code: "// & 369 "RIGPW should only be used for testing.") 370 CALL get_qs_env(qs_env=qs_env, lri_env=lri_env) 371 CALL lri_env_basis("RI", qs_env, lri_env, qs_kind_set) 372 END IF 373 374 ! more kpoint stuff 375 kpoints%para_env => para_env 376 CALL cp_para_env_retain(para_env) 377 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env) 378 kpoints%blacs_env_all => blacs_env 379 CALL cp_blacs_env_retain(blacs_env) 380 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 381 IF (do_kpoints) THEN 382 CALL kpoint_env_initialize(kpoints) 383 CALL kpoint_initialize_mos(kpoints, qs_env%mos) 384 CALL get_qs_env(qs_env=qs_env, wf_history=wf_history) 385 CALL wfi_create_for_kp(wf_history) 386 END IF 387 388 do_hfx = .FALSE. 389 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF") 390 CALL section_vals_get(hfx_section, explicit=do_hfx) 391 CALL get_qs_env(qs_env, dft_control=dft_control) 392 IF (do_hfx) THEN 393 ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization) 394 natom = SIZE(particle_set) 395 CALL hfx_create(qs_env%x_data, para_env, hfx_section, natom, atomic_kind_set, & 396 qs_kind_set, dft_control, my_cell) 397 END IF 398 399 mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION") 400 CALL section_vals_get(mp2_section, explicit=mp2_present) 401 IF (mp2_present) THEN 402 CALL mp2_env_create(qs_env%mp2_env) 403 CALL read_mp2_section(qs_env%input, qs_env%mp2_env) 404 ! create the EXX section if necessary 405 do_exx = .FALSE. 406 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") 407 CALL section_vals_get(hfx_section, explicit=do_exx) 408 IF (do_exx) THEN 409 ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization) 410 natom = SIZE(particle_set) 411 412 ! do_exx in call of hfx_create decides whether to go without ADMM (do_exx=.TRUE.) or with 413 ! ADMM (do_exx=.FALSE.) 414 CALL section_vals_val_get(mp2_section, "RI_RPA%ADMM", l_val=do_admm_rpa) 415 416 CALL hfx_create(qs_env%mp2_env%ri_rpa%x_data, para_env, hfx_section, natom, atomic_kind_set, & 417 qs_kind_set, dft_control, my_cell, do_exx=(.NOT. do_admm_rpa)) 418 419 END IF 420 END IF 421 422 IF (dft_control%qs_control%do_kg) THEN 423 CALL cite_reference(Iannuzzi2006) 424 CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input) 425 END IF 426 427 et_coupling_section => section_vals_get_subs_vals(qs_env%input, & 428 "PROPERTIES%ET_COUPLING") 429 CALL section_vals_get(et_coupling_section, explicit=do_et) 430 IF (do_et) CALL et_coupling_create(qs_env%et_coupling) 431 432 transport_section => section_vals_get_subs_vals(qs_env%input, "DFT%TRANSPORT") 433 CALL section_vals_get(transport_section, explicit=qs_env%do_transport) 434 IF (qs_env%do_transport) THEN 435 CALL transport_env_create(qs_env) 436 END IF 437 438 IF (dft_control%qs_control%do_ls_scf) THEN 439 CALL ls_scf_create(qs_env) 440 ENDIF 441 442 IF (dft_control%qs_control%do_almo_scf) THEN 443 CALL almo_scf_env_create(qs_env) 444 ENDIF 445 446 ! see if we have atomic relativistic corrections 447 CALL get_qs_env(qs_env, rel_control=rel_control) 448 IF (rel_control%rel_method /= rel_none) THEN 449 IF (rel_control%rel_transformation == rel_trans_atom) THEN 450 nkind = SIZE(atomic_kind_set) 451 DO ikind = 1, nkind 452 NULLIFY (rtmat) 453 CALL calculate_atomic_relkin(atomic_kind_set(ikind), qs_kind_set(ikind), rel_control, rtmat) 454 IF (ASSOCIATED(rtmat)) CALL set_qs_kind(qs_kind_set(ikind), reltmat=rtmat) 455 END DO 456 END IF 457 END IF 458 459 CALL qs_subsys_release(subsys) 460 CALL qs_ks_release(ks_env) 461 462 END SUBROUTINE qs_init 463 464! ************************************************************************************************** 465!> \brief Initialize the qs environment (subsys) 466!> \param qs_env ... 467!> \param para_env ... 468!> \param subsys ... 469!> \param cell ... 470!> \param cell_ref ... 471!> \param use_ref_cell ... 472!> \param subsys_section ... 473!> \param silent ... 474!> \author Creation (22.05.2000,MK) 475! ************************************************************************************************** 476 SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell, subsys_section, & 477 silent) 478 479 TYPE(qs_environment_type), POINTER :: qs_env 480 TYPE(cp_para_env_type), POINTER :: para_env 481 TYPE(qs_subsys_type), POINTER :: subsys 482 TYPE(cell_type), POINTER :: cell, cell_ref 483 LOGICAL, INTENT(in) :: use_ref_cell 484 TYPE(section_vals_type), POINTER :: subsys_section 485 LOGICAL, INTENT(in), OPTIONAL :: silent 486 487 CHARACTER(len=*), PARAMETER :: routineN = 'qs_init_subsys', routineP = moduleN//':'//routineN 488 489 CHARACTER(len=2) :: element_symbol 490 INTEGER :: handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, maxlppl, & 491 maxlppnl, method_id, multiplicity, my_ival, n_ao, n_ao_aux_fit, n_mo_add, natom, & 492 nelectron, ngauss, nkind, output_unit, tnadd_method 493 INTEGER, DIMENSION(2) :: n_mo, nelectron_spin 494 LOGICAL :: all_potential_present, be_silent, & 495 do_kpoints, e1terms, has_unit_metric, & 496 lribas, orb_gradient, was_present 497 REAL(dp) :: alpha, ccore, ewald_rcut, maxocc, & 498 verlet_skin, zeff_correction 499 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 500 TYPE(cp_logger_type), POINTER :: logger 501 TYPE(dft_control_type), POINTER :: dft_control 502 TYPE(dftb_control_type), POINTER :: dftb_control 503 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 504 TYPE(ewald_environment_type), POINTER :: ewald_env 505 TYPE(ewald_pw_type), POINTER :: ewald_pw 506 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env 507 TYPE(gapw_control_type), POINTER :: gapw_control 508 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, & 509 ri_hfx_basis, ri_xas_basis, & 510 tmp_basis_set 511 TYPE(lri_environment_type), POINTER :: lri_env 512 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 513 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 514 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 515 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 516 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 517 TYPE(qs_control_type), POINTER :: qs_control 518 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & 519 POINTER :: dftb_potential 520 TYPE(qs_dispersion_type), POINTER :: dispersion_env 521 TYPE(qs_energy_type), POINTER :: energy 522 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 523 TYPE(qs_gcp_type), POINTER :: gcp_env 524 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 525 TYPE(qs_kind_type), POINTER :: qs_kind 526 TYPE(qs_ks_env_type), POINTER :: ks_env 527 TYPE(qs_wf_history_type), POINTER :: wf_history 528 TYPE(rel_control_type), POINTER :: rel_control 529 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 530 TYPE(scf_control_type), POINTER :: scf_control 531 TYPE(se_taper_type), POINTER :: se_taper 532 TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, & 533 ewald_section, lri_section, nl_section, poisson_section, pp_section, print_section, & 534 qs_section, se_section, tddfpt_section, xc_section, xtb_section 535 TYPE(semi_empirical_control_type), POINTER :: se_control 536 TYPE(semi_empirical_si_type), POINTER :: se_store_int_env 537 TYPE(xtb_control_type), POINTER :: xtb_control 538 539 CALL timeset(routineN, handle) 540 NULLIFY (logger) 541 logger => cp_get_default_logger() 542 output_unit = cp_logger_get_default_io_unit(logger) 543 was_present = .FALSE. 544 545 be_silent = .FALSE. 546 IF (PRESENT(silent)) be_silent = silent 547 548 ! Initialise the Quickstep environment 549 NULLIFY (mos, se_taper, mos_aux_fit) 550 NULLIFY (dft_control) 551 NULLIFY (energy) 552 NULLIFY (force) 553 NULLIFY (local_molecules) 554 NULLIFY (local_particles) 555 NULLIFY (scf_control) 556 NULLIFY (dft_section) 557 NULLIFY (et_coupling_section) 558 NULLIFY (ks_env) 559 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") 560 qs_section => section_vals_get_subs_vals(dft_section, "QS") 561 et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING") 562 ! reimplemented TDDFPT 563 tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT") 564 565 CALL qs_subsys_get(subsys, particle_set=particle_set, & 566 qs_kind_set=qs_kind_set, & 567 atomic_kind_set=atomic_kind_set, & 568 molecule_set=molecule_set, & 569 molecule_kind_set=molecule_kind_set) 570 571 ! *** Read the input section with the DFT control parameters *** 572 CALL read_dft_control(dft_control, dft_section) 573 574 ! original implementation of TDDFPT 575 IF (dft_control%do_tddfpt_calculation) THEN 576 CALL read_tddfpt_control(dft_control%tddfpt_control, dft_section) 577 END IF 578 579 ! *** Print the Quickstep program banner (copyright and version number) *** 580 IF (.NOT. be_silent) THEN 581 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PROGRAM_BANNER", extension=".Log") 582 CALL section_vals_val_get(qs_section, "METHOD", i_val=method_id) 583 SELECT CASE (method_id) 584 CASE DEFAULT 585 CALL qs_header(iw) 586 CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, & 587 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl) 588 CALL se_header(iw) 589 CASE (do_method_dftb) 590 CALL dftb_header(iw) 591 CASE (do_method_xtb) 592 CALL xtb_header(iw) 593 END SELECT 594 CALL cp_print_key_finished_output(iw, logger, dft_section, & 595 "PRINT%PROGRAM_BANNER") 596 END IF 597 598 ! *** Read the input section with the Quickstep control parameters *** 599 CALL read_qs_section(dft_control%qs_control, qs_section) 600 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 601 IF (do_kpoints) THEN 602 ! reset some of the settings for wfn extrapolation for kpoints 603 SELECT CASE (dft_control%qs_control%wf_interpolation_method_nr) 604 CASE (wfi_linear_wf_method_nr, wfi_linear_ps_method_nr, wfi_ps_method_nr, wfi_aspc_nr) 605 dft_control%qs_control%wf_interpolation_method_nr = wfi_use_guess_method_nr 606 END SELECT 607 END IF 608 609 ! ******* check if any kind of electron transfer calculation has to be performed 610 CALL section_vals_val_get(et_coupling_section, "TYPE_OF_CONSTRAINT", i_val=my_ival) 611 dft_control%qs_control%et_coupling_calc = .FALSE. 612 IF (my_ival == do_et_ddapc) THEN 613 et_ddapc_section => section_vals_get_subs_vals(et_coupling_section, "DDAPC_RESTRAINT_A") 614 dft_control%qs_control%et_coupling_calc = .TRUE. 615 dft_control%qs_control%ddapc_restraint = .TRUE. 616 CALL read_ddapc_section(dft_control%qs_control, ddapc_restraint_section=et_ddapc_section) 617 ENDIF 618 619 CALL read_mgrid_section(dft_control%qs_control, dft_section) 620 621 ! reimplemented TDDFPT 622 CALL read_tddfpt2_control(dft_control%tddfpt2_control, tddfpt_section, dft_control%qs_control) 623 624 ! Create relativistic control section 625 CALL rel_c_create(rel_control) 626 CALL rel_c_read_parameters(rel_control, dft_section) 627 CALL set_qs_env(qs_env, rel_control=rel_control) 628 CALL rel_c_release(rel_control) 629 630 ! *** Read DFTB parameter files *** 631 IF (dft_control%qs_control%method_id == do_method_dftb) THEN 632 NULLIFY (ewald_env, ewald_pw, dftb_potential) 633 dftb_control => dft_control%qs_control%dftb_control 634 CALL qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, & 635 subsys_section=subsys_section, para_env=para_env) 636 CALL set_qs_env(qs_env, dftb_potential=dftb_potential) 637 ! check for Ewald 638 IF (dftb_control%do_ewald) THEN 639 CALL ewald_env_create(ewald_env, para_env) 640 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON") 641 CALL ewald_env_set(ewald_env, poisson_section=poisson_section) 642 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD") 643 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION") 644 CALL get_qs_kind_set(qs_kind_set, basis_rcut=ewald_rcut) 645 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat) 646 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section) 647 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw) 648 CALL ewald_env_release(ewald_env) 649 CALL ewald_pw_release(ewald_pw) 650 END IF 651 ELSEIF (dft_control%qs_control%method_id == do_method_xtb) THEN 652 ! *** Read xTB parameter file *** 653 xtb_control => dft_control%qs_control%xtb_control 654 NULLIFY (ewald_env, ewald_pw) 655 CALL get_qs_env(qs_env, nkind=nkind) 656 DO ikind = 1, nkind 657 qs_kind => qs_kind_set(ikind) 658 ! Setup proper xTB parameters 659 CPASSERT(.NOT. ASSOCIATED(qs_kind%xtb_parameter)) 660 CALL allocate_xtb_atom_param(qs_kind%xtb_parameter) 661 ! Set default parameters 662 CALL get_qs_kind(qs_kind, element_symbol=element_symbol) 663 CALL xtb_parameters_init(qs_kind%xtb_parameter, element_symbol, & 664 xtb_control%parameter_file_path, xtb_control%parameter_file_name, & 665 para_env) 666 ! Read specific parameters from input 667 xtb_section => section_vals_get_subs_vals(qs_section, "xTB") 668 CALL xtb_parameters_read(qs_kind%xtb_parameter, element_symbol, xtb_section) 669 ! set dependent parameters 670 CALL xtb_parameters_set(qs_kind%xtb_parameter) 671 ! Generate basis set 672 NULLIFY (tmp_basis_set) 673 IF (qs_kind%xtb_parameter%z == 1) THEN 674 ! special case hydrogen 675 ngauss = xtb_control%h_sto_ng 676 ELSE 677 ngauss = xtb_control%sto_ng 678 END IF 679 CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss) 680 CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB") 681 ! potential 682 zeff_correction = 0.0_dp 683 CALL init_potential(qs_kind%all_potential, itype="BARE", & 684 zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction) 685 CALL get_potential(qs_kind%all_potential, alpha_core_charge=alpha) 686 ccore = qs_kind%xtb_parameter%zeff*SQRT((alpha/pi)**3) 687 CALL set_potential(qs_kind%all_potential, ccore_charge=ccore) 688 qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff - zeff_correction 689 END DO 690 ! 691 ! check for Ewald 692 IF (xtb_control%do_ewald) THEN 693 CALL ewald_env_create(ewald_env, para_env) 694 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON") 695 CALL ewald_env_set(ewald_env, poisson_section=poisson_section) 696 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD") 697 print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION") 698 CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat) 699 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section) 700 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw) 701 CALL ewald_env_release(ewald_env) 702 CALL ewald_pw_release(ewald_pw) 703 END IF 704 END IF 705 706 ! lri or ri env initialization 707 lri_section => section_vals_get_subs_vals(qs_section, "LRIGPW") 708 IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. & 709 dft_control%qs_control%lri_optbas .OR. & 710 dft_control%qs_control%method_id == do_method_rigpw) THEN 711 CALL lri_env_init(lri_env, lri_section) 712 CALL set_qs_env(qs_env, lri_env=lri_env) 713 END IF 714 715 ! *** Check basis and fill in missing parts *** 716 CALL check_qs_kind_set(qs_kind_set, dft_control, subsys_section=subsys_section) 717 718 ! *** Check that no all-electron potential is present if GPW or GAPW_XC 719 CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present) 720 IF ((dft_control%qs_control%method == "GPW") .OR. & 721 (dft_control%qs_control%method == "GAPW_XC") .OR. & 722 (dft_control%qs_control%method == "OFGPW")) THEN 723 IF (all_potential_present) THEN 724 CPABORT("All-el calculations with GPW, GAPW_XC, and OFGPW are not implemented ") 725 END IF 726 END IF 727 728 ! DFT+U 729 CALL get_qs_kind_set(qs_kind_set, dft_plus_u_atom_present=dft_control%dft_plus_u) 730 731 IF (dft_control%do_admm) THEN 732 ! Check if ADMM basis is available, auto-generate if needed 733 CALL get_qs_env(qs_env, nkind=nkind) 734 DO ikind = 1, nkind 735 NULLIFY (aux_fit_basis) 736 qs_kind => qs_kind_set(ikind) 737 CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT") 738 IF (.NOT. (ASSOCIATED(aux_fit_basis))) THEN 739 ! AUX_FIT basis set is not yet loaded 740 CALL cp_warn(__LOCATION__, "Automatic Generation of AUX_FIT basis. "// & 741 "This is experimental code.") 742 ! Generate a default basis 743 CALL create_aux_fit_basis_set(aux_fit_basis, qs_kind, dft_control%auto_basis_aux_fit) 744 CALL add_basis_set_to_container(qs_kind%basis_sets, aux_fit_basis, "AUX_FIT") 745 END IF 746 END DO 747 END IF 748 749 lribas = .FALSE. 750 e1terms = .FALSE. 751 IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN 752 lribas = .TRUE. 753 CALL get_qs_env(qs_env, lri_env=lri_env) 754 e1terms = lri_env%exact_1c_terms 755 END IF 756 IF (dft_control%qs_control%do_kg) THEN 757 CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method) 758 IF (tnadd_method == kg_tnadd_embed_ri) lribas = .TRUE. 759 END IF 760 IF (lribas) THEN 761 ! Check if LRI_AUX basis is available, auto-generate if needed 762 CALL get_qs_env(qs_env, nkind=nkind) 763 DO ikind = 1, nkind 764 NULLIFY (lri_aux_basis) 765 qs_kind => qs_kind_set(ikind) 766 CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX") 767 IF (.NOT. (ASSOCIATED(lri_aux_basis))) THEN 768 ! LRI_AUX basis set is not yet loaded 769 CALL cp_warn(__LOCATION__, "Automatic Generation of LRI_AUX basis. "// & 770 "This is experimental code.") 771 ! Generate a default basis 772 CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, dft_control%auto_basis_lri_aux, e1terms) 773 CALL add_basis_set_to_container(qs_kind%basis_sets, lri_aux_basis, "LRI_AUX") 774 END IF 775 END DO 776 END IF 777 778 IF (dft_control%qs_control%method_id == do_method_rigpw) THEN 779 ! Check if RI_HXC basis is available, auto-generate if needed 780 CALL get_qs_env(qs_env, nkind=nkind) 781 DO ikind = 1, nkind 782 NULLIFY (ri_hfx_basis) 783 qs_kind => qs_kind_set(ikind) 784 CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HXC") 785 IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN 786 ! RI_HXC basis set is not yet loaded 787 CALL cp_warn(__LOCATION__, "Automatic Generation of RI_HXC basis. "// & 788 "This is experimental code.") 789 ! Generate a default basis 790 CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hxc) 791 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HXC") 792 END IF 793 END DO 794 END IF 795 796 IF (dft_control%do_xas_tdp_calculation) THEN 797 ! Check if RI_XAS basis is given, auto-generate if not 798 CALL get_qs_env(qs_env, nkind=nkind) 799 DO ikind = 1, nkind 800 NULLIFY (ri_xas_basis) 801 qs_kind => qs_kind_set(ikind) 802 CALL get_qs_kind(qs_kind, basis_Set=ri_xas_basis, basis_type="RI_XAS") 803 IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN 804 CALL cp_warn(__LOCATION__, "Automatic Generation of RI_XAS basis. "// & 805 "This is experimental code.") 806 ! Generate a default basis 807 CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas) 808 CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS") 809 END IF 810 END DO 811 END IF 812 813 ! *** Initialize the spherical harmonics and *** 814 ! *** the orbital transformation matrices *** 815 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl) 816 817 lmax_sphere = dft_control%qs_control%gapw_control%lmax_sphere 818 IF (lmax_sphere .LT. 0) THEN 819 lmax_sphere = 2*maxlgto 820 dft_control%qs_control%gapw_control%lmax_sphere = lmax_sphere 821 END IF 822 IF (dft_control%qs_control%method == "LRIGPW" .OR. dft_control%qs_control%lri_optbas) THEN 823 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="LRI_AUX") 824 !take maxlgto from lri basis if larger (usually) 825 maxlgto = MAX(maxlgto, maxlgto_lri) 826 ELSE IF (dft_control%qs_control%method == "RIGPW") THEN 827 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC") 828 maxlgto = MAX(maxlgto, maxlgto_lri) 829 END IF 830 IF (dft_control%do_xas_tdp_calculation) THEN 831 !done as a precaution 832 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS") 833 maxlgto = MAX(maxlgto, maxlgto_lri) 834 END IF 835 maxl = MAX(2*maxlgto, maxlppl, maxlppnl, lmax_sphere) + 1 836 837 CALL init_orbital_pointers(maxl) 838 839 CALL init_spherical_harmonics(maxl, 0) 840 841 ! *** Initialise the qs_kind_set *** 842 CALL init_qs_kind_set(qs_kind_set) 843 844 ! *** Initialise GAPW soft basis and projectors 845 IF (dft_control%qs_control%method == "GAPW" .OR. & 846 dft_control%qs_control%method == "GAPW_XC") THEN 847 qs_control => dft_control%qs_control 848 CALL init_gapw_basis_set(qs_kind_set, qs_control, qs_env%input) 849 ENDIF 850 851 ! *** Initialize the pretabulation for the calculation of the *** 852 ! *** incomplete Gamma function F_n(t) after McMurchie-Davidson *** 853 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto) 854 maxl = MAX(3*maxlgto + 1, 0) 855 CALL init_md_ftable(maxl) 856 857 ! *** Initialize the atomic interaction radii *** 858 CALL init_interaction_radii(dft_control%qs_control, atomic_kind_set, qs_kind_set) 859 ! 860 IF (dft_control%qs_control%method_id == do_method_xtb) THEN 861 ! cutoff radius 862 DO ikind = 1, nkind 863 qs_kind => qs_kind_set(ikind) 864 CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set) 865 CALL get_gto_basis_set(tmp_basis_set, kind_radius=qs_kind%xtb_parameter%rcut) 866 END DO 867 END IF 868 869 IF (.NOT. be_silent) THEN 870 CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section) 871 CALL write_pgf_orb_radii("aux", atomic_kind_set, qs_kind_set, subsys_section) 872 CALL write_pgf_orb_radii("lri", atomic_kind_set, qs_kind_set, subsys_section) 873 CALL write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section) 874 CALL write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section) 875 CALL write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section) 876 CALL write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section) 877 END IF 878 879 ! *** Distribute molecules and atoms using the new data structures *** 880 CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, & 881 particle_set=particle_set, & 882 local_particles=local_particles, & 883 molecule_kind_set=molecule_kind_set, & 884 molecule_set=molecule_set, & 885 local_molecules=local_molecules, & 886 force_env_section=qs_env%input) 887 888 ! *** SCF parameters *** 889 CALL scf_c_create(scf_control) 890 CALL scf_c_read_parameters(scf_control, dft_section) 891 892 ! *** Allocate the data structure for Quickstep energies *** 893 CALL allocate_qs_energy(energy) 894 895 ! check for orthogonal basis 896 has_unit_metric = .FALSE. 897 IF (dft_control%qs_control%semi_empirical) THEN 898 IF (dft_control%qs_control%se_control%orthogonal_basis) has_unit_metric = .TRUE. 899 END IF 900 IF (dft_control%qs_control%dftb) THEN 901 IF (dft_control%qs_control%dftb_control%orthogonal_basis) has_unit_metric = .TRUE. 902 END IF 903 CALL set_qs_env(qs_env, has_unit_metric=has_unit_metric) 904 905 ! *** Activate the interpolation *** 906 CALL wfi_create(wf_history, & 907 interpolation_method_nr= & 908 dft_control%qs_control%wf_interpolation_method_nr, & 909 extrapolation_order=dft_control%qs_control%wf_extrapolation_order, & 910 has_unit_metric=has_unit_metric) 911 912 ! *** Set the current Quickstep environment *** 913 CALL set_qs_env(qs_env=qs_env, & 914 scf_control=scf_control, & 915 wf_history=wf_history) 916 917 CALL qs_subsys_set(subsys, & 918 cell_ref=cell_ref, & 919 use_ref_cell=use_ref_cell, & 920 energy=energy, & 921 force=force) 922 923 CALL get_qs_env(qs_env, ks_env=ks_env) 924 CALL set_ks_env(ks_env, dft_control=dft_control) 925 926 CALL qs_subsys_set(subsys, local_molecules=local_molecules, & 927 local_particles=local_particles, cell=cell) 928 929 CALL distribution_1d_release(local_particles) 930 CALL distribution_1d_release(local_molecules) 931 CALL scf_c_release(scf_control) 932 CALL wfi_release(wf_history) 933 CALL dft_control_release(dft_control) 934 935 CALL get_qs_env(qs_env=qs_env, & 936 atomic_kind_set=atomic_kind_set, & 937 dft_control=dft_control, & 938 scf_control=scf_control) 939 940 ! decide what conditions need mo_derivs 941 ! right now, this only appears to be OT 942 IF (dft_control%qs_control%do_ls_scf .OR. & 943 dft_control%qs_control%do_almo_scf) THEN 944 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.) 945 ELSE 946 IF (scf_control%use_ot) THEN 947 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.TRUE.) 948 ELSE 949 CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.) 950 ENDIF 951 ENDIF 952 953 ! XXXXXXX this is backwards XXXXXXXX 954 dft_control%smear = scf_control%smear%do_smear 955 956 ! Periodic efield needs equal occupation and orbital gradients 957 IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)) THEN 958 IF (dft_control%apply_period_efield) THEN 959 CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient) 960 IF (.NOT. orb_gradient) THEN 961 CALL cp_abort(__LOCATION__, "Periodic Efield needs orbital gradient and direct optimization."// & 962 " Use the OT optimization method.") 963 END IF 964 IF (dft_control%smear) THEN 965 CALL cp_abort(__LOCATION__, "Periodic Efield needs equal occupation numbers."// & 966 " Smearing option is not possible.") 967 END IF 968 END IF 969 END IF 970 971 ! Initialize the GAPW local densities and potentials 972 IF (dft_control%qs_control%method_id == do_method_gapw .OR. & 973 dft_control%qs_control%method_id == do_method_gapw_xc) THEN 974 ! *** Allocate and initialize the set of atomic densities *** 975 gapw_control => dft_control%qs_control%gapw_control 976 CALL init_rho_atom(qs_env, gapw_control) 977 IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN 978 CALL get_qs_env(qs_env=qs_env, natom=natom) 979 ! *** Allocate and initialize the compensation density rho0 *** 980 CALL init_rho0(qs_env, gapw_control) 981 ! *** Allocate and Initialize the local coulomb term *** 982 CALL init_coulomb_local(qs_env%hartree_local, natom) 983 END IF 984 ! NLCC 985 CALL init_gapw_nlcc(qs_kind_set) 986 ELSE IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN 987 ! allocate local ri environment 988 ! nothing to do here? 989 ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN 990 ! allocate ri environment 991 ! nothing to do here? 992 ELSE IF (dft_control%qs_control%semi_empirical) THEN 993 NULLIFY (se_store_int_env, se_nddo_mpole, se_nonbond_env) 994 natom = SIZE(particle_set) 995 se_section => section_vals_get_subs_vals(qs_section, "SE") 996 se_control => dft_control%qs_control%se_control 997 998 ! Make the cutoff radii choice a bit smarter 999 CALL se_cutoff_compatible(se_control, se_section, cell, output_unit) 1000 1001 SELECT CASE (dft_control%qs_control%method_id) 1002 CASE DEFAULT 1003 CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pm3, & 1004 do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl) 1005 ! Neighbor lists have to be MAX(interaction range, orbital range) 1006 ! set new kind radius 1007 CALL init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, subsys_section) 1008 END SELECT 1009 ! Initialize to zero the max multipole to treat in the EWALD scheme.. 1010 se_control%max_multipole = do_multipole_none 1011 ! check for Ewald 1012 IF (se_control%do_ewald .OR. se_control%do_ewald_gks) THEN 1013 CALL ewald_env_create(ewald_env, para_env) 1014 poisson_section => section_vals_get_subs_vals(dft_section, "POISSON") 1015 CALL ewald_env_set(ewald_env, poisson_section=poisson_section) 1016 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD") 1017 print_section => section_vals_get_subs_vals(qs_env%input, & 1018 "PRINT%GRID_INFORMATION") 1019 CALL read_ewald_section(ewald_env, ewald_section) 1020 ! Create ewald grids 1021 CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, & 1022 print_section=print_section) 1023 ! Initialize ewald grids 1024 CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat) 1025 ! Setup the nonbond environment (real space part of Ewald) 1026 CALL ewald_env_get(ewald_env, rcut=ewald_rcut) 1027 ! Setup the maximum level of multipoles to be treated in the periodic SE scheme 1028 IF (se_control%do_ewald) THEN 1029 CALL ewald_env_get(ewald_env, max_multipole=se_control%max_multipole) 1030 ENDIF 1031 CALL section_vals_val_get(se_section, "NEIGHBOR_LISTS%VERLET_SKIN", & 1032 r_val=verlet_skin) 1033 CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, & 1034 do_nonbonded=.TRUE., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, & 1035 ei_scale14=0.0_dp, vdw_scale14=0.0_dp, shift_cutoff=.FALSE.) 1036 ! Create and Setup NDDO multipole environment 1037 CALL nddo_mpole_setup(se_nddo_mpole, natom) 1038 CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw, & 1039 se_nonbond_env=se_nonbond_env, se_nddo_mpole=se_nddo_mpole) 1040 CALL ewald_env_release(ewald_env) 1041 CALL ewald_pw_release(ewald_pw) 1042 ! Handle the residual integral part 1/R^3 1043 CALL semi_empirical_expns3_setup(qs_kind_set, se_control, & 1044 dft_control%qs_control%method_id) 1045 END IF 1046 ! Taper function 1047 CALL se_taper_create(se_taper, se_control%integral_screening, se_control%do_ewald, & 1048 se_control%taper_cou, se_control%range_cou, & 1049 se_control%taper_exc, se_control%range_exc, & 1050 se_control%taper_scr, se_control%range_scr, & 1051 se_control%taper_lrc, se_control%range_lrc) 1052 CALL set_qs_env(qs_env, se_taper=se_taper) 1053 ! Store integral environment 1054 CALL semi_empirical_si_create(se_store_int_env, se_section) 1055 CALL set_qs_env(qs_env, se_store_int_env=se_store_int_env) 1056 ENDIF 1057 1058 ! Initialize possible dispersion parameters 1059 IF (dft_control%qs_control%method_id == do_method_gpw .OR. & 1060 dft_control%qs_control%method_id == do_method_gapw .OR. & 1061 dft_control%qs_control%method_id == do_method_gapw_xc .OR. & 1062 dft_control%qs_control%method_id == do_method_lrigpw .OR. & 1063 dft_control%qs_control%method_id == do_method_rigpw .OR. & 1064 dft_control%qs_control%method_id == do_method_ofgpw) THEN 1065 ALLOCATE (dispersion_env) 1066 NULLIFY (xc_section) 1067 xc_section => section_vals_get_subs_vals(dft_section, "XC") 1068 CALL qs_dispersion_env_set(dispersion_env, xc_section) 1069 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN 1070 NULLIFY (pp_section) 1071 pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL") 1072 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env) 1073 ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN 1074 NULLIFY (nl_section) 1075 nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL") 1076 CALL qs_dispersion_nonloc_init(dispersion_env, para_env) 1077 END IF 1078 CALL set_qs_env(qs_env, dispersion_env=dispersion_env) 1079 ELSE IF (dft_control%qs_control%method_id == do_method_dftb) THEN 1080 ALLOCATE (dispersion_env) 1081 ! set general defaults 1082 dispersion_env%doabc = .FALSE. 1083 dispersion_env%c9cnst = .FALSE. 1084 dispersion_env%lrc = .FALSE. 1085 dispersion_env%srb = .FALSE. 1086 dispersion_env%verbose = .FALSE. 1087 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, & 1088 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist) 1089 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, & 1090 dispersion_env%d2y_dx2, dispersion_env%dftd_section) 1091 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn) 1092 IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3) THEN 1093 dispersion_env%type = xc_vdw_fun_pairpot 1094 dispersion_env%pp_type = vdw_pairpot_dftd3 1095 dispersion_env%eps_cn = dftb_control%epscn 1096 dispersion_env%s6 = dftb_control%sd3(1) 1097 dispersion_env%sr6 = dftb_control%sd3(2) 1098 dispersion_env%s8 = dftb_control%sd3(3) 1099 dispersion_env%domol = .FALSE. 1100 dispersion_env%kgc8 = 0._dp 1101 dispersion_env%rc_disp = dftb_control%rcdisp 1102 dispersion_env%exp_pre = 0._dp 1103 dispersion_env%scaling = 0._dp 1104 dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file 1105 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env) 1106 ELSE 1107 dispersion_env%type = xc_vdw_fun_none 1108 END IF 1109 CALL set_qs_env(qs_env, dispersion_env=dispersion_env) 1110 ELSE IF (dft_control%qs_control%method_id == do_method_xtb) THEN 1111 ALLOCATE (dispersion_env) 1112 ! set general defaults 1113 dispersion_env%doabc = .FALSE. 1114 dispersion_env%c9cnst = .FALSE. 1115 dispersion_env%lrc = .FALSE. 1116 dispersion_env%srb = .FALSE. 1117 dispersion_env%verbose = .FALSE. 1118 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, & 1119 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist) 1120 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, & 1121 dispersion_env%d2y_dx2, dispersion_env%dftd_section) 1122 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn) 1123 dispersion_env%type = xc_vdw_fun_pairpot 1124 dispersion_env%pp_type = vdw_pairpot_dftd3bj 1125 dispersion_env%eps_cn = xtb_control%epscn 1126 dispersion_env%s6 = xtb_control%s6 1127 dispersion_env%s8 = xtb_control%s8 1128 dispersion_env%a1 = xtb_control%a1 1129 dispersion_env%a2 = xtb_control%a2 1130 dispersion_env%domol = .FALSE. 1131 dispersion_env%kgc8 = 0._dp 1132 dispersion_env%rc_disp = xtb_control%rcdisp 1133 dispersion_env%exp_pre = 0._dp 1134 dispersion_env%scaling = 0._dp 1135 dispersion_env%parameter_file_name = xtb_control%dispersion_parameter_file 1136 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env) 1137 CALL set_qs_env(qs_env, dispersion_env=dispersion_env) 1138 ELSE IF (dft_control%qs_control%semi_empirical) THEN 1139 ALLOCATE (dispersion_env) 1140 ! set general defaults 1141 dispersion_env%doabc = .FALSE. 1142 dispersion_env%c9cnst = .FALSE. 1143 dispersion_env%lrc = .FALSE. 1144 dispersion_env%srb = .FALSE. 1145 dispersion_env%verbose = .FALSE. 1146 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, & 1147 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist) 1148 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, & 1149 dispersion_env%d2y_dx2, dispersion_env%dftd_section) 1150 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn) 1151 IF (se_control%dispersion) THEN 1152 dispersion_env%type = xc_vdw_fun_pairpot 1153 dispersion_env%pp_type = vdw_pairpot_dftd3 1154 dispersion_env%eps_cn = se_control%epscn 1155 dispersion_env%s6 = se_control%sd3(1) 1156 dispersion_env%sr6 = se_control%sd3(2) 1157 dispersion_env%s8 = se_control%sd3(3) 1158 dispersion_env%domol = .FALSE. 1159 dispersion_env%kgc8 = 0._dp 1160 dispersion_env%rc_disp = se_control%rcdisp 1161 dispersion_env%exp_pre = 0._dp 1162 dispersion_env%scaling = 0._dp 1163 dispersion_env%parameter_file_name = se_control%dispersion_parameter_file 1164 CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env) 1165 ELSE 1166 dispersion_env%type = xc_vdw_fun_none 1167 END IF 1168 CALL set_qs_env(qs_env, dispersion_env=dispersion_env) 1169 END IF 1170 1171 ! Initialize possible geomertical counterpoise correction potential 1172 IF (dft_control%qs_control%method_id == do_method_gpw .OR. & 1173 dft_control%qs_control%method_id == do_method_gapw .OR. & 1174 dft_control%qs_control%method_id == do_method_gapw_xc .OR. & 1175 dft_control%qs_control%method_id == do_method_lrigpw .OR. & 1176 dft_control%qs_control%method_id == do_method_rigpw .OR. & 1177 dft_control%qs_control%method_id == do_method_ofgpw) THEN 1178 ALLOCATE (gcp_env) 1179 NULLIFY (xc_section) 1180 xc_section => section_vals_get_subs_vals(dft_section, "XC") 1181 CALL qs_gcp_env_set(gcp_env, xc_section) 1182 CALL qs_gcp_init(qs_env, gcp_env) 1183 CALL set_qs_env(qs_env, gcp_env=gcp_env) 1184 END IF 1185 1186 ! *** Allocate the MO data types *** 1187 CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron) 1188 1189 ! the total number of electrons 1190 nelectron = nelectron - dft_control%charge 1191 1192 IF (dft_control%multiplicity == 0) THEN 1193 IF (MODULO(nelectron, 2) == 0) THEN 1194 dft_control%multiplicity = 1 1195 ELSE 1196 dft_control%multiplicity = 2 1197 END IF 1198 END IF 1199 1200 multiplicity = dft_control%multiplicity 1201 1202 IF ((dft_control%nspins < 1) .OR. (dft_control%nspins > 2)) THEN 1203 CPABORT("nspins should be 1 or 2 for the time being ...") 1204 END IF 1205 1206 IF ((MODULO(nelectron, 2) /= 0) .AND. (dft_control%nspins == 1)) THEN 1207 IF (.NOT. dft_control%qs_control%ofgpw .AND. .NOT. dft_control%smear) THEN 1208 CPABORT("Use the LSD option for an odd number of electrons") 1209 END IF 1210 END IF 1211 1212 ! The transition potential method to calculate XAS needs LSD 1213 IF (dft_control%do_xas_calculation) THEN 1214 IF (dft_control%nspins == 1) THEN 1215 CPABORT("Use the LSD option for XAS with transition potential") 1216 END IF 1217 END IF 1218 1219 ! assigning the number of states per spin initial version, not yet very 1220 ! general. Should work for an even number of electrons and a single 1221 ! additional electron this set of options that requires full matrices, 1222 ! however, makes things a bit ugly right now.... we try to make a 1223 ! distinction between the number of electrons per spin and the number of 1224 ! MOs per spin this should allow the use of fractional occupations later 1225 ! on 1226 IF (dft_control%qs_control%ofgpw) THEN 1227 1228 IF (dft_control%nspins == 1) THEN 1229 maxocc = nelectron 1230 nelectron_spin(1) = nelectron 1231 nelectron_spin(2) = 0 1232 n_mo(1) = 1 1233 n_mo(2) = 0 1234 ELSE 1235 IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN 1236 CPABORT("LSD: try to use a different multiplicity") 1237 END IF 1238 nelectron_spin(1) = (nelectron + multiplicity - 1)/2 1239 nelectron_spin(2) = (nelectron - multiplicity + 1)/2 1240 IF (nelectron_spin(1) < 0) THEN 1241 CPABORT("LSD: too few electrons for this multiplicity") 1242 END IF 1243 maxocc = MAXVAL(nelectron_spin) 1244 n_mo(1) = MIN(nelectron_spin(1), 1) 1245 n_mo(2) = MIN(nelectron_spin(2), 1) 1246 END IF 1247 1248 ELSE 1249 1250 IF (dft_control%nspins == 1) THEN 1251 maxocc = 2.0_dp 1252 nelectron_spin(1) = nelectron 1253 nelectron_spin(2) = 0 1254 IF (MODULO(nelectron, 2) == 0) THEN 1255 n_mo(1) = nelectron/2 1256 ELSE 1257 n_mo(1) = INT(nelectron/2._dp) + 1 1258 END IF 1259 n_mo(2) = 0 1260 ELSE 1261 maxocc = 1.0_dp 1262 1263 ! The simplist spin distribution is written here. Special cases will 1264 ! need additional user input 1265 IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN 1266 CPABORT("LSD: try to use a different multiplicity") 1267 END IF 1268 1269 nelectron_spin(1) = (nelectron + multiplicity - 1)/2 1270 nelectron_spin(2) = (nelectron - multiplicity + 1)/2 1271 1272 IF (nelectron_spin(2) < 0) THEN 1273 CPABORT("LSD: too few electrons for this multiplicity") 1274 END IF 1275 1276 n_mo(1) = nelectron_spin(1) 1277 n_mo(2) = nelectron_spin(2) 1278 1279 END IF 1280 1281 END IF 1282 1283 ! store the number of electrons once an for all 1284 CALL qs_subsys_set(subsys, & 1285 nelectron_total=nelectron, & 1286 nelectron_spin=nelectron_spin) 1287 1288 ! Check and set number of added (unoccupied) MOs 1289 IF (dft_control%nspins == 2) THEN 1290 IF (scf_control%added_mos(2) > 0) THEN 1291 n_mo_add = scf_control%added_mos(2) 1292 ELSE 1293 n_mo_add = scf_control%added_mos(1) 1294 END IF 1295 IF (n_mo_add > n_ao - n_mo(2)) & 1296 CPWARN("More added MOs requested for beta spin than available.") 1297 scf_control%added_mos(2) = MIN(n_mo_add, n_ao - n_mo(2)) 1298 n_mo(2) = n_mo(2) + scf_control%added_mos(2) 1299 END IF 1300 1301 ! proceed alpha orbitals after the beta orbitals; this is essential to avoid 1302 ! reduction in the number of available unoccupied molecular orbitals. 1303 ! E.g. n_ao = 10, nelectrons = 10, multiplicity = 3 implies n_mo(1) = 6, n_mo(2) = 4; 1304 ! added_mos(1:2) = (6,undef) should increase the number of molecular orbitals as 1305 ! n_mo(1) = min(n_ao, n_mo(1) + added_mos(1)) = 10, n_mo(2) = 10. 1306 ! However, if we try to proceed alpha orbitals first, this leads us n_mo(1:2) = (10,8) 1307 ! due to the following assignment instruction above: 1308 ! IF (scf_control%added_mos(2) > 0) THEN ... ELSE; n_mo_add = scf_control%added_mos(1); END IF 1309 IF (scf_control%added_mos(1) > n_ao - n_mo(1)) & 1310 CALL cp_warn(__LOCATION__, & 1311 "More added MOs requested than available. "// & 1312 "The full set of unoccupied MOs will be used.") 1313 scf_control%added_mos(1) = MIN(scf_control%added_mos(1), n_ao - n_mo(1)) 1314 n_mo(1) = n_mo(1) + scf_control%added_mos(1) 1315 1316 IF (dft_control%nspins == 2) THEN 1317 IF (n_mo(2) > n_mo(1)) & 1318 CALL cp_warn(__LOCATION__, & 1319 "More beta than alpha MOs requested. "// & 1320 "The number of beta MOs will be reduced to the number alpha MOs.") 1321 n_mo(2) = MIN(n_mo(1), n_mo(2)) 1322 CPASSERT(n_mo(1) >= nelectron_spin(1)) 1323 CPASSERT(n_mo(2) >= nelectron_spin(2)) 1324 END IF 1325 1326 ! kpoints 1327 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 1328 IF (do_kpoints .AND. dft_control%nspins == 2) THEN 1329 ! we need equal number of calculated states 1330 IF (n_mo(2) /= n_mo(1)) & 1331 CALL cp_warn(__LOCATION__, & 1332 "Kpoints: Different number of MOs requested. "// & 1333 "The number of beta MOs will be set to the number alpha MOs.") 1334 n_mo(2) = n_mo(1) 1335 CPASSERT(n_mo(1) >= nelectron_spin(1)) 1336 CPASSERT(n_mo(2) >= nelectron_spin(2)) 1337 END IF 1338 1339 ! Compatibility checks for smearing 1340 IF (scf_control%smear%do_smear) THEN 1341 IF (scf_control%added_mos(1) == 0) THEN 1342 CPABORT("Extra MOs (ADDED_MOS) are required for smearing") 1343 END IF 1344 END IF 1345 1346 ! *** Some options require that all MOs are computed ... *** 1347 IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & 1348 "PRINT%MO/CARTESIAN"), & 1349 cp_p_file) .OR. & 1350 (scf_control%level_shift /= 0.0_dp) .OR. & 1351 (scf_control%diagonalization%eps_jacobi /= 0.0_dp) .OR. & 1352 (dft_control%roks .AND. (.NOT. scf_control%use_ot))) THEN 1353 n_mo(:) = n_ao 1354 END IF 1355 1356 ! Compatibility checks for ROKS 1357 IF (dft_control%roks .AND. (.NOT. scf_control%use_ot)) THEN 1358 IF (scf_control%roks_scheme == general_roks) & 1359 CPWARN("General ROKS scheme is not yet tested!") 1360 1361 IF (scf_control%smear%do_smear) THEN 1362 CALL cp_abort(__LOCATION__, & 1363 "The options ROKS and SMEAR are not compatible. "// & 1364 "Try UKS instead of ROKS") 1365 END IF 1366 END IF 1367 1368 ! in principle the restricted calculation could be performed 1369 ! using just one set of MOs and special casing most of the code 1370 ! right now we'll just take care of what is effectively an additional constraint 1371 ! at as few places as possible, just duplicating the beta orbitals 1372 IF (dft_control%restricted .AND. (output_unit > 0)) THEN 1373 ! it is really not yet tested till the end ! Joost 1374 WRITE (output_unit, *) "" 1375 WRITE (output_unit, *) " **************************************" 1376 WRITE (output_unit, *) " restricted calculation cutting corners" 1377 WRITE (output_unit, *) " experimental feature, check code " 1378 WRITE (output_unit, *) " **************************************" 1379 ENDIF 1380 1381 ! no point in allocating these things here ? 1382 IF (dft_control%qs_control%do_ls_scf) THEN 1383 NULLIFY (mos) 1384 ELSE 1385 ALLOCATE (mos(dft_control%nspins)) 1386 DO ispin = 1, dft_control%nspins 1387 NULLIFY (mos(ispin)%mo_set) 1388 CALL allocate_mo_set(mo_set=mos(ispin)%mo_set, & 1389 nao=n_ao, & 1390 nmo=n_mo(ispin), & 1391 nelectron=nelectron_spin(ispin), & 1392 n_el_f=REAL(nelectron_spin(ispin), dp), & 1393 maxocc=maxocc, & 1394 flexible_electron_count=dft_control%relax_multiplicity) 1395 END DO 1396 END IF 1397 1398 CALL set_qs_env(qs_env, mos=mos) 1399 1400 ! If we use auxiliary density matrix methods , set mo_set_aux_fit 1401 IF (dft_control%do_admm) THEN 1402 ! 1403 CALL get_qs_kind_set(qs_kind_set, nelectron=nelectron, nsgf=n_ao_aux_fit, & 1404 basis_type="AUX_FIT") 1405 ALLOCATE (mos_aux_fit(dft_control%nspins)) 1406 1407 DO ispin = 1, dft_control%nspins 1408 NULLIFY (mos_aux_fit(ispin)%mo_set) 1409 CALL allocate_mo_set(mo_set=mos_aux_fit(ispin)%mo_set, & 1410 nao=n_ao_aux_fit, & 1411 nmo=n_mo(ispin), & 1412 nelectron=nelectron_spin(ispin), & 1413 n_el_f=REAL(nelectron_spin(ispin), dp), & 1414 maxocc=maxocc, & 1415 flexible_electron_count=dft_control%relax_multiplicity) 1416 END DO 1417 CALL set_qs_env(qs_env, mos_aux_fit=mos_aux_fit) 1418 END IF 1419 1420 IF (.NOT. be_silent) THEN 1421 ! Print the DFT control parameters 1422 CALL write_dft_control(dft_control, dft_section) 1423 1424 ! Print the vdW control parameters 1425 IF (dft_control%qs_control%method_id == do_method_gpw .OR. & 1426 dft_control%qs_control%method_id == do_method_gapw .OR. & 1427 dft_control%qs_control%method_id == do_method_gapw_xc .OR. & 1428 dft_control%qs_control%method_id == do_method_lrigpw .OR. & 1429 dft_control%qs_control%method_id == do_method_rigpw .OR. & 1430 dft_control%qs_control%method_id == do_method_dftb .OR. & 1431 dft_control%qs_control%method_id == do_method_xtb .OR. & 1432 dft_control%qs_control%method_id == do_method_ofgpw) THEN 1433 CALL get_qs_env(qs_env, dispersion_env=dispersion_env) 1434 CALL qs_write_dispersion(qs_env, dispersion_env) 1435 END IF 1436 1437 ! Print the Quickstep control parameters 1438 CALL write_qs_control(dft_control%qs_control, dft_section) 1439 1440 ! Print XES/XAS control parameters 1441 IF (dft_control%do_xas_calculation) THEN 1442 CALL cite_reference(Iannuzzi2007) 1443 !CALL write_xas_control(dft_control%xas_control,dft_section) 1444 END IF 1445 1446 ! Print the unnormalized basis set information (input data) 1447 CALL write_gto_basis_sets(qs_kind_set, subsys_section) 1448 1449 ! Print the atomic kind set 1450 CALL write_qs_kind_set(qs_kind_set, subsys_section) 1451 1452 ! Print the molecule kind set 1453 CALL write_molecule_kind_set(molecule_kind_set, subsys_section) 1454 1455 ! Print the total number of kinds, atoms, basis functions etc. 1456 CALL write_total_numbers(qs_kind_set, particle_set, qs_env%input) 1457 1458 ! Print the atomic coordinates 1459 CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="QUICKSTEP") 1460 1461 ! Print the interatomic distances 1462 CALL write_particle_distances(particle_set, cell, subsys_section) 1463 1464 ! Print the requested structure data 1465 CALL write_structure_data(particle_set, cell, subsys_section) 1466 1467 ! Print symmetry information 1468 CALL write_symmetry(particle_set, cell, subsys_section) 1469 1470 ! Print the SCF parameters 1471 IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. & 1472 (.NOT. dft_control%qs_control%do_almo_scf)) THEN 1473 CALL scf_c_write_parameters(scf_control, dft_section) 1474 ENDIF 1475 ENDIF 1476 1477 ! Sets up pw_env, qs_charges, mpools ... 1478 CALL qs_env_setup(qs_env) 1479 1480 ! Allocate and Initialie rho0 soft on the global grid 1481 IF (dft_control%qs_control%method == "GAPW") THEN 1482 CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole) 1483 CALL rho0_s_grid_create(qs_env, rho0_mpole) 1484 END IF 1485 1486 IF (output_unit > 0) CALL m_flush(output_unit) 1487 CALL timestop(handle) 1488 1489 END SUBROUTINE qs_init_subsys 1490 1491! ************************************************************************************************** 1492!> \brief Write the total number of kinds, atoms, etc. to the logical unit 1493!> number lunit. 1494!> \param qs_kind_set ... 1495!> \param particle_set ... 1496!> \param force_env_section ... 1497!> \author Creation (06.10.2000) 1498! ************************************************************************************************** 1499 SUBROUTINE write_total_numbers(qs_kind_set, particle_set, force_env_section) 1500 1501 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1502 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1503 TYPE(section_vals_type), POINTER :: force_env_section 1504 1505 CHARACTER(len=*), PARAMETER :: routineN = 'write_total_numbers', & 1506 routineP = moduleN//':'//routineN 1507 1508 INTEGER :: maxlgto, maxlppl, maxlppnl, natom, ncgf, & 1509 nkind, npgf, nset, nsgf, nshell, & 1510 output_unit 1511 TYPE(cp_logger_type), POINTER :: logger 1512 1513 NULLIFY (logger) 1514 logger => cp_get_default_logger() 1515 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%TOTAL_NUMBERS", & 1516 extension=".Log") 1517 1518 IF (output_unit > 0) THEN 1519 natom = SIZE(particle_set) 1520 nkind = SIZE(qs_kind_set) 1521 1522 CALL get_qs_kind_set(qs_kind_set, & 1523 maxlgto=maxlgto, & 1524 ncgf=ncgf, & 1525 npgf=npgf, & 1526 nset=nset, & 1527 nsgf=nsgf, & 1528 nshell=nshell, & 1529 maxlppl=maxlppl, & 1530 maxlppnl=maxlppnl) 1531 1532 WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") & 1533 "TOTAL NUMBERS AND MAXIMUM NUMBERS" 1534 1535 IF (nset + npgf + ncgf > 0) THEN 1536 WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") & 1537 "Total number of", & 1538 "- Atomic kinds: ", nkind, & 1539 "- Atoms: ", natom, & 1540 "- Shell sets: ", nset, & 1541 "- Shells: ", nshell, & 1542 "- Primitive Cartesian functions: ", npgf, & 1543 "- Cartesian basis functions: ", ncgf, & 1544 "- Spherical basis functions: ", nsgf 1545 ELSE IF (nshell + nsgf > 0) THEN 1546 WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") & 1547 "Total number of", & 1548 "- Atomic kinds: ", nkind, & 1549 "- Atoms: ", natom, & 1550 "- Shells: ", nshell, & 1551 "- Spherical basis functions: ", nsgf 1552 ELSE 1553 WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") & 1554 "Total number of", & 1555 "- Atomic kinds: ", nkind, & 1556 "- Atoms: ", natom 1557 END IF 1558 1559 IF ((maxlppl > -1) .AND. (maxlppnl > -1)) THEN 1560 WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") & 1561 "Maximum angular momentum of the", & 1562 "- Orbital basis functions: ", maxlgto, & 1563 "- Local part of the GTH pseudopotential: ", maxlppl, & 1564 "- Non-local part of the GTH pseudopotential: ", maxlppnl 1565 ELSEIF (maxlppl > -1) THEN 1566 WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") & 1567 "Maximum angular momentum of the", & 1568 "- Orbital basis functions: ", maxlgto, & 1569 "- Local part of the GTH pseudopotential: ", maxlppl 1570 ELSE 1571 WRITE (UNIT=output_unit, FMT="(/,T3,A,T75,I6)") & 1572 "Maximum angular momentum of the orbital basis functions: ", maxlgto 1573 END IF 1574 1575 ! LRI_AUX BASIS 1576 CALL get_qs_kind_set(qs_kind_set, & 1577 maxlgto=maxlgto, & 1578 ncgf=ncgf, & 1579 npgf=npgf, & 1580 nset=nset, & 1581 nsgf=nsgf, & 1582 nshell=nshell, & 1583 basis_type="LRI_AUX") 1584 IF (nset + npgf + ncgf > 0) THEN 1585 WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") & 1586 "LRI_AUX Basis: ", & 1587 "Total number of", & 1588 "- Shell sets: ", nset, & 1589 "- Shells: ", nshell, & 1590 "- Primitive Cartesian functions: ", npgf, & 1591 "- Cartesian basis functions: ", ncgf, & 1592 "- Spherical basis functions: ", nsgf 1593 WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") & 1594 " Maximum angular momentum ", maxlgto 1595 END IF 1596 1597 ! RI_HXC BASIS 1598 CALL get_qs_kind_set(qs_kind_set, & 1599 maxlgto=maxlgto, & 1600 ncgf=ncgf, & 1601 npgf=npgf, & 1602 nset=nset, & 1603 nsgf=nsgf, & 1604 nshell=nshell, & 1605 basis_type="RI_HXC") 1606 IF (nset + npgf + ncgf > 0) THEN 1607 WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") & 1608 "RI_HXC Basis: ", & 1609 "Total number of", & 1610 "- Shell sets: ", nset, & 1611 "- Shells: ", nshell, & 1612 "- Primitive Cartesian functions: ", npgf, & 1613 "- Cartesian basis functions: ", ncgf, & 1614 "- Spherical basis functions: ", nsgf 1615 WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") & 1616 " Maximum angular momentum ", maxlgto 1617 END IF 1618 1619 ! AUX_FIT BASIS 1620 CALL get_qs_kind_set(qs_kind_set, & 1621 maxlgto=maxlgto, & 1622 ncgf=ncgf, & 1623 npgf=npgf, & 1624 nset=nset, & 1625 nsgf=nsgf, & 1626 nshell=nshell, & 1627 basis_type="AUX_FIT") 1628 IF (nset + npgf + ncgf > 0) THEN 1629 WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") & 1630 "AUX_FIT ADMM-Basis: ", & 1631 "Total number of", & 1632 "- Shell sets: ", nset, & 1633 "- Shells: ", nshell, & 1634 "- Primitive Cartesian functions: ", npgf, & 1635 "- Cartesian basis functions: ", ncgf, & 1636 "- Spherical basis functions: ", nsgf 1637 WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") & 1638 " Maximum angular momentum ", maxlgto 1639 END IF 1640 1641 END IF 1642 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, & 1643 "PRINT%TOTAL_NUMBERS") 1644 1645 END SUBROUTINE write_total_numbers 1646 1647END MODULE qs_environment 1648