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 : MODULE atom_energy
9 : USE atom_admm_methods, ONLY: atom_admm
10 : USE atom_electronic_structure, ONLY: calculate_atom
11 : USE atom_fit, ONLY: atom_fit_density,&
12 : atom_fit_kgpot
13 : USE atom_grb, ONLY: atom_grb_construction
14 : USE atom_operators, ONLY: atom_int_release,&
15 : atom_int_setup,&
16 : atom_ppint_release,&
17 : atom_ppint_setup,&
18 : atom_relint_release,&
19 : atom_relint_setup
20 : USE atom_output, ONLY: atom_print_basis,&
21 : atom_print_info,&
22 : atom_print_method,&
23 : atom_print_orbitals,&
24 : atom_print_potential,&
25 : atom_write_pseudo_param
26 : USE atom_sgp, ONLY: atom_sgp_construction
27 : USE atom_types, ONLY: &
28 : atom_basis_type, atom_gthpot_type, atom_integrals, atom_optimization_type, atom_orbitals, &
29 : atom_p_type, atom_potential_type, atom_state, atom_type, create_atom_orbs, &
30 : create_atom_type, gth_pseudo, init_atom_basis, init_atom_potential, lmat, &
31 : read_atom_opt_section, release_atom_basis, release_atom_potential, release_atom_type, &
32 : set_atom
33 : USE atom_utils, ONLY: &
34 : atom_completeness, atom_condnumber, atom_consistent_method, atom_core_density, &
35 : atom_density, atom_local_potential, atom_read_external_density, atom_read_external_vxc, &
36 : atom_set_occupation, atom_write_zmp_restart, get_maxl_occ, get_maxn_occ
37 : USE atom_xc, ONLY: calculate_atom_ext_vxc,&
38 : calculate_atom_zmp
39 : USE cp_log_handling, ONLY: cp_get_default_logger,&
40 : cp_logger_type
41 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
42 : cp_print_key_unit_nr
43 : USE input_constants, ONLY: do_analytic
44 : USE input_section_types, ONLY: section_vals_get,&
45 : section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: default_string_length,&
49 : dp
50 : USE mathconstants, ONLY: dfac,&
51 : gamma1,&
52 : pi,&
53 : rootpi
54 : USE periodic_table, ONLY: nelem,&
55 : ptable
56 : USE physcon, ONLY: bohr
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 : PRIVATE
61 : PUBLIC :: atom_energy_opt
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_energy'
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Compute the atomic energy.
69 : !> \param atom_section ATOM input section
70 : !> \par History
71 : !> * 11.2016 geometrical response basis set [Juerg Hutter]
72 : !> * 07.2016 ADMM analysis [Juerg Hutter]
73 : !> * 02.2014 non-additive kinetic energy term [Juerg Hutter]
74 : !> * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
75 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
76 : !> * 05.2009 electronic density fit [Juerg Hutter]
77 : !> * 04.2009 refactored and renamed to atom_energy_opt() [Juerg Hutter]
78 : !> * 08.2008 created as atom_energy() [Juerg Hutter]
79 : ! **************************************************************************************************
80 324 : SUBROUTINE atom_energy_opt(atom_section)
81 : TYPE(section_vals_type), POINTER :: atom_section
82 :
83 : CHARACTER(len=*), PARAMETER :: routineN = 'atom_energy_opt'
84 :
85 : CHARACTER(LEN=2) :: elem
86 : CHARACTER(LEN=default_string_length) :: filename
87 : CHARACTER(LEN=default_string_length), &
88 324 : DIMENSION(:), POINTER :: tmpstringlist
89 : INTEGER :: do_eric, do_erie, handle, i, im, in, iw, k, maxl, mb, method, mo, n_meth, n_rep, &
90 : nder, nr_gh, num_gau, num_gto, num_pol, reltyp, zcore, zval, zz
91 : INTEGER, DIMENSION(0:lmat) :: maxn
92 324 : INTEGER, DIMENSION(:), POINTER :: cn
93 : LOGICAL :: dm, do_gh, do_zmp, doread, eri_c, eri_e, &
94 : graph, had_ae, had_pp, lcomp, lcond, &
95 : pp_calc, read_vxc
96 : REAL(KIND=dp) :: crad, delta, lambda
97 324 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ext_density, ext_vxc
98 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: pocc
99 : TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
100 : TYPE(atom_integrals), POINTER :: ae_int, pp_int
101 : TYPE(atom_optimization_type) :: optimization
102 : TYPE(atom_orbitals), POINTER :: orbitals
103 324 : TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
104 : TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
105 : TYPE(atom_state), POINTER :: state
106 : TYPE(cp_logger_type), POINTER :: logger
107 : TYPE(section_vals_type), POINTER :: admm_section, basis_section, external_vxc_section, &
108 : method_section, opt_section, potential_section, powell_section, sgp_section, xc_section, &
109 : zmp_restart_section, zmp_section
110 :
111 324 : CALL timeset(routineN, handle)
112 :
113 : ! What atom do we calculate
114 324 : CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
115 324 : CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
116 324 : zz = 0
117 11574 : DO i = 1, nelem
118 11574 : IF (ptable(i)%symbol == elem) THEN
119 : zz = i
120 : EXIT
121 : END IF
122 : END DO
123 324 : IF (zz /= 1) zval = zz
124 :
125 : ! read and set up inofrmation on the basis sets
126 11988 : ALLOCATE (ae_basis, pp_basis)
127 324 : basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
128 324 : NULLIFY (ae_basis%grid)
129 324 : CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
130 324 : NULLIFY (pp_basis%grid)
131 324 : basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
132 324 : CALL init_atom_basis(pp_basis, basis_section, zval, "PP")
133 :
134 : ! print general and basis set information
135 324 : logger => cp_get_default_logger()
136 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
137 324 : IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw)
138 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
139 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
140 324 : IF (iw > 0) THEN
141 0 : CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
142 0 : CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
143 : END IF
144 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
145 :
146 : ! read and setup information on the pseudopotential
147 324 : NULLIFY (potential_section)
148 324 : potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
149 3403620 : ALLOCATE (ae_pot, p_pot)
150 324 : CALL init_atom_potential(p_pot, potential_section, zval)
151 324 : CALL init_atom_potential(ae_pot, potential_section, -1)
152 :
153 : ! if the ERI's are calculated analytically, we have to precalculate them
154 324 : eri_c = .FALSE.
155 324 : CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
156 324 : IF (do_eric == do_analytic) eri_c = .TRUE.
157 324 : eri_e = .FALSE.
158 324 : CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
159 324 : IF (do_erie == do_analytic) eri_e = .TRUE.
160 324 : CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
161 324 : CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
162 :
163 : ! information on the states to be calculated
164 324 : CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
165 324 : maxn = 0
166 324 : CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
167 684 : DO in = 1, MIN(SIZE(cn), 4)
168 684 : maxn(in - 1) = cn(in)
169 : END DO
170 2268 : DO in = 0, lmat
171 1944 : maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
172 2268 : maxn(in) = MIN(maxn(in), pp_basis%nbas(in))
173 : END DO
174 :
175 : ! read optimization section
176 324 : opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
177 324 : CALL read_atom_opt_section(optimization, opt_section)
178 :
179 324 : had_ae = .FALSE.
180 324 : had_pp = .FALSE.
181 :
182 : ! Check for the total number of electron configurations to be calculated
183 324 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
184 : ! Check for the total number of method types to be calculated
185 324 : method_section => section_vals_get_subs_vals(atom_section, "METHOD")
186 324 : CALL section_vals_get(method_section, n_repetition=n_meth)
187 :
188 : ! integrals
189 137700 : ALLOCATE (ae_int, pp_int)
190 :
191 1950 : ALLOCATE (atom_info(n_rep, n_meth))
192 :
193 654 : DO in = 1, n_rep
194 984 : DO im = 1, n_meth
195 :
196 330 : NULLIFY (atom_info(in, im)%atom)
197 330 : CALL create_atom_type(atom_info(in, im)%atom)
198 :
199 330 : atom_info(in, im)%atom%optimization = optimization
200 :
201 330 : atom_info(in, im)%atom%z = zval
202 330 : xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
203 330 : atom_info(in, im)%atom%xc_section => xc_section
204 :
205 : ! ZMP Reading input sections if they are found initialize everything
206 : do_zmp = .FALSE.
207 : doread = .FALSE.
208 : read_vxc = .FALSE.
209 :
210 330 : zmp_section => section_vals_get_subs_vals(method_section, "ZMP")
211 330 : CALL section_vals_get(zmp_section, explicit=do_zmp)
212 330 : atom_info(in, im)%atom%do_zmp = do_zmp
213 330 : CALL section_vals_val_get(zmp_section, "FILE_DENSITY", c_val=filename)
214 330 : atom_info(in, im)%atom%ext_file = filename
215 : CALL section_vals_val_get(zmp_section, "GRID_TOL", &
216 330 : r_val=atom_info(in, im)%atom%zmpgrid_tol)
217 330 : CALL section_vals_val_get(zmp_section, "LAMBDA", r_val=lambda)
218 330 : atom_info(in, im)%atom%lambda = lambda
219 330 : CALL section_vals_val_get(zmp_section, "DM", l_val=dm)
220 330 : atom_info(in, im)%atom%dm = dm
221 :
222 330 : zmp_restart_section => section_vals_get_subs_vals(zmp_section, "RESTART")
223 330 : CALL section_vals_get(zmp_restart_section, explicit=doread)
224 330 : atom_info(in, im)%atom%doread = doread
225 330 : CALL section_vals_val_get(zmp_restart_section, "FILE_RESTART", c_val=filename)
226 330 : atom_info(in, im)%atom%zmp_restart_file = filename
227 :
228 : ! ZMP Reading external vxc section, if found initialize
229 330 : external_vxc_section => section_vals_get_subs_vals(method_section, "EXTERNAL_VXC")
230 330 : CALL section_vals_get(external_vxc_section, explicit=read_vxc)
231 330 : atom_info(in, im)%atom%read_vxc = read_vxc
232 330 : CALL section_vals_val_get(external_vxc_section, "FILE_VXC", c_val=filename)
233 330 : atom_info(in, im)%atom%ext_vxc_file = filename
234 : CALL section_vals_val_get(external_vxc_section, "GRID_TOL", &
235 330 : r_val=atom_info(in, im)%atom%zmpvxcgrid_tol)
236 :
237 119790 : ALLOCATE (state)
238 :
239 : ! get the electronic configuration
240 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
241 330 : c_vals=tmpstringlist)
242 :
243 : ! set occupations
244 330 : CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
245 330 : state%maxl_occ = get_maxl_occ(state%occ)
246 2310 : state%maxn_occ = get_maxn_occ(state%occ)
247 :
248 : ! set number of states to be calculated
249 330 : state%maxl_calc = MAX(maxl, state%maxl_occ)
250 330 : state%maxl_calc = MIN(lmat, state%maxl_calc)
251 2310 : state%maxn_calc = 0
252 1524 : DO k = 0, state%maxl_calc
253 1524 : state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k))
254 : END DO
255 :
256 : ! is there a pseudo potential
257 1104 : pp_calc = ANY(INDEX(tmpstringlist(1:), "CORE") /= 0)
258 330 : IF (pp_calc) THEN
259 : ! get and set the core occupations
260 72 : CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
261 72 : CALL atom_set_occupation(tmpstringlist, state%core, pocc)
262 5112 : zcore = zval - NINT(SUM(state%core))
263 72 : CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
264 : ELSE
265 18318 : state%core = 0._dp
266 258 : CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
267 : END IF
268 :
269 330 : CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
270 330 : CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
271 330 : CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)
272 :
273 330 : IF (atom_consistent_method(method, state%multiplicity)) THEN
274 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
275 330 : CALL atom_print_method(atom_info(in, im)%atom, iw)
276 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
277 :
278 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
279 330 : IF (pp_calc) THEN
280 72 : IF (iw > 0) CALL atom_print_potential(p_pot, iw)
281 : ELSE
282 258 : IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
283 : END IF
284 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
285 : END IF
286 :
287 : ! calculate integrals
288 330 : IF (pp_calc) THEN
289 : ! general integrals
290 : CALL atom_int_setup(pp_int, pp_basis, &
291 72 : potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e)
292 : ! potential
293 72 : CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot)
294 : !
295 72 : NULLIFY (pp_int%tzora, pp_int%hdkh)
296 : !
297 72 : CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot)
298 936 : state%maxn_calc(:) = MIN(state%maxn_calc(:), pp_basis%nbas(:))
299 504 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
300 : had_pp = .TRUE.
301 : ELSE
302 : ! general integrals
303 : CALL atom_int_setup(ae_int, ae_basis, potential=ae_pot, &
304 258 : eri_coulomb=eri_c, eri_exchange=eri_e)
305 : ! potential
306 258 : CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
307 : ! relativistic correction terms
308 258 : CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
309 : !
310 258 : CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
311 3354 : state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
312 1806 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
313 : had_ae = .TRUE.
314 : END IF
315 :
316 330 : CALL set_atom(atom_info(in, im)%atom, state=state)
317 :
318 : CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
319 330 : exchange_integral_type=do_erie)
320 330 : atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
321 330 : atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
322 :
323 330 : NULLIFY (orbitals)
324 2310 : mo = MAXVAL(state%maxn_calc)
325 2310 : mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
326 330 : CALL create_atom_orbs(orbitals, mb, mo)
327 330 : CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
328 :
329 2310 : IF (atom_consistent_method(method, state%multiplicity)) THEN
330 : !Calculate the electronic structure
331 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
332 330 : CALL calculate_atom(atom_info(in, im)%atom, iw)
333 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
334 :
335 : ! ZMP If we have the external density do zmp
336 330 : IF (atom_info(in, im)%atom%do_zmp) THEN
337 0 : CALL atom_write_zmp_restart(atom_info(in, im)%atom)
338 :
339 0 : ALLOCATE (ext_density(atom_info(in, im)%atom%basis%grid%nr))
340 0 : ext_density = 0._dp
341 0 : CALL atom_read_external_density(ext_density, atom_info(in, im)%atom, iw)
342 : CALL calculate_atom_zmp(ext_density=ext_density, &
343 : atom=atom_info(in, im)%atom, &
344 0 : lprint=.TRUE.)
345 0 : DEALLOCATE (ext_density)
346 : END IF
347 : ! ZMP If we have the external v_xc calculate KS quantities
348 330 : IF (atom_info(in, im)%atom%read_vxc) THEN
349 0 : ALLOCATE (ext_vxc(atom_info(in, im)%atom%basis%grid%nr))
350 0 : ext_vxc = 0._dp
351 0 : CALL atom_read_external_vxc(ext_vxc, atom_info(in, im)%atom, iw)
352 : CALL calculate_atom_ext_vxc(vxc=ext_vxc, &
353 : atom=atom_info(in, im)%atom, &
354 0 : lprint=.TRUE.)
355 0 : DEALLOCATE (ext_vxc)
356 : END IF
357 :
358 : ! Print out the orbitals if requested
359 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ORBITALS", extension=".log")
360 330 : CALL section_vals_val_get(atom_section, "PRINT%ORBITALS%XMGRACE", l_val=graph)
361 330 : IF (iw > 0) THEN
362 0 : CALL atom_print_orbitals(atom_info(in, im)%atom, iw, xmgrace=graph)
363 : END IF
364 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ORBITALS")
365 :
366 : ! perform a fit of the total electronic density
367 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
368 330 : IF (iw > 0) THEN
369 0 : CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
370 0 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
371 0 : CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
372 : END IF
373 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_DENSITY")
374 :
375 : ! Optimize a local potential for the non-additive kinetic energy term in KG
376 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
377 330 : IF (iw > 0) THEN
378 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
379 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
380 1 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
381 1 : CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
382 : END IF
383 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")
384 :
385 : ! generate a response basis
386 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
387 330 : IF (iw > 0) THEN
388 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
389 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
390 1 : CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
391 : END IF
392 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%RESPONSE_BASIS")
393 :
394 : ! generate a UPF file
395 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%UPF_FILE", extension=".upf", &
396 330 : file_position="REWIND")
397 330 : IF (iw > 0) THEN
398 2 : CALL atom_write_upf(atom_info(in, im)%atom, iw)
399 : END IF
400 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")
401 : END IF
402 :
403 : END DO
404 : END DO
405 :
406 : ! generate a geometrical response basis (GRB)
407 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
408 324 : IF (iw > 0) THEN
409 2 : CALL atom_grb_construction(atom_info, atom_section, iw)
410 : END IF
411 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
412 :
413 : ! Analyze basis sets
414 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
415 324 : IF (iw > 0) THEN
416 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
417 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
418 3 : crad = ptable(zval)%covalent_radius*bohr
419 3 : IF (had_ae) THEN
420 2 : IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
421 2 : IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
422 : END IF
423 3 : IF (had_pp) THEN
424 1 : IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
425 1 : IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
426 : END IF
427 : END IF
428 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
429 :
430 : ! Analyse ADMM approximation
431 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
432 324 : admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
433 324 : IF (iw > 0) THEN
434 1 : CALL atom_admm(atom_info, admm_section, iw)
435 : END IF
436 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")
437 :
438 : ! Generate a separable form of the pseudopotential using Gaussian functions
439 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
440 324 : sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
441 324 : IF (iw > 0) THEN
442 5 : CALL atom_sgp_construction(atom_info, sgp_section, iw)
443 : END IF
444 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
445 :
446 : ! clean up
447 324 : IF (had_ae) THEN
448 256 : CALL atom_int_release(ae_int)
449 256 : CALL atom_ppint_release(ae_int)
450 256 : CALL atom_relint_release(ae_int)
451 : END IF
452 324 : IF (had_pp) THEN
453 70 : CALL atom_int_release(pp_int)
454 70 : CALL atom_ppint_release(pp_int)
455 70 : CALL atom_relint_release(pp_int)
456 : END IF
457 324 : CALL release_atom_basis(ae_basis)
458 324 : CALL release_atom_basis(pp_basis)
459 :
460 324 : CALL release_atom_potential(p_pot)
461 324 : CALL release_atom_potential(ae_pot)
462 :
463 654 : DO in = 1, n_rep
464 984 : DO im = 1, n_meth
465 660 : CALL release_atom_type(atom_info(in, im)%atom)
466 : END DO
467 : END DO
468 324 : DEALLOCATE (atom_info)
469 :
470 324 : DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
471 :
472 324 : CALL timestop(handle)
473 :
474 2268 : END SUBROUTINE atom_energy_opt
475 :
476 : ! **************************************************************************************************
477 : !> \brief Calculate response basis contraction coefficients.
478 : !> \param atom information about the atomic kind
479 : !> \param delta variation of charge used in finite difference derivative calculation
480 : !> \param nder number of wavefunction derivatives to calculate
481 : !> \param iw output file unit
482 : !> \par History
483 : !> * 10.2011 Gram-Schmidt orthogonalization [Joost VandeVondele]
484 : !> * 05.2009 created [Juerg Hutter]
485 : ! **************************************************************************************************
486 1 : SUBROUTINE atom_response_basis(atom, delta, nder, iw)
487 :
488 : TYPE(atom_type), POINTER :: atom
489 : REAL(KIND=dp), INTENT(IN) :: delta
490 : INTEGER, INTENT(IN) :: nder, iw
491 :
492 : INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
493 : s1, s2
494 : REAL(KIND=dp) :: dene, emax, expzet, fhomo, o, prefac, &
495 : zeta
496 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
497 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
498 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
499 1 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp
500 : TYPE(atom_state), POINTER :: state
501 :
502 1 : WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
503 :
504 1 : state => atom%state
505 1 : ovlp => atom%integrals%ovlp
506 :
507 : ! find HOMO
508 1 : lhomo = -1
509 1 : nhomo = -1
510 1 : emax = -HUGE(1._dp)
511 4 : DO l = 0, state%maxl_occ
512 6 : DO i = 1, state%maxn_occ(l)
513 5 : IF (atom%orbitals%ener(i, l) > emax) THEN
514 1 : lhomo = l
515 1 : nhomo = i
516 1 : emax = atom%orbitals%ener(i, l)
517 1 : fhomo = state%occupation(l, i)
518 : END IF
519 : END DO
520 : END DO
521 :
522 1 : s1 = SIZE(atom%orbitals%wfn, 1)
523 1 : s2 = SIZE(atom%orbitals%wfn, 2)
524 6 : ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
525 7 : s2 = MAXVAL(state%maxn_occ) + nder
526 5 : ALLOCATE (rbasis(s1, s2, 0:lmat))
527 91 : rbasis = 0._dp
528 :
529 4 : DO ider = -nder, nder
530 3 : dene = REAL(ider, KIND=dp)*delta
531 3 : CPASSERT(fhomo > ABS(dene))
532 3 : state%occupation(lhomo, nhomo) = fhomo + dene
533 3 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
534 147 : wfn(:, :, :, ider) = atom%orbitals%wfn
535 4 : state%occupation(lhomo, nhomo) = fhomo
536 : END DO
537 :
538 4 : DO l = 0, state%maxl_occ
539 : ! occupied states
540 6 : DO i = 1, MAX(state%maxn_occ(l), 1)
541 24 : rbasis(:, i, l) = wfn(:, i, l, 0)
542 : END DO
543 : ! differentiation
544 6 : DO ider = 1, nder
545 3 : i = MAX(state%maxn_occ(l), 1)
546 3 : SELECT CASE (ider)
547 : CASE (1)
548 21 : rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
549 : CASE (2)
550 0 : rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
551 : CASE (3)
552 : rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
553 0 : + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
554 : CASE DEFAULT
555 3 : CPABORT("")
556 : END SELECT
557 : END DO
558 :
559 : ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
560 3 : n = state%maxn_occ(l) + nder
561 3 : m = atom%basis%nbas(l)
562 8 : DO i = 1, n
563 7 : DO j = 1, i - 1
564 192 : o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
565 19 : rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
566 : END DO
567 480 : o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
568 38 : rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
569 : END DO
570 :
571 : ! check
572 12 : ALLOCATE (amat(n, n))
573 578 : amat(1:n, 1:n) = MATMUL(TRANSPOSE(rbasis(1:m, 1:n, l)), MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
574 8 : DO i = 1, n
575 8 : amat(i, i) = amat(i, i) - 1._dp
576 : END DO
577 17 : IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
578 0 : WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
579 : END IF
580 3 : DEALLOCATE (amat)
581 :
582 : ! Quickstep normalization
583 3 : WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
584 :
585 3 : WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
586 6 : n - nder, " valence ", nder, " response"
587 3 : expzet = 0.25_dp*REAL(2*l + 3, dp)
588 3 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
589 22 : DO i = 1, m
590 18 : zeta = (2._dp*atom%basis%am(i, l))**expzet
591 51 : WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
592 : END DO
593 :
594 : END DO
595 :
596 1 : DEALLOCATE (wfn, rbasis)
597 :
598 1 : WRITE (iw, '(" ",79("*"))')
599 :
600 2 : END SUBROUTINE atom_response_basis
601 :
602 : ! **************************************************************************************************
603 : !> \brief Write pseudo-potential in Quantum Espresso UPF format.
604 : !> \param atom information about the atomic kind
605 : !> \param iw output file unit
606 : !> \par History
607 : !> * 09.2012 created [Juerg Hutter]
608 : ! **************************************************************************************************
609 2 : SUBROUTINE atom_write_upf(atom, iw)
610 :
611 : TYPE(atom_type), POINTER :: atom
612 : INTEGER, INTENT(IN) :: iw
613 :
614 : CHARACTER(LEN=default_string_length) :: string
615 : INTEGER :: i, ibeta, j, k, l, lmax, nbeta, nr, &
616 : nwfc, nwfn
617 : LOGICAL :: up
618 : REAL(KIND=dp) :: pf, rl, rmax
619 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
620 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
621 : TYPE(atom_gthpot_type), POINTER :: pot
622 :
623 2 : IF (.NOT. atom%pp_calc) RETURN
624 2 : IF (atom%potential%ppot_type /= gth_pseudo) RETURN
625 2 : pot => atom%potential%gth_pot
626 2 : CPASSERT(.NOT. pot%lsdpot)
627 :
628 2 : WRITE (iw, '(A)') '<UPF version="2.0.1">'
629 :
630 2 : WRITE (iw, '(T4,A)') '<PP_INFO>'
631 2 : WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
632 2 : WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
633 2 : CALL atom_write_pseudo_param(pot, iw)
634 2 : WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
635 2 : WRITE (iw, '(T4,A)') '</PP_INFO>'
636 2 : WRITE (iw, '(T4,A)') '<PP_HEADER'
637 2 : WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
638 2 : WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
639 2 : WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
640 2 : WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
641 2 : CALL compose(string, "element", cval=pot%symbol)
642 2 : WRITE (iw, '(T8,A)') TRIM(string)
643 2 : WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
644 2 : WRITE (iw, '(T8,A)') 'relativistic="no"'
645 2 : WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
646 2 : WRITE (iw, '(T8,A)') 'is_paw="F"'
647 2 : WRITE (iw, '(T8,A)') 'is_coulomb="F"'
648 2 : WRITE (iw, '(T8,A)') 'has_so="F"'
649 2 : WRITE (iw, '(T8,A)') 'has_wfc="F"'
650 2 : WRITE (iw, '(T8,A)') 'has_gipaw="F"'
651 2 : WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
652 2 : IF (pot%nlcc) THEN
653 1 : WRITE (iw, '(T8,A)') 'core_correction="T"'
654 : ELSE
655 1 : WRITE (iw, '(T8,A)') 'core_correction="F"'
656 : END IF
657 2 : WRITE (iw, '(T8,A)') 'functional="DFT"'
658 2 : CALL compose(string, "z_valence", rval=pot%zion)
659 2 : WRITE (iw, '(T8,A)') TRIM(string)
660 2 : CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
661 2 : WRITE (iw, '(T8,A)') TRIM(string)
662 2 : WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
663 2 : WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
664 2 : lmax = -1
665 14 : DO l = 0, lmat
666 14 : IF (pot%nl(l) > 0) lmax = l
667 : END DO
668 2 : CALL compose(string, "l_max", ival=lmax)
669 2 : WRITE (iw, '(T8,A)') TRIM(string)
670 2 : WRITE (iw, '(T8,A)') 'l_max_rho="0"'
671 2 : WRITE (iw, '(T8,A)') 'l_local="-3"'
672 2 : nr = atom%basis%grid%nr
673 2 : CALL compose(string, "mesh_size", ival=nr)
674 2 : WRITE (iw, '(T8,A)') TRIM(string)
675 14 : nwfc = SUM(atom%state%maxn_occ)
676 2 : CALL compose(string, "number_of_wfc", ival=nwfc)
677 2 : WRITE (iw, '(T8,A)') TRIM(string)
678 14 : nbeta = SUM(pot%nl)
679 2 : CALL compose(string, "number_of_proj", ival=nbeta)
680 2 : WRITE (iw, '(T8,A)') TRIM(string)//'/>'
681 :
682 : ! Mesh
683 2 : up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
684 2 : WRITE (iw, '(T4,A)') '<PP_MESH'
685 2 : WRITE (iw, '(T8,A)') 'dx="1.E+00"'
686 2 : CALL compose(string, "mesh", ival=nr)
687 2 : WRITE (iw, '(T8,A)') TRIM(string)
688 2 : WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
689 804 : rmax = MAXVAL(atom%basis%grid%rad)
690 2 : CALL compose(string, "rmax", rval=rmax)
691 2 : WRITE (iw, '(T8,A)') TRIM(string)
692 2 : WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
693 2 : WRITE (iw, '(T8,A)') '<PP_R type="real"'
694 2 : CALL compose(string, "size", ival=nr)
695 2 : WRITE (iw, '(T12,A)') TRIM(string)
696 2 : WRITE (iw, '(T12,A)') 'columns="4">'
697 2 : IF (up) THEN
698 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
699 : ELSE
700 802 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
701 : END IF
702 2 : WRITE (iw, '(T8,A)') '</PP_R>'
703 2 : WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
704 2 : CALL compose(string, "size", ival=nr)
705 2 : WRITE (iw, '(T12,A)') TRIM(string)
706 2 : WRITE (iw, '(T12,A)') 'columns="4">'
707 2 : IF (up) THEN
708 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
709 : ELSE
710 802 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
711 : END IF
712 2 : WRITE (iw, '(T8,A)') '</PP_RAB>'
713 2 : WRITE (iw, '(T4,A)') '</PP_MESH>'
714 :
715 : ! NLCC
716 2 : IF (pot%nlcc) THEN
717 1 : WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
718 1 : CALL compose(string, "size", ival=nr)
719 1 : WRITE (iw, '(T8,A)') TRIM(string)
720 1 : WRITE (iw, '(T8,A)') 'columns="4">'
721 3 : ALLOCATE (corden(nr))
722 401 : corden(:) = 0.0_dp
723 1 : CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
724 1 : IF (up) THEN
725 0 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
726 : ELSE
727 1 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
728 : END IF
729 1 : DEALLOCATE (corden)
730 1 : WRITE (iw, '(T4,A)') '</PP_NLCC>'
731 : END IF
732 :
733 : ! local PP
734 2 : WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
735 2 : CALL compose(string, "size", ival=nr)
736 2 : WRITE (iw, '(T8,A)') TRIM(string)
737 2 : WRITE (iw, '(T8,A)') 'columns="4">'
738 6 : ALLOCATE (locpot(nr))
739 802 : locpot(:) = 0.0_dp
740 2 : CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
741 2 : IF (up) THEN
742 0 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
743 : ELSE
744 802 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
745 : END IF
746 2 : DEALLOCATE (locpot)
747 2 : WRITE (iw, '(T4,A)') '</PP_LOCAL>'
748 :
749 : ! nonlocal PP
750 2 : WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
751 8 : ALLOCATE (rp(nr), ef(nr), beta(nr))
752 2 : ibeta = 0
753 14 : DO l = 0, lmat
754 12 : IF (pot%nl(l) == 0) CYCLE
755 3 : rl = pot%rcnl(l)
756 1203 : rp(:) = atom%basis%grid%rad
757 1203 : ef(:) = EXP(-0.5_dp*rp*rp/(rl*rl))
758 8 : DO i = 1, pot%nl(l)
759 3 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
760 3 : j = l + 2*i - 1
761 3 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
762 1203 : beta(:) = pf*rp**(l + 2*i - 2)*ef
763 3 : ibeta = ibeta + 1
764 3 : CALL compose(string, "<PP_BETA", counter=ibeta)
765 3 : WRITE (iw, '(T8,A)') TRIM(string)
766 3 : CALL compose(string, "angular_momentum", ival=l)
767 3 : WRITE (iw, '(T12,A)') TRIM(string)
768 3 : WRITE (iw, '(T12,A)') 'type="real"'
769 3 : CALL compose(string, "size", ival=nr)
770 3 : WRITE (iw, '(T12,A)') TRIM(string)
771 3 : WRITE (iw, '(T12,A)') 'columns="4">'
772 1203 : beta(:) = beta*rp
773 3 : IF (up) THEN
774 0 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
775 : ELSE
776 1203 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
777 : END IF
778 3 : CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
779 15 : WRITE (iw, '(T8,A)') TRIM(string)
780 : END DO
781 : END DO
782 2 : DEALLOCATE (rp, ef, beta)
783 : ! nonlocal PP matrix elements
784 8 : ALLOCATE (dij(nbeta, nbeta))
785 10 : dij = 0._dp
786 14 : DO l = 0, lmat
787 12 : IF (pot%nl(l) == 0) CYCLE
788 4 : ibeta = SUM(pot%nl(0:l - 1)) + 1
789 3 : i = ibeta + pot%nl(l) - 1
790 11 : dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
791 : END DO
792 2 : WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
793 2 : WRITE (iw, '(T12,A)') 'columns="4">'
794 10 : WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
795 2 : WRITE (iw, '(T8,A)') '</PP_DIJ>'
796 2 : DEALLOCATE (dij)
797 2 : WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
798 :
799 : ! atomic wavefunctions
800 4 : ALLOCATE (beta(nr))
801 2 : WRITE (iw, '(T4,A)') '<PP_PSWFC>'
802 2 : nwfn = 0
803 14 : DO l = 0, lmat
804 134 : DO i = 1, 10
805 120 : IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
806 4 : nwfn = nwfn + 1
807 4 : CALL compose(string, "<PP_CHI", counter=nwfn)
808 4 : WRITE (iw, '(T8,A)') TRIM(string)
809 4 : CALL compose(string, "l", ival=l)
810 4 : WRITE (iw, '(T12,A)') TRIM(string)
811 4 : CALL compose(string, "occupation", rval=atom%state%occupation(l, i))
812 4 : WRITE (iw, '(T12,A)') TRIM(string)
813 4 : CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
814 4 : WRITE (iw, '(T12,A)') TRIM(string)
815 4 : WRITE (iw, '(T12,A)') 'columns="4">'
816 1604 : beta = 0._dp
817 75 : DO k = 1, atom%basis%nbas(l)
818 28475 : beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
819 : END DO
820 1604 : beta(:) = beta*atom%basis%grid%rad
821 4 : IF (up) THEN
822 0 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=1, nr)
823 : ELSE
824 1604 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=nr, 1, -1)
825 : END IF
826 4 : CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
827 132 : WRITE (iw, '(T8,A)') TRIM(string)
828 : END DO
829 : END DO
830 2 : WRITE (iw, '(T4,A)') '</PP_PSWFC>'
831 2 : DEALLOCATE (beta)
832 :
833 : ! atomic charge
834 4 : ALLOCATE (dens(nr))
835 2 : WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
836 2 : CALL compose(string, "size", ival=nr)
837 2 : WRITE (iw, '(T8,A)') TRIM(string)
838 2 : WRITE (iw, '(T8,A)') 'columns="4">'
839 : CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
840 2 : "RHO", atom%basis%grid%rad)
841 2 : IF (up) THEN
842 0 : WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
843 : ELSE
844 802 : WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
845 : END IF
846 2 : WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
847 2 : DEALLOCATE (dens)
848 :
849 2 : WRITE (iw, '(A)') '</UPF>'
850 :
851 2 : END SUBROUTINE atom_write_upf
852 :
853 : ! **************************************************************************************************
854 : !> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
855 : !> \param string composed string
856 : !> \param tag attribute tag
857 : !> \param counter counter
858 : !> \param rval real variable
859 : !> \param ival integer variable
860 : !> \param cval string variable
861 : !> \param isfinal close the current XML element if this is the last attribute
862 : !> \par History
863 : !> * 09.2012 created [Juerg Hutter]
864 : ! **************************************************************************************************
865 59 : SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
866 : CHARACTER(len=*), INTENT(OUT) :: string
867 : CHARACTER(len=*), INTENT(IN) :: tag
868 : INTEGER, INTENT(IN), OPTIONAL :: counter
869 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rval
870 : INTEGER, INTENT(IN), OPTIONAL :: ival
871 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
872 : LOGICAL, INTENT(IN), OPTIONAL :: isfinal
873 :
874 : CHARACTER(LEN=default_string_length) :: str
875 : LOGICAL :: fin
876 :
877 59 : IF (PRESENT(counter)) THEN
878 14 : WRITE (str, "(I12)") counter
879 45 : ELSEIF (PRESENT(rval)) THEN
880 14 : WRITE (str, "(G18.8)") rval
881 31 : ELSEIF (PRESENT(ival)) THEN
882 29 : WRITE (str, "(I12)") ival
883 2 : ELSEIF (PRESENT(cval)) THEN
884 2 : WRITE (str, "(A)") TRIM(ADJUSTL(cval))
885 : ELSE
886 0 : WRITE (str, "(A)") ""
887 : END IF
888 59 : fin = .FALSE.
889 59 : IF (PRESENT(isfinal)) fin = isfinal
890 59 : IF (PRESENT(counter)) THEN
891 14 : IF (fin) THEN
892 7 : WRITE (string, "(A,A1,A,A1)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str)), '>'
893 : ELSE
894 7 : WRITE (string, "(A,A1,A)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str))
895 : END IF
896 : ELSE
897 45 : IF (fin) THEN
898 0 : WRITE (string, "(A,A2,A,A2)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '>"'
899 : ELSE
900 45 : WRITE (string, "(A,A2,A,A1)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '"'
901 : END IF
902 : END IF
903 59 : END SUBROUTINE compose
904 :
905 10 : END MODULE atom_energy
|
