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