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 181 : 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 181 : ADJUSTL(TRIM(info)), TRIM(ptable(zval)%name), TRIM(ptable(zval)%symbol), zval
75 :
76 181 : 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 2283 : 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 2283 : WRITE (iw, '(/,T2,A)') "Electronic structure"
98 162093 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of core electrons", SUM(state%core)
99 162093 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of valence electrons", SUM(state%occ)
100 162093 : WRITE (iw, '(T5,A,T71,F10.2)') "Total number of electrons", SUM(state%occ + state%core)
101 4520 : SELECT CASE (state%multiplicity)
102 : CASE (-1)
103 2237 : 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 11 : 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 2283 : WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "septet"
122 : CASE DEFAULT
123 : END SELECT
124 :
125 2283 : mlo = get_maxl_occ(state%occ)
126 2283 : mlc = get_maxl_occ(state%core)
127 2283 : mm = get_maxn_occ(state%core)
128 :
129 2283 : IF (state%multiplicity == -1) THEN
130 5789 : DO l = 0, MAX(mlo, mlc)
131 3552 : mo = state%maxn_occ(l)
132 41309 : IF (SUM(state%core(l, :)) == 0) THEN
133 5818 : WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occ(l, j), j=1, mo)
134 : ELSE
135 1107 : mc = mm(l)
136 2504 : CPASSERT(SUM(state%occ(l, 1:mc)) == 0)
137 2504 : WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (state%core(l, j), j=1, mc)
138 2189 : 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 46 : WRITE (iw, '(T5,A)') "Alpha Electrons"
143 124 : DO l = 0, MAX(mlo, mlc)
144 78 : mo = state%maxn_occ(l)
145 904 : IF (SUM(state%core(l, :)) == 0) THEN
146 106 : 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 46 : WRITE (iw, '(T5,A)') "Beta Electrons"
154 124 : DO l = 0, MAX(mlo, mlc)
155 78 : mo = state%maxn_occ(l)
156 904 : IF (SUM(state%core(l, :)) == 0) THEN
157 106 : 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 2283 : WRITE (iw, *)
166 :
167 2283 : 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 4566 : 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 2283 : WRITE (iw, '(/,A,T36,A,T61,F20.12)') " Energy components [Hartree]", &
187 4566 : " Total Energy ::", atom%energy%etot
188 2283 : WRITE (iw, '(T36,A,T61,F20.12)') " Band Energy ::", atom%energy%eband
189 2283 : WRITE (iw, '(T36,A,T61,F20.12)') " Kinetic Energy ::", atom%energy%ekin
190 2283 : WRITE (iw, '(T36,A,T61,F20.12)') "Potential Energy ::", atom%energy%epot
191 2283 : IF (atom%energy%ekin /= 0.0_dp) THEN
192 2270 : WRITE (iw, '(T36,A,T61,F20.12)') " Virial (-V/T) ::", -atom%energy%epot/atom%energy%ekin
193 : END IF
194 2283 : WRITE (iw, '(T36,A,T61,F20.12)') " Core Energy ::", atom%energy%ecore
195 2283 : IF (atom%energy%exc /= 0._dp) &
196 2252 : WRITE (iw, '(T36,A,T61,F20.12)') " XC Energy ::", atom%energy%exc
197 2283 : WRITE (iw, '(T36,A,T61,F20.12)') " Coulomb Energy ::", atom%energy%ecoulomb
198 2283 : IF (atom%energy%eexchange /= 0._dp) &
199 46 : WRITE (iw, '(T34,A,T61,F20.12)') "HF Exchange Energy ::", atom%energy%eexchange
200 2283 : IF (atom%potential%ppot_type /= NO_PSEUDO) THEN
201 1881 : WRITE (iw, '(T20,A,T61,F20.12)') " Total Pseudopotential Energy ::", atom%energy%epseudo
202 1881 : WRITE (iw, '(T20,A,T61,F20.12)') " Local Pseudopotential Energy ::", atom%energy%eploc
203 1881 : IF (atom%energy%elsd /= 0._dp) &
204 0 : WRITE (iw, '(T20,A,T61,F20.12)') " Local Spin-potential Energy ::", atom%energy%elsd
205 1881 : WRITE (iw, '(T20,A,T61,F20.12)') " Nonlocal Pseudopotential Energy ::", atom%energy%epnl
206 : END IF
207 2283 : IF (atom%potential%confinement) THEN
208 1874 : WRITE (iw, '(T36,A,T61,F20.12)') " Confinement ::", atom%energy%econfinement
209 : END IF
210 :
211 2283 : IF (atom%state%multiplicity == -1) THEN
212 2237 : WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T49,A,T71,A,/)') " Orbital energies", &
213 4474 : "State", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
214 5899 : DO l = 0, atom%state%maxl_calc
215 3662 : n = atom%state%maxn_calc(l)
216 8232 : DO i = 1, n
217 : WRITE (iw, '(T23,I2,T30,I1,T36,F10.3,T46,F15.6,T66,F15.6)') &
218 8232 : i, l, atom%state%occupation(l, i), atom%orbitals%ener(i, l), atom%orbitals%ener(i, l)*evolt
219 : END DO
220 5899 : IF (n > 0) WRITE (iw, *)
221 : END DO
222 : ELSE
223 46 : WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T42,A,T55,A,T71,A,/)') " Orbital energies", &
224 92 : "State", "Spin", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
225 164 : DO l = 0, atom%state%maxl_calc
226 118 : n = atom%state%maxn_calc(l)
227 210 : DO i = 1, n
228 : WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
229 210 : i, "alpha", l, atom%state%occa(l, i), atom%orbitals%enera(i, l), atom%orbitals%enera(i, l)*evolt
230 : END DO
231 210 : DO i = 1, n
232 : WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
233 210 : i, " beta", l, atom%state%occb(l, i), atom%orbitals%enerb(i, l), atom%orbitals%enerb(i, l)*evolt
234 : END DO
235 164 : IF (n > 0) WRITE (iw, *)
236 : END DO
237 : END IF
238 :
239 2283 : CALL get_rho0(atom, drho)
240 2283 : WRITE (iw, '(/,A,T66,F15.6)') " Total Electron Density at R=0: ", drho
241 :
242 2283 : 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 10142 : 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 10142 : IF (iter == 1) THEN
282 : WRITE (iw, '(/," ",79("*"),/,T19,A,T38,A,T70,A,/," ",79("*"))') &
283 2283 : "Iteration", "Convergence", "Energy [au]"
284 : END IF
285 10142 : WRITE (iw, '(T20,i8,T34,G14.6,T61,F20.12)') iter, deps, etot
286 :
287 10142 : 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 368 : 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 368 : 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 368 : NULLIFY (xc_fun, xc_fun_section, xc_section)
520 :
521 368 : meth = atom%method_type
522 :
523 368 : xc_section => atom%xc_section
524 368 : 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 330 : 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 368 : myfun = xc_none
538 : END SELECT
539 :
540 0 : SELECT CASE (meth)
541 : CASE DEFAULT
542 0 : CPABORT("")
543 : CASE (do_rks_atom)
544 296 : 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 330 : IF (iw > 0) WRITE (iw, fmt="(/,' METHOD | Restricted Open-Shell Kohn-Sham Calculation')")
553 : END SELECT
554 :
555 : ! zmp
556 368 : 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 368 : 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 368 : IF (atom%pp_calc) THEN
574 80 : 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 368 : lsd = (meth == do_uks_atom)
601 :
602 368 : IF (myfun /= xc_none) THEN
603 328 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", c_val=tmpStr)
604 328 : IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| ROUTINE=',a)") TRIM(tmpStr)
605 328 : CALL xc_functionals_expand(xc_fun_section, xc_section)
606 328 : IF (iw > 0) THEN
607 164 : ifun = 0
608 233 : DO
609 397 : ifun = ifun + 1
610 397 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
611 397 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
612 233 : 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 230 : ALLOCATE (CHARACTER(LEN=20*default_string_length) :: reference)
616 : END IF
617 233 : CALL xc_functional_get_info(xc_fun, lsd=lsd, reference=reference, shortform=shortform)
618 : WRITE (iw, fmt="(' FUNCTIONAL| ',a,':')") &
619 233 : TRIM(xc_fun%section%name)
620 633 : DO il = 1, LEN_TRIM(reference), 67
621 633 : WRITE (iw, fmt="(' FUNCTIONAL| ',a67)") reference(il:)
622 : END DO
623 397 : 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 368 : 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 0 : IF (potential%gth_pot%soc) THEN
692 0 : WRITE (iw, '(T10,A)') " Spin-orbit coupling parameters "
693 0 : DO l = 1, lmat
694 0 : IF (potential%gth_pot%nl(l) > 0) THEN
695 0 : WRITE (iw, '(T10,A,T76,I5)') " Angular momentum ", l
696 0 : WRITE (iw, '(T10,A,T66,F15.6)') " Rcnl ", potential%gth_pot%rcnl(l)
697 0 : WRITE (iw, '(T10,A,T76,I5)') " Nl ", potential%gth_pot%nl(l)
698 0 : WRITE (pline, '(5F12.6)') (potential%gth_pot%knl(1, j, l), j=1, potential%gth_pot%nl(l))
699 0 : WRITE (iw, '(T10,A,T21,A60)') " Hnl ", ADJUSTR(pline)
700 0 : DO i = 2, potential%gth_pot%nl(l)
701 0 : WRITE (pline, '(T21,5F12.6)') (potential%gth_pot%knl(i, j, l), j=i, potential%gth_pot%nl(l))
702 0 : WRITE (iw, '(T21,A60)') ADJUSTR(pline)
703 : END DO
704 : END IF
705 : END DO
706 : END IF
707 : CASE (upf_pseudo)
708 0 : WRITE (iw, '(/," ",29("*"),A,29("*"))') " UPF Pseudopotential "
709 0 : DO k = 1, potential%upf_pot%maxinfo
710 0 : WRITE (iw, '(A80)') potential%upf_pot%info(k)
711 : END DO
712 : CASE (sgp_pseudo)
713 0 : WRITE (iw, '(/," ",29("*"),A,29("*"))') " SGP Pseudopotential "
714 0 : WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%sgp_pot%zion
715 : CASE (ecp_pseudo)
716 3 : WRITE (iw, '(/," ",26("*"),A,27("*"))') " Effective Core Potential "
717 3 : WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%ecp_pot%zion
718 6 : DO k = 1, potential%ecp_pot%nloc
719 6 : IF (k == 1) THEN
720 3 : WRITE (iw, '(T10,A,T40,I3,T49,2F16.8)') " Local Potential ", potential%ecp_pot%nrloc(k), &
721 6 : potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
722 : ELSE
723 0 : WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrloc(k), &
724 0 : potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
725 : END IF
726 : END DO
727 14 : DO l = 0, potential%ecp_pot%lmax
728 11 : WRITE (iw, '(T10,A,I3)') " ECP l-value ", l
729 29 : DO k = 1, potential%ecp_pot%npot(l)
730 15 : WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrpot(k, l), &
731 41 : potential%ecp_pot%bpot(k, l), potential%ecp_pot%apot(k, l)
732 : END DO
733 : END DO
734 : CASE DEFAULT
735 3 : CPABORT("")
736 : END SELECT
737 3 : IF (potential%confinement) THEN
738 0 : IF (potential%conf_type == poly_conf) THEN
739 : WRITE (iw, '(/,T10,A,T51,F12.6," * (R /",F6.2,")**",F6.2)') &
740 0 : " Confinement Potential ", potential%acon, potential%rcon, potential%scon
741 0 : ELSE IF (potential%conf_type == barrier_conf) THEN
742 0 : WRITE (iw, '(/,T10,A)') " Confinement Potential s*F[(r-ron)/w] "
743 0 : WRITE (iw, '(T57,A,F12.6,A)') "s =", potential%acon, " Ha"
744 0 : WRITE (iw, '(T57,A,F12.6,A)') "w =", potential%rcon, " Bohr"
745 0 : WRITE (iw, '(T57,A,F12.6,A)') "ron =", potential%scon, " Bohr"
746 : ELSE
747 0 : CPABORT("")
748 : END IF
749 : ELSE
750 3 : WRITE (iw, '(/,T10,A)') " No Confinement Potential is applied "
751 : END IF
752 3 : WRITE (iw, '(" ",79("*"))')
753 :
754 3 : END SUBROUTINE atom_print_potential
755 :
756 : ! **************************************************************************************************
757 : !> \brief Print GTH pseudo-potential parameters.
758 : !> \param gthpot pseudo-potential
759 : !> \param iunit output file unit
760 : !> \par History
761 : !> * 09.2012 created [Juerg Hutter]
762 : !> \note The pseudo-potential is written into the 'iunit' file unit or as the file 'GTH-PARAMETER'
763 : !> inside the current working directory if the I/O unit is not given explicitly.
764 : ! **************************************************************************************************
765 38 : SUBROUTINE atom_write_pseudo_param(gthpot, iunit)
766 : TYPE(atom_gthpot_type), INTENT(INOUT) :: gthpot
767 : INTEGER, INTENT(IN), OPTIONAL :: iunit
768 :
769 : INTEGER :: i, iw, j, k, n
770 :
771 38 : IF (PRESENT(iunit)) THEN
772 3 : iw = iunit
773 : ELSE
774 35 : CALL open_file(file_name="GTH-PARAMETER", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
775 : END IF
776 38 : WRITE (iw, '(A)') TRIM(ADJUSTL(gthpot%symbol))//" "//TRIM(ADJUSTL(gthpot%pname))
777 38 : WRITE (iw, '(4I5)') gthpot%econf(0:3)
778 114 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
779 38 : IF (gthpot%lpotextended) THEN
780 0 : WRITE (iw, '(A,I5)') " LPOT", gthpot%nexp_lpot
781 0 : DO i = 1, gthpot%nexp_lpot
782 0 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
783 0 : (gthpot%cval_lpot(j, i), j=1, gthpot%nct_lpot(i))
784 : END DO
785 : END IF
786 38 : IF (gthpot%lsdpot) THEN
787 0 : WRITE (iw, '(A,I5)') " LSD ", gthpot%nexp_lsd
788 0 : DO i = 1, gthpot%nexp_lsd
789 0 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
790 0 : (gthpot%cval_lsd(j, i), j=1, gthpot%nct_lsd(i))
791 : END DO
792 : END IF
793 38 : IF (gthpot%nlcc) THEN
794 24 : WRITE (iw, '(A,I5)') " NLCC ", gthpot%nexp_nlcc
795 48 : DO i = 1, gthpot%nexp_nlcc
796 24 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
797 96 : (gthpot%cval_nlcc(j, i)*4.0_dp*pi, j=1, gthpot%nct_nlcc(i))
798 : END DO
799 : END IF
800 38 : n = 0
801 205 : DO i = lmat, 0, -1
802 205 : IF (gthpot%nl(i) > 0) THEN
803 35 : n = i + 1
804 35 : EXIT
805 : END IF
806 : END DO
807 38 : WRITE (iw, '(I8)') n
808 99 : DO i = 0, n - 1
809 127 : WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
810 38 : SELECT CASE (gthpot%nl(i))
811 : CASE (2)
812 1 : WRITE (iw, '(T49,F20.14)') gthpot%hnl(2, 2, i)
813 : CASE (3)
814 2 : WRITE (iw, '(T49,2F20.14)') gthpot%hnl(2, 2, i), gthpot%hnl(2, 3, i)
815 2 : WRITE (iw, '(T69,F20.14)') gthpot%hnl(3, 3, i)
816 : CASE DEFAULT
817 119 : DO j = 2, gthpot%nl(i)
818 58 : WRITE (iw, '(T29,5F20.14)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
819 : END DO
820 : END SELECT
821 : END DO
822 38 : IF (gthpot%soc) THEN
823 3 : DO i = 1, n - 1
824 7 : WRITE (iw, '(T29,5F20.14)') (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
825 1 : SELECT CASE (gthpot%nl(i))
826 : CASE (2)
827 1 : WRITE (iw, '(T49,F20.14)') gthpot%knl(2, 2, i)
828 : CASE (3)
829 1 : WRITE (iw, '(T49,2F20.14)') gthpot%knl(2, 2, i), gthpot%knl(2, 3, i)
830 1 : WRITE (iw, '(T69,F20.14)') gthpot%knl(3, 3, i)
831 : CASE DEFAULT
832 2 : DO j = 2, gthpot%nl(i)
833 0 : WRITE (iw, '(T29,5F20.14)') (gthpot%knl(j, k, i), k=j, gthpot%nl(i))
834 : END DO
835 : END SELECT
836 : END DO
837 : END IF
838 38 : IF (.NOT. PRESENT(iunit)) CALL close_file(unit_number=iw)
839 :
840 38 : END SUBROUTINE atom_write_pseudo_param
841 :
842 : ! **************************************************************************************************
843 : !> \brief Print atomic orbitals.
844 : !> \param atom information about the atomic kind
845 : !> \param iw output file unit
846 : !> \param xmgrace ...
847 : !> \par History
848 : !> * 04.2013 created [Juerg Hutter]
849 : ! **************************************************************************************************
850 0 : SUBROUTINE atom_print_orbitals(atom, iw, xmgrace)
851 : TYPE(atom_type), POINTER :: atom
852 : INTEGER, INTENT(IN) :: iw
853 : LOGICAL, INTENT(IN), OPTIONAL :: xmgrace
854 :
855 : CHARACTER(LEN=40) :: fnbody
856 : INTEGER :: z
857 : LOGICAL :: graph
858 :
859 0 : SELECT CASE (atom%method_type)
860 : CASE DEFAULT
861 0 : CPABORT("")
862 : CASE (do_rks_atom)
863 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
864 : CASE (do_uks_atom)
865 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
866 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
867 : CASE (do_rhf_atom)
868 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
869 : CASE (do_uhf_atom)
870 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
871 0 : CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
872 : CASE (do_rohf_atom)
873 0 : CPABORT("")
874 : END SELECT
875 :
876 0 : graph = .FALSE.
877 0 : IF (PRESENT(xmgrace)) graph = xmgrace
878 0 : IF (graph .AND. iw > 0) THEN
879 0 : z = atom%z
880 0 : fnbody = TRIM(ptable(z)%symbol)//"_PPorbital"
881 0 : SELECT CASE (atom%method_type)
882 : CASE DEFAULT
883 0 : CPABORT("")
884 : CASE (do_rks_atom)
885 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
886 : CASE (do_uks_atom)
887 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
888 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
889 : CASE (do_rhf_atom)
890 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
891 : CASE (do_uhf_atom)
892 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
893 0 : CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
894 : CASE (do_rohf_atom)
895 0 : CPABORT("")
896 : END SELECT
897 : END IF
898 :
899 0 : END SUBROUTINE atom_print_orbitals
900 :
901 : ! **************************************************************************************************
902 : !> \brief Print atomic orbitals of the given spin.
903 : !> \param atom information about the atomic kind
904 : !> \param wfn atomic orbitals
905 : !> \param description description string
906 : !> \param iw output file unit
907 : !> \par History
908 : !> * 04.2013 created [Juerg Hutter]
909 : ! **************************************************************************************************
910 0 : SUBROUTINE atom_print_orbitals_helper(atom, wfn, description, iw)
911 : TYPE(atom_type), POINTER :: atom
912 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: wfn
913 : CHARACTER(len=*), INTENT(IN) :: description
914 : INTEGER, INTENT(IN) :: iw
915 :
916 : INTEGER :: b, l, maxl, nb, nv, v
917 :
918 0 : WRITE (iw, '(/,A,A,A)') " Atomic orbital expansion coefficients [", description, "]"
919 :
920 0 : maxl = atom%state%maxl_calc
921 0 : DO l = 0, maxl
922 :
923 0 : nb = atom%basis%nbas(l)
924 0 : nv = atom%state%maxn_calc(l)
925 0 : IF (nb > 0 .AND. nv > 0) THEN
926 0 : nv = MIN(nv, SIZE(wfn, 2))
927 0 : DO v = 1, nv
928 0 : WRITE (iw, '(/," ORBITAL L = ",I1," State = ",I3)') l, v
929 0 : DO b = 1, nb
930 0 : WRITE (iw, '(" ",ES23.15)') wfn(b, v, l)
931 : END DO
932 : END DO
933 : END IF
934 : END DO
935 0 : END SUBROUTINE atom_print_orbitals_helper
936 :
937 : ! **************************************************************************************************
938 : !> \brief Print atomic orbitals of the given spin.
939 : !> \param atom information about the atomic kind
940 : !> \param wfn atomic orbitals
941 : !> \param fnbody body of file name
942 : !> \par History
943 : !> * 02.2025 created [Juerg Hutter]
944 : ! **************************************************************************************************
945 0 : SUBROUTINE atom_orbitals_grace(atom, wfn, fnbody)
946 : TYPE(atom_type), POINTER :: atom
947 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: wfn
948 : CHARACTER(len=*), INTENT(IN) :: fnbody
949 :
950 : CHARACTER(LEN=1), DIMENSION(0:8) :: lname
951 : CHARACTER(LEN=1), DIMENSION(1:9) :: wnum
952 : CHARACTER(LEN=40) :: fname, legend
953 : INTEGER :: b, i, iw, l, m, maxl, nb, nv, v
954 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gdata, wfnr
955 : REAL(KIND=dp), DIMENSION(4) :: world_coord
956 :
957 0 : lname = ['s', 'p', 'd', 'f', 'g', 'h', 'j', 'k', 'l']
958 0 : wnum = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
959 0 : m = atom%basis%grid%nr
960 0 : maxl = atom%state%maxl_calc
961 0 : DO l = 0, maxl
962 0 : fname = TRIM(fnbody)//"_"//lname(l)//".agr"
963 0 : nb = atom%basis%nbas(l)
964 0 : nv = atom%state%maxn_calc(l)
965 0 : IF (nb > 0 .AND. nv > 0) THEN
966 0 : CALL open_file(file_name=fname, file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
967 0 : nv = MIN(nv, SIZE(wfn, 2))
968 0 : ALLOCATE (wfnr(m, nv))
969 0 : wfnr = 0.0_dp
970 0 : DO v = 1, nv
971 0 : DO b = 1, nb
972 0 : wfnr(:, v) = wfnr(:, v) + wfn(b, v, l)*atom%basis%bf(:, b, l)
973 : END DO
974 : END DO
975 0 : world_coord(1) = 0.0_dp
976 0 : world_coord(2) = MINVAL(wfnr) - 0.5_dp
977 0 : world_coord(3) = 15.0_dp
978 0 : world_coord(4) = MAXVAL(wfnr) + 0.5_dp
979 : !
980 0 : CALL xm_write_defaults(iw)
981 0 : CALL xm_write_frameport(iw)
982 : CALL xm_write_frame(iw, world_coord, &
983 : title="PP Radial Wavefunction", &
984 : subtitle=lname(l)//"-Quantum Number", &
985 : xlabel="Radius [Bohr]", &
986 0 : ylabel="")
987 0 : DO i = 0, nv - 1
988 0 : legend = "WFN "//wnum(i + 1)
989 0 : CALL xm_graph_info(iw, i, 2.5_dp, legend)
990 : END DO
991 0 : ALLOCATE (gdata(m, 2))
992 0 : gdata(1:m, 1) = atom%basis%grid%rad(1:m)
993 0 : DO i = 0, nv - 1
994 0 : gdata(1:m, 2) = wfnr(1:m, i + 1)
995 0 : CALL xm_graph_data(iw, i, gdata)
996 : END DO
997 0 : DEALLOCATE (gdata, wfnr)
998 0 : CALL close_file(iw)
999 : END IF
1000 : END DO
1001 0 : END SUBROUTINE atom_orbitals_grace
1002 :
1003 : END MODULE atom_output
|
