(git:badb799)
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 gamma1,&
55 pi,&
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 ! perform a fit of the total electronic density
407 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
408 IF (iw > 0) THEN
409 CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
410 powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
411 CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
412 END IF
413 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_DENSITY")
414
415 ! Optimize a local potential for the non-additive kinetic energy term in KG
416 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
417 IF (iw > 0) THEN
418 CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
419 CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
420 powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
421 CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
422 END IF
423 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")
424
425 ! generate a response basis
426 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
427 IF (iw > 0) THEN
428 CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
429 CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
430 CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
431 END IF
432 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%RESPONSE_BASIS")
433
434 ! generate a UPF file
435 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%UPF_FILE", extension=".upf", &
436 file_position="REWIND")
437 IF (iw > 0) THEN
438 CALL atom_write_upf(atom_info(in, im)%atom, iw, xcfstr)
439 END IF
440 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")
441 END IF
442
443 END DO
444 END DO
445
446 ! generate a geometrical response basis (GRB)
447 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
448 IF (iw > 0) THEN
449 CALL atom_grb_construction(atom_info, atom_section, iw)
450 END IF
451 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
452
453 ! Analyze basis sets
454 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
455 IF (iw > 0) THEN
456 CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
457 CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
458 crad = ptable(zval)%covalent_radius*bohr
459 IF (had_ae) THEN
460 IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
461 IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
462 END IF
463 IF (had_pp) THEN
464 IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
465 IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
466 END IF
467 END IF
468 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
469
470 ! Analyse ADMM approximation
471 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
472 admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
473 IF (iw > 0) THEN
474 CALL atom_admm(atom_info, admm_section, iw)
475 END IF
476 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")
477
478 ! Generate a separable form of the pseudopotential using Gaussian functions
479 iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
480 sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
481 IF (iw > 0) THEN
482 CALL atom_sgp_construction(atom_info, sgp_section, iw)
483 END IF
484 CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
485
486 ! clean up
487 IF (had_ae) THEN
488 CALL atom_int_release(ae_int)
489 CALL atom_ppint_release(ae_int)
490 CALL atom_relint_release(ae_int)
491 END IF
492 IF (had_pp) THEN
493 CALL atom_int_release(pp_int)
494 CALL atom_ppint_release(pp_int)
495 CALL atom_relint_release(pp_int)
496 END IF
497 CALL release_atom_basis(ae_basis)
498 CALL release_atom_basis(pp_basis)
499
500 CALL release_atom_potential(p_pot)
501 CALL release_atom_potential(ae_pot)
502
503 DO in = 1, n_rep
504 DO im = 1, n_meth
505 CALL release_atom_type(atom_info(in, im)%atom)
506 END DO
507 END DO
508 DEALLOCATE (atom_info)
509
510 DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
511
512 CALL timestop(handle)
513
514 END SUBROUTINE atom_energy_opt
515
516! **************************************************************************************************
517!> \brief Calculate response basis contraction coefficients.
518!> \param atom information about the atomic kind
519!> \param delta variation of charge used in finite difference derivative calculation
520!> \param nder number of wavefunction derivatives to calculate
521!> \param iw output file unit
522!> \par History
523!> * 10.2011 Gram-Schmidt orthogonalization [Joost VandeVondele]
524!> * 05.2009 created [Juerg Hutter]
525! **************************************************************************************************
526 SUBROUTINE atom_response_basis(atom, delta, nder, iw)
527
528 TYPE(atom_type), POINTER :: atom
529 REAL(kind=dp), INTENT(IN) :: delta
530 INTEGER, INTENT(IN) :: nder, iw
531
532 INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
533 s1, s2
534 REAL(kind=dp) :: dene, emax, expzet, fhomo, o, prefac, &
535 zeta
536 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
537 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
538 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
539 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp
540 TYPE(atom_state), POINTER :: state
541
542 WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
543
544 state => atom%state
545 ovlp => atom%integrals%ovlp
546
547 ! find HOMO
548 lhomo = -1
549 nhomo = -1
550 emax = -huge(1._dp)
551 DO l = 0, state%maxl_occ
552 DO i = 1, state%maxn_occ(l)
553 IF (atom%orbitals%ener(i, l) > emax) THEN
554 lhomo = l
555 nhomo = i
556 emax = atom%orbitals%ener(i, l)
557 fhomo = state%occupation(l, i)
558 END IF
559 END DO
560 END DO
561
562 s1 = SIZE(atom%orbitals%wfn, 1)
563 s2 = SIZE(atom%orbitals%wfn, 2)
564 ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
565 s2 = maxval(state%maxn_occ) + nder
566 ALLOCATE (rbasis(s1, s2, 0:lmat))
567 rbasis = 0._dp
568
569 DO ider = -nder, nder
570 dene = real(ider, kind=dp)*delta
571 cpassert(fhomo > abs(dene))
572 state%occupation(lhomo, nhomo) = fhomo + dene
573 CALL calculate_atom(atom, iw=0, noguess=.true.)
574 wfn(:, :, :, ider) = atom%orbitals%wfn
575 state%occupation(lhomo, nhomo) = fhomo
576 END DO
577
578 DO l = 0, state%maxl_occ
579 ! occupied states
580 DO i = 1, max(state%maxn_occ(l), 1)
581 rbasis(:, i, l) = wfn(:, i, l, 0)
582 END DO
583 ! differentiation
584 DO ider = 1, nder
585 i = max(state%maxn_occ(l), 1)
586 SELECT CASE (ider)
587 CASE (1)
588 rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
589 CASE (2)
590 rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
591 CASE (3)
592 rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
593 + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
594 CASE DEFAULT
595 cpabort("")
596 END SELECT
597 END DO
598
599 ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
600 n = state%maxn_occ(l) + nder
601 m = atom%basis%nbas(l)
602 DO i = 1, n
603 DO j = 1, i - 1
604 o = dot_product(rbasis(1:m, j, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
605 rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
606 END DO
607 o = dot_product(rbasis(1:m, i, 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)/sqrt(o)
609 END DO
610
611 ! check
612 ALLOCATE (amat(n, n))
613 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)))
614 DO i = 1, n
615 amat(i, i) = amat(i, i) - 1._dp
616 END DO
617 IF (maxval(abs(amat)) > 1.e-12) THEN
618 WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", maxval(abs(amat))
619 END IF
620 DEALLOCATE (amat)
621
622 ! Quickstep normalization
623 WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
624
625 WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
626 n - nder, " valence ", nder, " response"
627 expzet = 0.25_dp*real(2*l + 3, dp)
628 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
629 DO i = 1, m
630 zeta = (2._dp*atom%basis%am(i, l))**expzet
631 WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
632 END DO
633
634 END DO
635
636 DEALLOCATE (wfn, rbasis)
637
638 WRITE (iw, '(" ",79("*"))')
639
640 END SUBROUTINE atom_response_basis
641
642! **************************************************************************************************
643!> \brief Write pseudo-potential in Quantum Espresso UPF format.
644!> \param atom information about the atomic kind
645!> \param iw output file unit
646!> \param xcfstr ...
647!> \par History
648!> * 09.2012 created [Juerg Hutter]
649! **************************************************************************************************
650 SUBROUTINE atom_write_upf(atom, iw, xcfstr)
651
652 TYPE(atom_type), POINTER :: atom
653 INTEGER, INTENT(IN) :: iw
654 CHARACTER(LEN=*), INTENT(IN) :: xcfstr
655
656 CHARACTER(LEN=default_string_length) :: string
657 INTEGER :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
658 nbeta, nr, nwfc, nwfn
659 LOGICAL :: soc, up
660 REAL(kind=dp) :: occa, occb, pf, rl, rlp, rmax, vnlm, vnlp
661 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
662 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
663 TYPE(atom_gthpot_type), POINTER :: pot
664
665 IF (.NOT. atom%pp_calc) RETURN
666 IF (atom%potential%ppot_type /= gth_pseudo) RETURN
667 pot => atom%potential%gth_pot
668 cpassert(.NOT. pot%lsdpot)
669 soc = pot%soc
670
671 WRITE (iw, '(A)') '<UPF version="2.0.1">'
672
673 WRITE (iw, '(T4,A)') '<PP_INFO>'
674 WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
675 WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
676 CALL atom_write_pseudo_param(pot, iw)
677 WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
678 WRITE (iw, '(T4,A)') '</PP_INFO>'
679 WRITE (iw, '(T4,A)') '<PP_HEADER'
680 WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
681 WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
682 WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
683 WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
684 CALL compose(string, "element", cval=pot%symbol)
685 WRITE (iw, '(T8,A)') trim(string)
686 WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
687 WRITE (iw, '(T8,A)') 'relativistic="no"'
688 WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
689 WRITE (iw, '(T8,A)') 'is_paw="F"'
690 WRITE (iw, '(T8,A)') 'is_coulomb="F"'
691 IF (soc) THEN
692 WRITE (iw, '(T8,A)') 'has_so="T"'
693 ELSE
694 WRITE (iw, '(T8,A)') 'has_so="F"'
695 END IF
696 WRITE (iw, '(T8,A)') 'has_wfc="F"'
697 WRITE (iw, '(T8,A)') 'has_gipaw="F"'
698 WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
699 IF (pot%nlcc) THEN
700 WRITE (iw, '(T8,A)') 'core_correction="T"'
701 ELSE
702 WRITE (iw, '(T8,A)') 'core_correction="F"'
703 END IF
704 WRITE (iw, '(T8,A)') 'functional="'//trim(adjustl(xcfstr))//'"'
705 CALL compose(string, "z_valence", rval=pot%zion)
706 WRITE (iw, '(T8,A)') trim(string)
707 CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
708 WRITE (iw, '(T8,A)') trim(string)
709 WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
710 WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
711 lmax = -1
712 DO l = 0, lmat
713 IF (pot%nl(l) > 0) lmax = l
714 END DO
715 CALL compose(string, "l_max", ival=lmax)
716 WRITE (iw, '(T8,A)') trim(string)
717 WRITE (iw, '(T8,A)') 'l_max_rho="0"'
718 WRITE (iw, '(T8,A)') 'l_local="-3"'
719 nr = atom%basis%grid%nr
720 CALL compose(string, "mesh_size", ival=nr)
721 WRITE (iw, '(T8,A)') trim(string)
722 nwfc = sum(atom%state%maxn_occ)
723 CALL compose(string, "number_of_wfc", ival=nwfc)
724 WRITE (iw, '(T8,A)') trim(string)
725 IF (soc) THEN
726 nbeta = pot%nl(0) + 2*sum(pot%nl(1:))
727 ELSE
728 nbeta = sum(pot%nl)
729 END IF
730 CALL compose(string, "number_of_proj", ival=nbeta)
731 WRITE (iw, '(T8,A)') trim(string)//'/>'
732
733 ! Mesh
734 up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
735 WRITE (iw, '(T4,A)') '<PP_MESH'
736 WRITE (iw, '(T8,A)') 'dx="1.E+00"'
737 CALL compose(string, "mesh", ival=nr)
738 WRITE (iw, '(T8,A)') trim(string)
739 WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
740 rmax = maxval(atom%basis%grid%rad)
741 CALL compose(string, "rmax", rval=rmax)
742 WRITE (iw, '(T8,A)') trim(string)
743 WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
744 WRITE (iw, '(T8,A)') '<PP_R type="real"'
745 CALL compose(string, "size", ival=nr)
746 WRITE (iw, '(T12,A)') trim(string)
747 WRITE (iw, '(T12,A)') 'columns="4">'
748 IF (up) THEN
749 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
750 ELSE
751 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
752 END IF
753 WRITE (iw, '(T8,A)') '</PP_R>'
754 WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
755 CALL compose(string, "size", ival=nr)
756 WRITE (iw, '(T12,A)') trim(string)
757 WRITE (iw, '(T12,A)') 'columns="4">'
758 IF (up) THEN
759 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
760 ELSE
761 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
762 END IF
763 WRITE (iw, '(T8,A)') '</PP_RAB>'
764 WRITE (iw, '(T4,A)') '</PP_MESH>'
765
766 ! NLCC
767 IF (pot%nlcc) THEN
768 WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
769 CALL compose(string, "size", ival=nr)
770 WRITE (iw, '(T8,A)') trim(string)
771 WRITE (iw, '(T8,A)') 'columns="4">'
772 ALLOCATE (corden(nr))
773 corden(:) = 0.0_dp
774 CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
775 IF (up) THEN
776 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
777 ELSE
778 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
779 END IF
780 DEALLOCATE (corden)
781 WRITE (iw, '(T4,A)') '</PP_NLCC>'
782 END IF
783
784 ! local PP
785 WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
786 CALL compose(string, "size", ival=nr)
787 WRITE (iw, '(T8,A)') trim(string)
788 WRITE (iw, '(T8,A)') 'columns="4">'
789 ALLOCATE (locpot(nr))
790 locpot(:) = 0.0_dp
791 CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
792 IF (up) THEN
793 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
794 ELSE
795 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
796 END IF
797 DEALLOCATE (locpot)
798 WRITE (iw, '(T4,A)') '</PP_LOCAL>'
799
800 ! nonlocal PP
801 WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
802 ALLOCATE (rp(nr), ef(nr), beta(nr))
803 ibeta = 0
804 DO l = 0, lmat
805 IF (pot%nl(l) == 0) cycle
806 rl = pot%rcnl(l)
807 rp(:) = atom%basis%grid%rad
808 ef(:) = exp(-0.5_dp*rp*rp/(rl*rl))
809 DO i = 1, pot%nl(l)
810 pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
811 j = l + 2*i - 1
812 pf = sqrt(2._dp)/(pf*sqrt(gamma1(j)))
813 beta(:) = pf*rp**(l + 2*i - 2)*ef
814 beta(:) = beta*rp
815 !
816 ibeta = ibeta + 1
817 CALL compose(string, "<PP_BETA", counter=ibeta)
818 WRITE (iw, '(T8,A)') trim(string)
819 CALL compose(string, "angular_momentum", ival=l)
820 WRITE (iw, '(T12,A)') trim(string)
821 WRITE (iw, '(T12,A)') 'type="real"'
822 CALL compose(string, "size", ival=nr)
823 WRITE (iw, '(T12,A)') trim(string)
824 WRITE (iw, '(T12,A)') 'columns="4">'
825 IF (up) THEN
826 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
827 ELSE
828 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
829 END IF
830 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
831 WRITE (iw, '(T8,A)') trim(string)
832 IF (soc .AND. l /= 0) THEN
833 ibeta = ibeta + 1
834 CALL compose(string, "<PP_BETA", counter=ibeta)
835 WRITE (iw, '(T8,A)') trim(string)
836 CALL compose(string, "angular_momentum", ival=l)
837 WRITE (iw, '(T12,A)') trim(string)
838 WRITE (iw, '(T12,A)') 'type="real"'
839 CALL compose(string, "size", ival=nr)
840 WRITE (iw, '(T12,A)') trim(string)
841 WRITE (iw, '(T12,A)') 'columns="4">'
842 IF (up) THEN
843 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
844 ELSE
845 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
846 END IF
847 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
848 WRITE (iw, '(T8,A)') trim(string)
849 END IF
850 END DO
851 END DO
852 DEALLOCATE (rp, ef, beta)
853 ! nonlocal PP matrix elements
854 ALLOCATE (dij(nbeta, nbeta))
855 dij = 0._dp
856 IF (soc) THEN
857 ! from reverse engineering we guess: hnl, knl, hnl, knl, ...
858 n0 = pot%nl(0)
859 DO l = 0, lmat
860 IF (pot%nl(l) == 0) cycle
861 IF (l == 0) THEN
862 ! no SO term for l=0
863 dij(1:n0, 1:n0) = pot%hnl(1:n0, 1:n0, 0)
864 ELSE
865 ibeta = n0 + 2*sum(pot%nl(1:l - 1))
866 DO i = 1, pot%nl(l)
867 ii = 2*(i - 1) + 1
868 DO k = 1, pot%nl(l)
869 kk = 2*(k - 1) + 1
870 vnlp = pot%hnl(i, k, l) + 0.5_dp*l*pot%knl(i, k, l)
871 vnlm = pot%hnl(i, k, l) - 0.5_dp*(l + 1)*pot%knl(i, k, l)
872 dij(ibeta + ii, ibeta + kk) = vnlm
873 dij(ibeta + ii + 1, ibeta + kk + 1) = vnlp
874 END DO
875 END DO
876 END IF
877 END DO
878 ELSE
879 DO l = 0, lmat
880 IF (pot%nl(l) == 0) cycle
881 ibeta = sum(pot%nl(0:l - 1)) + 1
882 i = ibeta + pot%nl(l) - 1
883 dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
884 END DO
885 END IF
886 WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
887 WRITE (iw, '(T12,A)') 'columns="4">'
888 WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
889 WRITE (iw, '(T8,A)') '</PP_DIJ>'
890 DEALLOCATE (dij)
891 WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
892
893 ! atomic wavefunctions
894 ALLOCATE (beta(nr))
895 WRITE (iw, '(T4,A)') '<PP_PSWFC>'
896 nwfn = 0
897 DO l = 0, lmat
898 DO i = 1, 10
899 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
900 occa = atom%state%occupation(l, i)
901 occb = 0.0_dp
902 IF (soc .AND. l /= 0) THEN
903 occb = floor(occa/2)
904 occa = occa - occb
905 END IF
906 nwfn = nwfn + 1
907 CALL compose(string, "<PP_CHI", counter=nwfn)
908 WRITE (iw, '(T8,A)') trim(string)
909 CALL compose(string, "l", ival=l)
910 WRITE (iw, '(T12,A)') trim(string)
911 CALL compose(string, "occupation", rval=occa)
912 WRITE (iw, '(T12,A)') trim(string)
913 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
914 WRITE (iw, '(T12,A)') trim(string)
915 WRITE (iw, '(T12,A)') 'columns="4">'
916 beta = 0._dp
917 DO k = 1, atom%basis%nbas(l)
918 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
919 END DO
920 beta(:) = beta*atom%basis%grid%rad
921 IF (up) THEN
922 WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=1, nr)
923 ELSE
924 WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=nr, 1, -1)
925 END IF
926 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
927 WRITE (iw, '(T8,A)') trim(string)
928 IF (soc .AND. l /= 0) THEN
929 nwfn = nwfn + 1
930 CALL compose(string, "<PP_CHI", counter=nwfn)
931 WRITE (iw, '(T8,A)') trim(string)
932 CALL compose(string, "l", ival=l)
933 WRITE (iw, '(T12,A)') trim(string)
934 CALL compose(string, "occupation", rval=occb)
935 WRITE (iw, '(T12,A)') trim(string)
936 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
937 WRITE (iw, '(T12,A)') trim(string)
938 WRITE (iw, '(T12,A)') 'columns="4">'
939 beta = 0._dp
940 DO k = 1, atom%basis%nbas(l)
941 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
942 END DO
943 beta(:) = beta*atom%basis%grid%rad
944 IF (up) THEN
945 WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=1, nr)
946 ELSE
947 WRITE (iw, '(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j), j=nr, 1, -1)
948 END IF
949 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
950 WRITE (iw, '(T8,A)') trim(string)
951 END IF
952 END DO
953 END DO
954 WRITE (iw, '(T4,A)') '</PP_PSWFC>'
955 DEALLOCATE (beta)
956
957 ! atomic charge
958 ALLOCATE (dens(nr))
959 WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
960 CALL compose(string, "size", ival=nr)
961 WRITE (iw, '(T8,A)') trim(string)
962 WRITE (iw, '(T8,A)') 'columns="4">'
963 CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
964 "RHO", atom%basis%grid%rad)
965 IF (up) THEN
966 WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
967 ELSE
968 WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
969 END IF
970 WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
971 DEALLOCATE (dens)
972
973 ! PP SOC information
974 IF (soc) THEN
975 WRITE (iw, '(T4,A)') '<PP_SPIN_ORB>'
976! <PP_RELWFC.1 index="1" els="6S" nn="1" lchi="0" jchi="5.000000000000e-1" oc="1.000000000000e0"/>
977 nwfn = 0
978 DO l = 0, lmat
979 DO i = 1, 10
980 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
981 occa = atom%state%occupation(l, i)
982 occb = 0.0_dp
983 IF (l /= 0) THEN
984 occb = floor(occa/2)
985 occa = occa - occb
986 END IF
987 nwfn = nwfn + 1
988 CALL compose(string, "<PP_RELWFC", counter=nwfn)
989 WRITE (iw, '(T8,A)') trim(string)
990 CALL compose(string, "index", ival=nwfn)
991 WRITE (iw, '(T12,A)') trim(string)
992 !!CALL compose(string, "els", cval=xxxx)
993 !!WRITE (iw, '(T12,A)') TRIM(string)
994 CALL compose(string, "nn", ival=nwfn)
995 WRITE (iw, '(T12,A)') trim(string)
996 CALL compose(string, "lchi", ival=l)
997 WRITE (iw, '(T12,A)') trim(string)
998 IF (l == 0) THEN
999 rlp = 0.5_dp
1000 ELSE
1001 rlp = l - 0.5_dp
1002 END IF
1003 CALL compose(string, "jchi", rval=rlp)
1004 WRITE (iw, '(T12,A)') trim(string)
1005 CALL compose(string, "oc", rval=occa)
1006 WRITE (iw, '(T12,A)') trim(string)
1007 WRITE (iw, '(T12,A)') '/>'
1008 IF (l /= 0) THEN
1009 nwfn = nwfn + 1
1010 CALL compose(string, "<PP_RELWFC", counter=nwfn)
1011 WRITE (iw, '(T8,A)') trim(string)
1012 CALL compose(string, "index", ival=nwfn)
1013 WRITE (iw, '(T12,A)') trim(string)
1014 !!CALL compose(string, "els", cval=xxxx)
1015 !!WRITE (iw, '(T12,A)') TRIM(string)
1016 CALL compose(string, "nn", ival=nwfn)
1017 WRITE (iw, '(T12,A)') trim(string)
1018 CALL compose(string, "lchi", ival=l)
1019 WRITE (iw, '(T12,A)') trim(string)
1020 rlp = l + 0.5_dp
1021 CALL compose(string, "jchi", rval=rlp)
1022 WRITE (iw, '(T12,A)') trim(string)
1023 CALL compose(string, "oc", rval=occa)
1024 WRITE (iw, '(T12,A)') trim(string)
1025 WRITE (iw, '(T12,A)') '/>'
1026 END IF
1027 END DO
1028 END DO
1029 ibeta = 0
1030 DO l = 0, lmat
1031 IF (pot%nl(l) == 0) cycle
1032 DO i = 1, pot%nl(l)
1033 ibeta = ibeta + 1
1034 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1035 WRITE (iw, '(T8,A)') trim(string)
1036 CALL compose(string, "index", ival=ibeta)
1037 WRITE (iw, '(T12,A)') trim(string)
1038 CALL compose(string, "lll", ival=l)
1039 WRITE (iw, '(T12,A)') trim(string)
1040 IF (l == 0) THEN
1041 rlp = 0.5_dp
1042 ELSE
1043 rlp = l - 0.5_dp
1044 END IF
1045 CALL compose(string, "jjj", rval=rlp)
1046 WRITE (iw, '(T12,A)') trim(string)
1047 WRITE (iw, '(T12,A)') '/>'
1048 IF (l > 0) THEN
1049 ibeta = ibeta + 1
1050 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1051 WRITE (iw, '(T8,A)') trim(string)
1052 CALL compose(string, "index", ival=ibeta)
1053 WRITE (iw, '(T12,A)') trim(string)
1054 CALL compose(string, "lll", ival=l)
1055 WRITE (iw, '(T12,A)') trim(string)
1056 rlp = l + 0.5_dp
1057 CALL compose(string, "jjj", rval=rlp)
1058 WRITE (iw, '(T12,A)') trim(string)
1059 WRITE (iw, '(T12,A)') '/>'
1060 END IF
1061 END DO
1062 END DO
1063 WRITE (iw, '(T4,A)') '</PP_SPIN_ORB>'
1064 END IF
1065
1066 WRITE (iw, '(A)') '</UPF>'
1067
1068 END SUBROUTINE atom_write_upf
1069
1070! **************************************************************************************************
1071!> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
1072!> \param string composed string
1073!> \param tag attribute tag
1074!> \param counter counter
1075!> \param rval real variable
1076!> \param ival integer variable
1077!> \param cval string variable
1078!> \param isfinal close the current XML element if this is the last attribute
1079!> \par History
1080!> * 09.2012 created [Juerg Hutter]
1081! **************************************************************************************************
1082 SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
1083 CHARACTER(len=*), INTENT(OUT) :: string
1084 CHARACTER(len=*), INTENT(IN) :: tag
1085 INTEGER, INTENT(IN), OPTIONAL :: counter
1086 REAL(kind=dp), INTENT(IN), OPTIONAL :: rval
1087 INTEGER, INTENT(IN), OPTIONAL :: ival
1088 CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
1089 LOGICAL, INTENT(IN), OPTIONAL :: isfinal
1090
1091 CHARACTER(LEN=default_string_length) :: str
1092 LOGICAL :: fin
1093
1094 IF (PRESENT(counter)) THEN
1095 WRITE (str, "(I12)") counter
1096 ELSEIF (PRESENT(rval)) THEN
1097 WRITE (str, "(G18.8)") rval
1098 ELSEIF (PRESENT(ival)) THEN
1099 WRITE (str, "(I12)") ival
1100 ELSEIF (PRESENT(cval)) THEN
1101 WRITE (str, "(A)") trim(adjustl(cval))
1102 ELSE
1103 WRITE (str, "(A)") ""
1104 END IF
1105 fin = .false.
1106 IF (PRESENT(isfinal)) fin = isfinal
1107 IF (PRESENT(counter)) THEN
1108 IF (fin) THEN
1109 WRITE (string, "(A,A1,A,A1)") trim(adjustl(tag)), '.', trim(adjustl(str)), '>'
1110 ELSE
1111 WRITE (string, "(A,A1,A)") trim(adjustl(tag)), '.', trim(adjustl(str))
1112 END IF
1113 ELSE
1114 IF (fin) THEN
1115 WRITE (string, "(A,A2,A,A2)") trim(adjustl(tag)), '="', trim(adjustl(str)), '>"'
1116 ELSE
1117 WRITE (string, "(A,A2,A,A1)") trim(adjustl(tag)), '="', trim(adjustl(str)), '"'
1118 END IF
1119 END IF
1120 END SUBROUTINE compose
1121
1122END 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:2398
atom_types::create_atom_type
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:947
atom_types::release_atom_type
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:971
atom_types::release_atom_potential
subroutine, public release_atom_potential(potential)
...
Definition atom_types.F:2524
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:1012
atom_types::init_atom_potential
subroutine, public init_atom_potential(potential, potential_section, zval)
...
Definition atom_types.F:2424
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:913
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1058
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::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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
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