(git:936074a)
Loading...
Searching...
No Matches
atom_electronic_structure.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8MODULE atom_electronic_structure
9 USE atom_optimization, ONLY: atom_history_init,&
10 atom_history_release,&
11 atom_history_type,&
12 atom_history_update,&
13 atom_opt_fmat
14 USE atom_output, ONLY: atom_print_energies,&
15 atom_print_iteration,&
16 atom_print_state,&
17 atom_print_zmp_iteration
18 USE atom_types, ONLY: &
19 atom_type, create_opgrid, create_opmat, ecp_pseudo, gth_pseudo, lmat, no_pseudo, &
20 opgrid_type, opmat_type, release_opgrid, release_opmat, set_atom, setup_hf_section, &
21 sgp_pseudo, upf_pseudo
22 USE atom_utils, ONLY: &
23 atom_denmat, atom_density, atom_read_external_density, atom_read_external_vxc, &
24 atom_read_zmp_restart, atom_solve, atom_trace, ceri_contract, coulomb_potential_analytic, &
25 coulomb_potential_numeric, eeri_contract, err_matrix, exchange_numeric, &
26 exchange_semi_analytic, numpot_matrix, slater_density, wigner_slater_functional
27 USE atom_xc, ONLY: calculate_atom_ext_vxc,&
28 calculate_atom_vxc_lda,&
29 calculate_atom_vxc_lsd,&
30 calculate_atom_zmp
31 USE input_constants, ONLY: &
32 do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
33 do_numeric, do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_semi_analytic, &
34 do_uhf_atom, do_uks_atom, do_zoramp_atom
35 USE input_section_types, ONLY: section_vals_get,&
36 section_vals_get_subs_vals,&
37 section_vals_type
38 USE kinds, ONLY: dp
39#include "./base/base_uses.f90"
40
41 IMPLICIT NONE
42 PRIVATE
43 PUBLIC :: calculate_atom
44
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_electronic_structure'
46
47CONTAINS
48
49! **************************************************************************************************
50!> \brief General routine to perform electronic structure atomic calculations.
51!> \param atom information about the atomic kind
52!> \param iw output file unit
53!> \param noguess skip initial guess
54!> \param converged whether SCF iterations have been converged
55!> \par History
56!> * 06.2017 disable XC input options [Juerg Hutter]
57!> * 11.2009 created from the subroutine calculate_atom() [Juerg Hutter]
58! **************************************************************************************************
59 SUBROUTINE calculate_atom(atom, iw, noguess, converged)
60 TYPE(atom_type), POINTER :: atom
61 INTEGER, INTENT(IN) :: iw
62 LOGICAL, INTENT(IN), OPTIONAL :: noguess
63 LOGICAL, INTENT(OUT), OPTIONAL :: converged
64
65 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom'
66
67 INTEGER :: handle, method
68 LOGICAL :: explicit
69 TYPE(section_vals_type), POINTER :: sub_section
70
71 CALL timeset(routinen, handle)
72
73 ! test for not supported methods and functionals
74 IF (ASSOCIATED(atom%xc_section)) THEN
75 !
76 sub_section => section_vals_get_subs_vals(atom%xc_section, "ADIABATIC_RESCALING")
77 CALL section_vals_get(sub_section, explicit=explicit)
78 IF (explicit) CALL cp_abort(__location__, "ADIABATIC_RESCALING not supported in ATOM code")
79 !
80 sub_section => section_vals_get_subs_vals(atom%xc_section, "VDW_POTENTIAL")
81 CALL section_vals_get(sub_section, explicit=explicit)
82 IF (explicit) CALL cp_abort(__location__, "VDW_POTENTIAL not supported in ATOM code")
83 !
84 sub_section => section_vals_get_subs_vals(atom%xc_section, "XC_POTENTIAL")
85 CALL section_vals_get(sub_section, explicit=explicit)
86 IF (explicit) CALL cp_abort(__location__, "XC_POTENTIAL not supported in ATOM code")
87 !
88 sub_section => section_vals_get_subs_vals(atom%xc_section, "WF_CORRELATION")
89 CALL section_vals_get(sub_section, explicit=explicit)
90 IF (explicit) CALL cp_abort(__location__, "WF_CORRELATION methods not supported in ATOM code")
91 !
92 END IF
93
94 method = atom%method_type
95 SELECT CASE (method)
96 CASE (do_rks_atom, do_rhf_atom)
97 CALL calculate_atom_restricted(atom, iw, noguess, converged)
98 CASE (do_uks_atom, do_uhf_atom)
99 CALL calculate_atom_unrestricted(atom, iw, noguess, converged)
100 CASE (do_rohf_atom)
101 cpabort("The ROHF method is not yet implemented in the ATOM code.")
102 CASE DEFAULT
103 cpabort("Invalid method specified for ATOM code. Check the code!")
104 END SELECT
105
106 CALL timestop(handle)
107
108 END SUBROUTINE calculate_atom
109
110! **************************************************************************************************
111!> \brief Perform restricted (closed shell) electronic structure atomic calculations.
112!> \param atom information about the atomic kind
113!> \param iw output file unit
114!> \param noguess skip initial guess
115!> \param converged whether SCF iterations have been converged
116!> \par History
117!> * 02.2016 support UPF files and ECP potentials [Juerg Hutter]
118!> * 09.2015 Polarized Atomic Orbitals (PAO) method [Ole Schuett]
119!> * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
120!> * 07.2013 scaled ZORA with model potential [Juerg Hutter]
121!> * 11.2009 split into three subroutines calculate_atom(), calculate_atom_restricted(),
122!> and calculate_atom_unrestricted() [Juerg Hutter]
123!> * 12.2008 refactored and renamed to calculate_atom() [Juerg Hutter]
124!> * 08.2008 created as atom_electronic_structure() [Juerg Hutter]
125! **************************************************************************************************
126 SUBROUTINE calculate_atom_restricted(atom, iw, noguess, converged)
127 TYPE(atom_type), POINTER :: atom
128 INTEGER, INTENT(IN) :: iw
129 LOGICAL, INTENT(IN), OPTIONAL :: noguess
130 LOGICAL, INTENT(OUT), OPTIONAL :: converged
131
132 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom_restricted'
133
134 INTEGER :: handle, i, iter, l, max_iter, method, &
135 reltyp
136 LOGICAL :: do_hfx, doguess, is_converged, need_vxc, &
137 need_x, need_xc, need_zmp
138 REAL(kind=dp) :: deps, eps_scf, hf_frac
139 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ext_density, ext_vxc, tmp_dens
140 TYPE(atom_history_type) :: history
141 TYPE(opgrid_type), POINTER :: cpot, density
142 TYPE(opmat_type), POINTER :: fmat, hcore, jmat, kmat, xcmat
143 TYPE(section_vals_type), POINTER :: xc_section
144
145 CALL timeset(routinen, handle)
146
147 IF (PRESENT(noguess)) THEN
148 doguess = .NOT. noguess
149 ELSE
150 doguess = .true.
151 END IF
152
153 CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
154
155 method = atom%method_type
156 max_iter = atom%optimization%max_iter
157 eps_scf = atom%optimization%eps_scf
158
159 SELECT CASE (method)
160 CASE DEFAULT
161 cpabort("")
162 CASE (do_rks_atom)
163 need_x = do_hfx
164 need_xc = .true.
165 CASE (do_uks_atom)
166 cpabort("")
167 CASE (do_rhf_atom)
168 need_x = .true.
169 need_xc = .false.
170 hf_frac = 1._dp
171 CASE (do_uhf_atom)
172 cpabort("")
173 CASE (do_rohf_atom)
174 need_x = .true.
175 need_xc = .false.
176 hf_frac = 1._dp
177 cpabort("")
178 END SELECT
179
180 ! ZMP starting to read external density for the zmp calculation
181 need_zmp = atom%do_zmp
182 IF (need_zmp) THEN
183 need_x = .false.
184 need_xc = .false.
185 ALLOCATE (ext_density(atom%basis%grid%nr))
186 ext_density = 0._dp
187 CALL atom_read_external_density(ext_density, atom, iw)
188 END IF
189
190 ! ZMP starting to read external vxc for the zmp calculation
191 need_vxc = atom%read_vxc
192
193 IF (need_vxc) THEN
194 need_x = .false.
195 need_xc = .false.
196 need_zmp = .false.
197 ALLOCATE (ext_vxc(atom%basis%grid%nr))
198 ext_vxc = 0._dp
199 CALL atom_read_external_vxc(ext_vxc, atom, iw)
200 END IF
201
202 ! check for relativistic method
203 reltyp = atom%relativistic
204
205 IF (iw > 0) CALL atom_print_state(atom%state, iw)
206
207 NULLIFY (hcore)
208 CALL create_opmat(hcore, atom%basis%nbas)
209 ! Pseudopotentials
210 SELECT CASE (atom%potential%ppot_type)
211 CASE DEFAULT
212 cpabort("")
213 CASE (no_pseudo)
214 SELECT CASE (reltyp)
215 CASE DEFAULT
216 cpabort("")
217 CASE (do_nonrel_atom)
218 hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
219 CASE (do_zoramp_atom, do_sczoramp_atom)
220 hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
221 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
222 hcore%op = atom%integrals%hdkh
223 END SELECT
224 CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
225 hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
226 END SELECT
227 ! add confinement potential (not included in relativistic transformations)
228 IF (atom%potential%confinement) THEN
229 hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
230 END IF
231
232 NULLIFY (fmat, jmat, kmat, xcmat)
233 CALL create_opmat(fmat, atom%basis%nbas)
234 CALL create_opmat(jmat, atom%basis%nbas)
235 CALL create_opmat(kmat, atom%basis%nbas)
236 CALL create_opmat(xcmat, atom%basis%nbas)
237
238 NULLIFY (density, cpot)
239 CALL create_opgrid(density, atom%basis%grid)
240 CALL create_opgrid(cpot, atom%basis%grid)
241
242 ! ZMP reading the file to restart
243 IF (atom%doread) CALL atom_read_zmp_restart(atom, doguess, iw)
244
245 IF (doguess) THEN
246 ! initial guess
247 ALLOCATE (tmp_dens(SIZE(density%op)))
248 CALL slater_density(density%op, tmp_dens, atom%z, atom%state, atom%basis%grid)
249 density%op = density%op + tmp_dens
250 DEALLOCATE (tmp_dens)
251 CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
252 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
253 CALL wigner_slater_functional(density%op, cpot%op)
254 CALL numpot_matrix(xcmat%op, cpot%op, atom%basis, 0)
255 fmat%op = hcore%op + jmat%op + xcmat%op
256 CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
257 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
258 END IF
259 CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
260 atom%state%maxl_occ, atom%state%maxn_occ)
261
262 ! wavefunction history
263 NULLIFY (history%dmat, history%hmat)
264 CALL atom_history_init(history, atom%optimization, fmat%op)
265 is_converged = .false.
266 iter = 0
267 DO !SCF Loop
268
269 ! Kinetic energy
270 atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
271
272 ! Band energy
273 atom%energy%eband = 0._dp
274 DO l = 0, lmat
275 DO i = 1, min(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
276 atom%energy%eband = atom%energy%eband + atom%orbitals%ener(i, l)*atom%state%occupation(l, i)
277 END DO
278 END DO
279
280 ! Pseudopotential energy
281 SELECT CASE (atom%potential%ppot_type)
282 CASE DEFAULT
283 cpabort("")
284 CASE (no_pseudo)
285 atom%energy%eploc = 0._dp
286 atom%energy%epnl = 0._dp
287 CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
288 atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
289 atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
290 END SELECT
291 atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
292
293 ! Core energy
294 atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
295
296 ! Confinement energy
297 IF (atom%potential%confinement) THEN
298 atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
299 ELSE
300 atom%energy%econfinement = 0._dp
301 END IF
302
303 ! Hartree Term
304 jmat%op = 0._dp
305 SELECT CASE (atom%coulomb_integral_type)
306 CASE DEFAULT
307 cpabort("")
308 CASE (do_analytic)
309 CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
310 CASE (do_semi_analytic)
311 CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
312 atom%state%maxl_occ)
313 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
314 CASE (do_numeric)
315 CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
316 CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
317 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
318 END SELECT
319 atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
320
321 ! Exchange Term
322 IF (need_x) THEN
323 kmat%op = 0._dp
324 SELECT CASE (atom%exchange_integral_type)
325 CASE DEFAULT
326 cpabort("")
327 CASE (do_analytic)
328 CALL eeri_contract(kmat%op, atom%integrals%eeri, atom%orbitals%pmat, atom%integrals%n)
329 CASE (do_semi_analytic)
330 CALL exchange_semi_analytic(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
331 CASE (do_numeric)
332 CALL exchange_numeric(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
333 END SELECT
334 atom%energy%eexchange = hf_frac*0.5_dp*atom_trace(kmat%op, atom%orbitals%pmat)
335 kmat%op = hf_frac*kmat%op
336 ELSE
337 kmat%op = 0._dp
338 atom%energy%eexchange = 0._dp
339 END IF
340
341 ! XC
342 IF (need_xc) THEN
343 xcmat%op = 0._dp
344 CALL calculate_atom_vxc_lda(xcmat, atom, xc_section)
345 ! ZMP added options for the zmp calculations, building external density and vxc potential
346 ELSEIF (need_zmp) THEN
347 xcmat%op = 0._dp
348 CALL calculate_atom_zmp(ext_density=ext_density, atom=atom, lprint=.false., xcmat=xcmat)
349 ELSEIF (need_vxc) THEN
350 xcmat%op = 0._dp
351 CALL calculate_atom_ext_vxc(vxc=ext_vxc, atom=atom, lprint=.false., xcmat=xcmat)
352 ELSE
353 xcmat%op = 0._dp
354 atom%energy%exc = 0._dp
355 END IF
356
357 ! Zero this contribution
358 atom%energy%elsd = 0._dp
359
360 ! Total energy
361 atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
362
363 ! Potential energy
364 atom%energy%epot = atom%energy%etot - atom%energy%ekin
365
366 ! Total HF/KS matrix
367 fmat%op = hcore%op + jmat%op + kmat%op + xcmat%op
368
369 ! calculate error matrix
370 CALL err_matrix(jmat%op, deps, fmat%op, atom%orbitals%pmat, atom%integrals%utrans, &
371 atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
372
373 iter = iter + 1
374
375 IF (iw > 0) THEN
376 IF (need_zmp) THEN
377 CALL atom_print_zmp_iteration(iter, deps, atom, iw)
378 ELSE
379 CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
380 END IF
381 END IF
382
383 IF (deps < eps_scf) THEN
384 is_converged = .true.
385 EXIT
386 END IF
387 IF (iter >= max_iter) THEN
388 IF (iw > 0) THEN
389 WRITE (iw, "(A)") " No convergence within maximum number of iterations "
390 END IF
391 EXIT
392 END IF
393
394 ! update history container and extrapolate KS matrix
395 CALL atom_history_update(history, atom%orbitals%pmat, fmat%op, jmat%op, atom%energy%etot, deps)
396 CALL atom_opt_fmat(fmat%op, history, deps)
397
398 ! Solve HF/KS equations
399 CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
400 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
401 CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
402 atom%state%maxl_occ, atom%state%maxn_occ)
403
404 END DO !SCF Loop
405
406 IF (iw > 0) THEN
407 CALL atom_print_energies(atom, iw)
408 END IF
409
410 CALL atom_history_release(history)
411
412 ! conserve fmat within atom_type
413 CALL set_atom(atom, fmat=fmat)
414 CALL release_opmat(jmat)
415 CALL release_opmat(kmat)
416 CALL release_opmat(xcmat)
417 CALL release_opmat(hcore)
418
419 CALL release_opgrid(density)
420 CALL release_opgrid(cpot)
421
422 ! ZMP deallocating ext_density ext_vxc
423 IF (need_zmp) DEALLOCATE (ext_density)
424 IF (need_vxc) DEALLOCATE (ext_vxc)
425
426 IF (PRESENT(converged)) THEN
427 converged = is_converged
428 END IF
429
430 CALL timestop(handle)
431
432 END SUBROUTINE calculate_atom_restricted
433
434! **************************************************************************************************
435!> \brief Perform unrestricted (spin-polarised) electronic structure atomic calculations.
436!> \param atom information about the atomic kind
437!> \param iw output file unit
438!> \param noguess skip initial guess
439!> \param converged whether SCF iterations have been converged
440!> \par History
441!> * identical to the subroutine calculate_atom_restricted()
442! **************************************************************************************************
443 SUBROUTINE calculate_atom_unrestricted(atom, iw, noguess, converged)
444 TYPE(atom_type), POINTER :: atom
445 INTEGER, INTENT(IN) :: iw
446 LOGICAL, INTENT(IN), OPTIONAL :: noguess
447 LOGICAL, INTENT(OUT), OPTIONAL :: converged
448
449 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom_unrestricted'
450
451 INTEGER :: handle, i, iter, k, l, max_iter, method, &
452 reltyp
453 LOGICAL :: do_hfx, doguess, is_converged, lsdpot, &
454 need_x, need_xc
455 REAL(kind=dp) :: deps, depsa, depsb, eps_scf, hf_frac, &
456 ne, nm
457 TYPE(atom_history_type) :: historya, historyb
458 TYPE(opgrid_type), POINTER :: cpot, density, rhoa, rhob
459 TYPE(opmat_type), POINTER :: fmata, fmatb, hcore, hlsd, jmat, kmata, &
460 kmatb, xcmata, xcmatb
461 TYPE(section_vals_type), POINTER :: xc_section
462
463 CALL timeset(routinen, handle)
464
465 IF (PRESENT(noguess)) THEN
466 doguess = .NOT. noguess
467 ELSE
468 doguess = .true.
469 END IF
470
471 CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
472
473 method = atom%method_type
474 max_iter = atom%optimization%max_iter
475 eps_scf = atom%optimization%eps_scf
476
477 SELECT CASE (method)
478 CASE DEFAULT
479 cpabort("")
480 CASE (do_rks_atom)
481 cpabort("")
482 CASE (do_uks_atom)
483 need_x = do_hfx
484 need_xc = .true.
485 CASE (do_rhf_atom)
486 cpabort("")
487 CASE (do_uhf_atom)
488 need_x = .true.
489 need_xc = .false.
490 hf_frac = 1._dp
491 CASE (do_rohf_atom)
492 need_x = .true.
493 need_xc = .false.
494 hf_frac = 1._dp
495 cpabort("")
496 END SELECT
497
498 ! set alpha and beta occupations
499 IF (sum(abs(atom%state%occa) + abs(atom%state%occb)) == 0.0_dp) THEN
500 DO l = 0, 3
501 nm = real((2*l + 1), kind=dp)
502 DO k = 1, 10
503 ne = atom%state%occupation(l, k)
504 IF (ne == 0._dp) THEN !empty shell
505 EXIT !assume there are no holes
506 ELSEIF (ne == 2._dp*nm) THEN !closed shell
507 atom%state%occa(l, k) = nm
508 atom%state%occb(l, k) = nm
509 ELSEIF (atom%state%multiplicity == -2) THEN !High spin case
510 atom%state%occa(l, k) = min(ne, nm)
511 atom%state%occb(l, k) = max(0._dp, ne - nm)
512 ELSE
513 atom%state%occa(l, k) = 0.5_dp*(ne + atom%state%multiplicity - 1._dp)
514 atom%state%occb(l, k) = ne - atom%state%occa(l, k)
515 END IF
516 END DO
517 END DO
518 END IF
519 ! check for relativistic method
520 reltyp = atom%relativistic
521
522 IF (iw > 0) CALL atom_print_state(atom%state, iw)
523
524 NULLIFY (hcore, hlsd)
525 CALL create_opmat(hcore, atom%basis%nbas)
526 CALL create_opmat(hlsd, atom%basis%nbas)
527 hlsd%op = 0._dp
528 ! Pseudopotentials
529 lsdpot = .false.
530 SELECT CASE (atom%potential%ppot_type)
531 CASE DEFAULT
532 cpabort("")
533 CASE (no_pseudo)
534 SELECT CASE (reltyp)
535 CASE DEFAULT
536 cpabort("")
537 CASE (do_nonrel_atom)
538 hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
539 CASE (do_zoramp_atom, do_sczoramp_atom)
540 hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
541 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
542 hcore%op = atom%integrals%hdkh
543 END SELECT
544 CASE (gth_pseudo)
545 hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
546 IF (atom%potential%gth_pot%lsdpot) THEN
547 lsdpot = .true.
548 hlsd%op = atom%integrals%clsd
549 END IF
550 CASE (upf_pseudo, sgp_pseudo, ecp_pseudo)
551 hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
552 END SELECT
553 ! add confinement potential (not included in relativistic transformations)
554 IF (atom%potential%confinement) THEN
555 hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
556 END IF
557
558 NULLIFY (fmata, fmatb, jmat, kmata, kmatb, xcmata, xcmatb)
559 CALL create_opmat(fmata, atom%basis%nbas)
560 CALL create_opmat(fmatb, atom%basis%nbas)
561 CALL create_opmat(jmat, atom%basis%nbas)
562 CALL create_opmat(kmata, atom%basis%nbas)
563 CALL create_opmat(kmatb, atom%basis%nbas)
564 CALL create_opmat(xcmata, atom%basis%nbas)
565 CALL create_opmat(xcmatb, atom%basis%nbas)
566
567 NULLIFY (density, rhoa, rhob, cpot)
568 CALL create_opgrid(density, atom%basis%grid)
569 CALL create_opgrid(rhoa, atom%basis%grid)
570 CALL create_opgrid(rhob, atom%basis%grid)
571 CALL create_opgrid(cpot, atom%basis%grid)
572
573 IF (doguess) THEN
574 ! initial guess
575 CALL slater_density(rhoa%op, rhob%op, atom%z, atom%state, atom%basis%grid)
576 density%op = rhoa%op + rhob%op
577 CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
578 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
579 ! alpha spin
580 density%op = 2._dp*rhoa%op
581 CALL wigner_slater_functional(density%op, cpot%op)
582 CALL numpot_matrix(xcmata%op, cpot%op, atom%basis, 0)
583 fmata%op = hcore%op + hlsd%op + jmat%op + xcmata%op
584 CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
585 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
586 ! beta spin
587 density%op = 2._dp*rhob%op
588 CALL wigner_slater_functional(density%op, cpot%op)
589 CALL numpot_matrix(xcmatb%op, cpot%op, atom%basis, 0)
590 fmatb%op = hcore%op - hlsd%op + jmat%op + xcmatb%op
591 CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
592 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
593 END IF
594 CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
595 atom%state%maxl_occ, atom%state%maxn_occ)
596 CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
597 atom%state%maxl_occ, atom%state%maxn_occ)
598 atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
599
600 ! wavefunction history
601 NULLIFY (historya%dmat, historya%hmat)
602 CALL atom_history_init(historya, atom%optimization, fmata%op)
603 NULLIFY (historyb%dmat, historyb%hmat)
604 CALL atom_history_init(historyb, atom%optimization, fmatb%op)
605
606 is_converged = .false.
607 iter = 0
608 DO !SCF Loop
609
610 ! Kinetic energy
611 atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
612
613 ! Band energy
614 atom%energy%eband = 0._dp
615 DO l = 0, 3
616 DO i = 1, min(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
617 atom%energy%eband = atom%energy%eband + atom%orbitals%enera(i, l)*atom%state%occa(l, i)
618 atom%energy%eband = atom%energy%eband + atom%orbitals%enerb(i, l)*atom%state%occb(l, i)
619 END DO
620 END DO
621
622 ! Pseudopotential energy
623 SELECT CASE (atom%potential%ppot_type)
624 CASE DEFAULT
625 cpabort("")
626 CASE (no_pseudo)
627 atom%energy%eploc = 0._dp
628 atom%energy%epnl = 0._dp
629 CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
630 atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
631 atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
632 END SELECT
633 atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
634
635 ! Core energy
636 atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
637
638 ! Confinement energy
639 IF (atom%potential%confinement) THEN
640 atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
641 ELSE
642 atom%energy%econfinement = 0._dp
643 END IF
644
645 ! Hartree Term
646 jmat%op = 0._dp
647 SELECT CASE (atom%coulomb_integral_type)
648 CASE DEFAULT
649 cpabort("")
650 CASE (do_analytic)
651 CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
652 CASE (do_semi_analytic)
653 CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
654 atom%state%maxl_occ)
655 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
656 CASE (do_numeric)
657 CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
658 CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
659 CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
660 END SELECT
661 atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
662
663 ! Exchange Term
664 IF (need_x) THEN
665 kmata%op = 0._dp
666 kmatb%op = 0._dp
667 SELECT CASE (atom%exchange_integral_type)
668 CASE DEFAULT
669 cpabort("")
670 CASE (do_analytic)
671 CALL eeri_contract(kmata%op, atom%integrals%eeri, atom%orbitals%pmata, atom%integrals%n)
672 CALL eeri_contract(kmatb%op, atom%integrals%eeri, atom%orbitals%pmatb, atom%integrals%n)
673 CASE (do_semi_analytic)
674 CALL exchange_semi_analytic(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
675 CALL exchange_semi_analytic(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
676 CASE (do_numeric)
677 CALL exchange_numeric(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
678 CALL exchange_numeric(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
679 END SELECT
680 atom%energy%eexchange = hf_frac*(atom_trace(kmata%op, atom%orbitals%pmata) + &
681 atom_trace(kmatb%op, atom%orbitals%pmatb))
682 kmata%op = 2._dp*hf_frac*kmata%op
683 kmatb%op = 2._dp*hf_frac*kmatb%op
684 ELSE
685 kmata%op = 0._dp
686 kmatb%op = 0._dp
687 atom%energy%eexchange = 0._dp
688 END IF
689
690 ! XC
691 IF (need_xc) THEN
692 xcmata%op = 0._dp
693 xcmatb%op = 0._dp
694 CALL calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
695 ELSE
696 xcmata%op = 0._dp
697 xcmatb%op = 0._dp
698 atom%energy%exc = 0._dp
699 END IF
700
701 IF (lsdpot) THEN
702 atom%energy%elsd = atom_trace(hlsd%op, atom%orbitals%pmata) - &
703 atom_trace(hlsd%op, atom%orbitals%pmatb)
704 atom%energy%epseudo = atom%energy%epseudo + atom%energy%elsd
705 atom%energy%ecore = atom%energy%ecore + atom%energy%elsd
706 ELSE
707 atom%energy%elsd = 0._dp
708 END IF
709
710 ! Total energy
711 atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
712
713 ! Potential energy
714 atom%energy%epot = atom%energy%etot - atom%energy%ekin
715
716 ! Total HF/KS matrix
717 fmata%op = hcore%op + hlsd%op + jmat%op + kmata%op + xcmata%op
718 fmatb%op = hcore%op - hlsd%op + jmat%op + kmatb%op + xcmatb%op
719
720 ! calculate error matrix
721 CALL err_matrix(xcmata%op, depsa, fmata%op, atom%orbitals%pmata, atom%integrals%utrans, &
722 atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
723 CALL err_matrix(xcmatb%op, depsb, fmatb%op, atom%orbitals%pmatb, atom%integrals%utrans, &
724 atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
725 deps = 2._dp*max(depsa, depsb)
726
727 iter = iter + 1
728
729 IF (iw > 0) THEN
730 CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
731 END IF
732
733 IF (deps < eps_scf) THEN
734 is_converged = .true.
735 EXIT
736 END IF
737 IF (iter >= max_iter) THEN
738 IF (iw > 0) THEN
739 WRITE (iw, "(A)") " No convergence within maximum number of iterations "
740 END IF
741 EXIT
742 END IF
743
744 ! update history container and extrapolate KS matrix
745 CALL atom_history_update(historya, atom%orbitals%pmata, fmata%op, xcmata%op, atom%energy%etot, deps)
746 CALL atom_history_update(historyb, atom%orbitals%pmatb, fmatb%op, xcmatb%op, atom%energy%etot, deps)
747 CALL atom_opt_fmat(fmata%op, historya, depsa)
748 CALL atom_opt_fmat(fmatb%op, historyb, depsb)
749
750 ! Solve HF/KS equations
751 CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
752 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
753 CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
754 atom%state%maxl_occ, atom%state%maxn_occ)
755 CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
756 atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
757 CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
758 atom%state%maxl_occ, atom%state%maxn_occ)
759 atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
760
761 END DO !SCF Loop
762
763 IF (iw > 0) THEN
764 CALL atom_print_energies(atom, iw)
765 END IF
766
767 CALL atom_history_release(historya)
768 CALL atom_history_release(historyb)
769
770 CALL release_opgrid(density)
771 CALL release_opgrid(rhoa)
772 CALL release_opgrid(rhob)
773 CALL release_opgrid(cpot)
774
775 ! conserve fmata as fmat within atom_type.
776 CALL set_atom(atom, fmat=fmata)
777 CALL release_opmat(fmatb)
778 CALL release_opmat(jmat)
779 CALL release_opmat(kmata)
780 CALL release_opmat(kmatb)
781 CALL release_opmat(xcmata)
782 CALL release_opmat(xcmatb)
783 CALL release_opmat(hlsd)
784 CALL release_opmat(hcore)
785
786 IF (PRESENT(converged)) THEN
787 converged = is_converged
788 END IF
789
790 CALL timestop(handle)
791
792 END SUBROUTINE calculate_atom_unrestricted
793
794END MODULE atom_electronic_structure
atom_electronic_structure
Definition atom_electronic_structure.F:8
atom_electronic_structure::calculate_atom
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
Definition atom_electronic_structure.F:60
atom_optimization
Optimizer for the atomic code.
Definition atom_optimization.F:11
atom_optimization::atom_opt_fmat
subroutine, public atom_opt_fmat(fmat, history, err)
Construct a Kohn-Sham matrix for the next iteration based on the historic data.
Definition atom_optimization.F:171
atom_optimization::atom_history_init
pure subroutine, public atom_history_init(history, optimization, matrix)
Initialise a circular buffer to keep Kohn-Sham and density matrices from previous iteration.
Definition atom_optimization.F:56
atom_optimization::atom_history_release
pure subroutine, public atom_history_release(history)
Release circular buffer to keep historic matrices.
Definition atom_optimization.F:132
atom_optimization::atom_history_update
pure subroutine, public atom_history_update(history, pmat, fmat, emat, energy, error)
Add matrices from the current iteration into the circular buffer.
Definition atom_optimization.F:103
atom_output
Routines that print various information about an atomic kind.
Definition atom_output.F:11
atom_output::atom_print_energies
subroutine, public atom_print_energies(atom, iw)
Print energy components.
Definition atom_output.F:180
atom_output::atom_print_iteration
subroutine, public atom_print_iteration(iter, deps, etot, iw)
Print convergence information.
Definition atom_output.F:277
atom_output::atom_print_state
subroutine, public atom_print_state(state, iw)
Print information about electronic state.
Definition atom_output.F:88
atom_output::atom_print_zmp_iteration
subroutine, public atom_print_zmp_iteration(iter, deps, atom, iw)
Printing of the atomic iterations when ZMP is active.
Definition atom_output.F:253
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::release_opmat
subroutine, public release_opmat(opmat)
...
Definition atom_types.F:1382
atom_types::release_opgrid
subroutine, public release_opgrid(opgrid)
...
Definition atom_types.F:1422
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:67
atom_types::set_atom
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
Definition atom_types.F:1023
atom_types::create_opmat
subroutine, public create_opmat(opmat, n, lmax)
...
Definition atom_types.F:1354
atom_types::setup_hf_section
subroutine, public setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
...
Definition atom_types.F:1184
atom_types::create_opgrid
subroutine, public create_opgrid(opgrid, grid)
...
Definition atom_types.F:1399
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::slater_density
subroutine, public slater_density(density1, density2, zcore, state, grid)
Calculate Slater density on a radial grid.
Definition atom_utils.F:2093
atom_utils::atom_denmat
pure subroutine, public atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
Calculate a density matrix using atomic orbitals.
Definition atom_utils.F:327
atom_utils::atom_read_external_vxc
subroutine, public atom_read_external_vxc(vxc, atom, iw)
ZMP subroutine to read external v_xc in the atomic code.
Definition atom_utils.F:607
atom_utils::eeri_contract
pure subroutine, public eeri_contract(kmat, erint, pmat, nsize)
Contract exchange Electron Repulsion Integrals.
Definition atom_utils.F:1973
atom_utils::exchange_semi_analytic
subroutine, public exchange_semi_analytic(kmat, state, occ, wfn, basis, hfx_pot)
Calculate the exchange potential semi-analytically.
Definition atom_utils.F:1307
atom_utils::ceri_contract
pure subroutine, public ceri_contract(jmat, erint, pmat, nsize, all_nu)
Contract Coulomb Electron Repulsion Integrals.
Definition atom_utils.F:1899
atom_utils::coulomb_potential_analytic
subroutine, public coulomb_potential_analytic(cpot, pmat, basis, grid, maxl)
Analytically compute the Coulomb potential on an atomic radial grid.
Definition atom_utils.F:1114
atom_utils::err_matrix
pure subroutine, public err_matrix(emat, demax, kmat, pmat, umat, upmat, nval, nbs)
Calculate the error matrix for each angular momentum.
Definition atom_utils.F:2050
atom_utils::atom_density
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition atom_utils.F:366
atom_utils::wigner_slater_functional
pure subroutine, public wigner_slater_functional(rho, vxc)
Calculate the functional derivative of the Wigner (correlation) - Slater (exchange) density functiona...
Definition atom_utils.F:2158
atom_utils::atom_read_zmp_restart
subroutine, public atom_read_zmp_restart(atom, doguess, iw)
ZMP subroutine to read external restart file.
Definition atom_utils.F:457
atom_utils::exchange_numeric
subroutine, public exchange_numeric(kmat, state, occ, wfn, basis, hfx_pot)
Calculate the exchange potential numerically.
Definition atom_utils.F:1219
atom_utils::atom_trace
pure real(kind=dp) function, public atom_trace(opmat, pmat)
Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
Definition atom_utils.F:1822
atom_utils::atom_read_external_density
subroutine, public atom_read_external_density(density, atom, iw)
ZMP subroutine to read external density from linear grid of density matrix.
Definition atom_utils.F:516
atom_utils::numpot_matrix
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
Definition atom_utils.F:1839
atom_utils::atom_solve
subroutine, public atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
Solve the generalised eigenproblem for every angular momentum.
Definition atom_utils.F:943
atom_utils::coulomb_potential_numeric
subroutine, public coulomb_potential_numeric(cpot, density, grid)
Numerically compute the Coulomb potential on an atomic radial grid.
Definition atom_utils.F:1076
atom_xc
routines that build the integrals of the Vxc potential calculated for the atomic code
Definition atom_xc.F:12
atom_xc::calculate_atom_vxc_lda
subroutine, public calculate_atom_vxc_lda(xcmat, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-restricted case.
Definition atom_xc.F:199
atom_xc::calculate_atom_ext_vxc
subroutine, public calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
ZMP subroutine for the scf external density from the external v_xc.
Definition atom_xc.F:148
atom_xc::calculate_atom_zmp
subroutine, public calculate_atom_zmp(ext_density, atom, lprint, xcmat)
ZMP subroutine for the calculation of the effective constraint potential in the atomic code.
Definition atom_xc.F:78
atom_xc::calculate_atom_vxc_lsd
subroutine, public calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-polarized case.
Definition atom_xc.F:404
atom
Definition atom.F:9
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_rhf_atom
integer, parameter, public do_rhf_atom
Definition input_constants.F:159
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition input_constants.F:159
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:149
input_constants::sgp_pseudo
integer, parameter, public sgp_pseudo
Definition input_constants.F:174
input_constants::do_dkh3_atom
integer, parameter, public do_dkh3_atom
Definition input_constants.F:165
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
input_constants::do_nonrel_atom
integer, parameter, public do_nonrel_atom
Definition input_constants.F:165
input_constants::do_dkh0_atom
integer, parameter, public do_dkh0_atom
Definition input_constants.F:165
input_constants::do_uhf_atom
integer, parameter, public do_uhf_atom
Definition input_constants.F:159
input_constants::upf_pseudo
integer, parameter, public upf_pseudo
Definition input_constants.F:174
input_constants::do_dkh2_atom
integer, parameter, public do_dkh2_atom
Definition input_constants.F:165
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition input_constants.F:159
input_constants::do_numeric
integer, parameter, public do_numeric
Definition input_constants.F:149
input_constants::do_zoramp_atom
integer, parameter, public do_zoramp_atom
Definition input_constants.F:165
input_constants::do_dkh1_atom
integer, parameter, public do_dkh1_atom
Definition input_constants.F:165
input_constants::do_rohf_atom
integer, parameter, public do_rohf_atom
Definition input_constants.F:159
input_constants::do_semi_analytic
integer, parameter, public do_semi_analytic
Definition input_constants.F:149
input_constants::do_sczoramp_atom
integer, parameter, public do_sczoramp_atom
Definition input_constants.F:165
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
atom_optimization::atom_history_type
Definition atom_optimization.F:30
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:293
atom_types::opgrid_type
Operator grids.
Definition atom_types.F:258
atom_types::opmat_type
Operator matrices.
Definition atom_types.F:251
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127