Line data Source code
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 :
8 : ! **************************************************************************************************
9 : !> \brief Routines that print various information about an atomic kind.
10 : ! **************************************************************************************************
11 : MODULE atom_output
12 : USE atom_types, ONLY: &
13 : atom_basis_type, atom_gthpot_type, atom_potential_type, atom_state, atom_type, cgto_basis, &
14 : ecp_pseudo, gth_pseudo, gto_basis, lmat, no_pseudo, num_basis, sgp_pseudo, sto_basis, &
15 : upf_pseudo
16 : USE atom_utils, ONLY: get_maxl_occ,&
17 : get_maxn_occ,&
18 : get_rho0
19 : USE cp_files, ONLY: close_file,&
20 : open_file
21 : USE input_constants, ONLY: &
22 : barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
23 : do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_uhf_atom, do_uks_atom, &
24 : do_zoramp_atom, poly_conf, xc_none
25 : USE input_cp2k_check, ONLY: xc_functionals_expand
26 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 : section_vals_get_subs_vals2,&
28 : section_vals_type,&
29 : section_vals_val_get
30 : USE kinds, ONLY: default_string_length,&
31 : dp
32 : USE mathconstants, ONLY: dfac,&
33 : pi,&
34 : rootpi
35 : USE periodic_table, ONLY: ptable
36 : USE physcon, ONLY: evolt
37 : USE xc_derivatives, ONLY: xc_functional_get_info
38 : USE xc_libxc, ONLY: libxc_check_existence_in_libxc,&
39 : libxc_get_reference_length
40 : USE xmgrace, ONLY: xm_graph_data,&
41 : xm_graph_info,&
42 : xm_write_defaults,&
43 : xm_write_frame,&
44 : xm_write_frameport
45 : #include "./base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_output'
52 :
53 : PUBLIC :: atom_print_state, atom_print_energies, atom_print_iteration, &
54 : atom_print_basis, atom_print_method, atom_print_info, atom_print_potential, &
55 : atom_print_basis_file, atom_write_pseudo_param, atom_print_orbitals, &
56 : atom_print_zmp_iteration
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Print an information string related to the atomic kind.
62 : !> \param zval atomic number
63 : !> \param info information string
64 : !> \param iw output file unit
65 : !> \par History
66 : !> * 09.2008 created [Juerg Hutter]
67 : ! **************************************************************************************************
68 180 : SUBROUTINE atom_print_info(zval, info, iw)
69 : INTEGER, INTENT(IN) :: zval
70 : CHARACTER(len=*), INTENT(IN) :: info
71 : INTEGER, INTENT(IN) :: iw
72 :
73 : WRITE (iw, '(/," ",A,T40,A," [",A,"]",T62,"Atomic number:",T78,I3,/)') &
74 180 : ADJUSTL(TRIM(info)), TRIM(ptable(zval)%name), TRIM(ptable(zval)%symbol), zval
75 :
76 180 : END SUBROUTINE atom_print_info
77 :
78 : ! **************************************************************************************************
79 : !> \brief Print information about electronic state.
80 : !> \param state electronic state
81 : !> \param iw output file unit
82 : !> \par History
83 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
84 : !> * 11.2009 print multiplicity [Juerg Hutter]
85 : !> * 08.2008 created [Juerg Hutter]
86 : ! **************************************************************************************************
87 2258 : SUBROUTINE atom_print_state(state, iw)
88 : TYPE(atom_state) :: state
89 : INTEGER, INTENT(IN) :: iw
90 :
91 : CHARACTER(LEN=1), DIMENSION(0:7), PARAMETER :: &
92 : label = (/"S", "P", "D", "F", "G", "H", "I", "K"/)
93 :
94 : INTEGER :: j, l, mc, mlc, mlo, mm(0:lmat), mo
95 :
96 : CPASSERT(lmat <= 7)
97 2258 : WRITE (iw, '(/,T2,A)') "Electronic structure"
98 160318 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of core electrons", SUM(state%core)
99 160318 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of valence electrons", SUM(state%occ)
100 160318 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of electrons", SUM(state%occ + state%core)
101 4471 : SELECT CASE (state%multiplicity)
102 : CASE (-1)
103 2213 : WRITE (iw, '(T5,A,T68,A)') "Multiplicity", "not specified"
104 : CASE (-2)
105 10 : WRITE (iw, '(T5,A,T72,A)') "Multiplicity", "high spin"
106 : CASE (-3)
107 0 : WRITE (iw, '(T5,A,T73,A)') "Multiplicity", "low spin"
108 : CASE (1)
109 22 : WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "singlet"
110 : CASE (2)
111 10 : WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "doublet"
112 : CASE (3)
113 3 : WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "triplet"
114 : CASE (4)
115 0 : WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quartet"
116 : CASE (5)
117 0 : WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quintet"
118 : CASE (6)
119 0 : WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "sextet"
120 : CASE (7)
121 2258 : WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "septet"
122 : CASE DEFAULT
123 : END SELECT
124 :
125 2258 : mlo = get_maxl_occ(state%occ)
126 2258 : mlc = get_maxl_occ(state%core)
127 2258 : mm = get_maxn_occ(state%core)
128 :
129 2258 : IF (state%multiplicity == -1) THEN
130 5722 : DO l = 0, MAX(mlo, mlc)
131 3509 : mo = state%maxn_occ(l)
132 40812 : IF (SUM(state%core(l, :)) == 0) THEN
133 5771 : WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occ(l, j), j=1, mo)
134 : ELSE
135 1086 : mc = mm(l)
136 2448 : CPASSERT(SUM(state%occ(l, 1:mc)) == 0)
137 2448 : WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (state%core(l, j), j=1, mc)
138 2141 : WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occ(l, j), j=mc + 1, mc + mo)
139 : END IF
140 : END DO
141 : ELSE
142 45 : WRITE (iw, '(T5,A)') "Alpha Electrons"
143 122 : DO l = 0, MAX(mlo, mlc)
144 77 : mo = state%maxn_occ(l)
145 892 : IF (SUM(state%core(l, :)) == 0) THEN
146 104 : WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occa(l, j), j=1, mo)
147 : ELSE
148 35 : mc = mm(l)
149 86 : WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
150 60 : WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occa(l, j), j=1, mo)
151 : END IF
152 : END DO
153 45 : WRITE (iw, '(T5,A)') "Beta Electrons"
154 122 : DO l = 0, MAX(mlo, mlc)
155 77 : mo = state%maxn_occ(l)
156 892 : IF (SUM(state%core(l, :)) == 0) THEN
157 104 : WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occb(l, j), j=1, mo)
158 : ELSE
159 35 : mc = mm(l)
160 86 : WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
161 60 : WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occb(l, j), j=1, mo)
162 : END IF
163 : END DO
164 : END IF
165 2258 : WRITE (iw, *)
166 :
167 2258 : END SUBROUTINE atom_print_state
168 :
169 : ! **************************************************************************************************
170 : !> \brief Print energy components.
171 : !> \param atom information about the atomic kind
172 : !> \param iw output file unit
173 : !> \par History
174 : !> * 05.2010 print virial coefficient [Juerg Hutter]
175 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
176 : !> * 09.2008 print orbital energies [Juerg Hutter]
177 : !> * 08.2008 created [Juerg Hutter]
178 : ! **************************************************************************************************
179 4516 : SUBROUTINE atom_print_energies(atom, iw)
180 : TYPE(atom_type) :: atom
181 : INTEGER, INTENT(IN) :: iw
182 :
183 : INTEGER :: i, l, n
184 : REAL(KIND=dp) :: drho
185 :
186 2258 : WRITE (iw, '(/,A,T36,A,T61,F20.12)') " Energy components [Hartree]", &
187 4516 : " Total Energy ::", atom%energy%etot
188 2258 : WRITE (iw, '(T36,A,T61,F20.12)') " Band Energy ::", atom%energy%eband
189 2258 : WRITE (iw, '(T36,A,T61,F20.12)') " Kinetic Energy ::", atom%energy%ekin
190 2258 : WRITE (iw, '(T36,A,T61,F20.12)') "Potential Energy ::", atom%energy%epot
191 2258 : IF (atom%energy%ekin /= 0.0_dp) THEN
192 2245 : WRITE (iw, '(T36,A,T61,F20.12)') " Virial (-V/T) ::", -atom%energy%epot/atom%energy%ekin
193 : END IF
194 2258 : WRITE (iw, '(T36,A,T61,F20.12)') " Core Energy ::", atom%energy%ecore
195 2258 : IF (atom%energy%exc /= 0._dp) &
196 2227 : WRITE (iw, '(T36,A,T61,F20.12)') " XC Energy ::", atom%energy%exc
197 2258 : WRITE (iw, '(T36,A,T61,F20.12)') " Coulomb Energy ::", atom%energy%ecoulomb
198 2258 : IF (atom%energy%eexchange /= 0._dp) &
199 46 : WRITE (iw, '(T34,A,T61,F20.12)') "HF Exchange Energy ::", atom%energy%eexchange
200 2258 : IF (atom%potential%ppot_type /= NO_PSEUDO) THEN
201 1868 : WRITE (iw, '(T20,A,T61,F20.12)') " Total Pseudopotential Energy ::", atom%energy%epseudo
202 1868 : WRITE (iw, '(T20,A,T61,F20.12)') " Local Pseudopotential Energy ::", atom%energy%eploc
203 1868 : IF (atom%energy%elsd /= 0._dp) &
204 0 : WRITE (iw, '(T20,A,T61,F20.12)') " Local Spin-potential Energy ::", atom%energy%elsd
205 1868 : WRITE (iw, '(T20,A,T61,F20.12)') " Nonlocal Pseudopotential Energy ::", atom%energy%epnl
206 : END IF
207 2258 : IF (atom%potential%confinement) THEN
208 1862 : WRITE (iw, '(T36,A,T61,F20.12)') " Confinement ::", atom%energy%econfinement
209 : END IF
210 :
211 2258 : IF (atom%state%multiplicity == -1) THEN
212 2213 : WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T49,A,T71,A,/)') " Orbital energies", &
213 4426 : "State", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
214 5832 : DO l = 0, atom%state%maxl_calc
215 3619 : n = atom%state%maxn_calc(l)
216 8137 : DO i = 1, n
217 : WRITE (iw, '(T23,I2,T30,I1,T36,F10.3,T46,F15.6,T66,F15.6)') &
218 8137 : i, l, atom%state%occupation(l, i), atom%orbitals%ener(i, l), atom%orbitals%ener(i, l)*evolt
219 : END DO
220 5832 : IF (n > 0) WRITE (iw, *)
221 : END DO
222 : ELSE
223 45 : WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T42,A,T55,A,T71,A,/)') " Orbital energies", &
224 90 : "State", "Spin", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
225 162 : DO l = 0, atom%state%maxl_calc
226 117 : n = atom%state%maxn_calc(l)
227 208 : DO i = 1, n
228 : WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
229 208 : i, "alpha", l, atom%state%occa(l, i), atom%orbitals%enera(i, l), atom%orbitals%enera(i, l)*evolt
230 : END DO
231 208 : DO i = 1, n
232 : WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
233 208 : i, " beta", l, atom%state%occb(l, i), atom%orbitals%enerb(i, l), atom%orbitals%enerb(i, l)*evolt
234 : END DO
235 162 : IF (n > 0) WRITE (iw, *)
236 : END DO
237 : END IF
238 :
239 2258 : CALL get_rho0(atom, drho)
240 2258 : WRITE (iw, '(/,A,T66,F15.6)') " Total Electron Density at R=0: ", drho
241 :
242 2258 : END SUBROUTINE atom_print_energies
243 :
244 : ! **************************************************************************************************
245 : !> \brief Printing of the atomic iterations when ZMP is active.
246 : !> \param iter current iteration number
247 : !> \param deps convergence
248 : !> \param atom intormation about the atomic kind
249 : !> \param iw output file unit
250 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
251 : ! **************************************************************************************************
252 0 : SUBROUTINE atom_print_zmp_iteration(iter, deps, atom, iw)
253 : INTEGER, INTENT(IN) :: iter
254 : REAL(dp), INTENT(IN) :: deps
255 : TYPE(atom_type), INTENT(IN) :: atom
256 : INTEGER, INTENT(IN) :: iw
257 :
258 0 : IF (iter == 1) THEN
259 : WRITE (iw, '(/," ",79("*"),/,T33,"Integral",T48,"Integral",/,T3,A,T16,A,T33,A,T46,A,T69,A/," ",79("*"))') &
260 0 : "Iteration", "Convergence", "rho diff.", "rho*v_xc[au]", "Energy[au]"
261 : END IF
262 0 : WRITE (iw, '(T3,I9,T15,G13.6,T30,G13.6,T46,G13.6,T61,F20.12)') iter, deps, atom%rho_diff_integral, &
263 0 : atom%energy%exc, atom%energy%etot
264 :
265 0 : END SUBROUTINE atom_print_zmp_iteration
266 :
267 : ! **************************************************************************************************
268 : !> \brief Print convergence information.
269 : !> \param iter current iteration number
270 : !> \param deps convergency
271 : !> \param etot total energy
272 : !> \param iw output file unit
273 : !> \par History
274 : !> * 08.2008 created [Juerg Hutter]
275 : ! **************************************************************************************************
276 10025 : SUBROUTINE atom_print_iteration(iter, deps, etot, iw)
277 : INTEGER, INTENT(IN) :: iter
278 : REAL(dp), INTENT(IN) :: deps, etot
279 : INTEGER, INTENT(IN) :: iw
280 :
281 10025 : IF (iter == 1) THEN
282 : WRITE (iw, '(/," ",79("*"),/,T19,A,T38,A,T70,A,/," ",79("*"))') &
283 2258 : "Iteration", "Convergence", "Energy [au]"
284 : END IF
285 10025 : WRITE (iw, '(T20,i8,T34,G14.6,T61,F20.12)') iter, deps, etot
286 :
287 10025 : END SUBROUTINE atom_print_iteration
288 :
289 : ! **************************************************************************************************
290 : !> \brief Print atomic basis set.
291 : !> \param atom_basis atomic basis set
292 : !> \param iw output file unit
293 : !> \param title header to print on top of the basis set
294 : !> \par History
295 : !> * 09.2008 created [Juerg Hutter]
296 : ! **************************************************************************************************
297 22 : SUBROUTINE atom_print_basis(atom_basis, iw, title)
298 : TYPE(atom_basis_type) :: atom_basis
299 : INTEGER, INTENT(IN) :: iw
300 : CHARACTER(len=*) :: title
301 :
302 : INTEGER :: i, j, l
303 :
304 22 : WRITE (iw, '(/,A)') TRIM(title)
305 39 : SELECT CASE (atom_basis%basis_type)
306 : CASE (GTO_BASIS)
307 17 : IF (atom_basis%geometrical) THEN
308 15 : WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
309 15 : WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
310 30 : " Proportionality factor: ", atom_basis%cval
311 : ELSE
312 2 : WRITE (iw, '(/," ",21("*"),A,21("*"))') " Uncontracted Gaussian Type Orbitals "
313 : END IF
314 119 : DO l = 0, lmat
315 119 : IF (atom_basis%nbas(l) > 0) THEN
316 14 : SELECT CASE (l)
317 : CASE DEFAULT
318 : WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
319 280 : "X Exponents: ", (i, atom_basis%am(i, l), i=1, atom_basis%nbas(l))
320 : CASE (0)
321 : WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
322 440 : "s Exponents: ", (i, atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
323 : CASE (1)
324 : WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
325 408 : "p Exponents: ", (i, atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
326 : CASE (2)
327 : WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
328 376 : "d Exponents: ", (i, atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
329 : CASE (3)
330 : WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
331 395 : "f Exponents: ", (i, atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
332 : END SELECT
333 : END IF
334 : END DO
335 17 : WRITE (iw, '(" ",79("*"))')
336 : CASE (CGTO_BASIS)
337 1 : WRITE (iw, '(/," ",22("*"),A,22("*"))') " Contracted Gaussian Type Orbitals "
338 7 : DO l = 0, lmat
339 7 : IF (atom_basis%nbas(l) > 0) THEN
340 2 : IF (l == 0) WRITE (iw, '(A)') " s Functions"
341 2 : IF (l == 1) WRITE (iw, '(A)') " p Functions"
342 2 : IF (l == 2) WRITE (iw, '(A)') " d Functions"
343 2 : IF (l == 3) WRITE (iw, '(A)') " f Functions"
344 2 : IF (l >= 3) WRITE (iw, '(A)') " x Functions"
345 11 : DO i = 1, atom_basis%nprim(l)
346 : WRITE (iw, '(F15.6,5(T21,6F10.6,/))') &
347 35 : atom_basis%am(i, l), (atom_basis%cm(i, j, l), j=1, atom_basis%nbas(l))
348 : END DO
349 : END IF
350 : END DO
351 1 : WRITE (iw, '(" ",79("*"))')
352 : CASE (STO_BASIS)
353 4 : WRITE (iw, '(/," ",28("*"),A,29("*"))') " Slater Type Orbitals "
354 28 : DO l = 0, lmat
355 39 : DO i = 1, atom_basis%nbas(l)
356 24 : SELECT CASE (l)
357 : CASE DEFAULT
358 0 : WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, l), "X Exponent :", atom_basis%as(i, l)
359 : CASE (0)
360 8 : WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 0), "S Exponent :", atom_basis%as(i, 0)
361 : CASE (1)
362 3 : WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 1), "P Exponent :", atom_basis%as(i, 1)
363 : CASE (2)
364 0 : WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 2), "D Exponent :", atom_basis%as(i, 2)
365 : CASE (3)
366 11 : WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 3), "F Exponent :", atom_basis%as(i, 3)
367 : END SELECT
368 : END DO
369 : END DO
370 4 : WRITE (iw, '(" ",79("*"))')
371 : CASE (NUM_BASIS)
372 0 : CPABORT("")
373 : CASE DEFAULT
374 22 : CPABORT("")
375 : END SELECT
376 :
377 22 : END SUBROUTINE atom_print_basis
378 :
379 : ! **************************************************************************************************
380 : !> \brief Print the optimized atomic basis set into a file.
381 : !> \param atom_basis atomic basis set
382 : !> \param wfn ...
383 : !> \par History
384 : !> * 11.2016 revised output format [Matthias Krack]
385 : !> * 11.2011 Slater basis functions [Juerg Hutter]
386 : !> * 03.2011 created [Juerg Hutter]
387 : !> \note The basis set is stored as the file 'OPT_BASIS' inside the current working directory.
388 : !> It may be a good idea, however, to specify the name of this file via some input section.
389 : ! **************************************************************************************************
390 5 : SUBROUTINE atom_print_basis_file(atom_basis, wfn)
391 : TYPE(atom_basis_type) :: atom_basis
392 : REAL(KIND=dp), DIMENSION(:, :, 0:), OPTIONAL :: wfn
393 :
394 : INTEGER :: i, im, iw, l
395 : REAL(KIND=dp) :: expzet, prefac, zeta
396 :
397 5 : CALL open_file(file_name="OPT_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
398 7 : SELECT CASE (atom_basis%basis_type)
399 : CASE (GTO_BASIS)
400 2 : IF (atom_basis%geometrical) THEN
401 0 : WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
402 0 : WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
403 0 : " Proportionality factor: ", atom_basis%cval
404 : ELSE
405 2 : WRITE (iw, '(T3,A)') "BASIS_TYPE GAUSSIAN"
406 : END IF
407 14 : DO l = 0, lmat
408 14 : IF (atom_basis%nbas(l) > 0) THEN
409 0 : SELECT CASE (l)
410 : CASE DEFAULT
411 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
412 0 : "X_EXPONENTS ", (atom_basis%am(i, l), i=1, atom_basis%nbas(l))
413 : CASE (0)
414 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
415 14 : "S_EXPONENTS ", (atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
416 : CASE (1)
417 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
418 14 : "P_EXPONENTS ", (atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
419 : CASE (2)
420 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
421 14 : "D_EXPONENTS ", (atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
422 : CASE (3)
423 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
424 6 : "F_EXPONENTS ", (atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
425 : END SELECT
426 : END IF
427 : END DO
428 : CASE (CGTO_BASIS)
429 0 : CPABORT("")
430 : CASE (STO_BASIS)
431 3 : WRITE (iw, '(T3,A)') "BASIS_TYPE SLATER"
432 21 : DO l = 0, lmat
433 21 : IF (atom_basis%nbas(l) > 0) THEN
434 0 : SELECT CASE (l)
435 : CASE DEFAULT
436 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
437 0 : "X_EXPONENTS ", (atom_basis%as(i, l), i=1, atom_basis%nbas(l))
438 : WRITE (iw, '(T3,A,60I3)') &
439 0 : "X_QUANTUM_NUMBERS ", (atom_basis%ns(i, l), i=1, atom_basis%nbas(l))
440 : CASE (0)
441 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
442 10 : "S_EXPONENTS ", (atom_basis%as(i, 0), i=1, atom_basis%nbas(0))
443 : WRITE (iw, '(T3,A,60I3)') &
444 10 : "S_QUANTUM_NUMBERS ", (atom_basis%ns(i, 0), i=1, atom_basis%nbas(0))
445 : CASE (1)
446 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
447 3 : "P_EXPONENTS ", (atom_basis%as(i, 1), i=1, atom_basis%nbas(1))
448 : WRITE (iw, '(T3,A,60I3)') &
449 3 : "P_QUANTUM_NUMBERS ", (atom_basis%ns(i, 1), i=1, atom_basis%nbas(1))
450 : CASE (2)
451 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
452 0 : "D_EXPONENTS ", (atom_basis%as(i, 2), i=1, atom_basis%nbas(2))
453 : WRITE (iw, '(T3,A,60I3)') &
454 0 : "D_QUANTUM_NUMBERS ", (atom_basis%ns(i, 2), i=1, atom_basis%nbas(2))
455 : CASE (3)
456 : WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
457 0 : "F_EXPONENTS ", (atom_basis%as(i, 3), i=1, atom_basis%nbas(3))
458 : WRITE (iw, '(T3,A,60I3)') &
459 4 : "F_QUANTUM_NUMBERS ", (atom_basis%ns(i, 3), i=1, atom_basis%nbas(3))
460 : END SELECT
461 : END IF
462 : END DO
463 : CASE (NUM_BASIS)
464 0 : CPABORT("")
465 : CASE DEFAULT
466 5 : CPABORT("")
467 : END SELECT
468 :
469 5 : IF (PRESENT(wfn)) THEN
470 7 : SELECT CASE (atom_basis%basis_type)
471 : CASE DEFAULT
472 : CASE (GTO_BASIS)
473 5 : IF (.NOT. atom_basis%geometrical) THEN
474 2 : WRITE (iw, '(/,T3,A)') "ORBITAL COEFFICENTS (Quickstep normalization)"
475 2 : im = MIN(6, SIZE(wfn, 2))
476 14 : DO l = 0, lmat
477 14 : IF (atom_basis%nbas(l) > 0) THEN
478 6 : WRITE (iw, '(T3,A,I3)') "L Quantum Number:", l
479 : ! Quickstep normalization
480 6 : expzet = 0.25_dp*REAL(2*l + 3, dp)
481 6 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
482 42 : DO i = 1, atom_basis%nbas(l)
483 36 : zeta = (2._dp*atom_basis%am(i, l))**expzet
484 78 : WRITE (iw, '(T5,F14.8,2x,6F12.8)') atom_basis%am(i, l), wfn(i, 1:im, l)*prefac/zeta
485 : END DO
486 : END IF
487 : END DO
488 : END IF
489 : END SELECT
490 : END IF
491 :
492 5 : CALL close_file(unit_number=iw)
493 :
494 5 : END SUBROUTINE atom_print_basis_file
495 :
496 : ! **************************************************************************************************
497 : !> \brief Print information about the electronic structure method in use.
498 : !> \param atom information about the atomic kind
499 : !> \param iw output file unit
500 : !> \par History
501 : !> * 09.2015 direct use of the LibXC Fortran interface [Andreas Gloess]
502 : !> * 10.2012 LibXC interface [Fabien Tran]
503 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
504 : !> * 04.2009 print geometrical Gaussian type orbitals [Juerg Hutter]
505 : !> * 09.2008 new subroutine's prototype; print relativistic methods [Juerg Hutter]
506 : !> * 09.2008 created [Juerg Hutter]
507 : ! **************************************************************************************************
508 366 : SUBROUTINE atom_print_method(atom, iw)
509 : TYPE(atom_type) :: atom
510 : INTEGER, INTENT(IN) :: iw
511 :
512 : CHARACTER(len=160) :: shortform
513 : CHARACTER(LEN=20) :: tmpStr
514 366 : CHARACTER(len=:), ALLOCATABLE :: reference
515 : INTEGER :: ifun, il, meth, myfun, reltyp
516 : LOGICAL :: lsd
517 : TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section, xc_section
518 :
519 366 : NULLIFY (xc_fun, xc_fun_section, xc_section)
520 :
521 366 : meth = atom%method_type
522 :
523 366 : xc_section => atom%xc_section
524 366 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
525 0 : SELECT CASE (meth)
526 : CASE DEFAULT
527 0 : CPABORT("")
528 : CASE (do_rks_atom)
529 328 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
530 : CASE (do_uks_atom)
531 34 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
532 : CASE (do_rhf_atom)
533 34 : myfun = xc_none
534 : CASE (do_uhf_atom)
535 4 : myfun = xc_none
536 : CASE (do_rohf_atom)
537 366 : myfun = xc_none
538 : END SELECT
539 :
540 0 : SELECT CASE (meth)
541 : CASE DEFAULT
542 0 : CPABORT("")
543 : CASE (do_rks_atom)
544 294 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Restricted Kohn-Sham Calculation')")
545 : CASE (do_uks_atom)
546 34 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Unrestricted Kohn-Sham Calculation')")
547 : CASE (do_rhf_atom)
548 34 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Restricted Hartree-Fock Calculation')")
549 : CASE (do_uhf_atom)
550 4 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Unrestricted Hartree-Fock Calculation')")
551 : CASE (do_rohf_atom)
552 328 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Restricted Open-Shell Kohn-Sham Calculation')")
553 : END SELECT
554 :
555 : ! zmp
556 366 : IF (atom%do_zmp) THEN
557 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Method on atomic radial density')")
558 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Lambda : ',F5.1)") atom%lambda
559 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Reading external density : ',A20)") atom%ext_file
560 0 : IF (atom%dm) THEN
561 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | The file is in the form of a density matrix')")
562 : ELSE
563 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | The file is in the form of a linear density')")
564 : END IF
565 0 : IF (atom%doread) THEN
566 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Restarting calculation from ',A20,' file if present')") atom%zmp_restart_file
567 : END IF
568 366 : ELSE IF (atom%read_vxc) THEN
569 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Calculating density from external V_xc')")
570 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Reading external v_xc file : ',A20)") atom%ext_vxc_file
571 : END IF
572 :
573 366 : IF (atom%pp_calc) THEN
574 78 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Nonrelativistic Calculation')")
575 : ELSE
576 288 : reltyp = atom%relativistic
577 :
578 0 : SELECT CASE (reltyp)
579 : CASE DEFAULT
580 0 : CPABORT("")
581 : CASE (do_nonrel_atom)
582 220 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Nonrelativistic Calculation')")
583 : CASE (do_zoramp_atom)
584 14 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using ZORA(MP)')")
585 : CASE (do_sczoramp_atom)
586 2 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using scaled ZORA(MP)')")
587 : CASE (do_dkh0_atom)
588 2 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using Douglas-Kroll 0th order')")
589 2 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using kietic energy scaling')")
590 : CASE (do_dkh1_atom)
591 2 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using Douglas-Kroll 1st order')")
592 2 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using Foldy-Wouthuysen transformation')")
593 : CASE (do_dkh2_atom)
594 16 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using Douglas-Kroll 2nd order')")
595 : CASE (do_dkh3_atom)
596 288 : IF (iw > 0) WRITE (iw, fmt="(' METHOD | Relativistic Calculation using Douglas-Kroll 3rd order')")
597 : END SELECT
598 : END IF
599 :
600 366 : lsd = (meth == do_uks_atom)
601 :
602 366 : IF (myfun /= xc_none) THEN
603 326 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", c_val=tmpStr)
604 326 : IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| ROUTINE=',a)") TRIM(tmpStr)
605 326 : CALL xc_functionals_expand(xc_fun_section, xc_section)
606 326 : IF (iw > 0) THEN
607 163 : ifun = 0
608 232 : DO
609 395 : ifun = ifun + 1
610 395 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
611 395 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
612 232 : IF (libxc_check_existence_in_libxc(xc_fun)) THEN
613 3 : ALLOCATE (CHARACTER(LEN=libxc_get_reference_length(xc_fun, lsd)) :: reference)
614 : ELSE
615 229 : ALLOCATE (CHARACTER(LEN=20*default_string_length) :: reference)
616 : END IF
617 232 : CALL xc_functional_get_info(xc_fun, lsd=lsd, reference=reference, shortform=shortform)
618 : WRITE (iw, fmt="(' FUNCTIONAL| ',a,':')") &
619 232 : TRIM(xc_fun%section%name)
620 630 : DO il = 1, LEN_TRIM(reference), 67
621 630 : WRITE (iw, fmt="(' FUNCTIONAL| ',a67)") reference(il:)
622 : END DO
623 395 : DEALLOCATE (reference)
624 : END DO
625 : END IF
626 : ELSE
627 40 : IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| NO EXCHANGE-CORRELATION FUNCTIONAL USED.')")
628 : END IF
629 :
630 366 : END SUBROUTINE atom_print_method
631 :
632 : ! **************************************************************************************************
633 : !> \brief Print information about the pseudo-potential.
634 : !> \param potential pseudo-potential
635 : !> \param iw output file unit
636 : !> \par History
637 : !> * 05.2017 SGP pseudo-potentials [Juerg Hutter]
638 : !> * 02.2016 pseudo-potential in Quantum Espresso UPF format [Juerg Hutter]
639 : !> * 01.2016 new confinement potential form [Juerg Hutter]
640 : !> * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
641 : !> * 05.2009 GTH pseudo-potential [Juerg Hutter]
642 : !> * 09.2008 created [Juerg Hutter]
643 : ! **************************************************************************************************
644 3 : SUBROUTINE atom_print_potential(potential, iw)
645 : TYPE(atom_potential_type) :: potential
646 : INTEGER, INTENT(IN) :: iw
647 :
648 : CHARACTER(len=60) :: pline
649 : INTEGER :: i, j, k, l
650 :
651 3 : SELECT CASE (potential%ppot_type)
652 : CASE (no_pseudo)
653 0 : WRITE (iw, '(/," ",28("*"),A,27("*"))') " All Electron Potential "
654 : CASE (gth_pseudo)
655 0 : WRITE (iw, '(/," ",29("*"),A,29("*"))') " GTH Pseudopotential "
656 0 : WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%gth_pot%zion
657 0 : WRITE (iw, '(T10,A,T66,F15.6)') " Rc ", potential%gth_pot%rc
658 0 : WRITE (pline, '(5F12.6)') (potential%gth_pot%cl(i), i=1, potential%gth_pot%ncl)
659 0 : WRITE (iw, '(T10,A,T21,A60)') " C1 C2 ... ", ADJUSTR(pline)
660 0 : IF (potential%gth_pot%lpotextended) THEN
661 0 : DO k = 1, potential%gth_pot%nexp_lpot
662 0 : WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LPot: rc=", potential%gth_pot%alpha_lpot(k), &
663 0 : "CX=", (potential%gth_pot%cval_lpot(i, k), i=1, potential%gth_pot%nct_lpot(k))
664 : END DO
665 : END IF
666 0 : IF (potential%gth_pot%nlcc) THEN
667 0 : DO k = 1, potential%gth_pot%nexp_nlcc
668 0 : WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_nlcc(k), &
669 0 : "CX=", (potential%gth_pot%cval_nlcc(i, k)*4.0_dp*pi, i=1, potential%gth_pot%nct_nlcc(k))
670 : END DO
671 : END IF
672 0 : IF (potential%gth_pot%lsdpot) THEN
673 0 : DO k = 1, potential%gth_pot%nexp_lsd
674 0 : WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_lsd(k), &
675 0 : "CX=", (potential%gth_pot%cval_lsd(i, k), i=1, potential%gth_pot%nct_lsd(k))
676 : END DO
677 : END IF
678 0 : DO l = 0, lmat
679 0 : IF (potential%gth_pot%nl(l) > 0) THEN
680 0 : WRITE (iw, '(T10,A,T76,I5)') " Angular momentum ", l
681 0 : WRITE (iw, '(T10,A,T66,F15.6)') " Rcnl ", potential%gth_pot%rcnl(l)
682 0 : WRITE (iw, '(T10,A,T76,I5)') " Nl ", potential%gth_pot%nl(l)
683 0 : WRITE (pline, '(5F12.6)') (potential%gth_pot%hnl(1, j, l), j=1, potential%gth_pot%nl(l))
684 0 : WRITE (iw, '(T10,A,T21,A60)') " Hnl ", ADJUSTR(pline)
685 0 : DO i = 2, potential%gth_pot%nl(l)
686 0 : WRITE (pline, '(T21,5F12.6)') (potential%gth_pot%hnl(i, j, l), j=i, potential%gth_pot%nl(l))
687 0 : WRITE (iw, '(T21,A60)') ADJUSTR(pline)
688 : END DO
689 : END IF
690 : END DO
691 : CASE (upf_pseudo)
692 0 : WRITE (iw, '(/," ",29("*"),A,29("*"))') " UPF Pseudopotential "
693 0 : DO k = 1, potential%upf_pot%maxinfo
694 0 : WRITE (iw, '(A80)') potential%upf_pot%info(k)
695 : END DO
696 : CASE (sgp_pseudo)
697 0 : WRITE (iw, '(/," ",29("*"),A,29("*"))') " SGP Pseudopotential "
698 0 : WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%sgp_pot%zion
699 : CASE (ecp_pseudo)
700 3 : WRITE (iw, '(/," ",26("*"),A,27("*"))') " Effective Core Potential "
701 3 : WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%ecp_pot%zion
702 6 : DO k = 1, potential%ecp_pot%nloc
703 6 : IF (k == 1) THEN
704 3 : WRITE (iw, '(T10,A,T40,I3,T49,2F16.8)') " Local Potential ", potential%ecp_pot%nrloc(k), &
705 6 : potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
706 : ELSE
707 0 : WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrloc(k), &
708 0 : potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
709 : END IF
710 : END DO
711 14 : DO l = 0, potential%ecp_pot%lmax
712 11 : WRITE (iw, '(T10,A,I3)') " ECP l-value ", l
713 29 : DO k = 1, potential%ecp_pot%npot(l)
714 15 : WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrpot(k, l), &
715 41 : potential%ecp_pot%bpot(k, l), potential%ecp_pot%apot(k, l)
716 : END DO
717 : END DO
718 : CASE DEFAULT
719 3 : CPABORT("")
720 : END SELECT
721 3 : IF (potential%confinement) THEN
722 0 : IF (potential%conf_type == poly_conf) THEN
723 : WRITE (iw, '(/,T10,A,T51,F12.6," * (R /",F6.2,")**",F6.2)') &
724 0 : " Confinement Potential ", potential%acon, potential%rcon, potential%scon
725 0 : ELSE IF (potential%conf_type == barrier_conf) THEN
726 0 : WRITE (iw, '(/,T10,A)') " Confinement Potential s*F[(r-ron)/w] "
727 0 : WRITE (iw, '(T57,A,F12.6,A)') "s =", potential%acon, " Ha"
728 0 : WRITE (iw, '(T57,A,F12.6,A)') "w =", potential%rcon, " Bohr"
729 0 : WRITE (iw, '(T57,A,F12.6,A)') "ron =", potential%scon, " Bohr"
730 : ELSE
731 0 : CPABORT("")
732 : END IF
733 : ELSE
734 3 : WRITE (iw, '(/,T10,A)') " No Confinement Potential is applied "
735 : END IF
736 3 : WRITE (iw, '(" ",79("*"))')
737 :
738 3 : END SUBROUTINE atom_print_potential
739 :
740 : ! **************************************************************************************************
741 : !> \brief Print GTH pseudo-potential parameters.
742 : !> \param gthpot pseudo-potential
743 : !> \param iunit output file unit
744 : !> \par History
745 : !> * 09.2012 created [Juerg Hutter]
746 : !> \note The pseudo-potential is written into the 'iunit' file unit or as the file 'GTH-PARAMETER'
747 : !> inside the current working directory if the I/O unit is not given explicitly.
748 : ! **************************************************************************************************
749 37 : SUBROUTINE atom_write_pseudo_param(gthpot, iunit)
750 : TYPE(atom_gthpot_type), INTENT(INOUT) :: gthpot
751 : INTEGER, INTENT(IN), OPTIONAL :: iunit
752 :
753 : INTEGER :: i, iw, j, k, n
754 :
755 37 : IF (PRESENT(iunit)) THEN
756 2 : iw = iunit
757 : ELSE
758 35 : CALL open_file(file_name="GTH-PARAMETER", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
759 : END IF
760 37 : WRITE (iw, '(A2,A)') gthpot%symbol, ADJUSTL(TRIM(gthpot%pname))
761 37 : WRITE (iw, '(4I5)') gthpot%econf(0:3)
762 111 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
763 37 : IF (gthpot%lpotextended) THEN
764 0 : WRITE (iw, '(A,I5)') " LPOT", gthpot%nexp_lpot
765 0 : DO i = 1, gthpot%nexp_lpot
766 0 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
767 0 : (gthpot%cval_lpot(j, i), j=1, gthpot%nct_lpot(i))
768 : END DO
769 : END IF
770 37 : IF (gthpot%lsdpot) THEN
771 0 : WRITE (iw, '(A,I5)') " LSD ", gthpot%nexp_lsd
772 0 : DO i = 1, gthpot%nexp_lsd
773 0 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
774 0 : (gthpot%cval_lsd(j, i), j=1, gthpot%nct_lsd(i))
775 : END DO
776 : END IF
777 37 : IF (gthpot%nlcc) THEN
778 24 : WRITE (iw, '(A,I5)') " NLCC ", gthpot%nexp_nlcc
779 48 : DO i = 1, gthpot%nexp_nlcc
780 24 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
781 96 : (gthpot%cval_nlcc(j, i)*4.0_dp*pi, j=1, gthpot%nct_nlcc(i))
782 : END DO
783 : END IF
784 37 : n = 0
785 201 : DO i = lmat, 0, -1
786 201 : IF (gthpot%nl(i) > 0) THEN
787 34 : n = i + 1
788 34 : EXIT
789 : END IF
790 : END DO
791 37 : WRITE (iw, '(I8)') n
792 95 : DO i = 0, n - 1
793 116 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
794 37 : SELECT CASE (gthpot%nl(i))
795 : CASE (2)
796 0 : WRITE (iw, '(T49,F20.14)') gthpot%hnl(2, 2, i)
797 : CASE (3)
798 0 : WRITE (iw, '(T49,2F20.14)') gthpot%hnl(2, 2, i), gthpot%hnl(2, 3, i)
799 0 : WRITE (iw, '(T69,F20.14)') gthpot%hnl(3, 3, i)
800 : CASE DEFAULT
801 58 : DO j = 2, gthpot%nl(i)
802 58 : WRITE (iw, '(T29,5F20.14)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
803 : END DO
804 : END SELECT
805 : END DO
806 37 : IF (.NOT. PRESENT(iunit)) CALL close_file(unit_number=iw)
807 :
808 37 : END SUBROUTINE atom_write_pseudo_param
809 :
810 : ! **************************************************************************************************
811 : !> \brief Print atomic orbitals.
812 : !> \param atom information about the atomic kind
813 : !> \param iw output file unit
814 : !> \param xmgrace ...
815 : !> \par History
816 : !> * 04.2013 created [Juerg Hutter]
817 : ! **************************************************************************************************
818 0 : SUBROUTINE atom_print_orbitals(atom, iw, xmgrace)
819 : TYPE(atom_type), POINTER :: atom
820 : INTEGER, INTENT(IN) :: iw
821 : LOGICAL, INTENT(IN), OPTIONAL :: xmgrace
822 :
823 : CHARACTER(LEN=40) :: fnbody
824 : INTEGER :: z
825 : LOGICAL :: graph
826 :
827 0 : SELECT CASE (atom%method_type)
828 : CASE DEFAULT
829 0 : CPABORT("")
830 : CASE (do_rks_atom)
831 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
832 : CASE (do_uks_atom)
833 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
834 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
835 : CASE (do_rhf_atom)
836 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
837 : CASE (do_uhf_atom)
838 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
839 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
840 : CASE (do_rohf_atom)
841 0 : CPABORT("")
842 : END SELECT
843 :
844 0 : graph = .FALSE.
845 0 : IF (PRESENT(xmgrace)) graph = xmgrace
846 0 : IF (graph .AND. iw > 0) THEN
847 0 : z = atom%z
848 0 : fnbody = TRIM(ptable(z)%symbol)//"_PPorbital"
849 0 : SELECT CASE (atom%method_type)
850 : CASE DEFAULT
851 0 : CPABORT("")
852 : CASE (do_rks_atom)
853 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
854 : CASE (do_uks_atom)
855 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
856 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
857 : CASE (do_rhf_atom)
858 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
859 : CASE (do_uhf_atom)
860 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
861 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
862 : CASE (do_rohf_atom)
863 0 : CPABORT("")
864 : END SELECT
865 : END IF
866 :
867 0 : END SUBROUTINE atom_print_orbitals
868 :
869 : ! **************************************************************************************************
870 : !> \brief Print atomic orbitals of the given spin.
871 : !> \param atom information about the atomic kind
872 : !> \param wfn atomic orbitals
873 : !> \param description description string
874 : !> \param iw output file unit
875 : !> \par History
876 : !> * 04.2013 created [Juerg Hutter]
877 : ! **************************************************************************************************
878 0 : SUBROUTINE atom_print_orbitals_helper(atom, wfn, description, iw)
879 : TYPE(atom_type), POINTER :: atom
880 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: wfn
881 : CHARACTER(len=*), INTENT(IN) :: description
882 : INTEGER, INTENT(IN) :: iw
883 :
884 : INTEGER :: b, l, maxl, nb, nv, v
885 :
886 0 : WRITE (iw, '(/,A,A,A)') " Atomic orbital expansion coefficients [", description, "]"
887 :
888 0 : maxl = atom%state%maxl_calc
889 0 : DO l = 0, maxl
890 :
891 0 : nb = atom%basis%nbas(l)
892 0 : nv = atom%state%maxn_calc(l)
893 0 : IF (nb > 0 .AND. nv > 0) THEN
894 0 : nv = MIN(nv, SIZE(wfn, 2))
895 0 : DO v = 1, nv
896 0 : WRITE (iw, '(/," ORBITAL L = ",I1," State = ",I3)') l, v
897 0 : DO b = 1, nb
898 0 : WRITE (iw, '(" ",ES23.15)') wfn(b, v, l)
899 : END DO
900 : END DO
901 : END IF
902 : END DO
903 0 : END SUBROUTINE atom_print_orbitals_helper
904 :
905 : ! **************************************************************************************************
906 : !> \brief Print atomic orbitals of the given spin.
907 : !> \param atom information about the atomic kind
908 : !> \param wfn atomic orbitals
909 : !> \param fnbody body of file name
910 : !> \par History
911 : !> * 02.2025 created [Juerg Hutter]
912 : ! **************************************************************************************************
913 0 : SUBROUTINE atom_orbitals_grace(atom, wfn, fnbody)
914 : TYPE(atom_type), POINTER :: atom
915 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: wfn
916 : CHARACTER(len=*), INTENT(IN) :: fnbody
917 :
918 : CHARACTER(LEN=1), DIMENSION(0:8) :: lname
919 : CHARACTER(LEN=1), DIMENSION(1:9) :: wnum
920 : CHARACTER(LEN=40) :: fname, legend
921 : INTEGER :: b, i, iw, l, m, maxl, nb, nv, v
922 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gdata, wfnr
923 : REAL(KIND=dp), DIMENSION(4) :: world_coord
924 :
925 0 : lname = ['s', 'p', 'd', 'f', 'g', 'h', 'j', 'k', 'l']
926 0 : wnum = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
927 0 : m = atom%basis%grid%nr
928 0 : maxl = atom%state%maxl_calc
929 0 : DO l = 0, maxl
930 0 : fname = TRIM(fnbody)//"_"//lname(l)//".agr"
931 0 : nb = atom%basis%nbas(l)
932 0 : nv = atom%state%maxn_calc(l)
933 0 : IF (nb > 0 .AND. nv > 0) THEN
934 0 : CALL open_file(file_name=fname, file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
935 0 : nv = MIN(nv, SIZE(wfn, 2))
936 0 : ALLOCATE (wfnr(m, nv))
937 0 : wfnr = 0.0_dp
938 0 : DO v = 1, nv
939 0 : DO b = 1, nb
940 0 : wfnr(:, v) = wfnr(:, v) + wfn(b, v, l)*atom%basis%bf(:, b, l)
941 : END DO
942 : END DO
943 0 : world_coord(1) = 0.0_dp
944 0 : world_coord(2) = MINVAL(wfnr) - 0.5_dp
945 0 : world_coord(3) = 15.0_dp
946 0 : world_coord(4) = MAXVAL(wfnr) + 0.5_dp
947 : !
948 0 : CALL xm_write_defaults(iw)
949 0 : CALL xm_write_frameport(iw)
950 : CALL xm_write_frame(iw, world_coord, &
951 : title="PP Radial Wavefunction", &
952 : subtitle=lname(l)//"-Quantum Number", &
953 : xlabel="Radius [Bohr]", &
954 0 : ylabel="")
955 0 : DO i = 0, nv - 1
956 0 : legend = "WFN "//wnum(i + 1)
957 0 : CALL xm_graph_info(iw, i, 2.5_dp, legend)
958 : END DO
959 0 : ALLOCATE (gdata(m, 2))
960 0 : gdata(1:m, 1) = atom%basis%grid%rad(1:m)
961 0 : DO i = 0, nv - 1
962 0 : gdata(1:m, 2) = wfnr(1:m, i + 1)
963 0 : CALL xm_graph_data(iw, i, gdata)
964 : END DO
965 0 : DEALLOCATE (gdata, wfnr)
966 0 : CALL close_file(iw)
967 : END IF
968 : END DO
969 0 : END SUBROUTINE atom_orbitals_grace
970 :
971 : END MODULE atom_output
|
