Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE atom_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 : had_ae, had_pp, lcomp, lcond, pp_calc, &
95 : 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 : IF (iw > 0) THEN
361 0 : CALL atom_print_orbitals(atom_info(in, im)%atom, iw)
362 : END IF
363 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ORBITALS")
364 :
365 : ! perform a fit of the total electronic density
366 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
367 330 : IF (iw > 0) THEN
368 0 : CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
369 0 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
370 0 : CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
371 : END IF
372 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_DENSITY")
373 :
374 : ! Optimize a local potential for the non-additive kinetic energy term in KG
375 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
376 330 : IF (iw > 0) THEN
377 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
378 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
379 1 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
380 1 : CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
381 : END IF
382 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")
383 :
384 : ! generate a response basis
385 330 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
386 330 : IF (iw > 0) THEN
387 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
388 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
389 1 : CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
390 : END IF
391 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%RESPONSE_BASIS")
392 :
393 : ! generate a UPF file
394 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%UPF_FILE", extension=".upf", &
395 330 : file_position="REWIND")
396 330 : IF (iw > 0) THEN
397 2 : CALL atom_write_upf(atom_info(in, im)%atom, iw)
398 : END IF
399 330 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")
400 : END IF
401 :
402 : END DO
403 : END DO
404 :
405 : ! generate a geometrical response basis (GRB)
406 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
407 324 : IF (iw > 0) THEN
408 2 : CALL atom_grb_construction(atom_info, atom_section, iw)
409 : END IF
410 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
411 :
412 : ! Analyze basis sets
413 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
414 324 : IF (iw > 0) THEN
415 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
416 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
417 3 : crad = ptable(zval)%covalent_radius*bohr
418 3 : IF (had_ae) THEN
419 2 : IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
420 2 : IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
421 : END IF
422 3 : IF (had_pp) THEN
423 1 : IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
424 1 : IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
425 : END IF
426 : END IF
427 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
428 :
429 : ! Analyse ADMM approximation
430 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
431 324 : admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
432 324 : IF (iw > 0) THEN
433 1 : CALL atom_admm(atom_info, admm_section, iw)
434 : END IF
435 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")
436 :
437 : ! Generate a separable form of the pseudopotential using Gaussian functions
438 324 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
439 324 : sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
440 324 : IF (iw > 0) THEN
441 5 : CALL atom_sgp_construction(atom_info, sgp_section, iw)
442 : END IF
443 324 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
444 :
445 : ! clean up
446 324 : IF (had_ae) THEN
447 256 : CALL atom_int_release(ae_int)
448 256 : CALL atom_ppint_release(ae_int)
449 256 : CALL atom_relint_release(ae_int)
450 : END IF
451 324 : IF (had_pp) THEN
452 70 : CALL atom_int_release(pp_int)
453 70 : CALL atom_ppint_release(pp_int)
454 70 : CALL atom_relint_release(pp_int)
455 : END IF
456 324 : CALL release_atom_basis(ae_basis)
457 324 : CALL release_atom_basis(pp_basis)
458 :
459 324 : CALL release_atom_potential(p_pot)
460 324 : CALL release_atom_potential(ae_pot)
461 :
462 654 : DO in = 1, n_rep
463 984 : DO im = 1, n_meth
464 660 : CALL release_atom_type(atom_info(in, im)%atom)
465 : END DO
466 : END DO
467 324 : DEALLOCATE (atom_info)
468 :
469 324 : DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
470 :
471 324 : CALL timestop(handle)
472 :
473 2268 : END SUBROUTINE atom_energy_opt
474 :
475 : ! **************************************************************************************************
476 : !> \brief Calculate response basis contraction coefficients.
477 : !> \param atom information about the atomic kind
478 : !> \param delta variation of charge used in finite difference derivative calculation
479 : !> \param nder number of wavefunction derivatives to calculate
480 : !> \param iw output file unit
481 : !> \par History
482 : !> * 10.2011 Gram-Schmidt orthogonalization [Joost VandeVondele]
483 : !> * 05.2009 created [Juerg Hutter]
484 : ! **************************************************************************************************
485 1 : SUBROUTINE atom_response_basis(atom, delta, nder, iw)
486 :
487 : TYPE(atom_type), POINTER :: atom
488 : REAL(KIND=dp), INTENT(IN) :: delta
489 : INTEGER, INTENT(IN) :: nder, iw
490 :
491 : INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
492 : s1, s2
493 : REAL(KIND=dp) :: dene, emax, expzet, fhomo, o, prefac, &
494 : zeta
495 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
496 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
497 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
498 1 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp
499 : TYPE(atom_state), POINTER :: state
500 :
501 1 : WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
502 :
503 1 : state => atom%state
504 1 : ovlp => atom%integrals%ovlp
505 :
506 : ! find HOMO
507 1 : lhomo = -1
508 1 : nhomo = -1
509 1 : emax = -HUGE(1._dp)
510 4 : DO l = 0, state%maxl_occ
511 6 : DO i = 1, state%maxn_occ(l)
512 5 : IF (atom%orbitals%ener(i, l) > emax) THEN
513 1 : lhomo = l
514 1 : nhomo = i
515 1 : emax = atom%orbitals%ener(i, l)
516 1 : fhomo = state%occupation(l, i)
517 : END IF
518 : END DO
519 : END DO
520 :
521 1 : s1 = SIZE(atom%orbitals%wfn, 1)
522 1 : s2 = SIZE(atom%orbitals%wfn, 2)
523 6 : ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
524 7 : s2 = MAXVAL(state%maxn_occ) + nder
525 5 : ALLOCATE (rbasis(s1, s2, 0:lmat))
526 91 : rbasis = 0._dp
527 :
528 4 : DO ider = -nder, nder
529 3 : dene = REAL(ider, KIND=dp)*delta
530 3 : CPASSERT(fhomo > ABS(dene))
531 3 : state%occupation(lhomo, nhomo) = fhomo + dene
532 3 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
533 147 : wfn(:, :, :, ider) = atom%orbitals%wfn
534 4 : state%occupation(lhomo, nhomo) = fhomo
535 : END DO
536 :
537 4 : DO l = 0, state%maxl_occ
538 : ! occupied states
539 6 : DO i = 1, MAX(state%maxn_occ(l), 1)
540 24 : rbasis(:, i, l) = wfn(:, i, l, 0)
541 : END DO
542 : ! differentiation
543 6 : DO ider = 1, nder
544 3 : i = MAX(state%maxn_occ(l), 1)
545 3 : SELECT CASE (ider)
546 : CASE (1)
547 21 : rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
548 : CASE (2)
549 0 : rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
550 : CASE (3)
551 : rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
552 0 : + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
553 : CASE DEFAULT
554 3 : CPABORT("")
555 : END SELECT
556 : END DO
557 :
558 : ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
559 3 : n = state%maxn_occ(l) + nder
560 3 : m = atom%basis%nbas(l)
561 8 : DO i = 1, n
562 7 : DO j = 1, i - 1
563 192 : o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
564 19 : rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
565 : END DO
566 480 : o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
567 38 : rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
568 : END DO
569 :
570 : ! check
571 12 : ALLOCATE (amat(n, n))
572 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)))
573 8 : DO i = 1, n
574 8 : amat(i, i) = amat(i, i) - 1._dp
575 : END DO
576 17 : IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
577 0 : WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
578 : END IF
579 3 : DEALLOCATE (amat)
580 :
581 : ! Quickstep normalization
582 3 : WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
583 :
584 3 : WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
585 6 : n - nder, " valence ", nder, " response"
586 3 : expzet = 0.25_dp*REAL(2*l + 3, dp)
587 3 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
588 22 : DO i = 1, m
589 18 : zeta = (2._dp*atom%basis%am(i, l))**expzet
590 51 : WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
591 : END DO
592 :
593 : END DO
594 :
595 1 : DEALLOCATE (wfn, rbasis)
596 :
597 1 : WRITE (iw, '(" ",79("*"))')
598 :
599 2 : END SUBROUTINE atom_response_basis
600 :
601 : ! **************************************************************************************************
602 : !> \brief Write pseudo-potential in Quantum Espresso UPF format.
603 : !> \param atom information about the atomic kind
604 : !> \param iw output file unit
605 : !> \par History
606 : !> * 09.2012 created [Juerg Hutter]
607 : ! **************************************************************************************************
608 2 : SUBROUTINE atom_write_upf(atom, iw)
609 :
610 : TYPE(atom_type), POINTER :: atom
611 : INTEGER, INTENT(IN) :: iw
612 :
613 : CHARACTER(LEN=default_string_length) :: string
614 : INTEGER :: i, ibeta, j, k, l, lmax, nbeta, nr, &
615 : nwfc, nwfn
616 : LOGICAL :: up
617 : REAL(KIND=dp) :: pf, rl, rmax
618 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
619 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
620 : TYPE(atom_gthpot_type), POINTER :: pot
621 :
622 2 : IF (.NOT. atom%pp_calc) RETURN
623 2 : IF (atom%potential%ppot_type /= gth_pseudo) RETURN
624 2 : pot => atom%potential%gth_pot
625 2 : CPASSERT(.NOT. pot%lsdpot)
626 :
627 2 : WRITE (iw, '(A)') '<UPF version="2.0.1">'
628 :
629 2 : WRITE (iw, '(T4,A)') '<PP_INFO>'
630 2 : WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
631 2 : WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
632 2 : CALL atom_write_pseudo_param(pot, iw)
633 2 : WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
634 2 : WRITE (iw, '(T4,A)') '</PP_INFO>'
635 2 : WRITE (iw, '(T4,A)') '<PP_HEADER'
636 2 : WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
637 2 : WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
638 2 : WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
639 2 : WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
640 2 : CALL compose(string, "element", cval=pot%symbol)
641 2 : WRITE (iw, '(T8,A)') TRIM(string)
642 2 : WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
643 2 : WRITE (iw, '(T8,A)') 'relativistic="no"'
644 2 : WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
645 2 : WRITE (iw, '(T8,A)') 'is_paw="F"'
646 2 : WRITE (iw, '(T8,A)') 'is_coulomb="F"'
647 2 : WRITE (iw, '(T8,A)') 'has_so="F"'
648 2 : WRITE (iw, '(T8,A)') 'has_wfc="F"'
649 2 : WRITE (iw, '(T8,A)') 'has_gipaw="F"'
650 2 : WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
651 2 : IF (pot%nlcc) THEN
652 1 : WRITE (iw, '(T8,A)') 'core_correction="T"'
653 : ELSE
654 1 : WRITE (iw, '(T8,A)') 'core_correction="F"'
655 : END IF
656 2 : WRITE (iw, '(T8,A)') 'functional="DFT"'
657 2 : CALL compose(string, "z_valence", rval=pot%zion)
658 2 : WRITE (iw, '(T8,A)') TRIM(string)
659 2 : CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
660 2 : WRITE (iw, '(T8,A)') TRIM(string)
661 2 : WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
662 2 : WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
663 2 : lmax = -1
664 14 : DO l = 0, lmat
665 14 : IF (pot%nl(l) > 0) lmax = l
666 : END DO
667 2 : CALL compose(string, "l_max", ival=lmax)
668 2 : WRITE (iw, '(T8,A)') TRIM(string)
669 2 : WRITE (iw, '(T8,A)') 'l_max_rho="0"'
670 2 : WRITE (iw, '(T8,A)') 'l_local="-3"'
671 2 : nr = atom%basis%grid%nr
672 2 : CALL compose(string, "mesh_size", ival=nr)
673 2 : WRITE (iw, '(T8,A)') TRIM(string)
674 14 : nwfc = SUM(atom%state%maxn_occ)
675 2 : CALL compose(string, "number_of_wfc", ival=nwfc)
676 2 : WRITE (iw, '(T8,A)') TRIM(string)
677 14 : nbeta = SUM(pot%nl)
678 2 : CALL compose(string, "number_of_proj", ival=nbeta)
679 2 : WRITE (iw, '(T8,A)') TRIM(string)//'/>'
680 :
681 : ! Mesh
682 2 : up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
683 2 : WRITE (iw, '(T4,A)') '<PP_MESH'
684 2 : WRITE (iw, '(T8,A)') 'dx="1.E+00"'
685 2 : CALL compose(string, "mesh", ival=nr)
686 2 : WRITE (iw, '(T8,A)') TRIM(string)
687 2 : WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
688 804 : rmax = MAXVAL(atom%basis%grid%rad)
689 2 : CALL compose(string, "rmax", rval=rmax)
690 2 : WRITE (iw, '(T8,A)') TRIM(string)
691 2 : WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
692 2 : WRITE (iw, '(T8,A)') '<PP_R type="real"'
693 2 : CALL compose(string, "size", ival=nr)
694 2 : WRITE (iw, '(T12,A)') TRIM(string)
695 2 : WRITE (iw, '(T12,A)') 'columns="4">'
696 2 : IF (up) THEN
697 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
698 : ELSE
699 802 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
700 : END IF
701 2 : WRITE (iw, '(T8,A)') '</PP_R>'
702 2 : WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
703 2 : CALL compose(string, "size", ival=nr)
704 2 : WRITE (iw, '(T12,A)') TRIM(string)
705 2 : WRITE (iw, '(T12,A)') 'columns="4">'
706 2 : IF (up) THEN
707 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
708 : ELSE
709 802 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
710 : END IF
711 2 : WRITE (iw, '(T8,A)') '</PP_RAB>'
712 2 : WRITE (iw, '(T4,A)') '</PP_MESH>'
713 :
714 : ! NLCC
715 2 : IF (pot%nlcc) THEN
716 1 : WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
717 1 : CALL compose(string, "size", ival=nr)
718 1 : WRITE (iw, '(T8,A)') TRIM(string)
719 1 : WRITE (iw, '(T8,A)') 'columns="4">'
720 3 : ALLOCATE (corden(nr))
721 401 : corden(:) = 0.0_dp
722 1 : CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
723 1 : IF (up) THEN
724 0 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
725 : ELSE
726 1 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
727 : END IF
728 1 : DEALLOCATE (corden)
729 1 : WRITE (iw, '(T4,A)') '</PP_NLCC>'
730 : END IF
731 :
732 : ! local PP
733 2 : WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
734 2 : CALL compose(string, "size", ival=nr)
735 2 : WRITE (iw, '(T8,A)') TRIM(string)
736 2 : WRITE (iw, '(T8,A)') 'columns="4">'
737 6 : ALLOCATE (locpot(nr))
738 802 : locpot(:) = 0.0_dp
739 2 : CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
740 2 : IF (up) THEN
741 0 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
742 : ELSE
743 802 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
744 : END IF
745 2 : DEALLOCATE (locpot)
746 2 : WRITE (iw, '(T4,A)') '</PP_LOCAL>'
747 :
748 : ! nonlocal PP
749 2 : WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
750 14 : ALLOCATE (rp(nr), ef(nr), beta(nr))
751 2 : ibeta = 0
752 14 : DO l = 0, lmat
753 12 : IF (pot%nl(l) == 0) CYCLE
754 3 : rl = pot%rcnl(l)
755 1203 : rp(:) = atom%basis%grid%rad
756 1203 : ef(:) = EXP(-0.5_dp*rp*rp/(rl*rl))
757 8 : DO i = 1, pot%nl(l)
758 3 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
759 3 : j = l + 2*i - 1
760 3 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
761 1203 : beta(:) = pf*rp**(l + 2*i - 2)*ef
762 3 : ibeta = ibeta + 1
763 3 : CALL compose(string, "<PP_BETA", counter=ibeta)
764 3 : WRITE (iw, '(T8,A)') TRIM(string)
765 3 : CALL compose(string, "angular_momentum", ival=l)
766 3 : WRITE (iw, '(T12,A)') TRIM(string)
767 3 : WRITE (iw, '(T12,A)') 'type="real"'
768 3 : CALL compose(string, "size", ival=nr)
769 3 : WRITE (iw, '(T12,A)') TRIM(string)
770 3 : WRITE (iw, '(T12,A)') 'columns="4">'
771 1203 : beta(:) = beta*rp
772 3 : IF (up) THEN
773 0 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
774 : ELSE
775 1203 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
776 : END IF
777 3 : CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
778 15 : WRITE (iw, '(T8,A)') TRIM(string)
779 : END DO
780 : END DO
781 2 : DEALLOCATE (rp, ef, beta)
782 : ! nonlocal PP matrix elements
783 8 : ALLOCATE (dij(nbeta, nbeta))
784 10 : dij = 0._dp
785 14 : DO l = 0, lmat
786 12 : IF (pot%nl(l) == 0) CYCLE
787 4 : ibeta = SUM(pot%nl(0:l - 1)) + 1
788 3 : i = ibeta + pot%nl(l) - 1
789 11 : dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
790 : END DO
791 2 : WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
792 2 : WRITE (iw, '(T12,A)') 'columns="4">'
793 10 : WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
794 2 : WRITE (iw, '(T8,A)') '</PP_DIJ>'
795 2 : DEALLOCATE (dij)
796 2 : WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
797 :
798 : ! atomic wavefunctions
799 6 : ALLOCATE (beta(nr))
800 2 : WRITE (iw, '(T4,A)') '<PP_PSWFC>'
801 2 : nwfn = 0
802 14 : DO l = 0, lmat
803 134 : DO i = 1, 10
804 120 : IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
805 4 : nwfn = nwfn + 1
806 4 : CALL compose(string, "<PP_CHI", counter=nwfn)
807 4 : WRITE (iw, '(T8,A)') TRIM(string)
808 4 : CALL compose(string, "l", ival=l)
809 4 : WRITE (iw, '(T12,A)') TRIM(string)
810 4 : CALL compose(string, "occupation", rval=atom%state%occupation(l, i))
811 4 : WRITE (iw, '(T12,A)') TRIM(string)
812 4 : CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
813 4 : WRITE (iw, '(T12,A)') TRIM(string)
814 4 : WRITE (iw, '(T12,A)') 'columns="4">'
815 1604 : beta = 0._dp
816 75 : DO k = 1, atom%basis%nbas(l)
817 28475 : beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
818 : END DO
819 1604 : beta(:) = beta*atom%basis%grid%rad
820 4 : IF (up) THEN
821 0 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=1, nr)
822 : ELSE
823 1604 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=nr, 1, -1)
824 : END IF
825 4 : CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
826 132 : WRITE (iw, '(T8,A)') TRIM(string)
827 : END DO
828 : END DO
829 2 : WRITE (iw, '(T4,A)') '</PP_PSWFC>'
830 2 : DEALLOCATE (beta)
831 :
832 : ! atomic charge
833 6 : ALLOCATE (dens(nr))
834 2 : WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
835 2 : CALL compose(string, "size", ival=nr)
836 2 : WRITE (iw, '(T8,A)') TRIM(string)
837 2 : WRITE (iw, '(T8,A)') 'columns="4">'
838 : CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
839 2 : "RHO", atom%basis%grid%rad)
840 2 : IF (up) THEN
841 0 : WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
842 : ELSE
843 802 : WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
844 : END IF
845 2 : WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
846 2 : DEALLOCATE (dens)
847 :
848 2 : WRITE (iw, '(A)') '</UPF>'
849 :
850 2 : END SUBROUTINE atom_write_upf
851 :
852 : ! **************************************************************************************************
853 : !> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
854 : !> \param string composed string
855 : !> \param tag attribute tag
856 : !> \param counter counter
857 : !> \param rval real variable
858 : !> \param ival integer variable
859 : !> \param cval string variable
860 : !> \param isfinal close the current XML element if this is the last attribute
861 : !> \par History
862 : !> * 09.2012 created [Juerg Hutter]
863 : ! **************************************************************************************************
864 59 : SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
865 : CHARACTER(len=*), INTENT(OUT) :: string
866 : CHARACTER(len=*), INTENT(IN) :: tag
867 : INTEGER, INTENT(IN), OPTIONAL :: counter
868 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rval
869 : INTEGER, INTENT(IN), OPTIONAL :: ival
870 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
871 : LOGICAL, INTENT(IN), OPTIONAL :: isfinal
872 :
873 : CHARACTER(LEN=default_string_length) :: str
874 : LOGICAL :: fin
875 :
876 59 : IF (PRESENT(counter)) THEN
877 14 : WRITE (str, "(I12)") counter
878 45 : ELSEIF (PRESENT(rval)) THEN
879 14 : WRITE (str, "(G18.8)") rval
880 31 : ELSEIF (PRESENT(ival)) THEN
881 29 : WRITE (str, "(I12)") ival
882 2 : ELSEIF (PRESENT(cval)) THEN
883 2 : WRITE (str, "(A)") TRIM(ADJUSTL(cval))
884 : ELSE
885 0 : WRITE (str, "(A)") ""
886 : END IF
887 59 : fin = .FALSE.
888 59 : IF (PRESENT(isfinal)) fin = isfinal
889 59 : IF (PRESENT(counter)) THEN
890 14 : IF (fin) THEN
891 7 : WRITE (string, "(A,A1,A,A1)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str)), '>'
892 : ELSE
893 7 : WRITE (string, "(A,A1,A)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str))
894 : END IF
895 : ELSE
896 45 : IF (fin) THEN
897 0 : WRITE (string, "(A,A2,A,A2)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '>"'
898 : ELSE
899 45 : WRITE (string, "(A,A2,A,A1)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '"'
900 : END IF
901 : END IF
902 59 : END SUBROUTINE compose
903 :
904 10 : END MODULE atom_energy
|
