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