(git:81d9eb2)
Loading...
Searching...
No Matches
atom_energy.F
Go to the documentation of this file.
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
8MODULE 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: &
44 do_analytic, xc_funct_b3lyp, xc_funct_beefvdw, xc_funct_blyp, xc_funct_bp, &
45 xc_funct_hcth120, xc_funct_no_shortcut, xc_funct_olyp, xc_funct_pade, xc_funct_pbe, &
46 xc_funct_pbe0, xc_funct_tpss, xc_none
47 USE input_section_types, ONLY: section_vals_get,&
48 section_vals_get_subs_vals,&
49 section_vals_type,&
50 section_vals_val_get
51 USE kinds, ONLY: default_string_length,&
52 dp
53 USE mathconstants, ONLY: dfac,&
54 fourpi,&
55 gamma1,&
56 rootpi
57 USE periodic_table, ONLY: nelem,&
58 ptable
59 USE physcon, ONLY: bohr
60#include "./base/base_uses.f90"
61
62 IMPLICIT NONE
63 PRIVATE
64 PUBLIC :: atom_energy_opt
65
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_energy'
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief Compute the atomic energy.
72!> \param atom_section ATOM input section
73!> \par History
74!> * 11.2016 geometrical response basis set [Juerg Hutter]
75!> * 07.2016 ADMM analysis [Juerg Hutter]
76!> * 02.2014 non-additive kinetic energy term [Juerg Hutter]
77!> * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
78!> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
79!> * 05.2009 electronic density fit [Juerg Hutter]
80!> * 04.2009 refactored and renamed to atom_energy_opt() [Juerg Hutter]
81!> * 08.2008 created as atom_energy() [Juerg Hutter]
82! **************************************************************************************************
83 SUBROUTINE atom_energy_opt(atom_section)
84 TYPE(section_vals_type), POINTER :: atom_section
85
86 CHARACTER(len=*), PARAMETER :: routinen = 'atom_energy_opt'
87
88 CHARACTER(LEN=2) :: elem
89 CHARACTER(LEN=default_string_length) :: filename, xcfstr
90 CHARACTER(LEN=default_string_length), &
91 DIMENSION(:), POINTER :: tmpstringlist
92 INTEGER :: do_eric, do_erie, handle, i, im, in, iw, k, maxl, mb, method, mo, n_meth, n_rep, &
93 nder, nr_gh, num_gau, num_gto, num_pol, reltyp, shortcut, zcore, zval, zz
94 INTEGER, DIMENSION(0:lmat) :: maxn
95 INTEGER, DIMENSION(:), POINTER :: cn
96 LOGICAL :: dm, do_gh, do_zmp, doread, eri_c, eri_e, &
97 explicit, graph, had_ae, had_pp, &
98 lcomp, lcond, pp_calc, read_vxc
99 REAL(kind=dp) :: crad, delta, lambda
100 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ext_density, ext_vxc
101 REAL(kind=dp), DIMENSION(0:lmat, 10) :: pocc
102 TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
103 TYPE(atom_integrals), POINTER :: ae_int, pp_int
104 TYPE(atom_optimization_type) :: optimization
105 TYPE(atom_orbitals), POINTER :: orbitals
106 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
107 TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
108 TYPE(atom_state), POINTER :: state
109 TYPE(cp_logger_type), POINTER :: logger
110 TYPE(section_vals_type), POINTER :: admm_section, basis_section, external_vxc_section, &
111 method_section, opt_section, potential_section, powell_section, sgp_section, xc_func, &
112 xc_section, zmp_restart_section, zmp_section
113
114 CALL timeset(routinen, handle)
115
116 ! What atom do we calculate
117 CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
118 CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
119 zz = 0
120 DO i = 1, nelem
121 IF (ptable(i)%symbol == elem) THEN
122 zz = i
123 EXIT
124 END IF
125 END DO
126 IF (zz /= 1) zval = zz
127
128 ! read and set up inofrmation on the basis sets
129 ALLOCATE (ae_basis, pp_basis)
130 basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
131 NULLIFY (ae_basis%grid)
132 CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
133 NULLIFY (pp_basis%grid)
134 basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
135 CALL init_atom_basis(pp_basis, basis_section, zval, "PP")
136
137 ! print general and basis set information
138 logger => cp_get_default_logger()
139 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
140 IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw)
141 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
142 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
143 IF (iw > 0) THEN
144 CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
145 CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
146 END IF
147 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
148
149 ! read and setup information on the pseudopotential
150 NULLIFY (potential_section)
151 potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
152 ALLOCATE (ae_pot, p_pot)
153 CALL init_atom_potential(p_pot, potential_section, zval)
154 CALL init_atom_potential(ae_pot, potential_section, -1)
155
156 ! if the ERI's are calculated analytically, we have to precalculate them
157 eri_c = .false.
158 CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
159 IF (do_eric == do_analytic) eri_c = .true.
160 eri_e = .false.
161 CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
162 IF (do_erie == do_analytic) eri_e = .true.
163 CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
164 CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
165
166 ! information on the states to be calculated
167 CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
168 maxn = 0
169 CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
170 DO in = 1, min(SIZE(cn), 4)
171 maxn(in - 1) = cn(in)
172 END DO
173 DO in = 0, lmat
174 maxn(in) = min(maxn(in), ae_basis%nbas(in))
175 maxn(in) = min(maxn(in), pp_basis%nbas(in))
176 END DO
177
178 ! read optimization section
179 opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
180 CALL read_atom_opt_section(optimization, opt_section)
181
182 had_ae = .false.
183 had_pp = .false.
184
185 ! Check for the total number of electron configurations to be calculated
186 CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
187 ! Check for the total number of method types to be calculated
188 method_section => section_vals_get_subs_vals(atom_section, "METHOD")
189 CALL section_vals_get(method_section, n_repetition=n_meth)
190
191 ! integrals
192 ALLOCATE (ae_int, pp_int)
193
194 ALLOCATE (atom_info(n_rep, n_meth))
195
196 DO in = 1, n_rep
197 DO im = 1, n_meth
198
199 NULLIFY (atom_info(in, im)%atom)
200 CALL create_atom_type(atom_info(in, im)%atom)
201
202 atom_info(in, im)%atom%optimization = optimization
203
204 atom_info(in, im)%atom%z = zval
205
206 xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
207 atom_info(in, im)%atom%xc_section => xc_section
208
209 ! Retrieve the XC functional string for upf output
210 ! This does not work within the routine atom_write_upf below for unknown reasons
211 xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
212 CALL section_vals_get(xc_func, explicit=explicit)
213 IF (explicit) THEN
214 CALL section_vals_val_get(xc_func, "_SECTION_PARAMETERS_", i_val=shortcut)
215 SELECT CASE (shortcut)
216 CASE (xc_funct_no_shortcut, xc_none)
217 xcfstr = "NONE"
218 CASE (xc_funct_pbe0)
219 xcfstr = "PBE0"
220 CASE (xc_funct_beefvdw)
221 xcfstr = "BEEFVDW"
222 CASE (xc_funct_b3lyp)
223 xcfstr = "B3LYP"
224 CASE (xc_funct_blyp)
225 xcfstr = "BLYP"
226 CASE (xc_funct_bp)
227 xcfstr = "BP"
228 CASE (xc_funct_pade)
229 xcfstr = "PADE"
230 CASE (xc_funct_pbe)
231 xcfstr = "PBE"
232 CASE (xc_funct_tpss)
233 xcfstr = "TPSS"
234 CASE (xc_funct_olyp)
235 xcfstr = "OLYP"
236 CASE (xc_funct_hcth120)
237 xcfstr = "HCTH120"
238 CASE DEFAULT
239 xcfstr = "DFT"
240 END SELECT
241 ELSE
242 xcfstr = "DFT"
243 END IF
244
245 ! ZMP Reading input sections if they are found initialize everything
246 do_zmp = .false.
247 doread = .false.
248 read_vxc = .false.
249
250 zmp_section => section_vals_get_subs_vals(method_section, "ZMP")
251 CALL section_vals_get(zmp_section, explicit=do_zmp)
252 atom_info(in, im)%atom%do_zmp = do_zmp
253 CALL section_vals_val_get(zmp_section, "FILE_DENSITY", c_val=filename)
254 atom_info(in, im)%atom%ext_file = filename
255 CALL section_vals_val_get(zmp_section, "GRID_TOL", &
256 r_val=atom_info(in, im)%atom%zmpgrid_tol)
257 CALL section_vals_val_get(zmp_section, "LAMBDA", r_val=lambda)
258 atom_info(in, im)%atom%lambda = lambda
259 CALL section_vals_val_get(zmp_section, "DM", l_val=dm)
260 atom_info(in, im)%atom%dm = dm
261
262 zmp_restart_section => section_vals_get_subs_vals(zmp_section, "RESTART")
263 CALL section_vals_get(zmp_restart_section, explicit=doread)
264 atom_info(in, im)%atom%doread = doread
265 CALL section_vals_val_get(zmp_restart_section, "FILE_RESTART", c_val=filename)
266 atom_info(in, im)%atom%zmp_restart_file = filename
267
268 ! ZMP Reading external vxc section, if found initialize
269 external_vxc_section => section_vals_get_subs_vals(method_section, "EXTERNAL_VXC")
270 CALL section_vals_get(external_vxc_section, explicit=read_vxc)
271 atom_info(in, im)%atom%read_vxc = read_vxc
272 CALL section_vals_val_get(external_vxc_section, "FILE_VXC", c_val=filename)
273 atom_info(in, im)%atom%ext_vxc_file = filename
274 CALL section_vals_val_get(external_vxc_section, "GRID_TOL", &
275 r_val=atom_info(in, im)%atom%zmpvxcgrid_tol)
276
277 ALLOCATE (state)
278
279 ! get the electronic configuration
280 CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
281 c_vals=tmpstringlist)
282
283 ! set occupations
284 CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
285 state%maxl_occ = get_maxl_occ(state%occ)
286 state%maxn_occ = get_maxn_occ(state%occ)
287
288 ! set number of states to be calculated
289 state%maxl_calc = max(maxl, state%maxl_occ)
290 state%maxl_calc = min(lmat, state%maxl_calc)
291 state%maxn_calc = 0
292 DO k = 0, state%maxl_calc
293 state%maxn_calc(k) = max(maxn(k), state%maxn_occ(k))
294 END DO
295
296 ! is there a pseudo potential
297 pp_calc = any(index(tmpstringlist(1:), "CORE") /= 0)
298 IF (pp_calc) THEN
299 ! get and set the core occupations
300 CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
301 CALL atom_set_occupation(tmpstringlist, state%core, pocc)
302 zcore = zval - nint(sum(state%core))
303 CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.true.)
304 ELSE
305 state%core = 0._dp
306 CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.false.)
307 END IF
308
309 CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
310 CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
311 CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)
312
313 IF (atom_consistent_method(method, state%multiplicity)) THEN
314 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
315 CALL atom_print_method(atom_info(in, im)%atom, iw)
316 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
317
318 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
319 IF (pp_calc) THEN
320 IF (iw > 0) CALL atom_print_potential(p_pot, iw)
321 ELSE
322 IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
323 END IF
324 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
325 END IF
326
327 ! calculate integrals
328 IF (pp_calc) THEN
329 ! general integrals
330 CALL atom_int_setup(pp_int, pp_basis, &
331 potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e)
332 ! potential
333 CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot)
334 !
335 NULLIFY (pp_int%tzora, pp_int%hdkh)
336 !
337 CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot)
338 state%maxn_calc(:) = min(state%maxn_calc(:), pp_basis%nbas(:))
339 cpassert(all(state%maxn_calc(:) >= state%maxn_occ))
340 had_pp = .true.
341 ELSE
342 ! general integrals
343 CALL atom_int_setup(ae_int, ae_basis, potential=ae_pot, &
344 eri_coulomb=eri_c, eri_exchange=eri_e)
345 ! potential
346 CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
347 ! relativistic correction terms
348 CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=real(zval, dp))
349 !
350 CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
351 state%maxn_calc(:) = min(state%maxn_calc(:), ae_basis%nbas(:))
352 cpassert(all(state%maxn_calc(:) >= state%maxn_occ))
353 had_ae = .true.
354 END IF
355
356 CALL set_atom(atom_info(in, im)%atom, state=state)
357
358 CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
359 exchange_integral_type=do_erie)
360 atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
361 atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
362
363 NULLIFY (orbitals)
364 mo = maxval(state%maxn_calc)
365 mb = maxval(atom_info(in, im)%atom%basis%nbas)
366 CALL create_atom_orbs(orbitals, mb, mo)
367 CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
368
369 IF (atom_consistent_method(method, state%multiplicity)) THEN
370 !Calculate the electronic structure
371 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
372 CALL calculate_atom(atom_info(in, im)%atom, iw)
373 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
374
375 ! ZMP If we have the external density do zmp
376 IF (atom_info(in, im)%atom%do_zmp) THEN
377 CALL atom_write_zmp_restart(atom_info(in, im)%atom)
378
379 ALLOCATE (ext_density(atom_info(in, im)%atom%basis%grid%nr))
380 ext_density = 0._dp
381 CALL atom_read_external_density(ext_density, atom_info(in, im)%atom, iw)
382 CALL calculate_atom_zmp(ext_density=ext_density, &
383 atom=atom_info(in, im)%atom, &
384 lprint=.true.)
385 DEALLOCATE (ext_density)
386 END IF
387 ! ZMP If we have the external v_xc calculate KS quantities
388 IF (atom_info(in, im)%atom%read_vxc) THEN
389 ALLOCATE (ext_vxc(atom_info(in, im)%atom%basis%grid%nr))
390 ext_vxc = 0._dp
391 CALL atom_read_external_vxc(ext_vxc, atom_info(in, im)%atom, iw)
392 CALL calculate_atom_ext_vxc(vxc=ext_vxc, &
393 atom=atom_info(in, im)%atom, &
394 lprint=.true.)
395 DEALLOCATE (ext_vxc)
396 END IF
397
398 ! Print out the orbitals if requested
399 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ORBITALS", extension=".log")
400 CALL section_vals_val_get(atom_section, "PRINT%ORBITALS%XMGRACE", l_val=graph)
401 IF (iw > 0) THEN
402 CALL atom_print_orbitals(atom_info(in, im)%atom, iw, xmgrace=graph)
403 END IF
404 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ORBITALS")
405
406 ! generate a UPF file
407 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%UPF_FILE", extension=".upf", &
408 file_position="REWIND")
409 IF (iw > 0) THEN
410 CALL atom_write_upf(atom_info(in, im)%atom, iw, xcfstr)
411 END IF
412 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")
413
414 ! perform a fit of the total electronic density
415 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
416 IF (iw > 0) THEN
417 CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
418 powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
419 CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
420 END IF
421 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_DENSITY")
422
423 ! Optimize a local potential for the non-additive kinetic energy term in KG
424 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
425 IF (iw > 0) THEN
426 CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
427 CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
428 powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
429 CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
430 END IF
431 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")
432
433 ! generate a response basis
434 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
435 IF (iw > 0) THEN
436 CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
437 CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
438 CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
439 END IF
440 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%RESPONSE_BASIS")
441
442 END IF
443
444 END DO
445 END DO
446
447 ! generate a geometrical response basis (GRB)
448 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
449 IF (iw > 0) THEN
450 CALL atom_grb_construction(atom_info, atom_section, iw)
451 END IF
452 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
453
454 ! Analyze basis sets
455 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
456 IF (iw > 0) THEN
457 CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
458 CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
459 crad = ptable(zval)%covalent_radius*bohr
460 IF (had_ae) THEN
461 IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
462 IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
463 END IF
464 IF (had_pp) THEN
465 IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
466 IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
467 END IF
468 END IF
469 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
470
471 ! Analyse ADMM approximation
472 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
473 admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
474 IF (iw > 0) THEN
475 CALL atom_admm(atom_info, admm_section, iw)
476 END IF
477 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")
478
479 ! Generate a separable form of the pseudopotential using Gaussian functions
480 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
481 sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
482 IF (iw > 0) THEN
483 CALL atom_sgp_construction(atom_info, sgp_section, iw)
484 END IF
485 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
486
487 ! clean up
488 IF (had_ae) THEN
489 CALL atom_int_release(ae_int)
490 CALL atom_ppint_release(ae_int)
491 CALL atom_relint_release(ae_int)
492 END IF
493 IF (had_pp) THEN
494 CALL atom_int_release(pp_int)
495 CALL atom_ppint_release(pp_int)
496 CALL atom_relint_release(pp_int)
497 END IF
498 CALL release_atom_basis(ae_basis)
499 CALL release_atom_basis(pp_basis)
500
501 CALL release_atom_potential(p_pot)
502 CALL release_atom_potential(ae_pot)
503
504 DO in = 1, n_rep
505 DO im = 1, n_meth
506 CALL release_atom_type(atom_info(in, im)%atom)
507 END DO
508 END DO
509 DEALLOCATE (atom_info)
510
511 DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
512
513 CALL timestop(handle)
514
515 END SUBROUTINE atom_energy_opt
516
517! **************************************************************************************************
518!> \brief Calculate response basis contraction coefficients.
519!> \param atom information about the atomic kind
520!> \param delta variation of charge used in finite difference derivative calculation
521!> \param nder number of wavefunction derivatives to calculate
522!> \param iw output file unit
523!> \par History
524!> * 10.2011 Gram-Schmidt orthogonalization [Joost VandeVondele]
525!> * 05.2009 created [Juerg Hutter]
526! **************************************************************************************************
527 SUBROUTINE atom_response_basis(atom, delta, nder, iw)
528
529 TYPE(atom_type), POINTER :: atom
530 REAL(kind=dp), INTENT(IN) :: delta
531 INTEGER, INTENT(IN) :: nder, iw
532
533 INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
534 s1, s2
535 REAL(kind=dp) :: dene, emax, expzet, fhomo, o, prefac, &
536 zeta
537 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
538 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
539 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
540 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp
541 TYPE(atom_state), POINTER :: state
542
543 WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
544
545 state => atom%state
546 ovlp => atom%integrals%ovlp
547
548 ! find HOMO
549 lhomo = -1
550 nhomo = -1
551 emax = -huge(1._dp)
552 DO l = 0, state%maxl_occ
553 DO i = 1, state%maxn_occ(l)
554 IF (atom%orbitals%ener(i, l) > emax) THEN
555 lhomo = l
556 nhomo = i
557 emax = atom%orbitals%ener(i, l)
558 fhomo = state%occupation(l, i)
559 END IF
560 END DO
561 END DO
562
563 s1 = SIZE(atom%orbitals%wfn, 1)
564 s2 = SIZE(atom%orbitals%wfn, 2)
565 ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
566 s2 = maxval(state%maxn_occ) + nder
567 ALLOCATE (rbasis(s1, s2, 0:lmat))
568 rbasis = 0._dp
569
570 DO ider = -nder, nder
571 dene = real(ider, kind=dp)*delta
572 cpassert(fhomo > abs(dene))
573 state%occupation(lhomo, nhomo) = fhomo + dene
574 CALL calculate_atom(atom, iw=0, noguess=.true.)
575 wfn(:, :, :, ider) = atom%orbitals%wfn
576 state%occupation(lhomo, nhomo) = fhomo
577 END DO
578 ! reset
579 CALL calculate_atom(atom, iw=0, noguess=.true.)
580
581 DO l = 0, state%maxl_occ
582 ! occupied states
583 DO i = 1, max(state%maxn_occ(l), 1)
584 rbasis(:, i, l) = wfn(:, i, l, 0)
585 END DO
586 ! differentiation
587 DO ider = 1, nder
588 i = max(state%maxn_occ(l), 1)
589 SELECT CASE (ider)
590 CASE (1)
591 rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
592 CASE (2)
593 rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
594 CASE (3)
595 rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
596 + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
597 CASE DEFAULT
598 cpabort("")
599 END SELECT
600 END DO
601
602 ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
603 n = state%maxn_occ(l) + nder
604 m = atom%basis%nbas(l)
605 DO i = 1, n
606 DO j = 1, i - 1
607 o = dot_product(rbasis(1:m, j, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
608 rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
609 END DO
610 o = dot_product(rbasis(1:m, i, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
611 rbasis(1:m, i, l) = rbasis(1:m, i, l)/sqrt(o)
612 END DO
613
614 ! check
615 ALLOCATE (amat(n, n))
616 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)))
617 DO i = 1, n
618 amat(i, i) = amat(i, i) - 1._dp
619 END DO
620 IF (maxval(abs(amat)) > 1.e-12) THEN
621 WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", maxval(abs(amat))
622 END IF
623 DEALLOCATE (amat)
624
625 ! Quickstep normalization
626 WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
627
628 WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
629 n - nder, " valence ", nder, " response"
630 expzet = 0.25_dp*real(2*l + 3, dp)
631 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
632 DO i = 1, m
633 zeta = (2._dp*atom%basis%am(i, l))**expzet
634 WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
635 END DO
636
637 END DO
638
639 DEALLOCATE (wfn, rbasis)
640
641 WRITE (iw, '(" ",79("*"))')
642
643 END SUBROUTINE atom_response_basis
644
645! **************************************************************************************************
646!> \brief Write pseudo-potential in Quantum Espresso UPF format.
647!> \param atom information about the atomic kind
648!> \param iw output file unit
649!> \param xcfstr ...
650!> \par History
651!> * 09.2012 created [Juerg Hutter]
652! **************************************************************************************************
653 SUBROUTINE atom_write_upf(atom, iw, xcfstr)
654
655 TYPE(atom_type), POINTER :: atom
656 INTEGER, INTENT(IN) :: iw
657 CHARACTER(LEN=*), INTENT(IN) :: xcfstr
658
659 CHARACTER(LEN=default_string_length) :: string
660 INTEGER :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
661 nbeta, nr, nwfc, nwfn
662 LOGICAL :: soc, up
663 REAL(kind=dp) :: occa, occb, pf, rhoatom, rl, rlp, rmax, &
664 vnlm, vnlp
665 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
666 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
667 TYPE(atom_gthpot_type), POINTER :: pot
668
669 IF (.NOT. atom%pp_calc) RETURN
670 IF (atom%potential%ppot_type /= gth_pseudo) RETURN
671 pot => atom%potential%gth_pot
672 cpassert(.NOT. pot%lsdpot)
673 soc = pot%soc
674
675 WRITE (iw, '(A)') '<UPF version="2.0.1">'
676
677 WRITE (iw, '(T4,A)') '<PP_INFO>'
678 WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
679 WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
680 CALL atom_write_pseudo_param(pot, iw)
681 WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
682 WRITE (iw, '(T4,A)') '</PP_INFO>'
683 WRITE (iw, '(T4,A)') '<PP_HEADER'
684 WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
685 WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
686 WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
687 WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
688 CALL compose(string, "element", cval=pot%symbol)
689 WRITE (iw, '(T8,A)') trim(string)
690 WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
691 WRITE (iw, '(T8,A)') 'relativistic="no"'
692 WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
693 WRITE (iw, '(T8,A)') 'is_paw="F"'
694 WRITE (iw, '(T8,A)') 'is_coulomb="F"'
695 IF (soc) THEN
696 WRITE (iw, '(T8,A)') 'has_so="T"'
697 ELSE
698 WRITE (iw, '(T8,A)') 'has_so="F"'
699 END IF
700 WRITE (iw, '(T8,A)') 'has_wfc="F"'
701 WRITE (iw, '(T8,A)') 'has_gipaw="F"'
702 WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
703 IF (pot%nlcc) THEN
704 WRITE (iw, '(T8,A)') 'core_correction="T"'
705 ELSE
706 WRITE (iw, '(T8,A)') 'core_correction="F"'
707 END IF
708 WRITE (iw, '(T8,A)') 'functional="'//trim(adjustl(xcfstr))//'"'
709 CALL compose(string, "z_valence", rval=pot%zion)
710 WRITE (iw, '(T8,A)') trim(string)
711 CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
712 WRITE (iw, '(T8,A)') trim(string)
713 WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
714 WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
715 lmax = -1
716 DO l = 0, lmat
717 IF (pot%nl(l) > 0) lmax = l
718 END DO
719 CALL compose(string, "l_max", ival=lmax)
720 WRITE (iw, '(T8,A)') trim(string)
721 WRITE (iw, '(T8,A)') 'l_max_rho="0"'
722 WRITE (iw, '(T8,A)') 'l_local="-3"'
723 nr = atom%basis%grid%nr
724 CALL compose(string, "mesh_size", ival=nr)
725 WRITE (iw, '(T8,A)') trim(string)
726 nwfc = sum(atom%state%maxn_occ)
727 CALL compose(string, "number_of_wfc", ival=nwfc)
728 WRITE (iw, '(T8,A)') trim(string)
729 IF (soc) THEN
730 nbeta = pot%nl(0) + 2*sum(pot%nl(1:))
731 ELSE
732 nbeta = sum(pot%nl)
733 END IF
734 CALL compose(string, "number_of_proj", ival=nbeta)
735 WRITE (iw, '(T8,A)') trim(string)//'/>'
736
737 ! Mesh
738 up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
739 WRITE (iw, '(T4,A)') '<PP_MESH'
740 WRITE (iw, '(T8,A)') 'dx="1.E+00"'
741 CALL compose(string, "mesh", ival=nr)
742 WRITE (iw, '(T8,A)') trim(string)
743 WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
744 rmax = maxval(atom%basis%grid%rad)
745 CALL compose(string, "rmax", rval=rmax)
746 WRITE (iw, '(T8,A)') trim(string)
747 WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
748 WRITE (iw, '(T8,A)') '<PP_R type="real"'
749 CALL compose(string, "size", ival=nr)
750 WRITE (iw, '(T12,A)') trim(string)
751 WRITE (iw, '(T12,A)') 'columns="4">'
752 IF (up) THEN
753 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
754 ELSE
755 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
756 END IF
757 WRITE (iw, '(T8,A)') '</PP_R>'
758 WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
759 CALL compose(string, "size", ival=nr)
760 WRITE (iw, '(T12,A)') trim(string)
761 WRITE (iw, '(T12,A)') 'columns="4">'
762 IF (up) THEN
763 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
764 ELSE
765 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
766 END IF
767 WRITE (iw, '(T8,A)') '</PP_RAB>'
768 WRITE (iw, '(T4,A)') '</PP_MESH>'
769
770 ! NLCC
771 IF (pot%nlcc) THEN
772 WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
773 CALL compose(string, "size", ival=nr)
774 WRITE (iw, '(T8,A)') trim(string)
775 WRITE (iw, '(T8,A)') 'columns="4">'
776 ALLOCATE (corden(nr))
777 corden(:) = 0.0_dp
778 CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
779 IF (up) THEN
780 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
781 ELSE
782 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
783 END IF
784 DEALLOCATE (corden)
785 WRITE (iw, '(T4,A)') '</PP_NLCC>'
786 END IF
787
788 ! local PP
789 WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
790 CALL compose(string, "size", ival=nr)
791 WRITE (iw, '(T8,A)') trim(string)
792 WRITE (iw, '(T8,A)') 'columns="4">'
793 ALLOCATE (locpot(nr))
794 locpot(:) = 0.0_dp
795 CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
796 IF (up) THEN
797 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
798 ELSE
799 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
800 END IF
801 DEALLOCATE (locpot)
802 WRITE (iw, '(T4,A)') '</PP_LOCAL>'
803
804 ! nonlocal PP
805 WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
806 ALLOCATE (rp(nr), ef(nr), beta(nr))
807 ibeta = 0
808 DO l = 0, lmat
809 IF (pot%nl(l) == 0) cycle
810 rl = pot%rcnl(l)
811 rp(:) = atom%basis%grid%rad
812 ef(:) = exp(-0.5_dp*rp*rp/(rl*rl))
813 DO i = 1, pot%nl(l)
814 pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
815 j = l + 2*i - 1
816 pf = sqrt(2._dp)/(pf*sqrt(gamma1(j)))
817 beta(:) = pf*rp**(l + 2*i - 2)*ef
818 beta(:) = beta*rp
819 !
820 ibeta = ibeta + 1
821 CALL compose(string, "<PP_BETA", counter=ibeta)
822 WRITE (iw, '(T8,A)') trim(string)
823 CALL compose(string, "angular_momentum", ival=l)
824 WRITE (iw, '(T12,A)') trim(string)
825 WRITE (iw, '(T12,A)') 'type="real"'
826 CALL compose(string, "size", ival=nr)
827 WRITE (iw, '(T12,A)') trim(string)
828 WRITE (iw, '(T12,A)') 'columns="4">'
829 IF (up) THEN
830 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
831 ELSE
832 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
833 END IF
834 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
835 WRITE (iw, '(T8,A)') trim(string)
836 IF (soc .AND. l /= 0) THEN
837 ibeta = ibeta + 1
838 CALL compose(string, "<PP_BETA", counter=ibeta)
839 WRITE (iw, '(T8,A)') trim(string)
840 CALL compose(string, "angular_momentum", ival=l)
841 WRITE (iw, '(T12,A)') trim(string)
842 WRITE (iw, '(T12,A)') 'type="real"'
843 CALL compose(string, "size", ival=nr)
844 WRITE (iw, '(T12,A)') trim(string)
845 WRITE (iw, '(T12,A)') 'columns="4">'
846 IF (up) THEN
847 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
848 ELSE
849 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
850 END IF
851 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
852 WRITE (iw, '(T8,A)') trim(string)
853 END IF
854 END DO
855 END DO
856 DEALLOCATE (rp, ef, beta)
857 ! nonlocal PP matrix elements
858 ALLOCATE (dij(nbeta, nbeta))
859 dij = 0._dp
860 IF (soc) THEN
861 ! from reverse engineering we guess: hnl, knl, hnl, knl, ...
862 n0 = pot%nl(0)
863 DO l = 0, lmat
864 IF (pot%nl(l) == 0) cycle
865 IF (l == 0) THEN
866 ! no SO term for l=0
867 dij(1:n0, 1:n0) = pot%hnl(1:n0, 1:n0, 0)
868 ELSE
869 ibeta = n0 + 2*sum(pot%nl(1:l - 1))
870 DO i = 1, pot%nl(l)
871 ii = 2*(i - 1) + 1
872 DO k = 1, pot%nl(l)
873 kk = 2*(k - 1) + 1
874 vnlp = pot%hnl(i, k, l) + 0.5_dp*l*pot%knl(i, k, l)
875 vnlm = pot%hnl(i, k, l) - 0.5_dp*(l + 1)*pot%knl(i, k, l)
876 dij(ibeta + ii, ibeta + kk) = vnlm
877 dij(ibeta + ii + 1, ibeta + kk + 1) = vnlp
878 END DO
879 END DO
880 END IF
881 END DO
882 ELSE
883 DO l = 0, lmat
884 IF (pot%nl(l) == 0) cycle
885 ibeta = sum(pot%nl(0:l - 1)) + 1
886 i = ibeta + pot%nl(l) - 1
887 dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
888 END DO
889 END IF
890 WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
891 WRITE (iw, '(T12,A)') 'columns="4">'
892 WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
893 WRITE (iw, '(T8,A)') '</PP_DIJ>'
894 DEALLOCATE (dij)
895 WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
896
897 ! atomic wavefunctions
898 ALLOCATE (beta(nr))
899 WRITE (iw, '(T4,A)') '<PP_PSWFC>'
900 nwfn = 0
901 DO l = 0, lmat
902 DO i = 1, 10
903 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
904 occa = atom%state%occupation(l, i)
905 occb = 0.0_dp
906 IF (soc .AND. l /= 0) THEN
907 occb = floor(occa/2)
908 occa = occa - occb
909 END IF
910 nwfn = nwfn + 1
911 CALL compose(string, "<PP_CHI", counter=nwfn)
912 WRITE (iw, '(T8,A)') trim(string)
913 CALL compose(string, "l", ival=l)
914 WRITE (iw, '(T12,A)') trim(string)
915 CALL compose(string, "occupation", rval=occa)
916 WRITE (iw, '(T12,A)') trim(string)
917 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
918 WRITE (iw, '(T12,A)') trim(string)
919 WRITE (iw, '(T12,A)') 'columns="4">'
920 beta = 0._dp
921 DO k = 1, atom%basis%nbas(l)
922 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
923 END DO
924 beta(:) = beta*atom%basis%grid%rad
925 IF (up) THEN
926 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
927 ELSE
928 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
929 END IF
930 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
931 WRITE (iw, '(T8,A)') trim(string)
932 IF (soc .AND. l /= 0) THEN
933 nwfn = nwfn + 1
934 CALL compose(string, "<PP_CHI", counter=nwfn)
935 WRITE (iw, '(T8,A)') trim(string)
936 CALL compose(string, "l", ival=l)
937 WRITE (iw, '(T12,A)') trim(string)
938 CALL compose(string, "occupation", rval=occb)
939 WRITE (iw, '(T12,A)') trim(string)
940 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
941 WRITE (iw, '(T12,A)') trim(string)
942 WRITE (iw, '(T12,A)') 'columns="4">'
943 beta = 0._dp
944 DO k = 1, atom%basis%nbas(l)
945 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
946 END DO
947 beta(:) = beta*atom%basis%grid%rad
948 IF (up) THEN
949 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
950 ELSE
951 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
952 END IF
953 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
954 WRITE (iw, '(T8,A)') trim(string)
955 END IF
956 END DO
957 END DO
958 WRITE (iw, '(T4,A)') '</PP_PSWFC>'
959 DEALLOCATE (beta)
960
961 ! atomic charge
962 ALLOCATE (dens(nr))
963 WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
964 CALL compose(string, "size", ival=nr)
965 WRITE (iw, '(T8,A)') trim(string)
966 WRITE (iw, '(T8,A)') 'columns="4">'
967 CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
968 "RHO", atom%basis%grid%rad)
969 IF (up) THEN
970 WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
971 ELSE
972 WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
973 END IF
974 WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
975
976 ! Check atomic electron count
977 rhoatom = sum(fourpi*atom%basis%grid%wr(:)*dens(:))
978 IF (abs(pot%zion - rhoatom) > 0.01_dp) THEN
979 CALL cp_warn(__location__, "Valence charge and electron count differ by more than 0.01")
980 END IF
981
982 DEALLOCATE (dens)
983 ! PP SOC information
984 IF (soc) THEN
985 WRITE (iw, '(T4,A)') '<PP_SPIN_ORB>'
986! <PP_RELWFC.1 index="1" els="6S" nn="1" lchi="0" jchi="5.000000000000e-1" oc="1.000000000000e0"/>
987 nwfn = 0
988 DO l = 0, lmat
989 DO i = 1, 10
990 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
991 occa = atom%state%occupation(l, i)
992 occb = 0.0_dp
993 IF (l /= 0) THEN
994 occb = floor(occa/2)
995 occa = occa - occb
996 END IF
997 nwfn = nwfn + 1
998 CALL compose(string, "<PP_RELWFC", counter=nwfn)
999 WRITE (iw, '(T8,A)') trim(string)
1000 CALL compose(string, "index", ival=nwfn)
1001 WRITE (iw, '(T12,A)') trim(string)
1002 !!CALL compose(string, "els", cval=xxxx)
1003 !!WRITE (iw, '(T12,A)') TRIM(string)
1004 CALL compose(string, "nn", ival=nwfn)
1005 WRITE (iw, '(T12,A)') trim(string)
1006 CALL compose(string, "lchi", ival=l)
1007 WRITE (iw, '(T12,A)') trim(string)
1008 IF (l == 0) THEN
1009 rlp = 0.5_dp
1010 ELSE
1011 rlp = l - 0.5_dp
1012 END IF
1013 CALL compose(string, "jchi", rval=rlp)
1014 WRITE (iw, '(T12,A)') trim(string)
1015 CALL compose(string, "oc", rval=occa)
1016 WRITE (iw, '(T12,A)') trim(string)
1017 WRITE (iw, '(T12,A)') '/>'
1018 IF (l /= 0) THEN
1019 nwfn = nwfn + 1
1020 CALL compose(string, "<PP_RELWFC", counter=nwfn)
1021 WRITE (iw, '(T8,A)') trim(string)
1022 CALL compose(string, "index", ival=nwfn)
1023 WRITE (iw, '(T12,A)') trim(string)
1024 !!CALL compose(string, "els", cval=xxxx)
1025 !!WRITE (iw, '(T12,A)') TRIM(string)
1026 CALL compose(string, "nn", ival=nwfn)
1027 WRITE (iw, '(T12,A)') trim(string)
1028 CALL compose(string, "lchi", ival=l)
1029 WRITE (iw, '(T12,A)') trim(string)
1030 rlp = l + 0.5_dp
1031 CALL compose(string, "jchi", rval=rlp)
1032 WRITE (iw, '(T12,A)') trim(string)
1033 CALL compose(string, "oc", rval=occa)
1034 WRITE (iw, '(T12,A)') trim(string)
1035 WRITE (iw, '(T12,A)') '/>'
1036 END IF
1037 END DO
1038 END DO
1039 ibeta = 0
1040 DO l = 0, lmat
1041 IF (pot%nl(l) == 0) cycle
1042 DO i = 1, pot%nl(l)
1043 ibeta = ibeta + 1
1044 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1045 WRITE (iw, '(T8,A)') trim(string)
1046 CALL compose(string, "index", ival=ibeta)
1047 WRITE (iw, '(T12,A)') trim(string)
1048 CALL compose(string, "lll", ival=l)
1049 WRITE (iw, '(T12,A)') trim(string)
1050 IF (l == 0) THEN
1051 rlp = 0.5_dp
1052 ELSE
1053 rlp = l - 0.5_dp
1054 END IF
1055 CALL compose(string, "jjj", rval=rlp)
1056 WRITE (iw, '(T12,A)') trim(string)
1057 WRITE (iw, '(T12,A)') '/>'
1058 IF (l > 0) THEN
1059 ibeta = ibeta + 1
1060 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1061 WRITE (iw, '(T8,A)') trim(string)
1062 CALL compose(string, "index", ival=ibeta)
1063 WRITE (iw, '(T12,A)') trim(string)
1064 CALL compose(string, "lll", ival=l)
1065 WRITE (iw, '(T12,A)') trim(string)
1066 rlp = l + 0.5_dp
1067 CALL compose(string, "jjj", rval=rlp)
1068 WRITE (iw, '(T12,A)') trim(string)
1069 WRITE (iw, '(T12,A)') '/>'
1070 END IF
1071 END DO
1072 END DO
1073 WRITE (iw, '(T4,A)') '</PP_SPIN_ORB>'
1074 END IF
1075
1076 WRITE (iw, '(A)') '</UPF>'
1077
1078 END SUBROUTINE atom_write_upf
1079
1080! **************************************************************************************************
1081!> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
1082!> \param string composed string
1083!> \param tag attribute tag
1084!> \param counter counter
1085!> \param rval real variable
1086!> \param ival integer variable
1087!> \param cval string variable
1088!> \param isfinal close the current XML element if this is the last attribute
1089!> \par History
1090!> * 09.2012 created [Juerg Hutter]
1091! **************************************************************************************************
1092 SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
1093 CHARACTER(len=*), INTENT(OUT) :: string
1094 CHARACTER(len=*), INTENT(IN) :: tag
1095 INTEGER, INTENT(IN), OPTIONAL :: counter
1096 REAL(kind=dp), INTENT(IN), OPTIONAL :: rval
1097 INTEGER, INTENT(IN), OPTIONAL :: ival
1098 CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
1099 LOGICAL, INTENT(IN), OPTIONAL :: isfinal
1100
1101 CHARACTER(LEN=default_string_length) :: str
1102 LOGICAL :: fin
1103
1104 IF (PRESENT(counter)) THEN
1105 WRITE (str, "(I12)") counter
1106 ELSEIF (PRESENT(rval)) THEN
1107 WRITE (str, "(G18.8)") rval
1108 ELSEIF (PRESENT(ival)) THEN
1109 WRITE (str, "(I12)") ival
1110 ELSEIF (PRESENT(cval)) THEN
1111 WRITE (str, "(A)") trim(adjustl(cval))
1112 ELSE
1113 WRITE (str, "(A)") ""
1114 END IF
1115 fin = .false.
1116 IF (PRESENT(isfinal)) fin = isfinal
1117 IF (PRESENT(counter)) THEN
1118 IF (fin) THEN
1119 WRITE (string, "(A,A1,A,A1)") trim(adjustl(tag)), '.', trim(adjustl(str)), '>'
1120 ELSE
1121 WRITE (string, "(A,A1,A)") trim(adjustl(tag)), '.', trim(adjustl(str))
1122 END IF
1123 ELSE
1124 IF (fin) THEN
1125 WRITE (string, "(A,A2,A,A2)") trim(adjustl(tag)), '="', trim(adjustl(str)), '>"'
1126 ELSE
1127 WRITE (string, "(A,A2,A,A1)") trim(adjustl(tag)), '="', trim(adjustl(str)), '"'
1128 END IF
1129 END IF
1130 END SUBROUTINE compose
1131
1132END MODULE atom_energy
graph
program graph
Program to Map on grid the hills spawned during a metadynamics run.
Definition graph.F:19
up
static GRID_HOST_DEVICE orbital up(const int i, const orbital a)
Increase i'th component of given orbital angular momentum.
Definition grid_common.h:136
atom_admm_methods
Definition atom_admm_methods.F:9
atom_admm_methods::atom_admm
subroutine, public atom_admm(atom_info, admm_section, iw)
Analysis of ADMM approximation to exact exchange.
Definition atom_admm_methods.F:61
atom_electronic_structure
Definition atom_electronic_structure.F:8
atom_electronic_structure::calculate_atom
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
Definition atom_electronic_structure.F:60
atom_energy
Definition atom_energy.F:8
atom_energy::atom_energy_opt
subroutine, public atom_energy_opt(atom_section)
Compute the atomic energy.
Definition atom_energy.F:84
atom_fit
routines that fit parameters for /from atomic calculations
Definition atom_fit.F:11
atom_fit::atom_fit_density
subroutine, public atom_fit_density(atom, num_gto, norder, iunit, agto, powell_section, results)
Fit the atomic electron density using a geometrical Gaussian basis set.
Definition atom_fit.F:78
atom_fit::atom_fit_kgpot
subroutine, public atom_fit_kgpot(atom, num_gau, num_pol, iunit, powell_section, results)
...
Definition atom_fit.F:1879
atom_grb
Definition atom_grb.F:8
atom_grb::atom_grb_construction
subroutine, public atom_grb_construction(atom_info, atom_section, iw)
Construct geometrical response basis set.
Definition atom_grb.F:77
atom_operators
Calculate the atomic operator matrices.
Definition atom_operators.F:15
atom_operators::atom_ppint_release
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
Definition atom_operators.F:1160
atom_operators::atom_int_setup
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
Definition atom_operators.F:63
atom_operators::atom_relint_setup
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
Definition atom_operators.F:817
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1181
atom_operators::atom_ppint_setup
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
Definition atom_operators.F:336
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1123
atom_output
Routines that print various information about an atomic kind.
Definition atom_output.F:11
atom_output::atom_print_orbitals
subroutine, public atom_print_orbitals(atom, iw, xmgrace)
Print atomic orbitals.
Definition atom_output.F:851
atom_output::atom_print_basis
subroutine, public atom_print_basis(atom_basis, iw, title)
Print atomic basis set.
Definition atom_output.F:298
atom_output::atom_write_pseudo_param
subroutine, public atom_write_pseudo_param(gthpot, iunit)
Print GTH pseudo-potential parameters.
Definition atom_output.F:766
atom_output::atom_print_method
subroutine, public atom_print_method(atom, iw)
Print information about the electronic structure method in use.
Definition atom_output.F:509
atom_output::atom_print_potential
subroutine, public atom_print_potential(potential, iw)
Print information about the pseudo-potential.
Definition atom_output.F:645
atom_output::atom_print_info
subroutine, public atom_print_info(zval, info, iw)
Print an information string related to the atomic kind.
Definition atom_output.F:69
atom_sgp
Definition atom_sgp.F:9
atom_sgp::atom_sgp_construction
subroutine, public atom_sgp_construction(atom_info, input_section, iw)
...
Definition atom_sgp.F:201
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::read_atom_opt_section
subroutine, public read_atom_opt_section(optimization, opt_section)
...
Definition atom_types.F:2409
atom_types::create_atom_type
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:958
atom_types::release_atom_type
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:982
atom_types::release_atom_potential
subroutine, public release_atom_potential(potential)
...
Definition atom_types.F:2535
atom_types::init_atom_basis
subroutine, public init_atom_basis(basis, basis_section, zval, btyp)
Initialize the basis for the atomic code.
Definition atom_types.F:378
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:67
atom_types::set_atom
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
Definition atom_types.F:1023
atom_types::init_atom_potential
subroutine, public init_atom_potential(potential, potential_section, zval)
...
Definition atom_types.F:2435
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:924
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1069
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::atom_condnumber
subroutine, public atom_condnumber(basis, crad, iw)
Print condition numbers of the given atomic basis set.
Definition atom_utils.F:2309
atom_utils::atom_local_potential
pure subroutine, public atom_local_potential(locpot, gthpot, rr)
...
Definition atom_utils.F:799
atom_utils::atom_read_external_vxc
subroutine, public atom_read_external_vxc(vxc, atom, iw)
ZMP subroutine to read external v_xc in the atomic code.
Definition atom_utils.F:607
atom_utils::atom_consistent_method
pure logical function, public atom_consistent_method(method, multiplicity)
Check that the atomic multiplicity is consistent with the electronic structure method.
Definition atom_utils.F:2188
atom_utils::atom_core_density
subroutine, public atom_core_density(corden, potential, typ, rr)
...
Definition atom_utils.F:693
atom_utils::get_maxn_occ
pure integer function, dimension(0:lmat), public get_maxn_occ(occupation)
Return the maximum principal quantum number of occupied orbitals.
Definition atom_utils.F:301
atom_utils::atom_density
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition atom_utils.F:366
atom_utils::atom_read_external_density
subroutine, public atom_read_external_density(density, atom, iw)
ZMP subroutine to read external density from linear grid of density matrix.
Definition atom_utils.F:516
atom_utils::atom_completeness
subroutine, public atom_completeness(basis, zv, iw)
Estimate completeness of the given atomic basis set.
Definition atom_utils.F:2344
atom_utils::atom_write_zmp_restart
subroutine, public atom_write_zmp_restart(atom)
ZMP subroutine to write external restart file.
Definition atom_utils.F:423
atom_utils::atom_set_occupation
subroutine, public atom_set_occupation(ostring, occupation, wfnocc, multiplicity)
Set occupation of atomic orbitals.
Definition atom_utils.F:103
atom_utils::get_maxl_occ
pure integer function, public get_maxl_occ(occupation)
Return the maximum orbital quantum number of occupied orbitals.
Definition atom_utils.F:281
atom_xc
routines that build the integrals of the Vxc potential calculated for the atomic code
Definition atom_xc.F:12
atom_xc::calculate_atom_ext_vxc
subroutine, public calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
ZMP subroutine for the scf external density from the external v_xc.
Definition atom_xc.F:148
atom_xc::calculate_atom_zmp
subroutine, public calculate_atom_zmp(ext_density, atom, lprint, xcmat)
ZMP subroutine for the calculation of the effective constraint potential in the atomic code.
Definition atom_xc.F:78
atom
Definition atom.F:9
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::xc_funct_tpss
integer, parameter, public xc_funct_tpss
Definition input_constants.F:553
input_constants::xc_funct_bp
integer, parameter, public xc_funct_bp
Definition input_constants.F:553
input_constants::xc_funct_olyp
integer, parameter, public xc_funct_olyp
Definition input_constants.F:553
input_constants::xc_funct_pbe
integer, parameter, public xc_funct_pbe
Definition input_constants.F:553
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:149
input_constants::xc_funct_no_shortcut
integer, parameter, public xc_funct_no_shortcut
Definition input_constants.F:553
input_constants::xc_funct_pbe0
integer, parameter, public xc_funct_pbe0
Definition input_constants.F:553
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::xc_funct_beefvdw
integer, parameter, public xc_funct_beefvdw
Definition input_constants.F:553
input_constants::xc_funct_pade
integer, parameter, public xc_funct_pade
Definition input_constants.F:553
input_constants::xc_funct_blyp
integer, parameter, public xc_funct_blyp
Definition input_constants.F:553
input_constants::xc_funct_b3lyp
integer, parameter, public xc_funct_b3lyp
Definition input_constants.F:553
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:553
input_constants::xc_funct_hcth120
integer, parameter, public xc_funct_hcth120
Definition input_constants.F:553
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::gamma1
real(kind=dp), dimension(0:maxfac), parameter, public gamma1
Definition mathconstants.F:95
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
periodic_table::nelem
integer, parameter, public nelem
Definition periodic_table.F:63
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
xmgrace
Routines to facilitate writing XMGRACE files.
Definition xmgrace.F:14
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_gthpot_type
Provides all information about a pseudopotential.
Definition atom_types.F:98
atom_types::atom_integrals
Definition atom_types.F:216
atom_types::atom_optimization_type
Information on optimization procedure.
Definition atom_types.F:283
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:237
atom_types::atom_p_type
Definition atom_types.F:323
atom_types::atom_potential_type
Definition atom_types.F:171
atom_types::atom_state
Provides all information on states and occupation.
Definition atom_types.F:198
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:293
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127