1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE atom_pseudo 8 USE atom_electronic_structure, ONLY: calculate_atom 9 USE atom_fit, ONLY: atom_fit_pseudo 10 USE atom_operators, ONLY: atom_int_release,& 11 atom_int_setup,& 12 atom_ppint_release,& 13 atom_ppint_setup,& 14 atom_relint_release,& 15 atom_relint_setup 16 USE atom_output, ONLY: atom_print_basis,& 17 atom_print_info,& 18 atom_print_method,& 19 atom_print_potential 20 USE atom_types, ONLY: & 21 atom_basis_type, atom_integrals, atom_optimization_type, atom_orbitals, atom_p_type, & 22 atom_potential_type, atom_state, create_atom_orbs, create_atom_type, init_atom_basis, & 23 init_atom_potential, lmat, read_atom_opt_section, release_atom_basis, & 24 release_atom_potential, release_atom_type, set_atom 25 USE atom_utils, ONLY: atom_consistent_method,& 26 atom_set_occupation,& 27 get_maxl_occ,& 28 get_maxn_occ 29 USE cp_log_handling, ONLY: cp_get_default_logger,& 30 cp_logger_type 31 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 32 cp_print_key_unit_nr 33 USE input_constants, ONLY: do_analytic,& 34 poly_conf 35 USE input_section_types, ONLY: section_vals_get,& 36 section_vals_get_subs_vals,& 37 section_vals_type,& 38 section_vals_val_get 39 USE kinds, ONLY: default_string_length,& 40 dp 41 USE periodic_table, ONLY: nelem,& 42 ptable 43 USE physcon, ONLY: bohr 44#include "./base/base_uses.f90" 45 46 IMPLICIT NONE 47 PRIVATE 48 PUBLIC :: atom_pseudo_opt 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_pseudo' 51 52! ************************************************************************************************** 53 54CONTAINS 55 56! ************************************************************************************************** 57 58! ************************************************************************************************** 59!> \brief ... 60!> \param atom_section ... 61! ************************************************************************************************** 62 SUBROUTINE atom_pseudo_opt(atom_section) 63 TYPE(section_vals_type), POINTER :: atom_section 64 65 CHARACTER(len=*), PARAMETER :: routineN = 'atom_pseudo_opt', & 66 routineP = moduleN//':'//routineN 67 68 CHARACTER(LEN=2) :: elem 69 CHARACTER(LEN=default_string_length), & 70 DIMENSION(:), POINTER :: tmpstringlist 71 INTEGER :: ads, do_eric, do_erie, handle, i, im, & 72 in, iw, k, l, maxl, mb, method, mo, & 73 n_meth, n_rep, nr_gh, reltyp, zcore, & 74 zval, zz 75 INTEGER, DIMENSION(0:lmat) :: maxn 76 INTEGER, DIMENSION(:), POINTER :: cn 77 LOGICAL :: do_gh, eri_c, eri_e, pp_calc 78 REAL(KIND=dp) :: ne, nm 79 REAL(KIND=dp), DIMENSION(0:lmat, 10) :: pocc 80 TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis 81 TYPE(atom_integrals), POINTER :: ae_int, pp_int 82 TYPE(atom_optimization_type) :: optimization 83 TYPE(atom_orbitals), POINTER :: orbitals 84 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs 85 TYPE(atom_potential_type), POINTER :: ae_pot, p_pot 86 TYPE(atom_state), POINTER :: state, statepp 87 TYPE(cp_logger_type), POINTER :: logger 88 TYPE(section_vals_type), POINTER :: basis_section, method_section, & 89 opt_section, potential_section, & 90 powell_section, xc_section 91 92 CALL timeset(routineN, handle) 93 94 ! What atom do we calculate 95 CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval) 96 CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem) 97 zz = 0 98 DO i = 1, nelem 99 IF (ptable(i)%symbol == elem) THEN 100 zz = i 101 EXIT 102 END IF 103 END DO 104 IF (zz /= 1) zval = zz 105 106 ! read and set up information on the basis sets 107 ALLOCATE (ae_basis, pp_basis) 108 basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS") 109 NULLIFY (ae_basis%grid) 110 CALL init_atom_basis(ae_basis, basis_section, zval, "AA") 111 NULLIFY (pp_basis%grid) 112 basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS") 113 CALL init_atom_basis(pp_basis, basis_section, zval, "AP") 114 115 ! print general and basis set information 116 logger => cp_get_default_logger() 117 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log") 118 IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw) 119 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER") 120 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log") 121 IF (iw > 0) THEN 122 CALL atom_print_basis(ae_basis, iw, " All Electron Basis") 123 CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis") 124 END IF 125 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET") 126 127 ! read and setup information on the pseudopotential 128 NULLIFY (potential_section) 129 potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL") 130 ALLOCATE (ae_pot, p_pot) 131 CALL init_atom_potential(p_pot, potential_section, zval) 132 CALL init_atom_potential(ae_pot, potential_section, -1) 133 IF (.NOT. p_pot%confinement .AND. .NOT. ae_pot%confinement) THEN 134 !set default confinement potential 135 p_pot%confinement = .TRUE. 136 p_pot%conf_type = poly_conf 137 p_pot%scon = 2.0_dp 138 p_pot%acon = 0.5_dp 139 ! this seems to be the default in the old code 140 p_pot%rcon = (2._dp*ptable(zval)%covalent_radius*bohr)**2 141 ae_pot%confinement = .TRUE. 142 ae_pot%conf_type = poly_conf 143 ae_pot%scon = 2.0_dp 144 ae_pot%acon = 0.5_dp 145 ! this seems to be the default in the old code 146 ae_pot%rcon = (2._dp*ptable(zval)%covalent_radius*bohr)**2 147 END IF 148 149 ! if the ERI's are calculated analytically, we have to precalculate them 150 eri_c = .FALSE. 151 CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric) 152 IF (do_eric == do_analytic) eri_c = .TRUE. 153 eri_e = .FALSE. 154 CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie) 155 IF (do_erie == do_analytic) eri_e = .TRUE. 156 CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh) 157 CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh) 158 159 ! information on the states to be calculated 160 CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl) 161 maxn = 0 162 CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn) 163 DO in = 1, MIN(SIZE(cn), 4) 164 maxn(in - 1) = cn(in) 165 END DO 166 DO in = 0, lmat 167 maxn(in) = MIN(maxn(in), ae_basis%nbas(in)) 168 END DO 169 170 ! read optimization section 171 opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION") 172 CALL read_atom_opt_section(optimization, opt_section) 173 174 ! Check for the total number of electron configurations to be calculated 175 CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep) 176 ! Check for the total number of method types to be calculated 177 method_section => section_vals_get_subs_vals(atom_section, "METHOD") 178 CALL section_vals_get(method_section, n_repetition=n_meth) 179 180 ! integrals 181 ALLOCATE (ae_int, pp_int) 182 183 ALLOCATE (atom_info(n_rep, n_meth), atom_refs(n_rep, n_meth)) 184 185 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log") 186 IF (iw > 0) THEN 187 WRITE (iw, '(/," ",79("*"))') 188 WRITE (iw, '(" ",26("*"),A,25("*"))') " Calculate Reference States " 189 WRITE (iw, '(" ",79("*"))') 190 END IF 191 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER") 192 193 DO in = 1, n_rep 194 DO im = 1, n_meth 195 196 NULLIFY (atom_info(in, im)%atom, atom_refs(in, im)%atom) 197 CALL create_atom_type(atom_info(in, im)%atom) 198 CALL create_atom_type(atom_refs(in, im)%atom) 199 200 atom_info(in, im)%atom%optimization = optimization 201 atom_refs(in, im)%atom%optimization = optimization 202 203 atom_info(in, im)%atom%z = zval 204 atom_refs(in, im)%atom%z = zval 205 xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im) 206 atom_info(in, im)%atom%xc_section => xc_section 207 atom_refs(in, im)%atom%xc_section => xc_section 208 209 ALLOCATE (state, statepp) 210 211 ! get the electronic configuration 212 CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, & 213 c_vals=tmpstringlist) 214 ! all electron configurations have to be with full core 215 pp_calc = INDEX(tmpstringlist(1), "CORE") /= 0 216 CPASSERT(.NOT. pp_calc) 217 218 ! set occupations 219 CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity) 220 state%maxl_occ = get_maxl_occ(state%occ) 221 state%maxn_occ = get_maxn_occ(state%occ) 222 ! set number of states to be calculated 223 state%maxl_calc = MAX(maxl, state%maxl_occ) 224 state%maxl_calc = MIN(lmat, state%maxl_calc) 225 state%maxn_calc = 0 226 DO k = 0, state%maxl_calc 227 ads = 2 228 IF (state%maxn_occ(k) == 0) ads = 1 229 state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k) + ads) 230 state%maxn_calc(k) = MIN(state%maxn_calc(k), ae_basis%nbas(k)) 231 END DO 232 state%core = 0._dp 233 CALL set_atom(atom_refs(in, im)%atom, zcore=zval, pp_calc=.FALSE.) 234 235 IF (state%multiplicity /= -1) THEN 236 ! set alpha and beta occupations 237 state%occa = 0._dp 238 state%occb = 0._dp 239 DO l = 0, lmat 240 nm = REAL((2*l + 1), KIND=dp) 241 DO k = 1, 10 242 ne = state%occupation(l, k) 243 IF (ne == 0._dp) THEN !empty shell 244 EXIT !assume there are no holes 245 ELSEIF (ne == 2._dp*nm) THEN !closed shell 246 state%occa(l, k) = nm 247 state%occb(l, k) = nm 248 ELSEIF (state%multiplicity == -2) THEN !High spin case 249 state%occa(l, k) = MIN(ne, nm) 250 state%occb(l, k) = MAX(0._dp, ne - nm) 251 ELSE 252 state%occa(l, k) = 0.5_dp*(ne + state%multiplicity - 1._dp) 253 state%occb(l, k) = ne - state%occa(l, k) 254 END IF 255 END DO 256 END DO 257 END IF 258 259 ! set occupations for pseudopotential calculation 260 CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist) 261 CALL atom_set_occupation(tmpstringlist, statepp%core, pocc) 262 zcore = zval - NINT(SUM(statepp%core)) 263 CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.) 264 265 statepp%occ = state%occ - statepp%core 266 statepp%occupation = 0._dp 267 DO l = 0, lmat 268 k = 0 269 DO i = 1, 10 270 IF (statepp%occ(l, i) /= 0._dp) THEN 271 k = k + 1 272 statepp%occupation(l, k) = state%occ(l, i) 273 IF (state%multiplicity /= -1) THEN 274 statepp%occa(l, k) = state%occa(l, i) - statepp%core(l, i)/2 275 statepp%occb(l, k) = state%occb(l, i) - statepp%core(l, i)/2 276 END IF 277 END IF 278 END DO 279 END DO 280 281 statepp%maxl_occ = get_maxl_occ(statepp%occ) 282 statepp%maxn_occ = get_maxn_occ(statepp%occ) 283 statepp%maxl_calc = state%maxl_calc 284 statepp%maxn_calc = 0 285 maxn = get_maxn_occ(statepp%core) 286 DO k = 0, statepp%maxl_calc 287 statepp%maxn_calc(k) = state%maxn_calc(k) - maxn(k) 288 statepp%maxn_calc(k) = MIN(statepp%maxn_calc(k), pp_basis%nbas(k)) 289 END DO 290 statepp%multiplicity = state%multiplicity 291 292 CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im) 293 CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im) 294 CALL set_atom(atom_info(in, im)%atom, method_type=method) 295 CALL set_atom(atom_refs(in, im)%atom, method_type=method, relativistic=reltyp) 296 297 ! calculate integrals: pseudopotential basis 298 ! general integrals 299 CALL atom_int_setup(pp_int, pp_basis, potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e) 300 ! 301 NULLIFY (pp_int%tzora, pp_int%hdkh) 302 ! potential 303 CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot) 304 ! 305 CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot) 306 statepp%maxn_calc(:) = MIN(statepp%maxn_calc(:), pp_basis%nbas(:)) 307 CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ)) 308 309 ! calculate integrals: all electron basis 310 ! general integrals 311 CALL atom_int_setup(ae_int, ae_basis, potential=ae_pot, & 312 eri_coulomb=eri_c, eri_exchange=eri_e) 313 ! potential 314 CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot) 315 ! relativistic correction terms 316 CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp)) 317 ! 318 CALL set_atom(atom_refs(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot) 319 state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:)) 320 CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ)) 321 322 CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, & 323 exchange_integral_type=do_erie) 324 CALL set_atom(atom_refs(in, im)%atom, coulomb_integral_type=do_eric, & 325 exchange_integral_type=do_erie) 326 atom_info(in, im)%atom%hfx_pot%do_gh = do_gh 327 atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh 328 atom_refs(in, im)%atom%hfx_pot%do_gh = do_gh 329 atom_refs(in, im)%atom%hfx_pot%nr_gh = nr_gh 330 331 CALL set_atom(atom_info(in, im)%atom, state=statepp) 332 NULLIFY (orbitals) 333 mo = MAXVAL(statepp%maxn_calc) 334 mb = MAXVAL(atom_info(in, im)%atom%basis%nbas) 335 CALL create_atom_orbs(orbitals, mb, mo) 336 CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals) 337 338 CALL set_atom(atom_refs(in, im)%atom, state=state) 339 NULLIFY (orbitals) 340 mo = MAXVAL(state%maxn_calc) 341 mb = MAXVAL(atom_refs(in, im)%atom%basis%nbas) 342 CALL create_atom_orbs(orbitals, mb, mo) 343 CALL set_atom(atom_refs(in, im)%atom, orbitals=orbitals) 344 345 IF (atom_consistent_method(atom_refs(in, im)%atom%method_type, atom_refs(in, im)%atom%state%multiplicity)) THEN 346 !Print method info 347 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log") 348 CALL atom_print_method(atom_refs(in, im)%atom, iw) 349 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO") 350 !Calculate the electronic structure 351 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log") 352 CALL calculate_atom(atom_refs(in, im)%atom, iw) 353 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO") 354 END IF 355 END DO 356 END DO 357 358 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_PSEUDO", extension=".log") 359 IF (iw > 0) THEN 360 WRITE (iw, '(/," ",79("*"))') 361 WRITE (iw, '(" ",21("*"),A,21("*"))') " Optimize Pseudopotential Parameters " 362 WRITE (iw, '(" ",79("*"))') 363 END IF 364 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_PSEUDO") 365 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log") 366 IF (iw > 0) THEN 367 CALL atom_print_potential(p_pot, iw) 368 END IF 369 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL") 370 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_PSEUDO", extension=".log") 371 IF (iw > 0) THEN 372 powell_section => section_vals_get_subs_vals(atom_section, "POWELL") 373 CALL atom_fit_pseudo(atom_info, atom_refs, p_pot, iw, powell_section) 374 END IF 375 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_PSEUDO") 376 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log") 377 IF (iw > 0) THEN 378 CALL atom_print_potential(p_pot, iw) 379 END IF 380 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL") 381 382 ! clean up 383 CALL atom_int_release(ae_int) 384 CALL atom_ppint_release(ae_int) 385 CALL atom_relint_release(ae_int) 386 387 CALL atom_int_release(pp_int) 388 CALL atom_ppint_release(pp_int) 389 CALL atom_relint_release(pp_int) 390 391 CALL release_atom_basis(ae_basis) 392 CALL release_atom_basis(pp_basis) 393 394 CALL release_atom_potential(p_pot) 395 CALL release_atom_potential(ae_pot) 396 397 DO in = 1, n_rep 398 DO im = 1, n_meth 399 CALL release_atom_type(atom_info(in, im)%atom) 400 CALL release_atom_type(atom_refs(in, im)%atom) 401 END DO 402 END DO 403 DEALLOCATE (atom_info, atom_refs) 404 405 DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int) 406 407 CALL timestop(handle) 408 409 END SUBROUTINE atom_pseudo_opt 410 411! ************************************************************************************************** 412 413END MODULE atom_pseudo 414