1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utilities to set up the control types 8! ************************************************************************************************** 9MODULE cp_control_utils 10 USE bibliography, ONLY: & 11 Andreussi2012, Dewar1977, Dewar1985, Elstner1998, Fattebert2002, Grimme2017, Hu2007, & 12 Krack2000, Lippert1997, Lippert1999, Porezag1995, Repasky2002, Rocha2006, Schenter2008, & 13 Seifert1996, Souza2002, Stengel2009, Stewart1989, Stewart2007, Thiel1992, Umari2002, & 14 VanVoorhis2015, VandeVondele2005a, VandeVondele2005b, Zhechkov2005, cite_reference 15 USE cp_control_types, ONLY: & 16 admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, & 17 dft_control_create, dft_control_type, efield_type, qs_control_type, sccs_control_create, & 18 tddfpt2_control_type, tddfpt_control_create, tddfpt_control_type, xtb_control_type 19 USE cp_files, ONLY: close_file,& 20 open_file 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type 23 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 24 cp_print_key_unit_nr 25 USE cp_units, ONLY: cp_unit_from_cp2k,& 26 cp_unit_to_cp2k 27 USE input_constants, ONLY: & 28 constant_env, custom_env, do_admm_basis_projection, do_admm_blocked_projection, & 29 do_admm_blocking_purify_full, do_admm_charge_constrained_projection, & 30 do_admm_exch_scaling_merlot, do_admm_purify_mcweeny, do_admm_purify_mo_diag, & 31 do_admm_purify_mo_no_diag, do_admm_purify_none, do_admm_purify_none_dm, & 32 do_ddapc_constraint, do_ddapc_restraint, do_method_am1, do_method_dftb, do_method_gapw, & 33 do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, & 34 do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, do_method_pm6fm, & 35 do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, do_pwgrid_ns_fullspace, & 36 do_pwgrid_ns_halfspace, do_pwgrid_spherical, do_s2_constraint, do_s2_restraint, & 37 do_se_is_kdso, do_se_is_kdso_d, do_se_is_slater, do_se_lr_ewald, do_se_lr_ewald_gks, & 38 do_se_lr_ewald_r3, do_se_lr_none, gaussian_env, numerical, ramp_env, & 39 real_time_propagation, sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, & 40 sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, sic_ad, sic_eo, & 41 sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, slater, & 42 tddfpt_excitations, use_mom_ref_user, xas_dip_len 43 USE input_cp2k_check, ONLY: xc_functionals_expand 44 USE input_cp2k_dft, ONLY: create_dft_section,& 45 create_qs_section 46 USE input_enumeration_types, ONLY: enum_i2c,& 47 enumeration_type 48 USE input_keyword_types, ONLY: keyword_get,& 49 keyword_type 50 USE input_section_types, ONLY: & 51 section_get_ival, section_get_keyword, section_release, section_type, section_vals_get, & 52 section_vals_get_subs_vals, section_vals_type, section_vals_val_get, section_vals_val_set 53 USE kinds, ONLY: default_path_length,& 54 default_string_length,& 55 dp 56 USE qs_cdft_utils, ONLY: read_cdft_control_section 57 USE string_utilities, ONLY: uppercase 58 USE util, ONLY: sort 59 USE xc_derivatives, ONLY: xc_functionals_get_needs 60 USE xc_input_constants, ONLY: xc_deriv_collocate 61 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 62 USE xc_write_output, ONLY: xc_write 63#include "./base/base_uses.f90" 64 65 IMPLICIT NONE 66 67 PRIVATE 68 69 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_control_utils' 70 71 PUBLIC :: read_dft_control, & 72 read_mgrid_section, & 73 read_qs_section, & 74 read_tddfpt_control, & 75 read_tddfpt2_control, & 76 write_dft_control, & 77 write_qs_control, & 78 read_ddapc_section 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief ... 83!> \param dft_control ... 84!> \param dft_section ... 85! ************************************************************************************************** 86 SUBROUTINE read_dft_control(dft_control, dft_section) 87 TYPE(dft_control_type), POINTER :: dft_control 88 TYPE(section_vals_type), POINTER :: dft_section 89 90 CHARACTER(len=*), PARAMETER :: routineN = 'read_dft_control', & 91 routineP = moduleN//':'//routineN 92 93 CHARACTER(len=default_path_length) :: basis_set_file_name, potential_file_name 94 CHARACTER(LEN=default_string_length), & 95 DIMENSION(:), POINTER :: tmpstringlist 96 INTEGER :: excitations, irep, isize, method_id, & 97 nrep, xc_deriv_method_id 98 LOGICAL :: do_ot, do_rtp, explicit, is_present, & 99 l_param, not_SE, was_present 100 REAL(KIND=dp) :: density_cut, gradient_cut, tau_cut 101 REAL(KIND=dp), DIMENSION(:), POINTER :: pol 102 TYPE(cp_logger_type), POINTER :: logger 103 TYPE(section_vals_type), POINTER :: sccs_section, scf_section, tmp_section, & 104 xc_fun_section, xc_section 105 TYPE(xc_rho_cflags_type) :: needs 106 107 was_present = .FALSE. 108 109 logger => cp_get_default_logger() 110 111 NULLIFY (tmp_section, xc_fun_section, xc_section) 112 CALL dft_control_create(dft_control) 113 ! determine wheather this is a semiempirical or DFTB run 114 ! --> (no XC section needs to be provided) 115 not_SE = .TRUE. 116 CALL section_vals_val_get(dft_section, "QS%METHOD", i_val=method_id) 117 SELECT CASE (method_id) 118 CASE (do_method_dftb, do_method_xtb, do_method_mndo, do_method_am1, do_method_pm3, do_method_pnnl, & 119 do_method_pm6, do_method_pm6fm, do_method_pdg, do_method_rm1, do_method_mndod) 120 not_SE = .FALSE. 121 END SELECT 122 ! Check for XC section and XC_FUNCTIONAL section 123 xc_section => section_vals_get_subs_vals(dft_section, "XC") 124 CALL section_vals_get(xc_section, explicit=is_present) 125 IF (.NOT. is_present .AND. not_SE) THEN 126 CPABORT("XC section missing.") 127 END IF 128 IF (is_present) THEN 129 CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut) 130 CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut) 131 CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut) 132 ! Perform numerical stability checks and possibly correct the issues 133 IF (density_cut <= EPSILON(0.0_dp)*100.0_dp) & 134 CALL cp_warn(__LOCATION__, & 135 "DENSITY_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// & 136 "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ") 137 density_cut = MAX(EPSILON(0.0_dp)*100.0_dp, density_cut) 138 IF (gradient_cut <= EPSILON(0.0_dp)*100.0_dp) & 139 CALL cp_warn(__LOCATION__, & 140 "GRADIENT_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// & 141 "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ") 142 gradient_cut = MAX(EPSILON(0.0_dp)*100.0_dp, gradient_cut) 143 IF (tau_cut <= EPSILON(0.0_dp)*100.0_dp) & 144 CALL cp_warn(__LOCATION__, & 145 "TAU_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// & 146 "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ") 147 tau_cut = MAX(EPSILON(0.0_dp)*100.0_dp, tau_cut) 148 CALL section_vals_val_set(xc_section, "density_cutoff", r_val=density_cut) 149 CALL section_vals_val_set(xc_section, "gradient_cutoff", r_val=gradient_cut) 150 CALL section_vals_val_set(xc_section, "tau_cutoff", r_val=tau_cut) 151 END IF 152 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 153 CALL section_vals_get(xc_fun_section, explicit=is_present) 154 IF (.NOT. is_present .AND. not_SE) THEN 155 CPABORT("XC_FUNCTIONAL section missing.") 156 END IF 157 scf_section => section_vals_get_subs_vals(dft_section, "SCF") 158 CALL section_vals_val_get(dft_section, "UKS", l_val=dft_control%uks) 159 CALL section_vals_val_get(dft_section, "ROKS", l_val=dft_control%roks) 160 IF (dft_control%uks .OR. dft_control%roks) THEN 161 dft_control%nspins = 2 162 ELSE 163 dft_control%nspins = 1 164 END IF 165 166 dft_control%lsd = (dft_control%nspins > 1) 167 needs = xc_functionals_get_needs(xc_fun_section, & 168 lsd=dft_control%lsd, & 169 add_basic_components=.TRUE.) 170 dft_control%use_kinetic_energy_density = (needs%tau_spin .OR. needs%tau) 171 172 xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV") 173 dft_control%drho_by_collocation = (needs%norm_drho .AND. (xc_deriv_method_id == xc_deriv_collocate)) 174 IF (dft_control%drho_by_collocation) THEN 175 CPABORT("derivatives by collocation not implemented") 176 END IF 177 178 ! Automatic auxiliary basis set generation 179 CALL section_vals_val_get(dft_section, "AUTO_BASIS", n_rep_val=nrep) 180 DO irep = 1, nrep 181 CALL section_vals_val_get(dft_section, "AUTO_BASIS", i_rep_val=irep, c_vals=tmpstringlist) 182 IF (SIZE(tmpstringlist) == 2) THEN 183 CALL uppercase(tmpstringlist(2)) 184 SELECT CASE (tmpstringlist(2)) 185 CASE ("X") 186 isize = -1 187 CASE ("SMALL") 188 isize = 0 189 CASE ("MEDIUM") 190 isize = 1 191 CASE ("LARGE") 192 isize = 2 193 CASE ("HUGE") 194 isize = 3 195 CASE DEFAULT 196 CPWARN("Unknown basis size in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1))) 197 END SELECT 198 ! 199 SELECT CASE (tmpstringlist(1)) 200 CASE ("X") 201 CASE ("RI_AUX") 202 dft_control%auto_basis_ri_aux = isize 203 CASE ("AUX_FIT") 204 dft_control%auto_basis_aux_fit = isize 205 CASE ("LRI_AUX") 206 dft_control%auto_basis_lri_aux = isize 207 CASE ("RI_HXC") 208 dft_control%auto_basis_ri_hxc = isize 209 CASE ("RI_XAS") 210 dft_control%auto_basis_ri_xas = isize 211 CASE DEFAULT 212 CPWARN("Unknown basis type in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1))) 213 END SELECT 214 ELSE 215 CALL cp_abort(__LOCATION__, & 216 "AUTO_BASIS keyword in &DFT section has a wrong number of arguments.") 217 END IF 218 END DO 219 220 !! check if we do wavefunction fitting 221 tmp_section => section_vals_get_subs_vals(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD") 222 CALL section_vals_get(tmp_section, explicit=is_present) 223 dft_control%do_admm = is_present 224 dft_control%do_admm_mo = .FALSE. 225 dft_control%do_admm_dm = .FALSE. 226 IF (is_present) THEN 227 do_ot = .FALSE. 228 CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=do_ot) 229 CALL admm_control_create(dft_control%admm_control) 230 231 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EPS_FILTER", & 232 r_val=dft_control%admm_control%eps_filter) 233 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", i_val=method_id) 234 dft_control%admm_control%purification_method = method_id 235 236 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%METHOD", i_val=method_id) 237 dft_control%admm_control%method = method_id 238 239 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_SCALING_MODEL", i_val=method_id) 240 dft_control%admm_control%scaling_model = method_id 241 242 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_CORRECTION_FUNC", i_val=method_id) 243 dft_control%admm_control%aux_exch_func = method_id 244 245 ! parameters for X functional 246 dft_control%admm_control%aux_exch_func_param = .FALSE. 247 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A1", explicit=explicit, & 248 r_val=dft_control%admm_control%aux_x_param(1)) 249 IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE. 250 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A2", explicit=explicit, & 251 r_val=dft_control%admm_control%aux_x_param(2)) 252 IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE. 253 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_GAMMA", explicit=explicit, & 254 r_val=dft_control%admm_control%aux_x_param(3)) 255 IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE. 256 257 CALL read_admm_block_list(dft_control%admm_control, dft_section) 258 259 ! In the case of charge-constrained projection (e.g. according to Merlot), 260 ! there is no purification needed and hence, do_admm_purify_none has to be set. 261 262 IF ((dft_control%admm_control%method == do_admm_blocking_purify_full .OR. & 263 dft_control%admm_control%method == do_admm_blocked_projection) & 264 .AND. dft_control%admm_control%scaling_model == do_admm_exch_scaling_merlot) THEN 265 CPABORT("ADMM: Blocking and Merlot scaling are mutually exclusive.") 266 END IF 267 268 IF (dft_control%admm_control%method == do_admm_charge_constrained_projection .AND. & 269 dft_control%admm_control%purification_method /= do_admm_purify_none) THEN 270 CALL cp_abort(__LOCATION__, & 271 "ADMM: In the case of METHOD=CHARGE_CONSTRAINED_PROJECTION, "// & 272 "ADMM_PURIFICATION_METHOD=NONE has to be set.") 273 END IF 274 275 IF (dft_control%admm_control%purification_method == do_admm_purify_mo_diag .OR. & 276 dft_control%admm_control%purification_method == do_admm_purify_mo_no_diag) THEN 277 IF (dft_control%admm_control%method /= do_admm_basis_projection) & 278 CPABORT("ADMM: Chosen purification requires BASIS_PROJECTION") 279 280 IF (.NOT. do_ot) CPABORT("ADMM: MO-based purification requires OT.") 281 END IF 282 283 IF (dft_control%admm_control%purification_method == do_admm_purify_none_dm .OR. & 284 dft_control%admm_control%purification_method == do_admm_purify_mcweeny) THEN 285 dft_control%do_admm_dm = .TRUE. 286 ELSE 287 dft_control%do_admm_mo = .TRUE. 288 ENDIF 289 END IF 290 291 ! Set restricted to true, if both OT and ROKS are requested 292 !MK in principle dft_control%restricted could be dropped completely like the 293 !MK input key by using only dft_control%roks now 294 CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=l_param) 295 dft_control%restricted = (dft_control%roks .AND. l_param) 296 297 CALL section_vals_val_get(dft_section, "CHARGE", i_val=dft_control%charge) 298 CALL section_vals_val_get(dft_section, "MULTIPLICITY", i_val=dft_control%multiplicity) 299 CALL section_vals_val_get(dft_section, "RELAX_MULTIPLICITY", r_val=dft_control%relax_multiplicity) 300 IF (dft_control%relax_multiplicity > 0.0_dp) THEN 301 IF (.NOT. dft_control%uks) & 302 CALL cp_abort(__LOCATION__, "The option RELAX_MULTIPLICITY is only valid for "// & 303 "unrestricted Kohn-Sham (UKS) calculations") 304 END IF 305 306 ! check for the presence of the low spin roks section 307 tmp_section => section_vals_get_subs_vals(dft_section, "LOW_SPIN_ROKS") 308 CALL section_vals_get(tmp_section, explicit=dft_control%low_spin_roks) 309 310 dft_control%sic_method_id = sic_none 311 dft_control%sic_scaling_a = 1.0_dp 312 dft_control%sic_scaling_b = 1.0_dp 313 314 ! DFT+U 315 dft_control%dft_plus_u = .FALSE. 316 CALL section_vals_val_get(dft_section, "PLUS_U_METHOD", i_val=method_id) 317 dft_control%plus_u_method_id = method_id 318 319 ! Smearing in use 320 dft_control%smear = .FALSE. 321 322 ! Surface dipole correction 323 dft_control%correct_surf_dip = .FALSE. 324 CALL section_vals_val_get(dft_section, "SURFACE_DIPOLE_CORRECTION", l_val=dft_control%correct_surf_dip) 325 CALL section_vals_val_get(dft_section, "SURF_DIP_DIR", i_val=dft_control%dir_surf_dip) 326 327 CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", & 328 c_val=basis_set_file_name) 329 CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", & 330 c_val=potential_file_name) 331 332 ! Read the input section 333 tmp_section => section_vals_get_subs_vals(dft_section, "sic") 334 CALL section_vals_val_get(tmp_section, "SIC_METHOD", & 335 i_val=dft_control%sic_method_id) 336 CALL section_vals_val_get(tmp_section, "ORBITAL_SET", & 337 i_val=dft_control%sic_list_id) 338 CALL section_vals_val_get(tmp_section, "SIC_SCALING_A", & 339 r_val=dft_control%sic_scaling_a) 340 CALL section_vals_val_get(tmp_section, "SIC_SCALING_B", & 341 r_val=dft_control%sic_scaling_b) 342 343 ! Determine if this is a TDDFPT run 344 CALL section_vals_val_get(dft_section, "EXCITATIONS", i_val=excitations) 345 dft_control%do_tddfpt_calculation = (excitations == tddfpt_excitations) 346 IF (dft_control%do_tddfpt_calculation) THEN 347 CALL tddfpt_control_create(dft_control%tddfpt_control) 348 END IF 349 350 do_rtp = .FALSE. 351 tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION") 352 CALL section_vals_get(tmp_section, explicit=is_present) 353 IF (is_present) THEN 354 CALL read_rtp_section(dft_control, tmp_section) 355 do_rtp = .TRUE. 356 END IF 357 358 ! Read the input section 359 tmp_section => section_vals_get_subs_vals(dft_section, "XAS") 360 CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_calculation) 361 IF (dft_control%do_xas_calculation) THEN 362 ! Override with section parameter 363 CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", & 364 l_val=dft_control%do_xas_calculation) 365 END IF 366 367 tmp_section => section_vals_get_subs_vals(dft_section, "XAS_TDP") 368 CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_tdp_calculation) 369 IF (dft_control%do_xas_tdp_calculation) THEN 370 ! Override with section parameter 371 CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", & 372 l_val=dft_control%do_xas_tdp_calculation) 373 END IF 374 375 ! Read the finite field input section 376 dft_control%apply_efield = .FALSE. 377 dft_control%apply_efield_field = .FALSE. !this is for RTP 378 tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD") 379 CALL section_vals_get(tmp_section, n_repetition=nrep, explicit=is_present) 380 IF (is_present) THEN 381 ALLOCATE (dft_control%efield_fields(nrep)) 382 CALL read_efield_sections(dft_control, tmp_section) 383 IF (do_rtp) THEN 384 dft_control%apply_efield_field = .TRUE. 385 ELSE 386 dft_control%apply_efield = .TRUE. 387 CPASSERT(nrep == 1) 388 END IF 389 END IF 390 391 ! Read the finite field input section for periodic fields 392 tmp_section => section_vals_get_subs_vals(dft_section, "PERIODIC_EFIELD") 393 CALL section_vals_get(tmp_section, explicit=dft_control%apply_period_efield) 394 IF (dft_control%apply_period_efield) THEN 395 ALLOCATE (dft_control%period_efield) 396 CALL section_vals_val_get(tmp_section, "POLARISATION", r_vals=pol) 397 dft_control%period_efield%polarisation(1:3) = pol(1:3) 398 CALL section_vals_val_get(tmp_section, "D_FILTER", r_vals=pol) 399 dft_control%period_efield%d_filter(1:3) = pol(1:3) 400 CALL section_vals_val_get(tmp_section, "INTENSITY", & 401 r_val=dft_control%period_efield%strength) 402 dft_control%period_efield%displacement_field = .FALSE. 403 CALL section_vals_val_get(tmp_section, "DISPLACEMENT_FIELD", & 404 l_val=dft_control%period_efield%displacement_field) 405 ! periodic fields don't work with RTP 406 CPASSERT(.NOT. do_rtp) 407 IF (dft_control%period_efield%displacement_field) THEN 408 CALL cite_reference(Stengel2009) 409 ELSE 410 CALL cite_reference(Souza2002) 411 CALL cite_reference(Umari2002) 412 END IF 413 END IF 414 415 ! Read the external potential input section 416 tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_POTENTIAL") 417 CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_potential) 418 419 !MK Read the SCCS input section if present 420 sccs_section => section_vals_get_subs_vals(dft_section, "SCCS") 421 CALL section_vals_get(sccs_section, explicit=is_present) 422 IF (is_present) THEN 423 ! Check section parameter if SCCS is activated 424 CALL section_vals_val_get(sccs_section, "_SECTION_PARAMETERS_", & 425 l_val=dft_control%do_sccs) 426 IF (dft_control%do_sccs) THEN 427 CALL sccs_control_create(dft_control%sccs_control) 428 CALL section_vals_val_get(sccs_section, "ALPHA", & 429 r_val=dft_control%sccs_control%alpha_solvent) 430 CALL section_vals_val_get(sccs_section, "BETA", & 431 r_val=dft_control%sccs_control%beta_solvent) 432 CALL section_vals_val_get(sccs_section, "DELTA_RHO", & 433 r_val=dft_control%sccs_control%delta_rho) 434 CALL section_vals_val_get(sccs_section, "DERIVATIVE_METHOD", & 435 i_val=dft_control%sccs_control%derivative_method) 436 CALL section_vals_val_get(sccs_section, "METHOD", & 437 i_val=dft_control%sccs_control%method_id) 438 CALL section_vals_val_get(sccs_section, "DIELECTRIC_CONSTANT", & 439 r_val=dft_control%sccs_control%epsilon_solvent) 440 CALL section_vals_val_get(sccs_section, "EPS_SCCS", & 441 r_val=dft_control%sccs_control%eps_sccs) 442 CALL section_vals_val_get(sccs_section, "EPS_SCF", & 443 r_val=dft_control%sccs_control%eps_scf) 444 CALL section_vals_val_get(sccs_section, "GAMMA", & 445 r_val=dft_control%sccs_control%gamma_solvent) 446 CALL section_vals_val_get(sccs_section, "MAX_ITER", & 447 i_val=dft_control%sccs_control%max_iter) 448 CALL section_vals_val_get(sccs_section, "MIXING", & 449 r_val=dft_control%sccs_control%mixing) 450 SELECT CASE (dft_control%sccs_control%method_id) 451 CASE (sccs_andreussi) 452 tmp_section => section_vals_get_subs_vals(sccs_section, "ANDREUSSI") 453 CALL section_vals_val_get(tmp_section, "RHO_MAX", & 454 r_val=dft_control%sccs_control%rho_max) 455 CALL section_vals_val_get(tmp_section, "RHO_MIN", & 456 r_val=dft_control%sccs_control%rho_min) 457 IF (dft_control%sccs_control%rho_max < dft_control%sccs_control%rho_min) THEN 458 CALL cp_abort(__LOCATION__, & 459 "The SCCS parameter RHO_MAX is smaller than RHO_MIN. "// & 460 "Please, check your input!") 461 END IF 462 CASE (sccs_fattebert_gygi) 463 tmp_section => section_vals_get_subs_vals(sccs_section, "FATTEBERT-GYGI") 464 CALL section_vals_val_get(tmp_section, "BETA", & 465 r_val=dft_control%sccs_control%beta) 466 IF (dft_control%sccs_control%beta < 0.5_dp) THEN 467 CALL cp_abort(__LOCATION__, & 468 "A value smaller than 0.5 for the SCCS parameter beta "// & 469 "causes numerical problems. Please, check your input!") 470 END IF 471 CALL section_vals_val_get(tmp_section, "RHO_ZERO", & 472 r_val=dft_control%sccs_control%rho_zero) 473 CASE DEFAULT 474 CPABORT("Invalid SCCS model specified. Please, check your input!") 475 END SELECT 476 CALL cite_reference(Fattebert2002) 477 CALL cite_reference(Andreussi2012) 478 END IF 479 END IF 480 481 ! ZMP added input sections 482 ! Read the external density input section 483 tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_DENSITY") 484 CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_density) 485 486 ! Read the external vxc input section 487 tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_VXC") 488 CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_vxc) 489 490 END SUBROUTINE read_dft_control 491 492! ************************************************************************************************** 493!> \brief ... 494!> \param qs_control ... 495!> \param dft_section ... 496! ************************************************************************************************** 497 SUBROUTINE read_mgrid_section(qs_control, dft_section) 498 499 TYPE(qs_control_type), INTENT(INOUT) :: qs_control 500 TYPE(section_vals_type), POINTER :: dft_section 501 502 CHARACTER(len=*), PARAMETER :: routineN = 'read_mgrid_section', & 503 routineP = moduleN//':'//routineN 504 505 INTEGER :: handle, igrid_level, ngrid_level 506 LOGICAL :: explicit, multigrid_set 507 REAL(dp) :: cutoff 508 REAL(dp), DIMENSION(:), POINTER :: cutofflist 509 TYPE(section_vals_type), POINTER :: mgrid_section 510 511 CALL timeset(routineN, handle) 512 513 NULLIFY (mgrid_section, cutofflist) 514 mgrid_section => section_vals_get_subs_vals(dft_section, "MGRID") 515 516 CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=ngrid_level) 517 CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set) 518 CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=cutoff) 519 CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", r_val=qs_control%progression_factor) 520 CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=qs_control%commensurate_mgrids) 521 CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=qs_control%realspace_mgrids) 522 CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=qs_control%relative_cutoff) 523 CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", & 524 l_val=qs_control%skip_load_balance_distributed) 525 526 ! For SE and DFTB possibly override with new defaults 527 IF (qs_control%semi_empirical .OR. qs_control%dftb .OR. qs_control%xtb) THEN 528 ngrid_level = 1 529 multigrid_set = .FALSE. 530 ! Override default cutoff value unless user specified an explicit argument.. 531 CALL section_vals_val_get(mgrid_section, "CUTOFF", explicit=explicit, r_val=cutoff) 532 IF (.NOT. explicit) cutoff = 1.0_dp 533 END IF 534 535 ALLOCATE (qs_control%e_cutoff(ngrid_level)) 536 qs_control%cutoff = cutoff 537 538 IF (multigrid_set) THEN 539 ! Read the values from input 540 IF (qs_control%commensurate_mgrids) THEN 541 CPABORT("Do not specify cutoffs for the commensurate grids (NYI)") 542 END IF 543 544 CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=cutofflist) 545 IF (ASSOCIATED(cutofflist)) THEN 546 IF (SIZE(cutofflist, 1) /= ngrid_level) THEN 547 CPABORT("Number of multi-grids requested and number of cutoff values do not match") 548 END IF 549 DO igrid_level = 1, ngrid_level 550 qs_control%e_cutoff(igrid_level) = cutofflist(igrid_level) 551 END DO 552 END IF 553 ! set cutoff to smallest value in multgrid available with >= cutoff 554 DO igrid_level = ngrid_level, 1, -1 555 IF (qs_control%cutoff <= qs_control%e_cutoff(igrid_level)) THEN 556 qs_control%cutoff = qs_control%e_cutoff(igrid_level) 557 EXIT 558 END IF 559 ! set largest grid value to cutoff 560 IF (igrid_level == 1) THEN 561 qs_control%cutoff = qs_control%e_cutoff(1) 562 END IF 563 END DO 564 ELSE 565 IF (qs_control%commensurate_mgrids) qs_control%progression_factor = 4.0_dp 566 qs_control%e_cutoff(1) = qs_control%cutoff 567 DO igrid_level = 2, ngrid_level 568 qs_control%e_cutoff(igrid_level) = qs_control%e_cutoff(igrid_level - 1)/ & 569 qs_control%progression_factor 570 END DO 571 END IF 572 ! check that multigrids are ordered 573 DO igrid_level = 2, ngrid_level 574 IF (qs_control%e_cutoff(igrid_level) > qs_control%e_cutoff(igrid_level - 1)) THEN 575 CPABORT("The cutoff values for the multi-grids are not ordered from large to small") 576 ELSE IF (qs_control%e_cutoff(igrid_level) == qs_control%e_cutoff(igrid_level - 1)) THEN 577 CPABORT("The same cutoff value was specified for two multi-grids") 578 END IF 579 END DO 580 CALL timestop(handle) 581 END SUBROUTINE read_mgrid_section 582 583! ************************************************************************************************** 584!> \brief ... 585!> \param qs_control ... 586!> \param qs_section ... 587! ************************************************************************************************** 588 SUBROUTINE read_qs_section(qs_control, qs_section) 589 590 TYPE(qs_control_type), INTENT(INOUT) :: qs_control 591 TYPE(section_vals_type), POINTER :: qs_section 592 593 CHARACTER(len=*), PARAMETER :: routineN = 'read_qs_section', & 594 routineP = moduleN//':'//routineN 595 596 CHARACTER(LEN=default_string_length), & 597 DIMENSION(:), POINTER :: clist 598 INTEGER :: handle, itmp, j, jj, k, n_rep, n_var, & 599 ngauss, nrep 600 INTEGER, DIMENSION(:), POINTER :: tmplist 601 LOGICAL :: explicit, was_present 602 REAL(dp) :: tmp, tmpsqrt, value 603 REAL(dp), POINTER :: scal(:) 604 TYPE(enumeration_type), POINTER :: enum 605 TYPE(keyword_type), POINTER :: keyword 606 TYPE(section_type), POINTER :: section 607 TYPE(section_vals_type), POINTER :: cdft_control_section, ddapc_restraint_section, & 608 dftb_parameter, dftb_section, lri_optbas_section, mull_section, s2_restraint_section, & 609 se_section, xtb_parameter, xtb_section 610 611 CALL timeset(routineN, handle) 612 613 was_present = .FALSE. 614 NULLIFY (mull_section, ddapc_restraint_section, s2_restraint_section, & 615 se_section, dftb_section, xtb_section, dftb_parameter, xtb_parameter, lri_optbas_section, & 616 cdft_control_section) 617 618 mull_section => section_vals_get_subs_vals(qs_section, "MULLIKEN_RESTRAINT") 619 ddapc_restraint_section => section_vals_get_subs_vals(qs_section, "DDAPC_RESTRAINT") 620 s2_restraint_section => section_vals_get_subs_vals(qs_section, "S2_RESTRAINT") 621 se_section => section_vals_get_subs_vals(qs_section, "SE") 622 dftb_section => section_vals_get_subs_vals(qs_section, "DFTB") 623 xtb_section => section_vals_get_subs_vals(qs_section, "xTB") 624 dftb_parameter => section_vals_get_subs_vals(dftb_section, "PARAMETER") 625 xtb_parameter => section_vals_get_subs_vals(xtb_section, "PARAMETER") 626 lri_optbas_section => section_vals_get_subs_vals(qs_section, "OPTIMIZE_LRI_BASIS") 627 cdft_control_section => section_vals_get_subs_vals(qs_section, "CDFT") 628 629 ! Setup all defaults values and overwrite input parameters 630 ! EPS_DEFAULT should set the target accuracy in the total energy (~per electron) or a closely related value 631 CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=value) 632 tmpsqrt = SQRT(value) ! a trick to work around a NAG 5.1 optimizer bug 633 634 ! random choice ? 635 qs_control%eps_core_charge = value/100.0_dp 636 ! correct if all Gaussians would have the same radius (overlap will be smaller than eps_pgf_orb**2). 637 ! Can be significantly in error if not... requires fully new screening/pairlist procedures 638 qs_control%eps_pgf_orb = tmpsqrt 639 qs_control%eps_kg_orb = qs_control%eps_pgf_orb 640 ! consistent since also a kind of overlap 641 qs_control%eps_ppnl = qs_control%eps_pgf_orb/100.0_dp 642 ! accuracy is basically set by the overlap, this sets an empirical shift 643 qs_control%eps_ppl = 1.0E-2_dp 644 ! 645 qs_control%gapw_control%eps_cpc = value 646 ! expexted error in the density 647 qs_control%eps_rho_gspace = value 648 qs_control%eps_rho_rspace = value 649 ! error in the gradient, can be the sqrt of the error in the energy, ignored if map_consistent 650 qs_control%eps_gvg_rspace = tmpsqrt 651 ! 652 CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", n_rep_val=n_rep) 653 IF (n_rep /= 0) THEN 654 CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", r_val=qs_control%eps_core_charge) 655 END IF 656 CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", n_rep_val=n_rep) 657 IF (n_rep /= 0) THEN 658 CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", r_val=qs_control%eps_gvg_rspace) 659 END IF 660 CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep) 661 IF (n_rep /= 0) THEN 662 CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=qs_control%eps_pgf_orb) 663 END IF 664 CALL section_vals_val_get(qs_section, "EPS_KG_ORB", n_rep_val=n_rep) 665 IF (n_rep /= 0) THEN 666 CALL section_vals_val_get(qs_section, "EPS_KG_ORB", r_val=tmp) 667 qs_control%eps_kg_orb = SQRT(tmp) 668 END IF 669 CALL section_vals_val_get(qs_section, "EPS_PPL", n_rep_val=n_rep) 670 IF (n_rep /= 0) THEN 671 CALL section_vals_val_get(qs_section, "EPS_PPL", r_val=qs_control%eps_ppl) 672 END IF 673 CALL section_vals_val_get(qs_section, "EPS_PPNL", n_rep_val=n_rep) 674 IF (n_rep /= 0) THEN 675 CALL section_vals_val_get(qs_section, "EPS_PPNL", r_val=qs_control%eps_ppnl) 676 END IF 677 CALL section_vals_val_get(qs_section, "EPS_RHO", n_rep_val=n_rep) 678 IF (n_rep /= 0) THEN 679 CALL section_vals_val_get(qs_section, "EPS_RHO", r_val=qs_control%eps_rho_gspace) 680 qs_control%eps_rho_rspace = qs_control%eps_rho_gspace 681 END IF 682 CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", n_rep_val=n_rep) 683 IF (n_rep /= 0) THEN 684 CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", r_val=qs_control%eps_rho_rspace) 685 END IF 686 CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", n_rep_val=n_rep) 687 IF (n_rep /= 0) THEN 688 CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", r_val=qs_control%eps_rho_gspace) 689 END IF 690 CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", n_rep_val=n_rep) 691 IF (n_rep /= 0) THEN 692 CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", r_val=qs_control%eps_filter_matrix) 693 END IF 694 CALL section_vals_val_get(qs_section, "EPS_CPC", n_rep_val=n_rep) 695 IF (n_rep /= 0) THEN 696 CALL section_vals_val_get(qs_section, "EPS_CPC", r_val=qs_control%gapw_control%eps_cpc) 697 END IF 698 699 CALL section_vals_val_get(qs_section, "EPSFIT", r_val=qs_control%gapw_control%eps_fit) 700 CALL section_vals_val_get(qs_section, "EPSISO", r_val=qs_control%gapw_control%eps_iso) 701 CALL section_vals_val_get(qs_section, "EPSSVD", r_val=qs_control%gapw_control%eps_svd) 702 CALL section_vals_val_get(qs_section, "EPSRHO0", r_val=qs_control%gapw_control%eps_Vrho0) 703 CALL section_vals_val_get(qs_section, "ALPHA0_HARD", r_val=qs_control%gapw_control%alpha0_hard) 704 qs_control%gapw_control%lrho1_eq_lrho0 = .FALSE. 705 qs_control%gapw_control%alpha0_hard_from_input = .FALSE. 706 IF (qs_control%gapw_control%alpha0_hard /= 0.0_dp) qs_control%gapw_control%alpha0_hard_from_input = .TRUE. 707 CALL section_vals_val_get(qs_section, "FORCE_PAW", l_val=qs_control%gapw_control%force_paw) 708 CALL section_vals_val_get(qs_section, "MAX_RAD_LOCAL", r_val=qs_control%gapw_control%max_rad_local) 709 710 CALL section_vals_val_get(qs_section, "LS_SCF", l_val=qs_control%do_ls_scf) 711 CALL section_vals_val_get(qs_section, "ALMO_SCF", l_val=qs_control%do_almo_scf) 712 CALL section_vals_val_get(qs_section, "KG_METHOD", l_val=qs_control%do_kg) 713 714 ! Logicals 715 CALL section_vals_val_get(qs_section, "MAP_CONSISTENT", l_val=qs_control%map_consistent) 716 CALL section_vals_val_get(qs_section, "REF_EMBED_SUBSYS", l_val=qs_control%ref_embed_subsys) 717 CALL section_vals_val_get(qs_section, "CLUSTER_EMBED_SUBSYS", l_val=qs_control%cluster_embed_subsys) 718 CALL section_vals_val_get(qs_section, "HIGH_LEVEL_EMBED_SUBSYS", l_val=qs_control%high_level_embed_subsys) 719 CALL section_vals_val_get(qs_section, "DFET_EMBEDDED", l_val=qs_control%dfet_embedded) 720 CALL section_vals_val_get(qs_section, "DMFET_EMBEDDED", l_val=qs_control%dmfet_embedded) 721 722 ! Integers gapw 723 CALL section_vals_val_get(qs_section, "LMAXN1", i_val=qs_control%gapw_control%lmax_sphere) 724 CALL section_vals_val_get(qs_section, "LMAXN0", i_val=qs_control%gapw_control%lmax_rho0) 725 CALL section_vals_val_get(qs_section, "LADDN0", i_val=qs_control%gapw_control%ladd_rho0) 726 CALL section_vals_val_get(qs_section, "QUADRATURE", i_val=qs_control%gapw_control%quadrature) 727 728 ! Integers grids 729 CALL section_vals_val_get(qs_section, "PW_GRID", i_val=itmp) 730 SELECT CASE (itmp) 731 CASE (do_pwgrid_spherical) 732 qs_control%pw_grid_opt%spherical = .TRUE. 733 qs_control%pw_grid_opt%fullspace = .FALSE. 734 CASE (do_pwgrid_ns_fullspace) 735 qs_control%pw_grid_opt%spherical = .FALSE. 736 qs_control%pw_grid_opt%fullspace = .TRUE. 737 CASE (do_pwgrid_ns_halfspace) 738 qs_control%pw_grid_opt%spherical = .FALSE. 739 qs_control%pw_grid_opt%fullspace = .FALSE. 740 END SELECT 741 742 ! Method for PPL calculation 743 CALL section_vals_val_get(qs_section, "CORE_PPL", i_val=itmp) 744 qs_control%do_ppl_method = itmp 745 746 CALL section_vals_val_get(qs_section, "PW_GRID_LAYOUT", i_vals=tmplist) 747 qs_control%pw_grid_opt%distribution_layout = tmplist 748 CALL section_vals_val_get(qs_section, "PW_GRID_BLOCKED", i_val=qs_control%pw_grid_opt%blocked) 749 750 !Integers extrapolation 751 CALL section_vals_val_get(qs_section, "EXTRAPOLATION", i_val=qs_control%wf_interpolation_method_nr) 752 CALL section_vals_val_get(qs_section, "EXTRAPOLATION_ORDER", i_val=qs_control%wf_extrapolation_order) 753 754 !Method 755 CALL section_vals_val_get(qs_section, "METHOD", i_val=qs_control%method_id) 756 NULLIFY (section, keyword, enum) 757 CALL create_qs_section(section) 758 keyword => section_get_keyword(section, "METHOD") 759 CALL keyword_get(keyword, enum=enum) 760 qs_control%method = enum_i2c(enum, qs_control%method_id) 761 CALL section_release(section) 762 qs_control%gapw = .FALSE. 763 qs_control%gapw_xc = .FALSE. 764 qs_control%gpw = .FALSE. 765 qs_control%pao = .FALSE. 766 qs_control%dftb = .FALSE. 767 qs_control%xtb = .FALSE. 768 qs_control%semi_empirical = .FALSE. 769 qs_control%ofgpw = .FALSE. 770 qs_control%lrigpw = .FALSE. 771 qs_control%rigpw = .FALSE. 772 SELECT CASE (qs_control%method_id) 773 CASE (do_method_gapw) 774 CALL cite_reference(Lippert1999) 775 CALL cite_reference(Krack2000) 776 qs_control%gapw = .TRUE. 777 CASE (do_method_gapw_xc) 778 qs_control%gapw_xc = .TRUE. 779 CASE (do_method_gpw) 780 CALL cite_reference(Lippert1997) 781 CALL cite_reference(VandeVondele2005a) 782 qs_control%gpw = .TRUE. 783 CASE (do_method_ofgpw) 784 qs_control%ofgpw = .TRUE. 785 CASE (do_method_lrigpw) 786 qs_control%lrigpw = .TRUE. 787 CASE (do_method_rigpw) 788 qs_control%rigpw = .TRUE. 789 CASE (do_method_dftb) 790 qs_control%dftb = .TRUE. 791 CALL cite_reference(Porezag1995) 792 CALL cite_reference(Seifert1996) 793 CASE (do_method_xtb) 794 qs_control%xtb = .TRUE. 795 CALL cite_reference(Grimme2017) 796 CASE (do_method_mndo) 797 CALL cite_reference(Dewar1977) 798 qs_control%semi_empirical = .TRUE. 799 CASE (do_method_am1) 800 CALL cite_reference(Dewar1985) 801 qs_control%semi_empirical = .TRUE. 802 CASE (do_method_pm3) 803 CALL cite_reference(Stewart1989) 804 qs_control%semi_empirical = .TRUE. 805 CASE (do_method_pnnl) 806 CALL cite_reference(Schenter2008) 807 qs_control%semi_empirical = .TRUE. 808 CASE (do_method_pm6) 809 CALL cite_reference(Stewart2007) 810 qs_control%semi_empirical = .TRUE. 811 CASE (do_method_pm6fm) 812 CALL cite_reference(VanVoorhis2015) 813 qs_control%semi_empirical = .TRUE. 814 CASE (do_method_pdg) 815 CALL cite_reference(Repasky2002) 816 qs_control%semi_empirical = .TRUE. 817 CASE (do_method_rm1) 818 CALL cite_reference(Rocha2006) 819 qs_control%semi_empirical = .TRUE. 820 CASE (do_method_mndod) 821 CALL cite_reference(Dewar1977) 822 CALL cite_reference(Thiel1992) 823 qs_control%semi_empirical = .TRUE. 824 END SELECT 825 826 CALL section_vals_get(mull_section, explicit=qs_control%mulliken_restraint) 827 828 IF (qs_control%mulliken_restraint) THEN 829 CALL section_vals_val_get(mull_section, "STRENGTH", r_val=qs_control%mulliken_restraint_control%strength) 830 CALL section_vals_val_get(mull_section, "TARGET", r_val=qs_control%mulliken_restraint_control%target) 831 CALL section_vals_val_get(mull_section, "ATOMS", n_rep_val=n_rep) 832 jj = 0 833 DO k = 1, n_rep 834 CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist) 835 jj = jj + SIZE(tmplist) 836 END DO 837 qs_control%mulliken_restraint_control%natoms = jj 838 IF (qs_control%mulliken_restraint_control%natoms < 1) & 839 CPABORT("Need at least 1 atom to use mulliken contraints") 840 ALLOCATE (qs_control%mulliken_restraint_control%atoms(qs_control%mulliken_restraint_control%natoms)) 841 jj = 0 842 DO k = 1, n_rep 843 CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist) 844 DO j = 1, SIZE(tmplist) 845 jj = jj + 1 846 qs_control%mulliken_restraint_control%atoms(jj) = tmplist(j) 847 END DO 848 END DO 849 ENDIF 850 CALL section_vals_get(ddapc_restraint_section, n_repetition=nrep, explicit=qs_control%ddapc_restraint) 851 IF (qs_control%ddapc_restraint) THEN 852 ALLOCATE (qs_control%ddapc_restraint_control(nrep)) 853 CALL read_ddapc_section(qs_control, qs_section=qs_section) 854 qs_control%ddapc_restraint_is_spin = .FALSE. 855 qs_control%ddapc_explicit_potential = .FALSE. 856 ENDIF 857 858 CALL section_vals_get(s2_restraint_section, explicit=qs_control%s2_restraint) 859 IF (qs_control%s2_restraint) THEN 860 CALL section_vals_val_get(s2_restraint_section, "STRENGTH", & 861 r_val=qs_control%s2_restraint_control%strength) 862 CALL section_vals_val_get(s2_restraint_section, "TARGET", & 863 r_val=qs_control%s2_restraint_control%target) 864 CALL section_vals_val_get(s2_restraint_section, "FUNCTIONAL_FORM", & 865 i_val=qs_control%s2_restraint_control%functional_form) 866 ENDIF 867 868 CALL section_vals_get(cdft_control_section, explicit=qs_control%cdft) 869 IF (qs_control%cdft) THEN 870 CALL read_cdft_control_section(qs_control, cdft_control_section) 871 ENDIF 872 873 ! Semi-empirical code 874 IF (qs_control%semi_empirical) THEN 875 CALL section_vals_val_get(se_section, "ORTHOGONAL_BASIS", & 876 l_val=qs_control%se_control%orthogonal_basis) 877 CALL section_vals_val_get(se_section, "DELTA", & 878 r_val=qs_control%se_control%delta) 879 CALL section_vals_val_get(se_section, "ANALYTICAL_GRADIENTS", & 880 l_val=qs_control%se_control%analytical_gradients) 881 CALL section_vals_val_get(se_section, "FORCE_KDSO-D_EXCHANGE", & 882 l_val=qs_control%se_control%force_kdsod_EX) 883 ! Integral Screening 884 CALL section_vals_val_get(se_section, "INTEGRAL_SCREENING", & 885 i_val=qs_control%se_control%integral_screening) 886 IF (qs_control%method_id == do_method_pnnl) THEN 887 IF (qs_control%se_control%integral_screening /= do_se_IS_slater) & 888 CALL cp_warn(__LOCATION__, & 889 "PNNL semi-empirical parameterization supports only the Slater type "// & 890 "integral scheme. Revert to Slater and continue the calculation.") 891 qs_control%se_control%integral_screening = do_se_IS_slater 892 END IF 893 ! Global Arrays variable 894 CALL section_vals_val_get(se_section, "GA%NCELLS", & 895 i_val=qs_control%se_control%ga_ncells) 896 ! Long-Range correction 897 CALL section_vals_val_get(se_section, "LR_CORRECTION%CUTOFF", & 898 r_val=qs_control%se_control%cutoff_lrc) 899 qs_control%se_control%taper_lrc = qs_control%se_control%cutoff_lrc 900 CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", & 901 explicit=explicit) 902 IF (explicit) THEN 903 CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", & 904 r_val=qs_control%se_control%taper_lrc) 905 END IF 906 CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_RANGE", & 907 r_val=qs_control%se_control%range_lrc) 908 ! Coulomb 909 CALL section_vals_val_get(se_section, "COULOMB%CUTOFF", & 910 r_val=qs_control%se_control%cutoff_cou) 911 qs_control%se_control%taper_cou = qs_control%se_control%cutoff_cou 912 CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", & 913 explicit=explicit) 914 IF (explicit) THEN 915 CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", & 916 r_val=qs_control%se_control%taper_cou) 917 END IF 918 CALL section_vals_val_get(se_section, "COULOMB%RC_RANGE", & 919 r_val=qs_control%se_control%range_cou) 920 ! Exchange 921 CALL section_vals_val_get(se_section, "EXCHANGE%CUTOFF", & 922 r_val=qs_control%se_control%cutoff_exc) 923 qs_control%se_control%taper_exc = qs_control%se_control%cutoff_exc 924 CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", & 925 explicit=explicit) 926 IF (explicit) THEN 927 CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", & 928 r_val=qs_control%se_control%taper_exc) 929 END IF 930 CALL section_vals_val_get(se_section, "EXCHANGE%RC_RANGE", & 931 r_val=qs_control%se_control%range_exc) 932 ! Screening (only if the integral scheme is of dumped type) 933 IF (qs_control%se_control%integral_screening == do_se_IS_kdso_d) THEN 934 CALL section_vals_val_get(se_section, "SCREENING%RC_TAPER", & 935 r_val=qs_control%se_control%taper_scr) 936 CALL section_vals_val_get(se_section, "SCREENING%RC_RANGE", & 937 r_val=qs_control%se_control%range_scr) 938 END IF 939 ! Periodic Type Calculation 940 CALL section_vals_val_get(se_section, "PERIODIC", & 941 i_val=qs_control%se_control%periodic_type) 942 SELECT CASE (qs_control%se_control%periodic_type) 943 CASE (do_se_lr_none) 944 qs_control%se_control%do_ewald = .FALSE. 945 qs_control%se_control%do_ewald_r3 = .FALSE. 946 qs_control%se_control%do_ewald_gks = .FALSE. 947 CASE (do_se_lr_ewald) 948 qs_control%se_control%do_ewald = .TRUE. 949 qs_control%se_control%do_ewald_r3 = .FALSE. 950 qs_control%se_control%do_ewald_gks = .FALSE. 951 CASE (do_se_lr_ewald_gks) 952 qs_control%se_control%do_ewald = .FALSE. 953 qs_control%se_control%do_ewald_r3 = .FALSE. 954 qs_control%se_control%do_ewald_gks = .TRUE. 955 IF (qs_control%method_id /= do_method_pnnl) & 956 CALL cp_abort(__LOCATION__, & 957 "A periodic semi-empirical calculation was requested with a long-range "// & 958 "summation on the single integral evaluation. This scheme is supported "// & 959 "only by the PNNL parameterization.") 960 CASE (do_se_lr_ewald_r3) 961 qs_control%se_control%do_ewald = .TRUE. 962 qs_control%se_control%do_ewald_r3 = .TRUE. 963 qs_control%se_control%do_ewald_gks = .FALSE. 964 IF (qs_control%se_control%integral_screening /= do_se_IS_kdso) & 965 CALL cp_abort(__LOCATION__, & 966 "A periodic semi-empirical calculation was requested with a long-range "// & 967 "summation for the slowly convergent part 1/R^3, which is not congruent "// & 968 "with the integral screening chosen. The only integral screening supported "// & 969 "by this periodic type calculation is the standard Klopman-Dewar-Sabelli-Ohno.") 970 END SELECT 971 972 ! dispersion pair potentials 973 CALL section_vals_val_get(se_section, "DISPERSION", & 974 l_val=qs_control%se_control%dispersion) 975 CALL section_vals_val_get(se_section, "DISPERSION_RADIUS", & 976 r_val=qs_control%se_control%rcdisp) 977 CALL section_vals_val_get(se_section, "COORDINATION_CUTOFF", & 978 r_val=qs_control%se_control%epscn) 979 CALL section_vals_val_get(se_section, "D3_SCALING", r_vals=scal) 980 qs_control%se_control%sd3(1) = scal(1) 981 qs_control%se_control%sd3(2) = scal(2) 982 qs_control%se_control%sd3(3) = scal(3) 983 CALL section_vals_val_get(se_section, "DISPERSION_PARAMETER_FILE", & 984 c_val=qs_control%se_control%dispersion_parameter_file) 985 986 ! Stop the execution for non-implemented features 987 IF (qs_control%se_control%periodic_type == do_se_lr_ewald_r3) THEN 988 CPABORT("EWALD_R3 not implemented yet!") 989 END IF 990 991 IF (qs_control%method_id == do_method_mndo .OR. & 992 qs_control%method_id == do_method_am1 .OR. & 993 qs_control%method_id == do_method_mndod .OR. & 994 qs_control%method_id == do_method_pdg .OR. & 995 qs_control%method_id == do_method_pm3 .OR. & 996 qs_control%method_id == do_method_pm6 .OR. & 997 qs_control%method_id == do_method_pm6fm .OR. & 998 qs_control%method_id == do_method_pnnl .OR. & 999 qs_control%method_id == do_method_rm1) THEN 1000 qs_control%se_control%orthogonal_basis = .TRUE. 1001 END IF 1002 END IF 1003 1004 ! DFTB code 1005 IF (qs_control%dftb) THEN 1006 CALL section_vals_val_get(dftb_section, "ORTHOGONAL_BASIS", & 1007 l_val=qs_control%dftb_control%orthogonal_basis) 1008 CALL section_vals_val_get(dftb_section, "SELF_CONSISTENT", & 1009 l_val=qs_control%dftb_control%self_consistent) 1010 CALL section_vals_val_get(dftb_section, "DISPERSION", & 1011 l_val=qs_control%dftb_control%dispersion) 1012 CALL section_vals_val_get(dftb_section, "DIAGONAL_DFTB3", & 1013 l_val=qs_control%dftb_control%dftb3_diagonal) 1014 CALL section_vals_val_get(dftb_section, "HB_SR_GAMMA", & 1015 l_val=qs_control%dftb_control%hb_sr_damp) 1016 CALL section_vals_val_get(dftb_section, "EPS_DISP", & 1017 r_val=qs_control%dftb_control%eps_disp) 1018 CALL section_vals_val_get(dftb_section, "DO_EWALD", & 1019 l_val=qs_control%dftb_control%do_ewald) 1020 CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_PATH", & 1021 c_val=qs_control%dftb_control%sk_file_path) 1022 CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_NAME", & 1023 c_val=qs_control%dftb_control%sk_file_list) 1024 CALL section_vals_val_get(dftb_parameter, "HB_SR_PARAM", & 1025 r_val=qs_control%dftb_control%hb_sr_para) 1026 CALL section_vals_val_get(dftb_parameter, "SK_FILE", n_rep_val=n_var) 1027 ALLOCATE (qs_control%dftb_control%sk_pair_list(3, n_var)) 1028 DO k = 1, n_var 1029 CALL section_vals_val_get(dftb_parameter, "SK_FILE", i_rep_val=k, & 1030 c_vals=clist) 1031 qs_control%dftb_control%sk_pair_list(1:3, k) = clist(1:3) 1032 END DO 1033 ! Dispersion type 1034 CALL section_vals_val_get(dftb_parameter, "DISPERSION_TYPE", & 1035 i_val=qs_control%dftb_control%dispersion_type) 1036 CALL section_vals_val_get(dftb_parameter, "UFF_FORCE_FIELD", & 1037 c_val=qs_control%dftb_control%uff_force_field) 1038 ! D3 Dispersion 1039 CALL section_vals_val_get(dftb_parameter, "DISPERSION_RADIUS", & 1040 r_val=qs_control%dftb_control%rcdisp) 1041 CALL section_vals_val_get(dftb_parameter, "COORDINATION_CUTOFF", & 1042 r_val=qs_control%dftb_control%epscn) 1043 CALL section_vals_val_get(dftb_parameter, "D3_SCALING", r_vals=scal) 1044 qs_control%dftb_control%sd3(1) = scal(1) 1045 qs_control%dftb_control%sd3(2) = scal(2) 1046 qs_control%dftb_control%sd3(3) = scal(3) 1047 CALL section_vals_val_get(dftb_parameter, "DISPERSION_PARAMETER_FILE", & 1048 c_val=qs_control%dftb_control%dispersion_parameter_file) 1049 1050 IF (qs_control%dftb_control%dispersion) CALL cite_reference(Zhechkov2005) 1051 IF (qs_control%dftb_control%self_consistent) CALL cite_reference(Elstner1998) 1052 IF (qs_control%dftb_control%hb_sr_damp) CALL cite_reference(Hu2007) 1053 END IF 1054 1055 ! xTB code 1056 IF (qs_control%xtb) THEN 1057 CALL section_vals_val_get(xtb_section, "DO_EWALD", & 1058 l_val=qs_control%xtb_control%do_ewald) 1059 CALL section_vals_val_get(xtb_section, "STO_NG", i_val=ngauss) 1060 qs_control%xtb_control%sto_ng = ngauss 1061 CALL section_vals_val_get(xtb_section, "HYDROGEN_STO_NG", i_val=ngauss) 1062 qs_control%xtb_control%h_sto_ng = ngauss 1063 CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_PATH", & 1064 c_val=qs_control%xtb_control%parameter_file_path) 1065 CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_NAME", & 1066 c_val=qs_control%xtb_control%parameter_file_name) 1067 ! D3 Dispersion 1068 CALL section_vals_val_get(xtb_parameter, "DISPERSION_RADIUS", & 1069 r_val=qs_control%xtb_control%rcdisp) 1070 CALL section_vals_val_get(xtb_parameter, "COORDINATION_CUTOFF", & 1071 r_val=qs_control%xtb_control%epscn) 1072 CALL section_vals_val_get(xtb_parameter, "D3BJ_SCALING", r_vals=scal) 1073 qs_control%xtb_control%s6 = scal(1) 1074 qs_control%xtb_control%s8 = scal(2) 1075 CALL section_vals_val_get(xtb_parameter, "D3BJ_PARAM", r_vals=scal) 1076 qs_control%xtb_control%a1 = scal(1) 1077 qs_control%xtb_control%a2 = scal(2) 1078 CALL section_vals_val_get(xtb_parameter, "DISPERSION_PARAMETER_FILE", & 1079 c_val=qs_control%xtb_control%dispersion_parameter_file) 1080 ! global parameters 1081 CALL section_vals_val_get(xtb_parameter, "HUCKEL_CONSTANTS", r_vals=scal) 1082 qs_control%xtb_control%ks = scal(1) 1083 qs_control%xtb_control%kp = scal(2) 1084 qs_control%xtb_control%kd = scal(3) 1085 qs_control%xtb_control%ksp = scal(4) 1086 qs_control%xtb_control%k2sh = scal(5) 1087 CALL section_vals_val_get(xtb_parameter, "COULOMB_CONSTANTS", r_vals=scal) 1088 qs_control%xtb_control%kg = scal(1) 1089 qs_control%xtb_control%kf = scal(2) 1090 CALL section_vals_val_get(xtb_parameter, "CN_CONSTANTS", r_vals=scal) 1091 qs_control%xtb_control%kcns = scal(1) 1092 qs_control%xtb_control%kcnp = scal(2) 1093 qs_control%xtb_control%kcnd = scal(3) 1094 CALL section_vals_val_get(xtb_parameter, "EN_CONSTANT", r_vals=scal) 1095 qs_control%xtb_control%ken = scal(1) 1096 ! XB 1097 CALL section_vals_val_get(xtb_section, "USE_HALOGEN_CORRECTION", & 1098 l_val=qs_control%xtb_control%xb_interaction) 1099 CALL section_vals_val_get(xtb_parameter, "HALOGEN_BINDING", r_vals=scal) 1100 qs_control%xtb_control%kxr = scal(1) 1101 qs_control%xtb_control%kx2 = scal(2) 1102 ! XB_radius 1103 CALL section_vals_val_get(xtb_parameter, "XB_RADIUS", r_val=qs_control%xtb_control%xb_radius) 1104 ! Kab 1105 CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", n_rep_val=n_rep) 1106 ! For debug purposes 1107 CALL section_vals_val_get(xtb_section, "COULOMB_INTERACTION", & 1108 l_val=qs_control%xtb_control%coulomb_interaction) 1109 CALL section_vals_val_get(xtb_section, "TB3_INTERACTION", & 1110 l_val=qs_control%xtb_control%tb3_interaction) 1111 ! Check for bad atomic charges 1112 CALL section_vals_val_get(xtb_section, "CHECK_ATOMIC_CHARGES", & 1113 l_val=qs_control%xtb_control%check_atomic_charges) 1114 IF (n_rep > 0) THEN 1115 ALLOCATE (qs_control%xtb_control%kab_param(3, n_rep)) 1116 DO j = 1, n_rep 1117 CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", i_rep_val=j, c_vals=clist) 1118 qs_control%xtb_control%kab_param(1, j) = clist(1) 1119 qs_control%xtb_control%kab_param(2, j) = clist(2) 1120 qs_control%xtb_control%kab_param(3, j) = clist(3) 1121 END DO 1122 END IF 1123 END IF 1124 1125 ! Optimize LRI basis set 1126 CALL section_vals_get(lri_optbas_section, explicit=qs_control%lri_optbas) 1127 1128 CALL timestop(handle) 1129 END SUBROUTINE read_qs_section 1130 1131! ************************************************************************************************** 1132!> \brief ... 1133!> \param t_control ... 1134!> \param dft_section ... 1135! ************************************************************************************************** 1136 SUBROUTINE read_tddfpt_control(t_control, dft_section) 1137 TYPE(tddfpt_control_type) :: t_control 1138 TYPE(section_vals_type), POINTER :: dft_section 1139 1140 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_tddfpt_control', & 1141 routineP = moduleN//':'//routineN 1142 1143 LOGICAL :: kenergy_den 1144 TYPE(section_vals_type), POINTER :: sic_section, t_section 1145 1146 kenergy_den = .FALSE. 1147 NULLIFY (sic_section, t_section) 1148 t_section => section_vals_get_subs_vals(dft_section, "TDDFPT") 1149 1150 CALL section_vals_val_get(t_section, "CONVERGENCE", r_val=t_control%tolerance) 1151 CALL section_vals_val_get(t_section, "NEV", i_val=t_control%n_ev) 1152 CALL section_vals_val_get(t_section, "MAX_KV", i_val=t_control%max_kv) 1153 CALL section_vals_val_get(t_section, "RESTARTS", i_val=t_control%n_restarts) 1154 CALL section_vals_val_get(t_section, "NREORTHO", i_val=t_control%n_reortho) 1155 CALL section_vals_val_get(t_section, "RES_ETYPE", i_val=t_control%res_etype) 1156 CALL section_vals_val_get(t_section, "DIAG_METHOD", i_val=t_control%diag_method) 1157 CALL section_vals_val_get(t_section, "KERNEL", l_val=t_control%do_kernel) 1158 CALL section_vals_val_get(t_section, "LSD_SINGLETS", l_val=t_control%lsd_singlets) 1159 CALL section_vals_val_get(t_section, "INVERT_S", l_val=t_control%invert_S) 1160 CALL section_vals_val_get(t_section, "PRECOND", l_val=t_control%precond) 1161 CALL section_vals_val_get(t_section, "OE_CORR", i_val=t_control%oe_corr) 1162 1163 t_control%use_kinetic_energy_density = .FALSE. 1164 sic_section => section_vals_get_subs_vals(t_section, "SIC") 1165 CALL section_vals_val_get(sic_section, "SIC_METHOD", i_val=t_control%sic_method_id) 1166 CALL section_vals_val_get(sic_section, "ORBITAL_SET", i_val=t_control%sic_list_id) 1167 CALL section_vals_val_get(sic_section, "SIC_SCALING_A", r_val=t_control%sic_scaling_a) 1168 CALL section_vals_val_get(sic_section, "SIC_SCALING_B", r_val=t_control%sic_scaling_b) 1169 1170 END SUBROUTINE read_tddfpt_control 1171 1172! ************************************************************************************************** 1173!> \brief Read TDDFPT-related input parameters. 1174!> \param t_control TDDFPT control parameters 1175!> \param t_section TDDFPT input section 1176!> \param qs_control Quickstep control parameters 1177! ************************************************************************************************** 1178 SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control) 1179 TYPE(tddfpt2_control_type), POINTER :: t_control 1180 TYPE(section_vals_type), POINTER :: t_section 1181 TYPE(qs_control_type), POINTER :: qs_control 1182 1183 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_tddfpt2_control', & 1184 routineP = moduleN//':'//routineN 1185 1186 INTEGER :: handle 1187 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds 1188 LOGICAL :: explicit, multigrid_set 1189 TYPE(section_vals_type), POINTER :: dipole_section, mgrid_section, xc_func, & 1190 xc_section 1191 1192 CALL timeset(routineN, handle) 1193 1194 CALL section_vals_val_get(t_section, "_SECTION_PARAMETERS_", l_val=t_control%enabled) 1195 1196 CALL section_vals_val_get(t_section, "NSTATES", i_val=t_control%nstates) 1197 CALL section_vals_val_get(t_section, "MAX_ITER", i_val=t_control%niters) 1198 CALL section_vals_val_get(t_section, "MAX_KV", i_val=t_control%nkvs) 1199 CALL section_vals_val_get(t_section, "NLUMO", i_val=t_control%nlumo) 1200 CALL section_vals_val_get(t_section, "NPROC_STATE", i_val=t_control%nprocs) 1201 1202 CALL section_vals_val_get(t_section, "CONVERGENCE", r_val=t_control%conv) 1203 CALL section_vals_val_get(t_section, "MIN_AMPLITUDE", r_val=t_control%min_excitation_amplitude) 1204 CALL section_vals_val_get(t_section, "ORTHOGONAL_EPS", r_val=t_control%orthogonal_eps) 1205 1206 CALL section_vals_val_get(t_section, "RESTART", l_val=t_control%is_restart) 1207 CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets) 1208 1209 IF (t_control%conv < 0) & 1210 t_control%conv = ABS(t_control%conv) 1211 1212 ! DIPOLE_MOMENTS subsection 1213 dipole_section => section_vals_get_subs_vals(t_section, "DIPOLE_MOMENTS") 1214 CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", i_val=t_control%dipole_form) 1215 CALL section_vals_val_get(dipole_section, "REFERENCE", i_val=t_control%dipole_reference) 1216 CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", explicit=explicit) 1217 IF (explicit) THEN 1218 CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", r_vals=t_control%dipole_ref_point) 1219 ELSE 1220 NULLIFY (t_control%dipole_ref_point) 1221 IF (t_control%dipole_form == xas_dip_len .AND. t_control%dipole_reference == use_mom_ref_user) THEN 1222 CPABORT("User-defined reference point should be given explicitly") 1223 END IF 1224 END IF 1225 1226 ! MGRID subsection 1227 mgrid_section => section_vals_get_subs_vals(t_section, "MGRID") 1228 CALL section_vals_get(mgrid_section, explicit=t_control%mgrid_is_explicit) 1229 1230 IF (t_control%mgrid_is_explicit) THEN 1231 CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=t_control%mgrid_ngrids, explicit=explicit) 1232 IF (.NOT. explicit) t_control%mgrid_ngrids = SIZE(qs_control%e_cutoff) 1233 1234 CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=t_control%mgrid_cutoff, explicit=explicit) 1235 IF (.NOT. explicit) t_control%mgrid_cutoff = qs_control%cutoff 1236 1237 CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", & 1238 r_val=t_control%mgrid_progression_factor, explicit=explicit) 1239 IF (explicit) THEN 1240 IF (t_control%mgrid_progression_factor <= 1.0_dp) & 1241 CALL cp_abort(__LOCATION__, & 1242 "Progression factor should be greater then 1.0 to ensure multi-grid ordering") 1243 ELSE 1244 t_control%mgrid_progression_factor = qs_control%progression_factor 1245 END IF 1246 1247 CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=t_control%mgrid_commensurate_mgrids, explicit=explicit) 1248 IF (.NOT. explicit) t_control%mgrid_commensurate_mgrids = qs_control%commensurate_mgrids 1249 IF (t_control%mgrid_commensurate_mgrids) THEN 1250 IF (explicit) THEN 1251 t_control%mgrid_progression_factor = 4.0_dp 1252 ELSE 1253 t_control%mgrid_progression_factor = qs_control%progression_factor 1254 END IF 1255 END IF 1256 1257 CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=t_control%mgrid_relative_cutoff, explicit=explicit) 1258 IF (.NOT. explicit) t_control%mgrid_relative_cutoff = qs_control%relative_cutoff 1259 1260 CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set, explicit=explicit) 1261 IF (.NOT. explicit) multigrid_set = .FALSE. 1262 IF (multigrid_set) THEN 1263 CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=t_control%mgrid_e_cutoff) 1264 ELSE 1265 NULLIFY (t_control%mgrid_e_cutoff) 1266 END IF 1267 1268 CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=t_control%mgrid_realspace_mgrids, explicit=explicit) 1269 IF (.NOT. explicit) t_control%mgrid_realspace_mgrids = qs_control%realspace_mgrids 1270 1271 CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", & 1272 l_val=t_control%mgrid_skip_load_balance, explicit=explicit) 1273 IF (.NOT. explicit) t_control%mgrid_skip_load_balance = qs_control%skip_load_balance_distributed 1274 1275 IF (ASSOCIATED(t_control%mgrid_e_cutoff)) THEN 1276 IF (SIZE(t_control%mgrid_e_cutoff) /= t_control%mgrid_ngrids) & 1277 CPABORT("Inconsistent values for number of multi-grids") 1278 1279 ! sort multi-grids in descending order according to their cutoff values 1280 t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff 1281 ALLOCATE (inds(t_control%mgrid_ngrids)) 1282 CALL sort(t_control%mgrid_e_cutoff, t_control%mgrid_ngrids, inds) 1283 DEALLOCATE (inds) 1284 t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff 1285 END IF 1286 END IF 1287 1288 ! expand XC subsection (if given explicitly) 1289 xc_section => section_vals_get_subs_vals(t_section, "XC") 1290 xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 1291 CALL section_vals_get(xc_func, explicit=explicit) 1292 IF (explicit) & 1293 CALL xc_functionals_expand(xc_func, xc_section) 1294 1295 CALL timestop(handle) 1296 END SUBROUTINE read_tddfpt2_control 1297 1298! ************************************************************************************************** 1299!> \brief Write the DFT control parameters to the output unit. 1300!> \param dft_control ... 1301!> \param dft_section ... 1302! ************************************************************************************************** 1303 SUBROUTINE write_dft_control(dft_control, dft_section) 1304 TYPE(dft_control_type), POINTER :: dft_control 1305 TYPE(section_vals_type), POINTER :: dft_section 1306 1307 CHARACTER(len=*), PARAMETER :: routineN = 'write_dft_control', & 1308 routineP = moduleN//':'//routineN 1309 1310 CHARACTER(LEN=20) :: tmpStr 1311 INTEGER :: handle, output_unit 1312 REAL(kind=dp) :: density_cut, density_smooth_cut_range, & 1313 gradient_cut, tau_cut 1314 TYPE(cp_logger_type), POINTER :: logger 1315 TYPE(enumeration_type), POINTER :: enum 1316 TYPE(keyword_type), POINTER :: keyword 1317 TYPE(section_type), POINTER :: section 1318 TYPE(section_vals_type), POINTER :: xc_section 1319 1320 IF (dft_control%qs_control%semi_empirical) RETURN 1321 IF (dft_control%qs_control%dftb) RETURN 1322 IF (dft_control%qs_control%xtb) THEN 1323 CALL write_xtb_control(dft_control%qs_control%xtb_control, dft_section) 1324 RETURN 1325 END IF 1326 CALL timeset(routineN, handle) 1327 1328 NULLIFY (logger) 1329 logger => cp_get_default_logger() 1330 1331 output_unit = cp_print_key_unit_nr(logger, dft_section, & 1332 "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log") 1333 1334 IF (output_unit > 0) THEN 1335 1336 xc_section => section_vals_get_subs_vals(dft_section, "XC") 1337 1338 IF (dft_control%uks) THEN 1339 WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A)") & 1340 "DFT| Spin unrestricted (spin-polarized) Kohn-Sham calculation", "UKS" 1341 ELSE IF (dft_control%roks) THEN 1342 WRITE (UNIT=output_unit, FMT="(/,T2,A,T77,A)") & 1343 "DFT| Spin restricted open Kohn-Sham calculation", "ROKS" 1344 ELSE 1345 WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A)") & 1346 "DFT| Spin restricted Kohn-Sham (RKS) calculation", "RKS" 1347 END IF 1348 1349 WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") & 1350 "DFT| Multiplicity", dft_control%multiplicity 1351 WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") & 1352 "DFT| Number of spin states", dft_control%nspins 1353 1354 WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") & 1355 "DFT| Charge", dft_control%charge 1356 1357 IF (dft_control%sic_method_id .NE. sic_none) CALL cite_reference(VandeVondele2005b) 1358 SELECT CASE (dft_control%sic_method_id) 1359 CASE (sic_none) 1360 tmpstr = "NO" 1361 CASE (sic_mauri_spz) 1362 tmpstr = "SPZ/MAURI SIC" 1363 CASE (sic_mauri_us) 1364 tmpstr = "US/MAURI SIC" 1365 CASE (sic_ad) 1366 tmpstr = "AD SIC" 1367 CASE (sic_eo) 1368 tmpstr = "Explicit Orbital SIC" 1369 CASE DEFAULT 1370 ! fix throughout the cp2k for this option 1371 CPABORT("SIC option unknown") 1372 END SELECT 1373 1374 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") & 1375 "DFT| Self-interaction correction (SIC)", ADJUSTR(TRIM(tmpstr)) 1376 1377 IF (dft_control%sic_method_id /= sic_none) THEN 1378 WRITE (UNIT=output_unit, FMT="(T2,A,T66,ES15.6)") & 1379 "DFT| SIC scaling parameter a", dft_control%sic_scaling_a, & 1380 "DFT| SIC scaling parameter b", dft_control%sic_scaling_b 1381 END IF 1382 1383 IF (dft_control%sic_method_id == sic_eo) THEN 1384 IF (dft_control%sic_list_id == sic_list_all) THEN 1385 WRITE (UNIT=output_unit, FMT="(T2,A,T66,A)") & 1386 "DFT| SIC orbitals", "ALL" 1387 ENDIF 1388 IF (dft_control%sic_list_id == sic_list_unpaired) THEN 1389 WRITE (UNIT=output_unit, FMT="(T2,A,T66,A)") & 1390 "DFT| SIC orbitals", "UNPAIRED" 1391 ENDIF 1392 END IF 1393 1394 CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut) 1395 CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut) 1396 CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut) 1397 CALL section_vals_val_get(xc_section, "density_smooth_cutoff_range", r_val=density_smooth_cut_range) 1398 1399 WRITE (UNIT=output_unit, FMT="(T2,A,T66,ES15.6)") & 1400 "DFT| Cutoffs: density ", density_cut, & 1401 "DFT| gradient", gradient_cut, & 1402 "DFT| tau ", tau_cut, & 1403 "DFT| cutoff_smoothing_range", density_smooth_cut_range 1404 CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", & 1405 c_val=tmpStr) 1406 WRITE (output_unit, '( A, T61, A )') & 1407 " DFT| XC density smoothing ", ADJUSTR(tmpStr) 1408 CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", & 1409 c_val=tmpStr) 1410 WRITE (output_unit, '( A, T61, A )') & 1411 " DFT| XC derivatives ", ADJUSTR(tmpStr) 1412 IF (dft_control%dft_plus_u) THEN 1413 NULLIFY (enum, keyword, section) 1414 CALL create_dft_section(section) 1415 keyword => section_get_keyword(section, "PLUS_U_METHOD") 1416 CALL keyword_get(keyword, enum=enum) 1417 WRITE (UNIT=output_unit, FMT="(/,T2,A,T41,A40)") & 1418 "DFT+U| Method", ADJUSTR(TRIM(enum_i2c(enum, dft_control%plus_u_method_id))) 1419 WRITE (UNIT=output_unit, FMT="(T2,A)") & 1420 "DFT+U| Check atomic kind information for details" 1421 CALL section_release(section) 1422 END IF 1423 1424 CALL xc_write(output_unit, xc_section, dft_control%lsd) 1425 1426 IF (dft_control%do_sccs) THEN 1427 IF (dft_control%qs_control%gapw) THEN 1428 CPABORT("SCCS is not yet implemented with GAPW") 1429 END IF 1430 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 1431 "SCCS| Self-consistent continuum solvation model" 1432 SELECT CASE (dft_control%sccs_control%method_id) 1433 CASE (sccs_andreussi) 1434 WRITE (UNIT=output_unit, FMT="(T2,A,/,(T2,A,T69,ES12.3))") & 1435 "SCCS| Dielectric function proposed by Andreussi et al.", & 1436 "SCCS| rho_max", dft_control%sccs_control%rho_max, & 1437 "SCCS| rho_min", dft_control%sccs_control%rho_min 1438 CASE (sccs_fattebert_gygi) 1439 WRITE (UNIT=output_unit, FMT="(T2,A,/,(T2,A,T69,ES12.3))") & 1440 "SCCS| Dielectric function proposed by Fattebert and Gygi", & 1441 "SCCS| beta", dft_control%sccs_control%beta, & 1442 "SCCS| rho_zero", dft_control%sccs_control%rho_zero 1443 CASE DEFAULT 1444 CPABORT("Invalid SCCS model specified. Please, check your input!") 1445 END SELECT 1446 SELECT CASE (dft_control%sccs_control%derivative_method) 1447 CASE (sccs_derivative_fft) 1448 WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") "SCCS| Numerical derivative calculation", & 1449 ADJUSTR("FFT") 1450 CASE (sccs_derivative_cd3) 1451 WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") "SCCS| Numerical derivative calculation", & 1452 ADJUSTR("3-point stencil central differences") 1453 CASE (sccs_derivative_cd5) 1454 WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") "SCCS| Numerical derivative calculation", & 1455 ADJUSTR("5-point stencil central differences") 1456 CASE (sccs_derivative_cd7) 1457 WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") "SCCS| Numerical derivative calculation", & 1458 ADJUSTR("7-point stencil central differences") 1459 CASE DEFAULT 1460 CALL cp_abort(__LOCATION__, & 1461 "Invalid derivative method specified for SCCS model. "// & 1462 "Please, check your input!") 1463 END SELECT 1464 WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.3)") & 1465 "SCCS| Repulsion parameter alpha [mN/m] = [dyn/cm]", & 1466 cp_unit_from_cp2k(dft_control%sccs_control%alpha_solvent, "mN/m") 1467 WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.3)") & 1468 "SCCS| Dispersion parameter beta [GPa]", & 1469 cp_unit_from_cp2k(dft_control%sccs_control%beta_solvent, "GPa") 1470 WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.3)") & 1471 "SCCS| Surface tension gamma [mN/m] = [dyn/cm]", & 1472 cp_unit_from_cp2k(dft_control%sccs_control%gamma_solvent, "mN/m") 1473 WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.3)") & 1474 "SCCS| Mixing parameter applied during the iteration cycle", & 1475 dft_control%sccs_control%mixing 1476 WRITE (UNIT=output_unit, FMT="(T2,A,T69,ES12.3)") & 1477 "SCCS| Tolerance for the convergence of the SCCS iteration cycle", & 1478 dft_control%sccs_control%eps_sccs 1479 WRITE (UNIT=output_unit, FMT="(T2,A,T69,I12)") & 1480 "SCCS| Maximum number of iteration steps", & 1481 dft_control%sccs_control%max_iter 1482 WRITE (UNIT=output_unit, FMT="(T2,A,T69,ES12.3)") & 1483 "SCCS| SCF convergence threshold for starting the SCCS iteration", & 1484 dft_control%sccs_control%eps_scf 1485 WRITE (UNIT=output_unit, FMT="(T2,A,T69,ES12.3)") & 1486 "SCCS| Numerical increment for the cavity surface calculation", & 1487 dft_control%sccs_control%delta_rho 1488 END IF 1489 1490 END IF 1491 1492 CALL cp_print_key_finished_output(output_unit, logger, dft_section, & 1493 "PRINT%DFT_CONTROL_PARAMETERS") 1494 1495 CALL timestop(handle) 1496 1497 END SUBROUTINE write_dft_control 1498 1499! ************************************************************************************************** 1500!> \brief Write the xTB control parameters to the output unit. 1501!> \param xtb_control ... 1502!> \param dft_section ... 1503! ************************************************************************************************** 1504 SUBROUTINE write_xtb_control(xtb_control, dft_section) 1505 TYPE(xtb_control_type), POINTER :: xtb_control 1506 TYPE(section_vals_type), POINTER :: dft_section 1507 1508 CHARACTER(len=*), PARAMETER :: routineN = 'write_xtb_control', & 1509 routineP = moduleN//':'//routineN 1510 1511 INTEGER :: handle, output_unit 1512 TYPE(cp_logger_type), POINTER :: logger 1513 1514 CALL timeset(routineN, handle) 1515 NULLIFY (logger) 1516 logger => cp_get_default_logger() 1517 1518 output_unit = cp_print_key_unit_nr(logger, dft_section, & 1519 "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log") 1520 1521 IF (output_unit > 0) THEN 1522 1523 WRITE (UNIT=output_unit, FMT="(/,T2,A,T31,A50)") & 1524 "xTB| Parameter file", ADJUSTR(TRIM(xtb_control%parameter_file_name)) 1525 WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & 1526 "xTB| Basis expansion STO-NG", xtb_control%sto_ng 1527 WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & 1528 "xTB| Basis expansion STO-NG for Hydrogen", xtb_control%h_sto_ng 1529 WRITE (UNIT=output_unit, FMT="(T2,A,T71,L10)") & 1530 "xTB| Halogen interaction potential", xtb_control%xb_interaction 1531 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & 1532 "xTB| Halogen interaction potential cutoff radius", xtb_control%xb_radius 1533 WRITE (UNIT=output_unit, FMT="(T2,A,T31,A50)") & 1534 "xTB| D3 Dispersion: Parameter file", ADJUSTR(TRIM(xtb_control%dispersion_parameter_file)) 1535 WRITE (UNIT=output_unit, FMT="(T2,A,T51,3F10.3)") & 1536 "xTB| Huckel constants ks kp kd", xtb_control%ks, xtb_control%kp, xtb_control%kd 1537 WRITE (UNIT=output_unit, FMT="(T2,A,T61,2F10.3)") & 1538 "xTB| Huckel constants ksp k2sh", xtb_control%ksp, xtb_control%k2sh 1539 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & 1540 "xTB| Mataga-Nishimoto exponent", xtb_control%kg 1541 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & 1542 "xTB| Repulsion potential exponent", xtb_control%kf 1543 WRITE (UNIT=output_unit, FMT="(T2,A,T51,3F10.3)") & 1544 "xTB| Coordination number scaling kcn(s) kcn(p) kcn(d)", & 1545 xtb_control%kcns, xtb_control%kcnp, xtb_control%kcnd 1546 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & 1547 "xTB| Electronegativity scaling", xtb_control%ken 1548 WRITE (UNIT=output_unit, FMT="(T2,A,T61,2F10.3)") & 1549 "xTB| Halogen potential scaling kxr kx2", xtb_control%kxr, xtb_control%kx2 1550 WRITE (UNIT=output_unit, FMT="(/)") 1551 1552 END IF 1553 1554 CALL cp_print_key_finished_output(output_unit, logger, dft_section, & 1555 "PRINT%DFT_CONTROL_PARAMETERS") 1556 1557 CALL timestop(handle) 1558 1559 END SUBROUTINE write_xtb_control 1560 1561! ************************************************************************************************** 1562!> \brief Purpose: Write the QS control parameters to the output unit. 1563!> \param qs_control ... 1564!> \param dft_section ... 1565! ************************************************************************************************** 1566 SUBROUTINE write_qs_control(qs_control, dft_section) 1567 TYPE(qs_control_type), INTENT(IN) :: qs_control 1568 TYPE(section_vals_type), POINTER :: dft_section 1569 1570 CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_control', & 1571 routineP = moduleN//':'//routineN 1572 1573 INTEGER :: handle, i, igrid_level, ngrid_level, & 1574 output_unit 1575 TYPE(cp_logger_type), POINTER :: logger 1576 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control 1577 TYPE(enumeration_type), POINTER :: enum 1578 TYPE(keyword_type), POINTER :: keyword 1579 TYPE(section_type), POINTER :: section 1580 1581 IF (qs_control%semi_empirical) RETURN 1582 IF (qs_control%dftb) RETURN 1583 IF (qs_control%xtb) RETURN 1584 CALL timeset(routineN, handle) 1585 NULLIFY (logger, enum, keyword, section) 1586 logger => cp_get_default_logger() 1587 1588 CALL create_qs_section(section) 1589 keyword => section_get_keyword(section, "QUADRATURE") 1590 CALL keyword_get(keyword, enum=enum) 1591 1592 output_unit = cp_print_key_unit_nr(logger, dft_section, & 1593 "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log") 1594 IF (output_unit > 0) THEN 1595 ngrid_level = SIZE(qs_control%e_cutoff) 1596 WRITE (UNIT=output_unit, FMT="(/,T2,A,T71,A)") & 1597 "QS| Method:", ADJUSTR(qs_control%method) 1598 IF (qs_control%pw_grid_opt%spherical) THEN 1599 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A)") & 1600 "QS| Density plane wave grid type", " SPHERICAL HALFSPACE" 1601 ELSE IF (qs_control%pw_grid_opt%fullspace) THEN 1602 WRITE (UNIT=output_unit, FMT="(T2,A,T57,A)") & 1603 "QS| Density plane wave grid type", " NON-SPHERICAL FULLSPACE" 1604 ELSE 1605 WRITE (UNIT=output_unit, FMT="(T2,A,T57,A)") & 1606 "QS| Density plane wave grid type", " NON-SPHERICAL HALFSPACE" 1607 END IF 1608 WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & 1609 "QS| Number of grid levels:", SIZE(qs_control%e_cutoff) 1610 IF (ngrid_level == 1) THEN 1611 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") & 1612 "QS| Density cutoff [a.u.]:", qs_control%e_cutoff(1) 1613 ELSE 1614 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") & 1615 "QS| Density cutoff [a.u.]:", qs_control%cutoff 1616 IF (qs_control%commensurate_mgrids) & 1617 WRITE (UNIT=output_unit, FMT="(T2,A)") "QS| Using commensurate multigrids" 1618 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") & 1619 "QS| Multi grid cutoff [a.u.]: 1) grid level", qs_control%e_cutoff(1) 1620 WRITE (UNIT=output_unit, FMT="(T2,A,I3,A,T71,F10.1)") & 1621 ("QS| ", igrid_level, ") grid level", & 1622 qs_control%e_cutoff(igrid_level), & 1623 igrid_level=2, SIZE(qs_control%e_cutoff)) 1624 END IF 1625 IF (qs_control%pao) THEN 1626 WRITE (UNIT=output_unit, FMT="(T2,A)") "QS| PAO active" 1627 END IF 1628 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") & 1629 "QS| Grid level progression factor:", qs_control%progression_factor 1630 WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") & 1631 "QS| Relative density cutoff [a.u.]:", qs_control%relative_cutoff 1632 IF (qs_control%map_consistent) THEN 1633 WRITE (UNIT=output_unit, FMT="(T2,A)") & 1634 "QS| Consistent realspace mapping and integration " 1635 ENDIF 1636 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1637 "QS| Interaction thresholds: eps_pgf_orb:", & 1638 qs_control%eps_pgf_orb, & 1639 "QS| eps_filter_matrix:", & 1640 qs_control%eps_filter_matrix, & 1641 "QS| eps_core_charge:", & 1642 qs_control%eps_core_charge, & 1643 "QS| eps_rho_gspace:", & 1644 qs_control%eps_rho_gspace, & 1645 "QS| eps_rho_rspace:", & 1646 qs_control%eps_rho_rspace, & 1647 "QS| eps_gvg_rspace:", & 1648 qs_control%eps_gvg_rspace, & 1649 "QS| eps_ppl:", & 1650 qs_control%eps_ppl, & 1651 "QS| eps_ppnl:", & 1652 qs_control%eps_ppnl 1653 IF (qs_control%gapw) THEN 1654 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1655 "QS| GAPW| eps_fit:", & 1656 qs_control%gapw_control%eps_fit, & 1657 "QS| GAPW| eps_iso:", & 1658 qs_control%gapw_control%eps_iso, & 1659 "QS| GAPW| eps_svd:", & 1660 qs_control%gapw_control%eps_svd, & 1661 "QS| GAPW| eps_cpc:", & 1662 qs_control%gapw_control%eps_cpc 1663 WRITE (UNIT=output_unit, FMT="(T2,A,T55,A30)") & 1664 "QS| GAPW| atom-r-grid: quadrature:", & 1665 enum_i2c(enum, qs_control%gapw_control%quadrature) 1666 WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & 1667 "QS| GAPW| atom-s-grid: max l :", & 1668 qs_control%gapw_control%lmax_sphere, & 1669 "QS| GAPW| max_l_rho0 :", & 1670 qs_control%gapw_control%lmax_rho0 1671 IF (qs_control%gapw_control%non_paw_atoms) THEN 1672 WRITE (UNIT=output_unit, FMT="(T2,A)") & 1673 "QS| GAPW| At least one kind is NOT PAW, i.e. it has only soft AO " 1674 END IF 1675 IF (qs_control%gapw_control%nopaw_as_gpw) THEN 1676 WRITE (UNIT=output_unit, FMT="(T2,A)") & 1677 "QS| GAPW| The NOT PAW atoms are treated fully GPW" 1678 END IF 1679 END IF 1680 IF (qs_control%gapw_xc) THEN 1681 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1682 "QS| GAPW_XC| eps_fit:", & 1683 qs_control%gapw_control%eps_fit, & 1684 "QS| GAPW_XC| eps_iso:", & 1685 qs_control%gapw_control%eps_iso, & 1686 "QS| GAPW_XC| eps_svd:", & 1687 qs_control%gapw_control%eps_svd 1688 WRITE (UNIT=output_unit, FMT="(T2,A,T55,A30)") & 1689 "QS| GAPW_XC|atom-r-grid: quadrature:", & 1690 enum_i2c(enum, qs_control%gapw_control%quadrature) 1691 WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & 1692 "QS| GAPW_XC| atom-s-grid: max l :", & 1693 qs_control%gapw_control%lmax_sphere 1694 END IF 1695 IF (qs_control%mulliken_restraint) THEN 1696 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1697 "QS| Mulliken restraint target", qs_control%mulliken_restraint_control%target 1698 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1699 "QS| Mulliken restraint strength", qs_control%mulliken_restraint_control%strength 1700 WRITE (UNIT=output_unit, FMT="(T2,A,T73,I8)") & 1701 "QS| Mulliken restraint atoms: ", qs_control%mulliken_restraint_control%natoms 1702 WRITE (UNIT=output_unit, FMT="(5I8)") qs_control%mulliken_restraint_control%atoms 1703 END IF 1704 IF (qs_control%ddapc_restraint) THEN 1705 DO i = 1, SIZE(qs_control%ddapc_restraint_control) 1706 ddapc_restraint_control => qs_control%ddapc_restraint_control(i)%ddapc_restraint_control 1707 IF (SIZE(qs_control%ddapc_restraint_control) .GT. 1) & 1708 WRITE (UNIT=output_unit, FMT="(T2,A,T3,I8)") & 1709 "QS| parameters for DDAPC restraint number", i 1710 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1711 "QS| ddapc restraint target", ddapc_restraint_control%target 1712 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1713 "QS| ddapc restraint strength", ddapc_restraint_control%strength 1714 WRITE (UNIT=output_unit, FMT="(T2,A,T73,I8)") & 1715 "QS| ddapc restraint atoms: ", ddapc_restraint_control%natoms 1716 WRITE (UNIT=output_unit, FMT="(5I8)") ddapc_restraint_control%atoms 1717 WRITE (UNIT=output_unit, FMT="(T2,A)") "Coefficients:" 1718 WRITE (UNIT=output_unit, FMT="(5F6.2)") ddapc_restraint_control%coeff 1719 SELECT CASE (ddapc_restraint_control%functional_form) 1720 CASE (do_ddapc_restraint) 1721 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") & 1722 "QS| ddapc restraint functional form :", "RESTRAINT" 1723 CASE (do_ddapc_constraint) 1724 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") & 1725 "QS| ddapc restraint functional form :", "CONSTRAINT" 1726 CASE DEFAULT 1727 CPABORT("Unknown ddapc restraint") 1728 END SELECT 1729 END DO 1730 END IF 1731 IF (qs_control%s2_restraint) THEN 1732 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1733 "QS| s2 restraint target", qs_control%s2_restraint_control%target 1734 WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") & 1735 "QS| s2 restraint strength", qs_control%s2_restraint_control%strength 1736 SELECT CASE (qs_control%s2_restraint_control%functional_form) 1737 CASE (do_s2_restraint) 1738 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") & 1739 "QS| s2 restraint functional form :", "RESTRAINT" 1740 CPABORT("Not yet implemented") 1741 CASE (do_s2_constraint) 1742 WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") & 1743 "QS| s2 restraint functional form :", "CONSTRAINT" 1744 CASE DEFAULT 1745 CPABORT("Unknown ddapc restraint") 1746 END SELECT 1747 END IF 1748 END IF 1749 CALL cp_print_key_finished_output(output_unit, logger, dft_section, & 1750 "PRINT%DFT_CONTROL_PARAMETERS") 1751 CALL section_release(section) 1752 1753 CALL timestop(handle) 1754 1755 END SUBROUTINE write_qs_control 1756 1757! ************************************************************************************************** 1758!> \brief reads the input parameters needed for ddapc. 1759!> \param qs_control ... 1760!> \param qs_section ... 1761!> \param ddapc_restraint_section ... 1762!> \author fschiff 1763!> \note 1764!> either reads DFT%QS%DDAPC_RESTRAINT or PROPERTIES%ET_coupling 1765!> if(qs_section is present the DFT part is read, if ddapc_restraint_section 1766!> is present ET_COUPLING is read. Avoid having both!!! 1767! ************************************************************************************************** 1768 SUBROUTINE read_ddapc_section(qs_control, qs_section, ddapc_restraint_section) 1769 1770 TYPE(qs_control_type), INTENT(INOUT) :: qs_control 1771 TYPE(section_vals_type), OPTIONAL, POINTER :: qs_section, ddapc_restraint_section 1772 1773 CHARACTER(len=*), PARAMETER :: routineN = 'read_ddapc_section', & 1774 routineP = moduleN//':'//routineN 1775 1776 INTEGER :: i, j, jj, k, n_rep 1777 INTEGER, DIMENSION(:), POINTER :: tmplist 1778 REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist 1779 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control 1780 TYPE(section_vals_type), POINTER :: ddapc_section 1781 1782 IF (PRESENT(ddapc_restraint_section)) THEN 1783 IF (ASSOCIATED(qs_control%ddapc_restraint_control)) THEN 1784 IF (SIZE(qs_control%ddapc_restraint_control) .GE. 2) & 1785 CPABORT("ET_COUPLING cannot be used in combination with a normal restraint") 1786 ELSE 1787 ddapc_section => ddapc_restraint_section 1788 ALLOCATE (qs_control%ddapc_restraint_control(1)) 1789 END IF 1790 END IF 1791 1792 IF (PRESENT(qs_section)) THEN 1793 NULLIFY (ddapc_section) 1794 ddapc_section => section_vals_get_subs_vals(qs_section, & 1795 "DDAPC_RESTRAINT") 1796 END IF 1797 1798 DO i = 1, SIZE(qs_control%ddapc_restraint_control) 1799 1800 NULLIFY (qs_control%ddapc_restraint_control(i)%ddapc_restraint_control) 1801 CALL ddapc_control_create(qs_control%ddapc_restraint_control(i)%ddapc_restraint_control) 1802 ddapc_restraint_control => qs_control%ddapc_restraint_control(i)%ddapc_restraint_control 1803 1804 CALL section_vals_val_get(ddapc_section, "STRENGTH", i_rep_section=i, & 1805 r_val=ddapc_restraint_control%strength) 1806 CALL section_vals_val_get(ddapc_section, "TARGET", i_rep_section=i, & 1807 r_val=ddapc_restraint_control%target) 1808 CALL section_vals_val_get(ddapc_section, "FUNCTIONAL_FORM", i_rep_section=i, & 1809 i_val=ddapc_restraint_control%functional_form) 1810 CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, & 1811 n_rep_val=n_rep) 1812 CALL section_vals_val_get(ddapc_section, "TYPE_OF_DENSITY", i_rep_section=i, & 1813 i_val=ddapc_restraint_control%density_type) 1814 1815 jj = 0 1816 DO k = 1, n_rep 1817 CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, & 1818 i_rep_val=k, i_vals=tmplist) 1819 DO j = 1, SIZE(tmplist) 1820 jj = jj + 1 1821 END DO 1822 END DO 1823 IF (jj < 1) CPABORT("Need at least 1 atom to use ddapc contraints") 1824 ddapc_restraint_control%natoms = jj 1825 IF (ASSOCIATED(ddapc_restraint_control%atoms)) & 1826 DEALLOCATE (ddapc_restraint_control%atoms) 1827 ALLOCATE (ddapc_restraint_control%atoms(ddapc_restraint_control%natoms)) 1828 jj = 0 1829 DO k = 1, n_rep 1830 CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, & 1831 i_rep_val=k, i_vals=tmplist) 1832 DO j = 1, SIZE(tmplist) 1833 jj = jj + 1 1834 ddapc_restraint_control%atoms(jj) = tmplist(j) 1835 END DO 1836 END DO 1837 1838 IF (ASSOCIATED(ddapc_restraint_control%coeff)) & 1839 DEALLOCATE (ddapc_restraint_control%coeff) 1840 ALLOCATE (ddapc_restraint_control%coeff(ddapc_restraint_control%natoms)) 1841 ddapc_restraint_control%coeff = 1.0_dp 1842 1843 CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, & 1844 n_rep_val=n_rep) 1845 jj = 0 1846 DO k = 1, n_rep 1847 CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, & 1848 i_rep_val=k, r_vals=rtmplist) 1849 DO j = 1, SIZE(rtmplist) 1850 jj = jj + 1 1851 IF (jj > ddapc_restraint_control%natoms) & 1852 CPABORT("Need the same number of coeff as there are atoms ") 1853 ddapc_restraint_control%coeff(jj) = rtmplist(j) 1854 END DO 1855 END DO 1856 IF (jj < ddapc_restraint_control%natoms .AND. jj .NE. 0) & 1857 CPABORT("Need no or the same number of coeff as there are atoms.") 1858 END DO 1859 k = 0 1860 DO i = 1, SIZE(qs_control%ddapc_restraint_control) 1861 IF (qs_control%ddapc_restraint_control(i)%ddapc_restraint_control%functional_form == & 1862 do_ddapc_constraint) k = k + 1 1863 END DO 1864 IF (k == 2) CALL cp_abort(__LOCATION__, & 1865 "Only a single constraint possible yet, try to use restraints instead ") 1866 1867 END SUBROUTINE read_ddapc_section 1868 1869! ************************************************************************************************** 1870!> \brief ... 1871!> \param dft_control ... 1872!> \param efield_section ... 1873! ************************************************************************************************** 1874 SUBROUTINE read_efield_sections(dft_control, efield_section) 1875 TYPE(dft_control_type), POINTER :: dft_control 1876 TYPE(section_vals_type), POINTER :: efield_section 1877 1878 CHARACTER(len=*), PARAMETER :: routineN = 'read_efield_sections', & 1879 routineP = moduleN//':'//routineN 1880 1881 CHARACTER(len=default_path_length) :: file_name 1882 INTEGER :: i, io, j, n, unit_nr 1883 REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_vals 1884 TYPE(efield_type), POINTER :: efield 1885 TYPE(section_vals_type), POINTER :: tmp_section 1886 1887 DO i = 1, SIZE(dft_control%efield_fields) 1888 NULLIFY (dft_control%efield_fields(i)%efield) 1889 ALLOCATE (dft_control%efield_fields(i)%efield) 1890 efield => dft_control%efield_fields(i)%efield 1891 NULLIFY (efield%envelop_i_vars, efield%envelop_r_vars) 1892 CALL section_vals_val_get(efield_section, "INTENSITY", i_rep_section=i, & 1893 r_val=efield%strength) 1894 1895 CALL section_vals_val_get(efield_section, "POLARISATION", i_rep_section=i, & 1896 r_vals=tmp_vals) 1897 ALLOCATE (efield%polarisation(SIZE(tmp_vals))) 1898 efield%polarisation = tmp_vals 1899 CALL section_vals_val_get(efield_section, "PHASE", i_rep_section=i, & 1900 r_val=efield%phase_offset) 1901 CALL section_vals_val_get(efield_section, "ENVELOP", i_rep_section=i, & 1902 i_val=efield%envelop_id) 1903 CALL section_vals_val_get(efield_section, "WAVELENGTH", i_rep_section=i, & 1904 r_val=efield%wavelength) 1905 1906 IF (efield%envelop_id == constant_env) THEN 1907 ALLOCATE (efield%envelop_i_vars(2)) 1908 tmp_section => section_vals_get_subs_vals(efield_section, "CONSTANT_ENV", i_rep_section=i) 1909 CALL section_vals_val_get(tmp_section, "START_STEP", & 1910 i_val=efield%envelop_i_vars(1)) 1911 CALL section_vals_val_get(tmp_section, "END_STEP", & 1912 i_val=efield%envelop_i_vars(2)) 1913 ELSE IF (efield%envelop_id == gaussian_env) THEN 1914 ALLOCATE (efield%envelop_r_vars(2)) 1915 tmp_section => section_vals_get_subs_vals(efield_section, "GAUSSIAN_ENV", i_rep_section=i) 1916 CALL section_vals_val_get(tmp_section, "T0", & 1917 r_val=efield%envelop_r_vars(1)) 1918 CALL section_vals_val_get(tmp_section, "SIGMA", & 1919 r_val=efield%envelop_r_vars(2)) 1920 ELSE IF (efield%envelop_id == ramp_env) THEN 1921 ALLOCATE (efield%envelop_i_vars(4)) 1922 tmp_section => section_vals_get_subs_vals(efield_section, "RAMP_ENV", i_rep_section=i) 1923 CALL section_vals_val_get(tmp_section, "START_STEP_IN", & 1924 i_val=efield%envelop_i_vars(1)) 1925 CALL section_vals_val_get(tmp_section, "END_STEP_IN", & 1926 i_val=efield%envelop_i_vars(2)) 1927 CALL section_vals_val_get(tmp_section, "START_STEP_OUT", & 1928 i_val=efield%envelop_i_vars(3)) 1929 CALL section_vals_val_get(tmp_section, "END_STEP_OUT", & 1930 i_val=efield%envelop_i_vars(4)) 1931 ELSE IF (efield%envelop_id == custom_env) THEN 1932 tmp_section => section_vals_get_subs_vals(efield_section, "CUSTOM_ENV", i_rep_section=i) 1933 CALL section_vals_val_get(tmp_section, "EFIELD_FILE_NAME", c_val=file_name) 1934 CALL open_file(file_name=TRIM(file_name), file_action="READ", file_status="OLD", unit_number=unit_nr) 1935 !Determine the number of lines in file 1936 n = 0 1937 DO WHILE (.TRUE.) 1938 READ (unit_nr, *, iostat=io) 1939 IF (io /= 0) EXIT 1940 n = n + 1 1941 END DO 1942 REWIND (unit_nr) 1943 ALLOCATE (efield%envelop_r_vars(n + 1)) 1944 !Store the timestep of the list in the first entry of the r_vars 1945 CALL section_vals_val_get(tmp_section, "TIMESTEP", r_val=efield%envelop_r_vars(1)) 1946 !Read the file 1947 DO j = 2, n + 1 1948 READ (unit_nr, *) efield%envelop_r_vars(j) 1949 efield%envelop_r_vars(j) = cp_unit_to_cp2k(efield%envelop_r_vars(j), "volt/m") 1950 END DO 1951 CALL close_file(unit_nr) 1952 END IF 1953 END DO 1954 END SUBROUTINE read_efield_sections 1955 1956! ************************************************************************************************** 1957!> \brief reads the input parameters needed real time propagation 1958!> \param dft_control ... 1959!> \param rtp_section ... 1960!> \author fschiff 1961! ************************************************************************************************** 1962 SUBROUTINE read_rtp_section(dft_control, rtp_section) 1963 1964 TYPE(dft_control_type), INTENT(INOUT) :: dft_control 1965 TYPE(section_vals_type), POINTER :: rtp_section 1966 1967 CHARACTER(len=*), PARAMETER :: routineN = 'read_rtp_section', & 1968 routineP = moduleN//':'//routineN 1969 1970 INTEGER, DIMENSION(:), POINTER :: tmp 1971 1972 ALLOCATE (dft_control%rtp_control) 1973 CALL section_vals_val_get(rtp_section, "MAX_ITER", & 1974 i_val=dft_control%rtp_control%max_iter) 1975 CALL section_vals_val_get(rtp_section, "MAT_EXP", & 1976 i_val=dft_control%rtp_control%mat_exp) 1977 CALL section_vals_val_get(rtp_section, "ASPC_ORDER", & 1978 i_val=dft_control%rtp_control%aspc_order) 1979 CALL section_vals_val_get(rtp_section, "EXP_ACCURACY", & 1980 r_val=dft_control%rtp_control%eps_exp) 1981 CALL section_vals_val_get(rtp_section, "PROPAGATOR", & 1982 i_val=dft_control%rtp_control%propagator) 1983 CALL section_vals_val_get(rtp_section, "EPS_ITER", & 1984 r_val=dft_control%rtp_control%eps_ener) 1985 CALL section_vals_val_get(rtp_section, "INITIAL_WFN", & 1986 i_val=dft_control%rtp_control%initial_wfn) 1987 CALL section_vals_val_get(rtp_section, "HFX_BALANCE_IN_CORE", & 1988 l_val=dft_control%rtp_control%hfx_redistribute) 1989 CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE", & 1990 l_val=dft_control%rtp_control%apply_delta_pulse) 1991 CALL section_vals_val_get(rtp_section, "PERIODIC", & 1992 l_val=dft_control%rtp_control%periodic) 1993 CALL section_vals_val_get(rtp_section, "DENSITY_PROPAGATION", & 1994 l_val=dft_control%rtp_control%linear_scaling) 1995 CALL section_vals_val_get(rtp_section, "MCWEENY_MAX_ITER", & 1996 i_val=dft_control%rtp_control%mcweeny_max_iter) 1997 CALL section_vals_val_get(rtp_section, "ACCURACY_REFINEMENT", & 1998 i_val=dft_control%rtp_control%acc_ref) 1999 CALL section_vals_val_get(rtp_section, "MCWEENY_EPS", & 2000 r_val=dft_control%rtp_control%mcweeny_eps) 2001 CALL section_vals_val_get(rtp_section, "DELTA_PULSE_SCALE", & 2002 r_val=dft_control%rtp_control%delta_pulse_scale) 2003 CALL section_vals_val_get(rtp_section, "DELTA_PULSE_DIRECTION", & 2004 i_vals=tmp) 2005 dft_control%rtp_control%delta_pulse_direction = tmp 2006 CALL section_vals_val_get(rtp_section, "SC_CHECK_START", & 2007 i_val=dft_control%rtp_control%sc_check_start) 2008 2009 END SUBROUTINE read_rtp_section 2010 2011! ************************************************************************************************** 2012!> \brief Parses the BLOCK_LIST keywords from the ADMM section 2013!> \param admm_control ... 2014!> \param dft_section ... 2015! ************************************************************************************************** 2016 SUBROUTINE read_admm_block_list(admm_control, dft_section) 2017 TYPE(admm_control_type), POINTER :: admm_control 2018 TYPE(section_vals_type), POINTER :: dft_section 2019 2020 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_admm_block_list', & 2021 routineP = moduleN//':'//routineN 2022 2023 INTEGER :: irep, list_size, n_rep 2024 INTEGER, DIMENSION(:), POINTER :: tmplist 2025 2026 NULLIFY (tmplist) 2027 2028 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", & 2029 n_rep_val=n_rep) 2030 2031 ALLOCATE (admm_control%blocks(n_rep)) 2032 2033 DO irep = 1, n_rep 2034 CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", & 2035 i_rep_val=irep, i_vals=tmplist) 2036 list_size = SIZE(tmplist) 2037 ALLOCATE (admm_control%blocks(irep)%list(list_size)) 2038 admm_control%blocks(irep)%list(:) = tmplist(:) 2039 END DO 2040 2041 END SUBROUTINE read_admm_block_list 2042 2043END MODULE cp_control_utils 2044