Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE 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 :
47 : CONTAINS
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 11282 : 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 11282 : CALL timeset(routineN, handle)
72 :
73 : ! test for not supported methods and functionals
74 11282 : IF (ASSOCIATED(atom%xc_section)) THEN
75 : !
76 2904 : sub_section => section_vals_get_subs_vals(atom%xc_section, "ADIABATIC_RESCALING")
77 2904 : CALL section_vals_get(sub_section, explicit=explicit)
78 2904 : IF (explicit) CALL cp_abort(__LOCATION__, "ADIABATIC_RESCALING not supported in ATOM code")
79 : !
80 2904 : sub_section => section_vals_get_subs_vals(atom%xc_section, "VDW_POTENTIAL")
81 2904 : CALL section_vals_get(sub_section, explicit=explicit)
82 2904 : IF (explicit) CALL cp_abort(__LOCATION__, "VDW_POTENTIAL not supported in ATOM code")
83 : !
84 2904 : sub_section => section_vals_get_subs_vals(atom%xc_section, "XC_POTENTIAL")
85 2904 : CALL section_vals_get(sub_section, explicit=explicit)
86 2904 : IF (explicit) CALL cp_abort(__LOCATION__, "XC_POTENTIAL not supported in ATOM code")
87 : !
88 2904 : sub_section => section_vals_get_subs_vals(atom%xc_section, "WF_CORRELATION")
89 2904 : CALL section_vals_get(sub_section, explicit=explicit)
90 2904 : IF (explicit) CALL cp_abort(__LOCATION__, "WF_CORRELATION methods not supported in ATOM code")
91 : !
92 : END IF
93 :
94 11282 : method = atom%method_type
95 11165 : SELECT CASE (method)
96 : CASE (do_rks_atom, do_rhf_atom)
97 11282 : CALL calculate_atom_restricted(atom, iw, noguess, converged)
98 : CASE (do_uks_atom, do_uhf_atom)
99 117 : CALL calculate_atom_unrestricted(atom, iw, noguess, converged)
100 : CASE (do_rohf_atom)
101 0 : CPABORT("")
102 : CASE DEFAULT
103 11282 : CPABORT("")
104 : END SELECT
105 :
106 11282 : CALL timestop(handle)
107 :
108 11282 : 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 11165 : 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 11165 : 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 11165 : CALL timeset(routineN, handle)
146 :
147 11165 : IF (PRESENT(noguess)) THEN
148 329 : doguess = .NOT. noguess
149 : ELSE
150 10836 : doguess = .TRUE.
151 : END IF
152 :
153 11165 : CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
154 :
155 11165 : method = atom%method_type
156 11165 : max_iter = atom%optimization%max_iter
157 11165 : eps_scf = atom%optimization%eps_scf
158 :
159 0 : SELECT CASE (method)
160 : CASE DEFAULT
161 0 : CPABORT("")
162 : CASE (do_rks_atom)
163 10022 : need_x = do_hfx
164 10022 : need_xc = .TRUE.
165 : CASE (do_uks_atom)
166 0 : CPABORT("")
167 : CASE (do_rhf_atom)
168 1143 : need_x = .TRUE.
169 1143 : need_xc = .FALSE.
170 1143 : hf_frac = 1._dp
171 : CASE (do_uhf_atom)
172 0 : CPABORT("")
173 : CASE (do_rohf_atom)
174 0 : need_x = .TRUE.
175 0 : need_xc = .FALSE.
176 0 : hf_frac = 1._dp
177 11165 : CPABORT("")
178 : END SELECT
179 :
180 : ! ZMP starting to read external density for the zmp calculation
181 11165 : need_zmp = atom%do_zmp
182 11165 : IF (need_zmp) THEN
183 0 : need_x = .FALSE.
184 0 : need_xc = .FALSE.
185 0 : ALLOCATE (ext_density(atom%basis%grid%nr))
186 0 : ext_density = 0._dp
187 0 : 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 11165 : need_vxc = atom%read_vxc
192 :
193 11165 : IF (need_vxc) THEN
194 0 : need_x = .FALSE.
195 0 : need_xc = .FALSE.
196 0 : need_zmp = .FALSE.
197 0 : ALLOCATE (ext_vxc(atom%basis%grid%nr))
198 0 : ext_vxc = 0._dp
199 0 : CALL atom_read_external_vxc(ext_vxc, atom, iw)
200 : END IF
201 :
202 : ! check for relativistic method
203 11165 : reltyp = atom%relativistic
204 :
205 11165 : IF (iw > 0) CALL atom_print_state(atom%state, iw)
206 :
207 11165 : NULLIFY (hcore)
208 11165 : CALL create_opmat(hcore, atom%basis%nbas)
209 : ! Pseudopotentials
210 11165 : SELECT CASE (atom%potential%ppot_type)
211 : CASE DEFAULT
212 0 : CPABORT("")
213 : CASE (NO_PSEUDO)
214 8913 : SELECT CASE (reltyp)
215 : CASE DEFAULT
216 0 : CPABORT("")
217 : CASE (do_nonrel_atom)
218 2547930 : hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
219 : CASE (do_zoramp_atom, do_sczoramp_atom)
220 242752 : 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 641396 : hcore%op = atom%integrals%hdkh
223 : END SELECT
224 : CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
225 2775569 : hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
226 : END SELECT
227 : ! add confinement potential (not included in relativistic transformations)
228 11165 : IF (atom%potential%confinement) THEN
229 3053407 : hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
230 : END IF
231 :
232 11165 : NULLIFY (fmat, jmat, kmat, xcmat)
233 11165 : CALL create_opmat(fmat, atom%basis%nbas)
234 11165 : CALL create_opmat(jmat, atom%basis%nbas)
235 11165 : CALL create_opmat(kmat, atom%basis%nbas)
236 11165 : CALL create_opmat(xcmat, atom%basis%nbas)
237 :
238 11165 : NULLIFY (density, cpot)
239 11165 : CALL create_opgrid(density, atom%basis%grid)
240 11165 : CALL create_opgrid(cpot, atom%basis%grid)
241 :
242 : ! ZMP reading the file to restart
243 11165 : IF (atom%doread) CALL atom_read_zmp_restart(atom, doguess, iw)
244 :
245 11165 : IF (doguess) THEN
246 : ! initial guess
247 32544 : ALLOCATE (tmp_dens(SIZE(density%op)))
248 10848 : CALL slater_density(density%op, tmp_dens, atom%z, atom%state, atom%basis%grid)
249 4352818 : density%op = density%op + tmp_dens
250 10848 : DEALLOCATE (tmp_dens)
251 10848 : CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
252 10848 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
253 10848 : CALL wigner_slater_functional(density%op, cpot%op)
254 10848 : CALL numpot_matrix(xcmat%op, cpot%op, atom%basis, 0)
255 5175360 : fmat%op = hcore%op + jmat%op + xcmat%op
256 : CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
257 10848 : 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 11165 : atom%state%maxl_occ, atom%state%maxn_occ)
261 :
262 : ! wavefunction history
263 11165 : NULLIFY (history%dmat, history%hmat)
264 11165 : CALL atom_history_init(history, atom%optimization, fmat%op)
265 11165 : is_converged = .FALSE.
266 11165 : iter = 0
267 64244 : DO !SCF Loop
268 :
269 : ! Kinetic energy
270 75409 : atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
271 :
272 : ! Band energy
273 75409 : atom%energy%eband = 0._dp
274 527863 : DO l = 0, lmat
275 1082509 : DO i = 1, MIN(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
276 1007100 : 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 75409 : SELECT CASE (atom%potential%ppot_type)
282 : CASE DEFAULT
283 0 : CPABORT("")
284 : CASE (no_pseudo)
285 24164 : atom%energy%eploc = 0._dp
286 24164 : atom%energy%epnl = 0._dp
287 : CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
288 51245 : atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
289 75409 : atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
290 : END SELECT
291 75409 : atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
292 :
293 : ! Core energy
294 75409 : atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
295 :
296 : ! Confinement energy
297 75409 : IF (atom%potential%confinement) THEN
298 51013 : atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
299 : ELSE
300 24396 : atom%energy%econfinement = 0._dp
301 : END IF
302 :
303 : ! Hartree Term
304 30672175 : jmat%op = 0._dp
305 75409 : SELECT CASE (atom%coulomb_integral_type)
306 : CASE DEFAULT
307 0 : CPABORT("")
308 : CASE (do_analytic)
309 1505 : 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 38 : atom%state%maxl_occ)
313 38 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
314 : CASE (do_numeric)
315 73866 : CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
316 73866 : CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
317 149275 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
318 : END SELECT
319 75409 : atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
320 :
321 : ! Exchange Term
322 75409 : IF (need_x) THEN
323 5230416 : kmat%op = 0._dp
324 18876 : SELECT CASE (atom%exchange_integral_type)
325 : CASE DEFAULT
326 0 : CPABORT("")
327 : CASE (do_analytic)
328 535 : CALL eeri_contract(kmat%op, atom%integrals%eeri, atom%orbitals%pmat, atom%integrals%n)
329 : CASE (do_semi_analytic)
330 191 : CALL exchange_semi_analytic(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
331 : CASE (do_numeric)
332 18876 : CALL exchange_numeric(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
333 : END SELECT
334 18876 : atom%energy%eexchange = hf_frac*0.5_dp*atom_trace(kmat%op, atom%orbitals%pmat)
335 5230416 : kmat%op = hf_frac*kmat%op
336 : ELSE
337 25441759 : kmat%op = 0._dp
338 56533 : atom%energy%eexchange = 0._dp
339 : END IF
340 :
341 : ! XC
342 75409 : IF (need_xc) THEN
343 26902719 : xcmat%op = 0._dp
344 57105 : CALL calculate_atom_vxc_lda(xcmat, atom, xc_section)
345 : ! ZMP added options for the zmp calculations, building external density and vxc potential
346 18304 : ELSEIF (need_zmp) THEN
347 0 : xcmat%op = 0._dp
348 0 : CALL calculate_atom_zmp(ext_density=ext_density, atom=atom, lprint=.FALSE., xcmat=xcmat)
349 18304 : ELSEIF (need_vxc) THEN
350 0 : xcmat%op = 0._dp
351 0 : CALL calculate_atom_ext_vxc(vxc=ext_vxc, atom=atom, lprint=.FALSE., xcmat=xcmat)
352 : ELSE
353 3769456 : xcmat%op = 0._dp
354 18304 : atom%energy%exc = 0._dp
355 : END IF
356 :
357 : ! Zero this contribution
358 75409 : atom%energy%elsd = 0._dp
359 :
360 : ! Total energy
361 75409 : atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
362 :
363 : ! Potential energy
364 75409 : atom%energy%epot = atom%energy%etot - atom%energy%ekin
365 :
366 : ! Total HF/KS matrix
367 61268941 : 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 75409 : atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
372 :
373 75409 : iter = iter + 1
374 :
375 75409 : IF (iw > 0) THEN
376 9462 : IF (need_zmp) THEN
377 0 : CALL atom_print_zmp_iteration(iter, deps, atom, iw)
378 : ELSE
379 9462 : CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
380 : END IF
381 : END IF
382 :
383 75409 : IF (deps < eps_scf) THEN
384 : is_converged = .TRUE.
385 : EXIT
386 : END IF
387 64244 : IF (iter >= max_iter) THEN
388 0 : IF (iw > 0) THEN
389 0 : 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 64244 : CALL atom_history_update(history, atom%orbitals%pmat, fmat%op, jmat%op, atom%energy%etot, deps)
396 64244 : 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 64244 : 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 64244 : atom%state%maxl_occ, atom%state%maxn_occ)
403 :
404 : END DO !SCF Loop
405 :
406 11165 : IF (iw > 0) THEN
407 2201 : CALL atom_print_energies(atom, iw)
408 : END IF
409 :
410 11165 : CALL atom_history_release(history)
411 :
412 : ! conserve fmat within atom_type
413 11165 : CALL set_atom(atom, fmat=fmat)
414 11165 : CALL release_opmat(jmat)
415 11165 : CALL release_opmat(kmat)
416 11165 : CALL release_opmat(xcmat)
417 11165 : CALL release_opmat(hcore)
418 :
419 11165 : CALL release_opgrid(density)
420 11165 : CALL release_opgrid(cpot)
421 :
422 : ! ZMP deallocating ext_density ext_vxc
423 11165 : IF (need_zmp) DEALLOCATE (ext_density)
424 11165 : IF (need_vxc) DEALLOCATE (ext_vxc)
425 :
426 11165 : IF (PRESENT(converged)) THEN
427 310 : converged = is_converged
428 : END IF
429 :
430 11165 : CALL timestop(handle)
431 :
432 22330 : 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 234 : 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 117 : CALL timeset(routineN, handle)
464 :
465 117 : IF (PRESENT(noguess)) THEN
466 22 : doguess = .NOT. noguess
467 : ELSE
468 : doguess = .TRUE.
469 : END IF
470 :
471 117 : CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
472 :
473 117 : method = atom%method_type
474 117 : max_iter = atom%optimization%max_iter
475 117 : eps_scf = atom%optimization%eps_scf
476 :
477 0 : SELECT CASE (method)
478 : CASE DEFAULT
479 0 : CPABORT("")
480 : CASE (do_rks_atom)
481 0 : CPABORT("")
482 : CASE (do_uks_atom)
483 113 : need_x = do_hfx
484 113 : need_xc = .TRUE.
485 : CASE (do_rhf_atom)
486 0 : CPABORT("")
487 : CASE (do_uhf_atom)
488 4 : need_x = .TRUE.
489 4 : need_xc = .FALSE.
490 4 : hf_frac = 1._dp
491 : CASE (do_rohf_atom)
492 0 : need_x = .TRUE.
493 0 : need_xc = .FALSE.
494 0 : hf_frac = 1._dp
495 117 : CPABORT("")
496 : END SELECT
497 :
498 : ! set alpha and beta occupations
499 8307 : IF (SUM(ABS(atom%state%occa) + ABS(atom%state%occb)) == 0.0_dp) THEN
500 240 : DO l = 0, 3
501 192 : nm = REAL((2*l + 1), KIND=dp)
502 342 : DO k = 1, 10
503 294 : ne = atom%state%occupation(l, k)
504 294 : IF (ne == 0._dp) THEN !empty shell
505 : EXIT !assume there are no holes
506 102 : ELSEIF (ne == 2._dp*nm) THEN !closed shell
507 48 : atom%state%occa(l, k) = nm
508 48 : atom%state%occb(l, k) = nm
509 54 : ELSEIF (atom%state%multiplicity == -2) THEN !High spin case
510 40 : atom%state%occa(l, k) = MIN(ne, nm)
511 40 : atom%state%occb(l, k) = MAX(0._dp, ne - nm)
512 : ELSE
513 14 : atom%state%occa(l, k) = 0.5_dp*(ne + atom%state%multiplicity - 1._dp)
514 14 : 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 117 : reltyp = atom%relativistic
521 :
522 117 : IF (iw > 0) CALL atom_print_state(atom%state, iw)
523 :
524 117 : NULLIFY (hcore, hlsd)
525 117 : CALL create_opmat(hcore, atom%basis%nbas)
526 117 : CALL create_opmat(hlsd, atom%basis%nbas)
527 275739 : hlsd%op = 0._dp
528 : ! Pseudopotentials
529 117 : lsdpot = .FALSE.
530 117 : SELECT CASE (atom%potential%ppot_type)
531 : CASE DEFAULT
532 0 : CPABORT("")
533 : CASE (no_pseudo)
534 93 : SELECT CASE (reltyp)
535 : CASE DEFAULT
536 0 : CPABORT("")
537 : CASE (do_nonrel_atom)
538 86856 : hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
539 : CASE (do_zoramp_atom, do_sczoramp_atom)
540 0 : 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 24 : hcore%op = atom%integrals%hdkh
543 : END SELECT
544 : CASE (gth_pseudo)
545 464505 : hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
546 93 : IF (atom%potential%gth_pot%lsdpot) THEN
547 0 : lsdpot = .TRUE.
548 0 : hlsd%op = atom%integrals%clsd
549 : END IF
550 : CASE (upf_pseudo, sgp_pseudo, ecp_pseudo)
551 117 : hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
552 : END SELECT
553 : ! add confinement potential (not included in relativistic transformations)
554 117 : IF (atom%potential%confinement) THEN
555 496413 : hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
556 : END IF
557 :
558 117 : NULLIFY (fmata, fmatb, jmat, kmata, kmatb, xcmata, xcmatb)
559 117 : CALL create_opmat(fmata, atom%basis%nbas)
560 117 : CALL create_opmat(fmatb, atom%basis%nbas)
561 117 : CALL create_opmat(jmat, atom%basis%nbas)
562 117 : CALL create_opmat(kmata, atom%basis%nbas)
563 117 : CALL create_opmat(kmatb, atom%basis%nbas)
564 117 : CALL create_opmat(xcmata, atom%basis%nbas)
565 117 : CALL create_opmat(xcmatb, atom%basis%nbas)
566 :
567 117 : NULLIFY (density, rhoa, rhob, cpot)
568 117 : CALL create_opgrid(density, atom%basis%grid)
569 117 : CALL create_opgrid(rhoa, atom%basis%grid)
570 117 : CALL create_opgrid(rhob, atom%basis%grid)
571 117 : CALL create_opgrid(cpot, atom%basis%grid)
572 :
573 117 : IF (doguess) THEN
574 : ! initial guess
575 96 : CALL slater_density(rhoa%op, rhob%op, atom%z, atom%state, atom%basis%grid)
576 78496 : density%op = rhoa%op + rhob%op
577 96 : CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
578 96 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
579 : ! alpha spin
580 78496 : density%op = 2._dp*rhoa%op
581 96 : CALL wigner_slater_functional(density%op, cpot%op)
582 96 : CALL numpot_matrix(xcmata%op, cpot%op, atom%basis, 0)
583 137808 : 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 96 : atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
586 : ! beta spin
587 78496 : density%op = 2._dp*rhob%op
588 96 : CALL wigner_slater_functional(density%op, cpot%op)
589 96 : CALL numpot_matrix(xcmatb%op, cpot%op, atom%basis, 0)
590 137808 : 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 96 : 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 117 : 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 117 : atom%state%maxl_occ, atom%state%maxn_occ)
598 275739 : atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
599 :
600 : ! wavefunction history
601 117 : NULLIFY (historya%dmat, historya%hmat)
602 117 : CALL atom_history_init(historya, atom%optimization, fmata%op)
603 117 : NULLIFY (historyb%dmat, historyb%hmat)
604 117 : CALL atom_history_init(historyb, atom%optimization, fmatb%op)
605 :
606 117 : is_converged = .FALSE.
607 117 : iter = 0
608 : DO !SCF Loop
609 :
610 : ! Kinetic energy
611 1351 : atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
612 :
613 : ! Band energy
614 1351 : atom%energy%eband = 0._dp
615 6755 : DO l = 0, 3
616 17231 : DO i = 1, MIN(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
617 10476 : atom%energy%eband = atom%energy%eband + atom%orbitals%enera(i, l)*atom%state%occa(l, i)
618 15880 : 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 1351 : SELECT CASE (atom%potential%ppot_type)
624 : CASE DEFAULT
625 0 : CPABORT("")
626 : CASE (no_pseudo)
627 148 : atom%energy%eploc = 0._dp
628 148 : atom%energy%epnl = 0._dp
629 : CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
630 1203 : atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
631 1351 : atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
632 : END SELECT
633 1351 : atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
634 :
635 : ! Core energy
636 1351 : atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
637 :
638 : ! Confinement energy
639 1351 : IF (atom%potential%confinement) THEN
640 747 : atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
641 : ELSE
642 604 : atom%energy%econfinement = 0._dp
643 : END IF
644 :
645 : ! Hartree Term
646 3696193 : jmat%op = 0._dp
647 1351 : SELECT CASE (atom%coulomb_integral_type)
648 : CASE DEFAULT
649 0 : CPABORT("")
650 : CASE (do_analytic)
651 10 : 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 0 : atom%state%maxl_occ)
655 0 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
656 : CASE (do_numeric)
657 1341 : CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
658 1341 : CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
659 2692 : CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
660 : END SELECT
661 1351 : atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
662 :
663 : ! Exchange Term
664 1351 : IF (need_x) THEN
665 295546 : kmata%op = 0._dp
666 295546 : kmatb%op = 0._dp
667 118 : SELECT CASE (atom%exchange_integral_type)
668 : CASE DEFAULT
669 0 : CPABORT("")
670 : CASE (do_analytic)
671 4 : CALL eeri_contract(kmata%op, atom%integrals%eeri, atom%orbitals%pmata, atom%integrals%n)
672 4 : CALL eeri_contract(kmatb%op, atom%integrals%eeri, atom%orbitals%pmatb, atom%integrals%n)
673 : CASE (do_semi_analytic)
674 0 : CALL exchange_semi_analytic(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
675 0 : CALL exchange_semi_analytic(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
676 : CASE (do_numeric)
677 114 : CALL exchange_numeric(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
678 232 : 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 118 : atom_trace(kmatb%op, atom%orbitals%pmatb))
682 295546 : kmata%op = 2._dp*hf_frac*kmata%op
683 295546 : kmatb%op = 2._dp*hf_frac*kmatb%op
684 : ELSE
685 3400647 : kmata%op = 0._dp
686 3400647 : kmatb%op = 0._dp
687 1233 : atom%energy%eexchange = 0._dp
688 : END IF
689 :
690 : ! XC
691 1351 : IF (need_xc) THEN
692 3558249 : xcmata%op = 0._dp
693 3558249 : xcmatb%op = 0._dp
694 1335 : CALL calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
695 : ELSE
696 137944 : xcmata%op = 0._dp
697 137944 : xcmatb%op = 0._dp
698 16 : atom%energy%exc = 0._dp
699 : END IF
700 :
701 1351 : IF (lsdpot) THEN
702 : atom%energy%elsd = atom_trace(hlsd%op, atom%orbitals%pmata) - &
703 0 : atom_trace(hlsd%op, atom%orbitals%pmatb)
704 0 : atom%energy%epseudo = atom%energy%epseudo + atom%energy%elsd
705 0 : atom%energy%ecore = atom%energy%ecore + atom%energy%elsd
706 : ELSE
707 1351 : atom%energy%elsd = 0._dp
708 : END IF
709 :
710 : ! Total energy
711 1351 : atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
712 :
713 : ! Potential energy
714 1351 : atom%energy%epot = atom%energy%etot - atom%energy%ekin
715 :
716 : ! Total HF/KS matrix
717 7391035 : fmata%op = hcore%op + hlsd%op + jmat%op + kmata%op + xcmata%op
718 7391035 : 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 1351 : 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 1351 : atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
725 1351 : deps = 2._dp*MAX(depsa, depsb)
726 :
727 1351 : iter = iter + 1
728 :
729 1351 : IF (iw > 0) THEN
730 504 : CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
731 : END IF
732 :
733 1351 : IF (deps < eps_scf) THEN
734 : is_converged = .TRUE.
735 : EXIT
736 : END IF
737 1234 : IF (iter >= max_iter) THEN
738 0 : IF (iw > 0) THEN
739 0 : 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 1234 : CALL atom_history_update(historya, atom%orbitals%pmata, fmata%op, xcmata%op, atom%energy%etot, deps)
746 1234 : CALL atom_history_update(historyb, atom%orbitals%pmatb, fmatb%op, xcmatb%op, atom%energy%etot, deps)
747 1234 : CALL atom_opt_fmat(fmata%op, historya, depsa)
748 1234 : 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 1234 : 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 1234 : atom%state%maxl_occ, atom%state%maxn_occ)
755 : CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
756 1234 : 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 1234 : atom%state%maxl_occ, atom%state%maxn_occ)
759 3420454 : atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
760 :
761 : END DO !SCF Loop
762 :
763 117 : IF (iw > 0) THEN
764 45 : CALL atom_print_energies(atom, iw)
765 : END IF
766 :
767 117 : CALL atom_history_release(historya)
768 117 : CALL atom_history_release(historyb)
769 :
770 117 : CALL release_opgrid(density)
771 117 : CALL release_opgrid(rhoa)
772 117 : CALL release_opgrid(rhob)
773 117 : CALL release_opgrid(cpot)
774 :
775 : ! conserve fmata as fmat within atom_type.
776 117 : CALL set_atom(atom, fmat=fmata)
777 117 : CALL release_opmat(fmatb)
778 117 : CALL release_opmat(jmat)
779 117 : CALL release_opmat(kmata)
780 117 : CALL release_opmat(kmatb)
781 117 : CALL release_opmat(xcmata)
782 117 : CALL release_opmat(xcmatb)
783 117 : CALL release_opmat(hlsd)
784 117 : CALL release_opmat(hcore)
785 :
786 117 : IF (PRESENT(converged)) THEN
787 22 : converged = is_converged
788 : END IF
789 :
790 117 : CALL timestop(handle)
791 :
792 117 : END SUBROUTINE calculate_atom_unrestricted
793 :
794 : END MODULE atom_electronic_structure
|
