(git:b310db4)
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-2026 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 REAL(kind=dp), PARAMETER :: error_threshold = 1.0e-12_dp
534
535 INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
536 s1, s2
537 REAL(kind=dp) :: dene, emax, expzet, fhomo, o, prefac, &
538 zeta
539 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
540 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
541 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
542 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp
543 TYPE(atom_state), POINTER :: state
544
545 WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
546
547 state => atom%state
548 ovlp => atom%integrals%ovlp
549
550 ! find HOMO
551 lhomo = -1
552 nhomo = -1
553 emax = -huge(1._dp)
554 DO l = 0, state%maxl_occ
555 DO i = 1, state%maxn_occ(l)
556 IF (atom%orbitals%ener(i, l) > emax) THEN
557 lhomo = l
558 nhomo = i
559 emax = atom%orbitals%ener(i, l)
560 fhomo = state%occupation(l, i)
561 END IF
562 END DO
563 END DO
564
565 s1 = SIZE(atom%orbitals%wfn, 1)
566 s2 = SIZE(atom%orbitals%wfn, 2)
567 ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
568 s2 = maxval(state%maxn_occ) + nder
569 ALLOCATE (rbasis(s1, s2, 0:lmat))
570 rbasis = 0._dp
571
572 DO ider = -nder, nder
573 dene = real(ider, kind=dp)*delta
574 cpassert(fhomo > abs(dene))
575 state%occupation(lhomo, nhomo) = fhomo + dene
576 CALL calculate_atom(atom, iw=0, noguess=.true.)
577 wfn(:, :, :, ider) = atom%orbitals%wfn
578 state%occupation(lhomo, nhomo) = fhomo
579 END DO
580 ! reset
581 CALL calculate_atom(atom, iw=0, noguess=.true.)
582
583 DO l = 0, state%maxl_occ
584 ! occupied states
585 DO i = 1, max(state%maxn_occ(l), 1)
586 rbasis(:, i, l) = wfn(:, i, l, 0)
587 END DO
588 ! differentiation
589 DO ider = 1, nder
590 i = max(state%maxn_occ(l), 1)
591 SELECT CASE (ider)
592 CASE (1)
593 rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
594 CASE (2)
595 rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
596 CASE (3)
597 rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
598 + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
599 CASE DEFAULT
600 cpabort("Only 1, 2, 3 are supported as the number of wavefunction derivatives")
601 END SELECT
602 END DO
603
604 ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
605 n = state%maxn_occ(l) + nder
606 m = atom%basis%nbas(l)
607 DO i = 1, n
608 DO j = 1, i - 1
609 o = dot_product(rbasis(1:m, j, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
610 rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
611 END DO
612 o = dot_product(rbasis(1:m, i, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
613 rbasis(1:m, i, l) = rbasis(1:m, i, l)/sqrt(o)
614 END DO
615
616 ! check
617 ALLOCATE (amat(n, n))
618 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)))
619 DO i = 1, n
620 amat(i, i) = amat(i, i) - 1._dp
621 END DO
622 IF (maxval(abs(amat)) > error_threshold) THEN
623 WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", maxval(abs(amat))
624 END IF
625 DEALLOCATE (amat)
626
627 ! Quickstep normalization
628 WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
629
630 WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
631 n - nder, " valence ", nder, " response"
632 expzet = 0.25_dp*real(2*l + 3, dp)
633 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
634 DO i = 1, m
635 zeta = (2._dp*atom%basis%am(i, l))**expzet
636 WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
637 END DO
638
639 END DO
640
641 DEALLOCATE (wfn, rbasis)
642
643 WRITE (iw, '(" ",79("*"))')
644
645 END SUBROUTINE atom_response_basis
646
647! **************************************************************************************************
648!> \brief Write pseudo-potential in Quantum Espresso UPF format.
649!> \param atom information about the atomic kind
650!> \param iw output file unit
651!> \param xcfstr ...
652!> \par History
653!> * 09.2012 created [Juerg Hutter]
654! **************************************************************************************************
655 SUBROUTINE atom_write_upf(atom, iw, xcfstr)
656
657 TYPE(atom_type), POINTER :: atom
658 INTEGER, INTENT(IN) :: iw
659 CHARACTER(LEN=*), INTENT(IN) :: xcfstr
660
661 CHARACTER(LEN=default_string_length) :: string
662 INTEGER :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
663 nbeta, nr, nwfc, nwfn
664 LOGICAL :: soc, up
665 REAL(kind=dp) :: occa, occb, pf, rhoatom, rl, rlp, rmax, &
666 vnlm, vnlp
667 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
668 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
669 TYPE(atom_gthpot_type), POINTER :: pot
670
671 IF (.NOT. atom%pp_calc) RETURN
672 IF (atom%potential%ppot_type /= gth_pseudo) RETURN
673 pot => atom%potential%gth_pot
674 cpassert(.NOT. pot%lsdpot)
675 soc = pot%soc
676
677 WRITE (iw, '(A)') '<UPF version="2.0.1">'
678
679 WRITE (iw, '(T4,A)') '<PP_INFO>'
680 WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
681 WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
682 CALL atom_write_pseudo_param(pot, iw)
683 WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
684 WRITE (iw, '(T4,A)') '</PP_INFO>'
685 WRITE (iw, '(T4,A)') '<PP_HEADER'
686 WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
687 WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
688 WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
689 WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
690 CALL compose(string, "element", cval=pot%symbol)
691 WRITE (iw, '(T8,A)') trim(string)
692 WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
693 WRITE (iw, '(T8,A)') 'relativistic="no"'
694 WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
695 WRITE (iw, '(T8,A)') 'is_paw="F"'
696 WRITE (iw, '(T8,A)') 'is_coulomb="F"'
697 IF (soc) THEN
698 WRITE (iw, '(T8,A)') 'has_so="T"'
699 ELSE
700 WRITE (iw, '(T8,A)') 'has_so="F"'
701 END IF
702 WRITE (iw, '(T8,A)') 'has_wfc="F"'
703 WRITE (iw, '(T8,A)') 'has_gipaw="F"'
704 WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
705 IF (pot%nlcc) THEN
706 WRITE (iw, '(T8,A)') 'core_correction="T"'
707 ELSE
708 WRITE (iw, '(T8,A)') 'core_correction="F"'
709 END IF
710 WRITE (iw, '(T8,A)') 'functional="'//trim(adjustl(xcfstr))//'"'
711 CALL compose(string, "z_valence", rval=pot%zion)
712 WRITE (iw, '(T8,A)') trim(string)
713 CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
714 WRITE (iw, '(T8,A)') trim(string)
715 WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
716 WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
717 lmax = -1
718 DO l = 0, lmat
719 IF (pot%nl(l) > 0) lmax = l
720 END DO
721 CALL compose(string, "l_max", ival=lmax)
722 WRITE (iw, '(T8,A)') trim(string)
723 WRITE (iw, '(T8,A)') 'l_max_rho="0"'
724 WRITE (iw, '(T8,A)') 'l_local="-3"'
725 nr = atom%basis%grid%nr
726 CALL compose(string, "mesh_size", ival=nr)
727 WRITE (iw, '(T8,A)') trim(string)
728 nwfc = sum(atom%state%maxn_occ)
729 CALL compose(string, "number_of_wfc", ival=nwfc)
730 WRITE (iw, '(T8,A)') trim(string)
731 IF (soc) THEN
732 nbeta = pot%nl(0) + 2*sum(pot%nl(1:))
733 ELSE
734 nbeta = sum(pot%nl)
735 END IF
736 CALL compose(string, "number_of_proj", ival=nbeta)
737 WRITE (iw, '(T8,A)') trim(string)//'/>'
738
739 ! Mesh
740 up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
741 WRITE (iw, '(T4,A)') '<PP_MESH'
742 WRITE (iw, '(T8,A)') 'dx="1.E+00"'
743 CALL compose(string, "mesh", ival=nr)
744 WRITE (iw, '(T8,A)') trim(string)
745 WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
746 rmax = maxval(atom%basis%grid%rad)
747 CALL compose(string, "rmax", rval=rmax)
748 WRITE (iw, '(T8,A)') trim(string)
749 WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
750 WRITE (iw, '(T8,A)') '<PP_R type="real"'
751 CALL compose(string, "size", ival=nr)
752 WRITE (iw, '(T12,A)') trim(string)
753 WRITE (iw, '(T12,A)') 'columns="4">'
754 IF (up) THEN
755 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
756 ELSE
757 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
758 END IF
759 WRITE (iw, '(T8,A)') '</PP_R>'
760 WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
761 CALL compose(string, "size", ival=nr)
762 WRITE (iw, '(T12,A)') trim(string)
763 WRITE (iw, '(T12,A)') 'columns="4">'
764 IF (up) THEN
765 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
766 ELSE
767 WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
768 END IF
769 WRITE (iw, '(T8,A)') '</PP_RAB>'
770 WRITE (iw, '(T4,A)') '</PP_MESH>'
771
772 ! NLCC
773 IF (pot%nlcc) THEN
774 WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
775 CALL compose(string, "size", ival=nr)
776 WRITE (iw, '(T8,A)') trim(string)
777 WRITE (iw, '(T8,A)') 'columns="4">'
778 ALLOCATE (corden(nr))
779 corden(:) = 0.0_dp
780 CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
781 IF (up) THEN
782 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
783 ELSE
784 WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
785 END IF
786 DEALLOCATE (corden)
787 WRITE (iw, '(T4,A)') '</PP_NLCC>'
788 END IF
789
790 ! local PP
791 WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
792 CALL compose(string, "size", ival=nr)
793 WRITE (iw, '(T8,A)') trim(string)
794 WRITE (iw, '(T8,A)') 'columns="4">'
795 ALLOCATE (locpot(nr))
796 locpot(:) = 0.0_dp
797 CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
798 IF (up) THEN
799 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
800 ELSE
801 WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
802 END IF
803 DEALLOCATE (locpot)
804 WRITE (iw, '(T4,A)') '</PP_LOCAL>'
805
806 ! nonlocal PP
807 WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
808 ALLOCATE (rp(nr), ef(nr), beta(nr))
809 ibeta = 0
810 DO l = 0, lmat
811 IF (pot%nl(l) == 0) cycle
812 rl = pot%rcnl(l)
813 rp(:) = atom%basis%grid%rad
814 ef(:) = exp(-0.5_dp*rp*rp/(rl*rl))
815 DO i = 1, pot%nl(l)
816 pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
817 j = l + 2*i - 1
818 pf = sqrt(2._dp)/(pf*sqrt(gamma1(j)))
819 beta(:) = pf*rp**(l + 2*i - 2)*ef
820 beta(:) = beta*rp
821 !
822 ibeta = ibeta + 1
823 CALL compose(string, "<PP_BETA", counter=ibeta)
824 WRITE (iw, '(T8,A)') trim(string)
825 CALL compose(string, "angular_momentum", ival=l)
826 WRITE (iw, '(T12,A)') trim(string)
827 WRITE (iw, '(T12,A)') 'type="real"'
828 CALL compose(string, "size", ival=nr)
829 WRITE (iw, '(T12,A)') trim(string)
830 WRITE (iw, '(T12,A)') 'columns="4">'
831 IF (up) THEN
832 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
833 ELSE
834 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
835 END IF
836 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
837 WRITE (iw, '(T8,A)') trim(string)
838 IF (soc .AND. l /= 0) THEN
839 ibeta = ibeta + 1
840 CALL compose(string, "<PP_BETA", counter=ibeta)
841 WRITE (iw, '(T8,A)') trim(string)
842 CALL compose(string, "angular_momentum", ival=l)
843 WRITE (iw, '(T12,A)') trim(string)
844 WRITE (iw, '(T12,A)') 'type="real"'
845 CALL compose(string, "size", ival=nr)
846 WRITE (iw, '(T12,A)') trim(string)
847 WRITE (iw, '(T12,A)') 'columns="4">'
848 IF (up) THEN
849 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
850 ELSE
851 WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
852 END IF
853 CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.true.)
854 WRITE (iw, '(T8,A)') trim(string)
855 END IF
856 END DO
857 END DO
858 DEALLOCATE (rp, ef, beta)
859 ! nonlocal PP matrix elements
860 ALLOCATE (dij(nbeta, nbeta))
861 dij = 0._dp
862 IF (soc) THEN
863 ! from reverse engineering we guess: hnl, knl, hnl, knl, ...
864 n0 = pot%nl(0)
865 DO l = 0, lmat
866 IF (pot%nl(l) == 0) cycle
867 IF (l == 0) THEN
868 ! no SO term for l=0
869 dij(1:n0, 1:n0) = pot%hnl(1:n0, 1:n0, 0)
870 ELSE
871 ibeta = n0 + 2*sum(pot%nl(1:l - 1))
872 DO i = 1, pot%nl(l)
873 ii = 2*(i - 1) + 1
874 DO k = 1, pot%nl(l)
875 kk = 2*(k - 1) + 1
876 vnlp = pot%hnl(i, k, l) + 0.5_dp*l*pot%knl(i, k, l)
877 vnlm = pot%hnl(i, k, l) - 0.5_dp*(l + 1)*pot%knl(i, k, l)
878 dij(ibeta + ii, ibeta + kk) = vnlm
879 dij(ibeta + ii + 1, ibeta + kk + 1) = vnlp
880 END DO
881 END DO
882 END IF
883 END DO
884 ELSE
885 DO l = 0, lmat
886 IF (pot%nl(l) == 0) cycle
887 ibeta = sum(pot%nl(0:l - 1)) + 1
888 i = ibeta + pot%nl(l) - 1
889 dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
890 END DO
891 END IF
892 WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
893 WRITE (iw, '(T12,A)') 'columns="4">'
894 WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
895 WRITE (iw, '(T8,A)') '</PP_DIJ>'
896 DEALLOCATE (dij)
897 WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
898
899 ! atomic wavefunctions
900 ALLOCATE (beta(nr))
901 WRITE (iw, '(T4,A)') '<PP_PSWFC>'
902 nwfn = 0
903 DO l = 0, lmat
904 DO i = 1, 10
905 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
906 occa = atom%state%occupation(l, i)
907 occb = 0.0_dp
908 IF (soc .AND. l /= 0) THEN
909 occb = floor(occa/2)
910 occa = occa - occb
911 END IF
912 nwfn = nwfn + 1
913 CALL compose(string, "<PP_CHI", counter=nwfn)
914 WRITE (iw, '(T8,A)') trim(string)
915 CALL compose(string, "l", ival=l)
916 WRITE (iw, '(T12,A)') trim(string)
917 CALL compose(string, "occupation", rval=occa)
918 WRITE (iw, '(T12,A)') trim(string)
919 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
920 WRITE (iw, '(T12,A)') trim(string)
921 WRITE (iw, '(T12,A)') 'columns="4">'
922 beta = 0._dp
923 DO k = 1, atom%basis%nbas(l)
924 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
925 END DO
926 beta(:) = beta*atom%basis%grid%rad
927 IF (up) THEN
928 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
929 ELSE
930 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
931 END IF
932 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
933 WRITE (iw, '(T8,A)') trim(string)
934 IF (soc .AND. l /= 0) THEN
935 nwfn = nwfn + 1
936 CALL compose(string, "<PP_CHI", counter=nwfn)
937 WRITE (iw, '(T8,A)') trim(string)
938 CALL compose(string, "l", ival=l)
939 WRITE (iw, '(T12,A)') trim(string)
940 CALL compose(string, "occupation", rval=occb)
941 WRITE (iw, '(T12,A)') trim(string)
942 CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
943 WRITE (iw, '(T12,A)') trim(string)
944 WRITE (iw, '(T12,A)') 'columns="4">'
945 beta = 0._dp
946 DO k = 1, atom%basis%nbas(l)
947 beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
948 END DO
949 beta(:) = beta*atom%basis%grid%rad
950 IF (up) THEN
951 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
952 ELSE
953 WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
954 END IF
955 CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.true.)
956 WRITE (iw, '(T8,A)') trim(string)
957 END IF
958 END DO
959 END DO
960 WRITE (iw, '(T4,A)') '</PP_PSWFC>'
961 DEALLOCATE (beta)
962
963 ! atomic charge
964 ALLOCATE (dens(nr))
965 WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
966 CALL compose(string, "size", ival=nr)
967 WRITE (iw, '(T8,A)') trim(string)
968 WRITE (iw, '(T8,A)') 'columns="4">'
969 CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
970 "RHO", atom%basis%grid%rad)
971 IF (up) THEN
972 WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
973 ELSE
974 WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
975 END IF
976 WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
977
978 ! Check atomic electron count
979 rhoatom = sum(fourpi*atom%basis%grid%wr(:)*dens(:))
980 IF (abs(pot%zion - rhoatom) > 0.01_dp) THEN
981 CALL cp_warn(__location__, "Valence charge and electron count differ by more than 0.01")
982 END IF
983
984 DEALLOCATE (dens)
985 ! PP SOC information
986 IF (soc) THEN
987 WRITE (iw, '(T4,A)') '<PP_SPIN_ORB>'
988! <PP_RELWFC.1 index="1" els="6S" nn="1" lchi="0" jchi="5.000000000000e-1" oc="1.000000000000e0"/>
989 nwfn = 0
990 DO l = 0, lmat
991 DO i = 1, 10
992 IF (abs(atom%state%occupation(l, i)) == 0._dp) cycle
993 occa = atom%state%occupation(l, i)
994 occb = 0.0_dp
995 IF (l /= 0) THEN
996 occb = floor(occa/2)
997 occa = occa - occb
998 END IF
999 nwfn = nwfn + 1
1000 CALL compose(string, "<PP_RELWFC", counter=nwfn)
1001 WRITE (iw, '(T8,A)') trim(string)
1002 CALL compose(string, "index", ival=nwfn)
1003 WRITE (iw, '(T12,A)') trim(string)
1004 !!CALL compose(string, "els", cval=xxxx)
1005 !!WRITE (iw, '(T12,A)') TRIM(string)
1006 CALL compose(string, "nn", ival=nwfn)
1007 WRITE (iw, '(T12,A)') trim(string)
1008 CALL compose(string, "lchi", ival=l)
1009 WRITE (iw, '(T12,A)') trim(string)
1010 IF (l == 0) THEN
1011 rlp = 0.5_dp
1012 ELSE
1013 rlp = l - 0.5_dp
1014 END IF
1015 CALL compose(string, "jchi", rval=rlp)
1016 WRITE (iw, '(T12,A)') trim(string)
1017 CALL compose(string, "oc", rval=occa)
1018 WRITE (iw, '(T12,A)') trim(string)
1019 WRITE (iw, '(T12,A)') '/>'
1020 IF (l /= 0) THEN
1021 nwfn = nwfn + 1
1022 CALL compose(string, "<PP_RELWFC", counter=nwfn)
1023 WRITE (iw, '(T8,A)') trim(string)
1024 CALL compose(string, "index", ival=nwfn)
1025 WRITE (iw, '(T12,A)') trim(string)
1026 !!CALL compose(string, "els", cval=xxxx)
1027 !!WRITE (iw, '(T12,A)') TRIM(string)
1028 CALL compose(string, "nn", ival=nwfn)
1029 WRITE (iw, '(T12,A)') trim(string)
1030 CALL compose(string, "lchi", ival=l)
1031 WRITE (iw, '(T12,A)') trim(string)
1032 rlp = l + 0.5_dp
1033 CALL compose(string, "jchi", rval=rlp)
1034 WRITE (iw, '(T12,A)') trim(string)
1035 CALL compose(string, "oc", rval=occa)
1036 WRITE (iw, '(T12,A)') trim(string)
1037 WRITE (iw, '(T12,A)') '/>'
1038 END IF
1039 END DO
1040 END DO
1041 ibeta = 0
1042 DO l = 0, lmat
1043 IF (pot%nl(l) == 0) cycle
1044 DO i = 1, pot%nl(l)
1045 ibeta = ibeta + 1
1046 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1047 WRITE (iw, '(T8,A)') trim(string)
1048 CALL compose(string, "index", ival=ibeta)
1049 WRITE (iw, '(T12,A)') trim(string)
1050 CALL compose(string, "lll", ival=l)
1051 WRITE (iw, '(T12,A)') trim(string)
1052 IF (l == 0) THEN
1053 rlp = 0.5_dp
1054 ELSE
1055 rlp = l - 0.5_dp
1056 END IF
1057 CALL compose(string, "jjj", rval=rlp)
1058 WRITE (iw, '(T12,A)') trim(string)
1059 WRITE (iw, '(T12,A)') '/>'
1060 IF (l > 0) THEN
1061 ibeta = ibeta + 1
1062 CALL compose(string, "<PP_RELBETA", counter=ibeta)
1063 WRITE (iw, '(T8,A)') trim(string)
1064 CALL compose(string, "index", ival=ibeta)
1065 WRITE (iw, '(T12,A)') trim(string)
1066 CALL compose(string, "lll", ival=l)
1067 WRITE (iw, '(T12,A)') trim(string)
1068 rlp = l + 0.5_dp
1069 CALL compose(string, "jjj", rval=rlp)
1070 WRITE (iw, '(T12,A)') trim(string)
1071 WRITE (iw, '(T12,A)') '/>'
1072 END IF
1073 END DO
1074 END DO
1075 WRITE (iw, '(T4,A)') '</PP_SPIN_ORB>'
1076 END IF
1077
1078 WRITE (iw, '(A)') '</UPF>'
1079
1080 END SUBROUTINE atom_write_upf
1081
1082! **************************************************************************************************
1083!> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
1084!> \param string composed string
1085!> \param tag attribute tag
1086!> \param counter counter
1087!> \param rval real variable
1088!> \param ival integer variable
1089!> \param cval string variable
1090!> \param isfinal close the current XML element if this is the last attribute
1091!> \par History
1092!> * 09.2012 created [Juerg Hutter]
1093! **************************************************************************************************
1094 SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
1095 CHARACTER(len=*), INTENT(OUT) :: string
1096 CHARACTER(len=*), INTENT(IN) :: tag
1097 INTEGER, INTENT(IN), OPTIONAL :: counter
1098 REAL(kind=dp), INTENT(IN), OPTIONAL :: rval
1099 INTEGER, INTENT(IN), OPTIONAL :: ival
1100 CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
1101 LOGICAL, INTENT(IN), OPTIONAL :: isfinal
1102
1103 CHARACTER(LEN=default_string_length) :: str
1104 LOGICAL :: fin
1105
1106 IF (PRESENT(counter)) THEN
1107 WRITE (str, "(I12)") counter
1108 ELSE IF (PRESENT(rval)) THEN
1109 WRITE (str, "(G18.8)") rval
1110 ELSE IF (PRESENT(ival)) THEN
1111 WRITE (str, "(I12)") ival
1112 ELSE IF (PRESENT(cval)) THEN
1113 WRITE (str, "(A)") trim(adjustl(cval))
1114 ELSE
1115 WRITE (str, "(A)") ""
1116 END IF
1117 fin = .false.
1118 IF (PRESENT(isfinal)) fin = isfinal
1119 IF (PRESENT(counter)) THEN
1120 IF (fin) THEN
1121 WRITE (string, "(A,A1,A,A1)") trim(adjustl(tag)), '.', trim(adjustl(str)), '>'
1122 ELSE
1123 WRITE (string, "(A,A1,A)") trim(adjustl(tag)), '.', trim(adjustl(str))
1124 END IF
1125 ELSE
1126 IF (fin) THEN
1127 WRITE (string, "(A,A2,A,A2)") trim(adjustl(tag)), '="', trim(adjustl(str)), '>"'
1128 ELSE
1129 WRITE (string, "(A,A2,A,A1)") trim(adjustl(tag)), '="', trim(adjustl(str)), '"'
1130 END IF
1131 END IF
1132 END SUBROUTINE compose
1133
1134END 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:142
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:1892
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:859
atom_output::atom_print_basis
subroutine, public atom_print_basis(atom_basis, iw, title)
Print atomic basis set.
Definition atom_output.F:301
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:514
atom_output::atom_write_pseudo_param
subroutine, public atom_write_pseudo_param(gthpot, iunit, fopt)
Print GTH pseudo-potential parameters.
Definition atom_output.F:769
atom_output::atom_print_potential
subroutine, public atom_print_potential(potential, iw)
Print information about the pseudo-potential.
Definition atom_output.F:647
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:204
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:2410
atom_types::create_atom_type
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:968
atom_types::release_atom_type
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:992
atom_types::release_atom_potential
subroutine, public release_atom_potential(potential)
...
Definition atom_types.F:2537
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:384
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:68
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:1033
atom_types::init_atom_potential
subroutine, public init_atom_potential(potential, potential_section, zval)
...
Definition atom_types.F:2436
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:934
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1079
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:2334
atom_utils::atom_local_potential
pure subroutine, public atom_local_potential(locpot, gthpot, rr)
...
Definition atom_utils.F:816
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:604
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:2213
atom_utils::atom_core_density
subroutine, public atom_core_density(corden, potential, typ, rr)
...
Definition atom_utils.F:689
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:300
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:365
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:515
atom_utils::atom_completeness
subroutine, public atom_completeness(basis, zv, iw)
Estimate completeness of the given atomic basis set.
Definition atom_utils.F:2369
atom_utils::atom_write_zmp_restart
subroutine, public atom_write_zmp_restart(atom)
ZMP subroutine to write external restart file.
Definition atom_utils.F:424
atom_utils::atom_set_occupation
subroutine, public atom_set_occupation(ostring, occupation, wfnocc, multiplicity)
Set occupation of atomic orbitals.
Definition atom_utils.F:102
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:280
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:146
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:77
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:871
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:1091
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::xc_funct_tpss
integer, parameter, public xc_funct_tpss
Definition input_constants.F:593
input_constants::xc_funct_bp
integer, parameter, public xc_funct_bp
Definition input_constants.F:593
input_constants::xc_funct_olyp
integer, parameter, public xc_funct_olyp
Definition input_constants.F:593
input_constants::xc_funct_pbe
integer, parameter, public xc_funct_pbe
Definition input_constants.F:593
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:160
input_constants::xc_funct_no_shortcut
integer, parameter, public xc_funct_no_shortcut
Definition input_constants.F:593
input_constants::xc_funct_pbe0
integer, parameter, public xc_funct_pbe0
Definition input_constants.F:593
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:185
input_constants::xc_funct_beefvdw
integer, parameter, public xc_funct_beefvdw
Definition input_constants.F:593
input_constants::xc_funct_pade
integer, parameter, public xc_funct_pade
Definition input_constants.F:593
input_constants::xc_funct_blyp
integer, parameter, public xc_funct_blyp
Definition input_constants.F:593
input_constants::xc_funct_b3lyp
integer, parameter, public xc_funct_b3lyp
Definition input_constants.F:593
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:593
input_constants::xc_funct_hcth120
integer, parameter, public xc_funct_hcth120
Definition input_constants.F:593
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:735
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:708
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:1052
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:79
atom_types::atom_gthpot_type
Provides all information about a pseudopotential.
Definition atom_types.F:99
atom_types::atom_integrals
Definition atom_types.F:217
atom_types::atom_optimization_type
Information on optimization procedure.
Definition atom_types.F:284
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:238
atom_types::atom_p_type
Definition atom_types.F:324
atom_types::atom_potential_type
Definition atom_types.F:172
atom_types::atom_state
Provides all information on states and occupation.
Definition atom_types.F:199
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:294
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