Line data Source code
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 :
8 : MODULE atom_energy
9 : USE atom_admm_methods, ONLY: atom_admm
10 : USE atom_electronic_structure, ONLY: calculate_atom
11 : USE atom_fit, ONLY: atom_fit_density,&
12 : atom_fit_kgpot
13 : USE atom_grb, ONLY: atom_grb_construction
14 : USE atom_operators, ONLY: atom_int_release,&
15 : atom_int_setup,&
16 : atom_ppint_release,&
17 : atom_ppint_setup,&
18 : atom_relint_release,&
19 : atom_relint_setup
20 : USE atom_output, ONLY: atom_print_basis,&
21 : atom_print_info,&
22 : atom_print_method,&
23 : atom_print_orbitals,&
24 : atom_print_potential,&
25 : atom_write_pseudo_param
26 : USE atom_sgp, ONLY: atom_sgp_construction
27 : USE atom_types, ONLY: &
28 : atom_basis_type, atom_gthpot_type, atom_integrals, atom_optimization_type, atom_orbitals, &
29 : atom_p_type, atom_potential_type, atom_state, atom_type, create_atom_orbs, &
30 : create_atom_type, gth_pseudo, init_atom_basis, init_atom_potential, lmat, &
31 : read_atom_opt_section, release_atom_basis, release_atom_potential, release_atom_type, &
32 : set_atom
33 : USE atom_utils, ONLY: &
34 : atom_completeness, atom_condnumber, atom_consistent_method, atom_core_density, &
35 : atom_density, atom_local_potential, atom_read_external_density, atom_read_external_vxc, &
36 : atom_set_occupation, atom_write_zmp_restart, get_maxl_occ, get_maxn_occ
37 : USE atom_xc, ONLY: calculate_atom_ext_vxc,&
38 : calculate_atom_zmp
39 : USE cp_log_handling, ONLY: cp_get_default_logger,&
40 : cp_logger_type
41 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
42 : cp_print_key_unit_nr
43 : USE input_constants, ONLY: &
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 :
68 : CONTAINS
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 326 : 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 326 : 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 326 : 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 326 : 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 326 : 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 326 : CALL timeset(routineN, handle)
115 :
116 : ! What atom do we calculate
117 326 : CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
118 326 : CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
119 326 : zz = 0
120 11732 : DO i = 1, nelem
121 11732 : IF (ptable(i)%symbol == elem) THEN
122 : zz = i
123 : EXIT
124 : END IF
125 : END DO
126 326 : IF (zz /= 1) zval = zz
127 :
128 : ! read and set up inofrmation on the basis sets
129 12062 : ALLOCATE (ae_basis, pp_basis)
130 326 : basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
131 326 : NULLIFY (ae_basis%grid)
132 326 : CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
133 326 : NULLIFY (pp_basis%grid)
134 326 : basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
135 326 : CALL init_atom_basis(pp_basis, basis_section, zval, "PP")
136 :
137 : ! print general and basis set information
138 326 : logger => cp_get_default_logger()
139 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
140 326 : IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw)
141 326 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
142 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
143 326 : IF (iw > 0) THEN
144 0 : CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
145 0 : CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
146 : END IF
147 326 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
148 :
149 : ! read and setup information on the pseudopotential
150 326 : NULLIFY (potential_section)
151 326 : potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
152 3506782 : ALLOCATE (ae_pot, p_pot)
153 326 : CALL init_atom_potential(p_pot, potential_section, zval)
154 326 : CALL init_atom_potential(ae_pot, potential_section, -1)
155 :
156 : ! if the ERI's are calculated analytically, we have to precalculate them
157 326 : eri_c = .FALSE.
158 326 : CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
159 326 : IF (do_eric == do_analytic) eri_c = .TRUE.
160 326 : eri_e = .FALSE.
161 326 : CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
162 326 : IF (do_erie == do_analytic) eri_e = .TRUE.
163 326 : CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
164 326 : CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
165 :
166 : ! information on the states to be calculated
167 326 : CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
168 326 : maxn = 0
169 326 : CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
170 688 : DO in = 1, MIN(SIZE(cn), 4)
171 688 : maxn(in - 1) = cn(in)
172 : END DO
173 2282 : DO in = 0, lmat
174 1956 : maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
175 2282 : maxn(in) = MIN(maxn(in), pp_basis%nbas(in))
176 : END DO
177 :
178 : ! read optimization section
179 326 : opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
180 326 : CALL read_atom_opt_section(optimization, opt_section)
181 :
182 326 : had_ae = .FALSE.
183 326 : had_pp = .FALSE.
184 :
185 : ! Check for the total number of electron configurations to be calculated
186 326 : 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 326 : method_section => section_vals_get_subs_vals(atom_section, "METHOD")
189 326 : CALL section_vals_get(method_section, n_repetition=n_meth)
190 :
191 : ! integrals
192 138550 : ALLOCATE (ae_int, pp_int)
193 :
194 1962 : ALLOCATE (atom_info(n_rep, n_meth))
195 :
196 658 : DO in = 1, n_rep
197 990 : DO im = 1, n_meth
198 :
199 332 : NULLIFY (atom_info(in, im)%atom)
200 332 : CALL create_atom_type(atom_info(in, im)%atom)
201 :
202 332 : atom_info(in, im)%atom%optimization = optimization
203 :
204 332 : atom_info(in, im)%atom%z = zval
205 :
206 332 : xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
207 332 : 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 332 : xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
212 332 : CALL section_vals_get(xc_func, explicit=explicit)
213 332 : IF (explicit) THEN
214 308 : CALL section_vals_val_get(xc_func, "_SECTION_PARAMETERS_", i_val=shortcut)
215 78 : SELECT CASE (shortcut)
216 : CASE (xc_funct_no_shortcut, xc_none)
217 78 : xcfstr = "NONE"
218 : CASE (xc_funct_pbe0)
219 26 : xcfstr = "PBE0"
220 : CASE (xc_funct_beefvdw)
221 0 : xcfstr = "BEEFVDW"
222 : CASE (xc_funct_b3lyp)
223 16 : xcfstr = "B3LYP"
224 : CASE (xc_funct_blyp)
225 58 : xcfstr = "BLYP"
226 : CASE (xc_funct_bp)
227 0 : xcfstr = "BP"
228 : CASE (xc_funct_pade)
229 62 : xcfstr = "PADE"
230 : CASE (xc_funct_pbe)
231 54 : xcfstr = "PBE"
232 : CASE (xc_funct_tpss)
233 10 : xcfstr = "TPSS"
234 : CASE (xc_funct_olyp)
235 2 : xcfstr = "OLYP"
236 : CASE (xc_funct_hcth120)
237 2 : xcfstr = "HCTH120"
238 : CASE DEFAULT
239 308 : xcfstr = "DFT"
240 : END SELECT
241 : ELSE
242 24 : 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 332 : zmp_section => section_vals_get_subs_vals(method_section, "ZMP")
251 332 : CALL section_vals_get(zmp_section, explicit=do_zmp)
252 332 : atom_info(in, im)%atom%do_zmp = do_zmp
253 332 : CALL section_vals_val_get(zmp_section, "FILE_DENSITY", c_val=filename)
254 332 : atom_info(in, im)%atom%ext_file = filename
255 : CALL section_vals_val_get(zmp_section, "GRID_TOL", &
256 332 : r_val=atom_info(in, im)%atom%zmpgrid_tol)
257 332 : CALL section_vals_val_get(zmp_section, "LAMBDA", r_val=lambda)
258 332 : atom_info(in, im)%atom%lambda = lambda
259 332 : CALL section_vals_val_get(zmp_section, "DM", l_val=dm)
260 332 : atom_info(in, im)%atom%dm = dm
261 :
262 332 : zmp_restart_section => section_vals_get_subs_vals(zmp_section, "RESTART")
263 332 : CALL section_vals_get(zmp_restart_section, explicit=doread)
264 332 : atom_info(in, im)%atom%doread = doread
265 332 : CALL section_vals_val_get(zmp_restart_section, "FILE_RESTART", c_val=filename)
266 332 : atom_info(in, im)%atom%zmp_restart_file = filename
267 :
268 : ! ZMP Reading external vxc section, if found initialize
269 332 : external_vxc_section => section_vals_get_subs_vals(method_section, "EXTERNAL_VXC")
270 332 : CALL section_vals_get(external_vxc_section, explicit=read_vxc)
271 332 : atom_info(in, im)%atom%read_vxc = read_vxc
272 332 : CALL section_vals_val_get(external_vxc_section, "FILE_VXC", c_val=filename)
273 332 : atom_info(in, im)%atom%ext_vxc_file = filename
274 : CALL section_vals_val_get(external_vxc_section, "GRID_TOL", &
275 332 : r_val=atom_info(in, im)%atom%zmpvxcgrid_tol)
276 :
277 120516 : ALLOCATE (state)
278 :
279 : ! get the electronic configuration
280 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
281 332 : c_vals=tmpstringlist)
282 :
283 : ! set occupations
284 332 : CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
285 332 : state%maxl_occ = get_maxl_occ(state%occ)
286 2324 : state%maxn_occ = get_maxn_occ(state%occ)
287 :
288 : ! set number of states to be calculated
289 332 : state%maxl_calc = MAX(maxl, state%maxl_occ)
290 332 : state%maxl_calc = MIN(lmat, state%maxl_calc)
291 2324 : state%maxn_calc = 0
292 1534 : DO k = 0, state%maxl_calc
293 1534 : state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k))
294 : END DO
295 :
296 : ! is there a pseudo potential
297 1106 : pp_calc = ANY(INDEX(tmpstringlist(1:), "CORE") /= 0)
298 332 : IF (pp_calc) THEN
299 : ! get and set the core occupations
300 74 : CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
301 74 : CALL atom_set_occupation(tmpstringlist, state%core, pocc)
302 5254 : zcore = zval - NINT(SUM(state%core))
303 74 : CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
304 : ELSE
305 18318 : state%core = 0._dp
306 258 : CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
307 : END IF
308 :
309 332 : CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
310 332 : CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
311 332 : CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)
312 :
313 332 : IF (atom_consistent_method(method, state%multiplicity)) THEN
314 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
315 332 : CALL atom_print_method(atom_info(in, im)%atom, iw)
316 332 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
317 :
318 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
319 332 : IF (pp_calc) THEN
320 74 : IF (iw > 0) CALL atom_print_potential(p_pot, iw)
321 : ELSE
322 258 : IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
323 : END IF
324 332 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
325 : END IF
326 :
327 : ! calculate integrals
328 332 : IF (pp_calc) THEN
329 : ! general integrals
330 : CALL atom_int_setup(pp_int, pp_basis, &
331 74 : potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e)
332 : ! potential
333 74 : CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot)
334 : !
335 74 : NULLIFY (pp_int%tzora, pp_int%hdkh)
336 : !
337 74 : CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot)
338 962 : state%maxn_calc(:) = MIN(state%maxn_calc(:), pp_basis%nbas(:))
339 518 : 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 258 : eri_coulomb=eri_c, eri_exchange=eri_e)
345 : ! potential
346 258 : CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
347 : ! relativistic correction terms
348 258 : CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
349 : !
350 258 : CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
351 3354 : state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
352 1806 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
353 : had_ae = .TRUE.
354 : END IF
355 :
356 332 : 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 332 : exchange_integral_type=do_erie)
360 332 : atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
361 332 : atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
362 :
363 332 : NULLIFY (orbitals)
364 2324 : mo = MAXVAL(state%maxn_calc)
365 2324 : mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
366 332 : CALL create_atom_orbs(orbitals, mb, mo)
367 332 : CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
368 :
369 2656 : IF (atom_consistent_method(method, state%multiplicity)) THEN
370 : !Calculate the electronic structure
371 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
372 332 : CALL calculate_atom(atom_info(in, im)%atom, iw)
373 332 : 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 332 : IF (atom_info(in, im)%atom%do_zmp) THEN
377 0 : CALL atom_write_zmp_restart(atom_info(in, im)%atom)
378 :
379 0 : ALLOCATE (ext_density(atom_info(in, im)%atom%basis%grid%nr))
380 0 : ext_density = 0._dp
381 0 : 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 0 : lprint=.TRUE.)
385 0 : DEALLOCATE (ext_density)
386 : END IF
387 : ! ZMP If we have the external v_xc calculate KS quantities
388 332 : IF (atom_info(in, im)%atom%read_vxc) THEN
389 0 : ALLOCATE (ext_vxc(atom_info(in, im)%atom%basis%grid%nr))
390 0 : ext_vxc = 0._dp
391 0 : 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 0 : lprint=.TRUE.)
395 0 : DEALLOCATE (ext_vxc)
396 : END IF
397 :
398 : ! Print out the orbitals if requested
399 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ORBITALS", extension=".log")
400 332 : CALL section_vals_val_get(atom_section, "PRINT%ORBITALS%XMGRACE", l_val=graph)
401 332 : IF (iw > 0) THEN
402 0 : CALL atom_print_orbitals(atom_info(in, im)%atom, iw, xmgrace=graph)
403 : END IF
404 332 : 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 332 : file_position="REWIND")
409 332 : IF (iw > 0) THEN
410 3 : CALL atom_write_upf(atom_info(in, im)%atom, iw, xcfstr)
411 : END IF
412 332 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%UPF_FILE")
413 :
414 : ! perform a fit of the total electronic density
415 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_DENSITY", extension=".log")
416 332 : IF (iw > 0) THEN
417 0 : CALL section_vals_val_get(atom_section, "PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto)
418 0 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
419 0 : CALL atom_fit_density(atom_info(in, im)%atom, num_gto, 0, iw, powell_section=powell_section)
420 : END IF
421 332 : 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 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_KGPOT", extension=".log")
425 332 : IF (iw > 0) THEN
426 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_GAUSSIAN", i_val=num_gau)
427 1 : CALL section_vals_val_get(atom_section, "PRINT%FIT_KGPOT%NUM_POLYNOM", i_val=num_pol)
428 1 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
429 1 : CALL atom_fit_kgpot(atom_info(in, im)%atom, num_gau, num_pol, iw, powell_section=powell_section)
430 : END IF
431 332 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_KGPOT")
432 :
433 : ! generate a response basis
434 332 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%RESPONSE_BASIS", extension=".log")
435 332 : IF (iw > 0) THEN
436 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta)
437 1 : CALL section_vals_val_get(atom_section, "PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder)
438 1 : CALL atom_response_basis(atom_info(in, im)%atom, delta, nder, iw)
439 : END IF
440 332 : 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 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
449 326 : IF (iw > 0) THEN
450 2 : CALL atom_grb_construction(atom_info, atom_section, iw)
451 : END IF
452 326 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
453 :
454 : ! Analyze basis sets
455 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
456 326 : IF (iw > 0) THEN
457 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
458 3 : CALL section_vals_val_get(atom_section, "PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
459 3 : crad = ptable(zval)%covalent_radius*bohr
460 3 : IF (had_ae) THEN
461 2 : IF (lcond) CALL atom_condnumber(ae_basis, crad, iw)
462 2 : IF (lcomp) CALL atom_completeness(ae_basis, zval, iw)
463 : END IF
464 3 : IF (had_pp) THEN
465 1 : IF (lcond) CALL atom_condnumber(pp_basis, crad, iw)
466 1 : IF (lcomp) CALL atom_completeness(pp_basis, zval, iw)
467 : END IF
468 : END IF
469 326 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
470 :
471 : ! Analyse ADMM approximation
472 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
473 326 : admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
474 326 : IF (iw > 0) THEN
475 1 : CALL atom_admm(atom_info, admm_section, iw)
476 : END IF
477 326 : 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 326 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO", extension=".log")
481 326 : sgp_section => section_vals_get_subs_vals(atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
482 326 : IF (iw > 0) THEN
483 5 : CALL atom_sgp_construction(atom_info, sgp_section, iw)
484 : END IF
485 326 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SEPARABLE_GAUSSIAN_PSEUDO")
486 :
487 : ! clean up
488 326 : IF (had_ae) THEN
489 256 : CALL atom_int_release(ae_int)
490 256 : CALL atom_ppint_release(ae_int)
491 256 : CALL atom_relint_release(ae_int)
492 : END IF
493 326 : IF (had_pp) THEN
494 72 : CALL atom_int_release(pp_int)
495 72 : CALL atom_ppint_release(pp_int)
496 72 : CALL atom_relint_release(pp_int)
497 : END IF
498 326 : CALL release_atom_basis(ae_basis)
499 326 : CALL release_atom_basis(pp_basis)
500 :
501 326 : CALL release_atom_potential(p_pot)
502 326 : CALL release_atom_potential(ae_pot)
503 :
504 658 : DO in = 1, n_rep
505 990 : DO im = 1, n_meth
506 664 : CALL release_atom_type(atom_info(in, im)%atom)
507 : END DO
508 : END DO
509 326 : DEALLOCATE (atom_info)
510 :
511 326 : DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
512 :
513 326 : CALL timestop(handle)
514 :
515 2282 : 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 1 : SUBROUTINE atom_response_basis(atom, delta, nder, iw)
528 :
529 : TYPE(atom_type), POINTER :: atom
530 : REAL(KIND=dp), INTENT(IN) :: delta
531 : INTEGER, INTENT(IN) :: nder, iw
532 :
533 : INTEGER :: i, ider, j, k, l, lhomo, m, n, nhomo, &
534 : s1, s2
535 : REAL(KIND=dp) :: dene, emax, expzet, fhomo, o, prefac, &
536 : zeta
537 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
538 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rbasis
539 1 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
540 1 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp
541 : TYPE(atom_state), POINTER :: state
542 :
543 1 : WRITE (iw, '(/," ",79("*"),/,T30,A,A/," ",79("*"))') "RESPONSE BASIS for ", ptable(atom%z)%symbol
544 :
545 1 : state => atom%state
546 1 : ovlp => atom%integrals%ovlp
547 :
548 : ! find HOMO
549 1 : lhomo = -1
550 1 : nhomo = -1
551 1 : emax = -HUGE(1._dp)
552 4 : DO l = 0, state%maxl_occ
553 6 : DO i = 1, state%maxn_occ(l)
554 5 : IF (atom%orbitals%ener(i, l) > emax) THEN
555 1 : lhomo = l
556 1 : nhomo = i
557 1 : emax = atom%orbitals%ener(i, l)
558 1 : fhomo = state%occupation(l, i)
559 : END IF
560 : END DO
561 : END DO
562 :
563 1 : s1 = SIZE(atom%orbitals%wfn, 1)
564 1 : s2 = SIZE(atom%orbitals%wfn, 2)
565 6 : ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
566 7 : s2 = MAXVAL(state%maxn_occ) + nder
567 5 : ALLOCATE (rbasis(s1, s2, 0:lmat))
568 91 : rbasis = 0._dp
569 :
570 4 : DO ider = -nder, nder
571 3 : dene = REAL(ider, KIND=dp)*delta
572 3 : CPASSERT(fhomo > ABS(dene))
573 3 : state%occupation(lhomo, nhomo) = fhomo + dene
574 3 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
575 147 : wfn(:, :, :, ider) = atom%orbitals%wfn
576 4 : state%occupation(lhomo, nhomo) = fhomo
577 : END DO
578 : ! reset
579 1 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
580 :
581 4 : DO l = 0, state%maxl_occ
582 : ! occupied states
583 6 : DO i = 1, MAX(state%maxn_occ(l), 1)
584 24 : rbasis(:, i, l) = wfn(:, i, l, 0)
585 : END DO
586 : ! differentiation
587 6 : DO ider = 1, nder
588 3 : i = MAX(state%maxn_occ(l), 1)
589 3 : SELECT CASE (ider)
590 : CASE (1)
591 21 : rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
592 : CASE (2)
593 0 : rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
594 : CASE (3)
595 : rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
596 0 : + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
597 : CASE DEFAULT
598 3 : CPABORT("")
599 : END SELECT
600 : END DO
601 :
602 : ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
603 3 : n = state%maxn_occ(l) + nder
604 3 : m = atom%basis%nbas(l)
605 8 : DO i = 1, n
606 7 : DO j = 1, i - 1
607 192 : o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
608 19 : rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
609 : END DO
610 480 : o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
611 38 : rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
612 : END DO
613 :
614 : ! check
615 12 : ALLOCATE (amat(n, n))
616 578 : amat(1:n, 1:n) = MATMUL(TRANSPOSE(rbasis(1:m, 1:n, l)), MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
617 8 : DO i = 1, n
618 8 : amat(i, i) = amat(i, i) - 1._dp
619 : END DO
620 17 : IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
621 0 : WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
622 : END IF
623 3 : DEALLOCATE (amat)
624 :
625 : ! Quickstep normalization
626 3 : WRITE (iw, '(/,A,T30,I3)') " Angular momentum :", l
627 :
628 3 : WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
629 6 : n - nder, " valence ", nder, " response"
630 3 : expzet = 0.25_dp*REAL(2*l + 3, dp)
631 3 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
632 22 : DO i = 1, m
633 18 : zeta = (2._dp*atom%basis%am(i, l))**expzet
634 51 : WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)
635 : END DO
636 :
637 : END DO
638 :
639 1 : DEALLOCATE (wfn, rbasis)
640 :
641 1 : WRITE (iw, '(" ",79("*"))')
642 :
643 2 : END SUBROUTINE atom_response_basis
644 :
645 : ! **************************************************************************************************
646 : !> \brief Write pseudo-potential in Quantum Espresso UPF format.
647 : !> \param atom information about the atomic kind
648 : !> \param iw output file unit
649 : !> \param xcfstr ...
650 : !> \par History
651 : !> * 09.2012 created [Juerg Hutter]
652 : ! **************************************************************************************************
653 3 : SUBROUTINE atom_write_upf(atom, iw, xcfstr)
654 :
655 : TYPE(atom_type), POINTER :: atom
656 : INTEGER, INTENT(IN) :: iw
657 : CHARACTER(LEN=*), INTENT(IN) :: xcfstr
658 :
659 : CHARACTER(LEN=default_string_length) :: string
660 : INTEGER :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
661 : nbeta, nr, nwfc, nwfn
662 : LOGICAL :: soc, up
663 : REAL(KIND=dp) :: occa, occb, pf, rhoatom, rl, rlp, rmax, &
664 : vnlm, vnlp
665 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
666 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
667 : TYPE(atom_gthpot_type), POINTER :: pot
668 :
669 3 : IF (.NOT. atom%pp_calc) RETURN
670 3 : IF (atom%potential%ppot_type /= gth_pseudo) RETURN
671 3 : pot => atom%potential%gth_pot
672 3 : CPASSERT(.NOT. pot%lsdpot)
673 3 : soc = pot%soc
674 :
675 3 : WRITE (iw, '(A)') '<UPF version="2.0.1">'
676 :
677 3 : WRITE (iw, '(T4,A)') '<PP_INFO>'
678 3 : WRITE (iw, '(T8,A)') 'Converted from CP2K GTH format'
679 3 : WRITE (iw, '(T8,A)') '<PP_INPUTFILE>'
680 3 : CALL atom_write_pseudo_param(pot, iw)
681 3 : WRITE (iw, '(T8,A)') '</PP_INPUTFILE>'
682 3 : WRITE (iw, '(T4,A)') '</PP_INFO>'
683 3 : WRITE (iw, '(T4,A)') '<PP_HEADER'
684 3 : WRITE (iw, '(T8,A)') 'generated="Generated in analytical, separable form"'
685 3 : WRITE (iw, '(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"'
686 3 : WRITE (iw, '(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"'
687 3 : WRITE (iw, '(T8,A)') 'comment="This file generated by CP2K ATOM code"'
688 3 : CALL compose(string, "element", cval=pot%symbol)
689 3 : WRITE (iw, '(T8,A)') TRIM(string)
690 3 : WRITE (iw, '(T8,A)') 'pseudo_type="NC"'
691 3 : WRITE (iw, '(T8,A)') 'relativistic="no"'
692 3 : WRITE (iw, '(T8,A)') 'is_ultrasoft="F"'
693 3 : WRITE (iw, '(T8,A)') 'is_paw="F"'
694 3 : WRITE (iw, '(T8,A)') 'is_coulomb="F"'
695 3 : IF (soc) THEN
696 1 : WRITE (iw, '(T8,A)') 'has_so="T"'
697 : ELSE
698 2 : WRITE (iw, '(T8,A)') 'has_so="F"'
699 : END IF
700 3 : WRITE (iw, '(T8,A)') 'has_wfc="F"'
701 3 : WRITE (iw, '(T8,A)') 'has_gipaw="F"'
702 3 : WRITE (iw, '(T8,A)') 'paw_as_gipaw="F"'
703 3 : IF (pot%nlcc) THEN
704 1 : WRITE (iw, '(T8,A)') 'core_correction="T"'
705 : ELSE
706 2 : WRITE (iw, '(T8,A)') 'core_correction="F"'
707 : END IF
708 3 : WRITE (iw, '(T8,A)') 'functional="'//TRIM(ADJUSTL(xcfstr))//'"'
709 3 : CALL compose(string, "z_valence", rval=pot%zion)
710 3 : WRITE (iw, '(T8,A)') TRIM(string)
711 3 : CALL compose(string, "total_psenergy", rval=2._dp*atom%energy%etot)
712 3 : WRITE (iw, '(T8,A)') TRIM(string)
713 3 : WRITE (iw, '(T8,A)') 'wfc_cutoff="0.0E+00"'
714 3 : WRITE (iw, '(T8,A)') 'rho_cutoff="0.0E+00"'
715 3 : lmax = -1
716 21 : DO l = 0, lmat
717 21 : IF (pot%nl(l) > 0) lmax = l
718 : END DO
719 3 : CALL compose(string, "l_max", ival=lmax)
720 3 : WRITE (iw, '(T8,A)') TRIM(string)
721 3 : WRITE (iw, '(T8,A)') 'l_max_rho="0"'
722 3 : WRITE (iw, '(T8,A)') 'l_local="-3"'
723 3 : nr = atom%basis%grid%nr
724 3 : CALL compose(string, "mesh_size", ival=nr)
725 3 : WRITE (iw, '(T8,A)') TRIM(string)
726 21 : nwfc = SUM(atom%state%maxn_occ)
727 3 : CALL compose(string, "number_of_wfc", ival=nwfc)
728 3 : WRITE (iw, '(T8,A)') TRIM(string)
729 3 : IF (soc) THEN
730 6 : nbeta = pot%nl(0) + 2*SUM(pot%nl(1:))
731 : ELSE
732 14 : nbeta = SUM(pot%nl)
733 : END IF
734 3 : CALL compose(string, "number_of_proj", ival=nbeta)
735 3 : WRITE (iw, '(T8,A)') TRIM(string)//'/>'
736 :
737 : ! Mesh
738 3 : up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr)
739 3 : WRITE (iw, '(T4,A)') '<PP_MESH'
740 3 : WRITE (iw, '(T8,A)') 'dx="1.E+00"'
741 3 : CALL compose(string, "mesh", ival=nr)
742 3 : WRITE (iw, '(T8,A)') TRIM(string)
743 3 : WRITE (iw, '(T8,A)') 'xmin="1.E+00"'
744 1206 : rmax = MAXVAL(atom%basis%grid%rad)
745 3 : CALL compose(string, "rmax", rval=rmax)
746 3 : WRITE (iw, '(T8,A)') TRIM(string)
747 3 : WRITE (iw, '(T8,A)') 'zmesh="1.E+00">'
748 3 : WRITE (iw, '(T8,A)') '<PP_R type="real"'
749 3 : CALL compose(string, "size", ival=nr)
750 3 : WRITE (iw, '(T12,A)') TRIM(string)
751 3 : WRITE (iw, '(T12,A)') 'columns="4">'
752 3 : IF (up) THEN
753 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=1, nr)
754 : ELSE
755 1203 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%rad(i), i=nr, 1, -1)
756 : END IF
757 3 : WRITE (iw, '(T8,A)') '</PP_R>'
758 3 : WRITE (iw, '(T8,A)') '<PP_RAB type="real"'
759 3 : CALL compose(string, "size", ival=nr)
760 3 : WRITE (iw, '(T12,A)') TRIM(string)
761 3 : WRITE (iw, '(T12,A)') 'columns="4">'
762 3 : IF (up) THEN
763 0 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=1, nr)
764 : ELSE
765 1203 : WRITE (iw, '(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i), i=nr, 1, -1)
766 : END IF
767 3 : WRITE (iw, '(T8,A)') '</PP_RAB>'
768 3 : WRITE (iw, '(T4,A)') '</PP_MESH>'
769 :
770 : ! NLCC
771 3 : IF (pot%nlcc) THEN
772 1 : WRITE (iw, '(T4,A)') '<PP_NLCC type="real"'
773 1 : CALL compose(string, "size", ival=nr)
774 1 : WRITE (iw, '(T8,A)') TRIM(string)
775 1 : WRITE (iw, '(T8,A)') 'columns="4">'
776 3 : ALLOCATE (corden(nr))
777 401 : corden(:) = 0.0_dp
778 1 : CALL atom_core_density(corden, atom%potential, "RHO", atom%basis%grid%rad)
779 1 : IF (up) THEN
780 0 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=1, nr)
781 : ELSE
782 1 : WRITE (iw, '(T8,4ES25.12E3)') (corden(i), i=nr, 1, -1)
783 : END IF
784 1 : DEALLOCATE (corden)
785 1 : WRITE (iw, '(T4,A)') '</PP_NLCC>'
786 : END IF
787 :
788 : ! local PP
789 3 : WRITE (iw, '(T4,A)') '<PP_LOCAL type="real"'
790 3 : CALL compose(string, "size", ival=nr)
791 3 : WRITE (iw, '(T8,A)') TRIM(string)
792 3 : WRITE (iw, '(T8,A)') 'columns="4">'
793 9 : ALLOCATE (locpot(nr))
794 1203 : locpot(:) = 0.0_dp
795 3 : CALL atom_local_potential(locpot, pot, atom%basis%grid%rad)
796 3 : IF (up) THEN
797 0 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=1, nr)
798 : ELSE
799 1203 : WRITE (iw, '(T8,4ES25.12E3)') (2.0_dp*locpot(i), i=nr, 1, -1)
800 : END IF
801 3 : DEALLOCATE (locpot)
802 3 : WRITE (iw, '(T4,A)') '</PP_LOCAL>'
803 :
804 : ! nonlocal PP
805 3 : WRITE (iw, '(T4,A)') '<PP_NONLOCAL>'
806 12 : ALLOCATE (rp(nr), ef(nr), beta(nr))
807 3 : ibeta = 0
808 21 : DO l = 0, lmat
809 18 : IF (pot%nl(l) == 0) CYCLE
810 6 : rl = pot%rcnl(l)
811 2406 : rp(:) = atom%basis%grid%rad
812 2406 : ef(:) = EXP(-0.5_dp*rp*rp/(rl*rl))
813 20 : DO i = 1, pot%nl(l)
814 11 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
815 11 : j = l + 2*i - 1
816 11 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
817 4411 : beta(:) = pf*rp**(l + 2*i - 2)*ef
818 4411 : beta(:) = beta*rp
819 : !
820 11 : ibeta = ibeta + 1
821 11 : CALL compose(string, "<PP_BETA", counter=ibeta)
822 11 : WRITE (iw, '(T8,A)') TRIM(string)
823 11 : CALL compose(string, "angular_momentum", ival=l)
824 11 : WRITE (iw, '(T12,A)') TRIM(string)
825 11 : WRITE (iw, '(T12,A)') 'type="real"'
826 11 : CALL compose(string, "size", ival=nr)
827 11 : WRITE (iw, '(T12,A)') TRIM(string)
828 11 : WRITE (iw, '(T12,A)') 'columns="4">'
829 11 : IF (up) THEN
830 0 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
831 : ELSE
832 4411 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
833 : END IF
834 11 : CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
835 11 : WRITE (iw, '(T8,A)') TRIM(string)
836 29 : IF (soc .AND. l /= 0) THEN
837 5 : ibeta = ibeta + 1
838 5 : CALL compose(string, "<PP_BETA", counter=ibeta)
839 5 : WRITE (iw, '(T8,A)') TRIM(string)
840 5 : CALL compose(string, "angular_momentum", ival=l)
841 5 : WRITE (iw, '(T12,A)') TRIM(string)
842 5 : WRITE (iw, '(T12,A)') 'type="real"'
843 5 : CALL compose(string, "size", ival=nr)
844 5 : WRITE (iw, '(T12,A)') TRIM(string)
845 5 : WRITE (iw, '(T12,A)') 'columns="4">'
846 5 : IF (up) THEN
847 0 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=1, nr)
848 : ELSE
849 2005 : WRITE (iw, '(T12,4ES25.12E3)') (2._dp*beta(j), j=nr, 1, -1)
850 : END IF
851 5 : CALL compose(string, "</PP_BETA", counter=ibeta, isfinal=.TRUE.)
852 5 : WRITE (iw, '(T8,A)') TRIM(string)
853 : END IF
854 : END DO
855 : END DO
856 3 : DEALLOCATE (rp, ef, beta)
857 : ! nonlocal PP matrix elements
858 12 : ALLOCATE (dij(nbeta, nbeta))
859 193 : dij = 0._dp
860 3 : IF (soc) THEN
861 : ! from reverse engineering we guess: hnl, knl, hnl, knl, ...
862 1 : n0 = pot%nl(0)
863 7 : DO l = 0, lmat
864 6 : IF (pot%nl(l) == 0) CYCLE
865 4 : IF (l == 0) THEN
866 : ! no SO term for l=0
867 13 : dij(1:n0, 1:n0) = pot%hnl(1:n0, 1:n0, 0)
868 : ELSE
869 3 : ibeta = n0 + 2*SUM(pot%nl(1:l - 1))
870 7 : DO i = 1, pot%nl(l)
871 5 : ii = 2*(i - 1) + 1
872 23 : DO k = 1, pot%nl(l)
873 13 : kk = 2*(k - 1) + 1
874 13 : vnlp = pot%hnl(i, k, l) + 0.5_dp*l*pot%knl(i, k, l)
875 13 : vnlm = pot%hnl(i, k, l) - 0.5_dp*(l + 1)*pot%knl(i, k, l)
876 13 : dij(ibeta + ii, ibeta + kk) = vnlm
877 18 : dij(ibeta + ii + 1, ibeta + kk + 1) = vnlp
878 : END DO
879 : END DO
880 : END IF
881 : END DO
882 : ELSE
883 14 : DO l = 0, lmat
884 12 : IF (pot%nl(l) == 0) CYCLE
885 4 : ibeta = SUM(pot%nl(0:l - 1)) + 1
886 3 : i = ibeta + pot%nl(l) - 1
887 11 : dij(ibeta:i, ibeta:i) = pot%hnl(1:pot%nl(l), 1:pot%nl(l), l)
888 : END DO
889 : END IF
890 3 : WRITE (iw, '(T8,A)') '<PP_DIJ type="real"'
891 3 : WRITE (iw, '(T12,A)') 'columns="4">'
892 193 : WRITE (iw, '(T12,4ES25.12E3)') ((0.5_dp*dij(i, j), j=1, nbeta), i=1, nbeta)
893 3 : WRITE (iw, '(T8,A)') '</PP_DIJ>'
894 3 : DEALLOCATE (dij)
895 3 : WRITE (iw, '(T4,A)') '</PP_NONLOCAL>'
896 :
897 : ! atomic wavefunctions
898 6 : ALLOCATE (beta(nr))
899 3 : WRITE (iw, '(T4,A)') '<PP_PSWFC>'
900 3 : nwfn = 0
901 21 : DO l = 0, lmat
902 201 : DO i = 1, 10
903 180 : IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
904 8 : occa = atom%state%occupation(l, i)
905 8 : occb = 0.0_dp
906 8 : IF (soc .AND. l /= 0) THEN
907 2 : occb = FLOOR(occa/2)
908 2 : occa = occa - occb
909 : END IF
910 8 : nwfn = nwfn + 1
911 8 : CALL compose(string, "<PP_CHI", counter=nwfn)
912 8 : WRITE (iw, '(T8,A)') TRIM(string)
913 8 : CALL compose(string, "l", ival=l)
914 8 : WRITE (iw, '(T12,A)') TRIM(string)
915 8 : CALL compose(string, "occupation", rval=occa)
916 8 : WRITE (iw, '(T12,A)') TRIM(string)
917 8 : CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
918 8 : WRITE (iw, '(T12,A)') TRIM(string)
919 8 : WRITE (iw, '(T12,A)') 'columns="4">'
920 3208 : beta = 0._dp
921 143 : DO k = 1, atom%basis%nbas(l)
922 54143 : beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
923 : END DO
924 3208 : beta(:) = beta*atom%basis%grid%rad
925 8 : IF (up) THEN
926 0 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
927 : ELSE
928 8 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
929 : END IF
930 8 : CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
931 8 : WRITE (iw, '(T8,A)') TRIM(string)
932 26 : IF (soc .AND. l /= 0) THEN
933 2 : nwfn = nwfn + 1
934 2 : CALL compose(string, "<PP_CHI", counter=nwfn)
935 2 : WRITE (iw, '(T8,A)') TRIM(string)
936 2 : CALL compose(string, "l", ival=l)
937 2 : WRITE (iw, '(T12,A)') TRIM(string)
938 2 : CALL compose(string, "occupation", rval=occb)
939 2 : WRITE (iw, '(T12,A)') TRIM(string)
940 2 : CALL compose(string, "pseudo_energy", rval=2._dp*atom%orbitals%ener(i, l))
941 2 : WRITE (iw, '(T12,A)') TRIM(string)
942 2 : WRITE (iw, '(T12,A)') 'columns="4">'
943 802 : beta = 0._dp
944 34 : DO k = 1, atom%basis%nbas(l)
945 12834 : beta(:) = beta(:) + atom%orbitals%wfn(k, i, l)*atom%basis%bf(:, k, l)
946 : END DO
947 802 : beta(:) = beta*atom%basis%grid%rad
948 2 : IF (up) THEN
949 0 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=1, nr)
950 : ELSE
951 2 : WRITE (iw, '(T12,4ES25.12E3)') (beta(j), j=nr, 1, -1)
952 : END IF
953 2 : CALL compose(string, "</PP_CHI", counter=nwfn, isfinal=.TRUE.)
954 2 : WRITE (iw, '(T8,A)') TRIM(string)
955 : END IF
956 : END DO
957 : END DO
958 3 : WRITE (iw, '(T4,A)') '</PP_PSWFC>'
959 3 : DEALLOCATE (beta)
960 :
961 : ! atomic charge
962 6 : ALLOCATE (dens(nr))
963 3 : WRITE (iw, '(T4,A)') '<PP_RHOATOM type="real"'
964 3 : CALL compose(string, "size", ival=nr)
965 3 : WRITE (iw, '(T8,A)') TRIM(string)
966 3 : WRITE (iw, '(T8,A)') 'columns="4">'
967 : CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
968 3 : "RHO", atom%basis%grid%rad)
969 3 : IF (up) THEN
970 0 : WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
971 : ELSE
972 1203 : WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
973 : END IF
974 3 : WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
975 :
976 : ! Check atomic electron count
977 1203 : rhoatom = SUM(fourpi*atom%basis%grid%wr(:)*dens(:))
978 3 : IF (ABS(pot%zion - rhoatom) > 0.01_dp) THEN
979 1 : CALL cp_warn(__LOCATION__, "Valence charge and electron count differ by more than 0.01")
980 : END IF
981 :
982 3 : DEALLOCATE (dens)
983 : ! PP SOC information
984 3 : IF (soc) THEN
985 1 : WRITE (iw, '(T4,A)') '<PP_SPIN_ORB>'
986 : ! <PP_RELWFC.1 index="1" els="6S" nn="1" lchi="0" jchi="5.000000000000e-1" oc="1.000000000000e0"/>
987 1 : nwfn = 0
988 7 : DO l = 0, lmat
989 67 : DO i = 1, 10
990 60 : IF (ABS(atom%state%occupation(l, i)) == 0._dp) CYCLE
991 4 : occa = atom%state%occupation(l, i)
992 4 : occb = 0.0_dp
993 4 : IF (l /= 0) THEN
994 2 : occb = FLOOR(occa/2)
995 2 : occa = occa - occb
996 : END IF
997 4 : nwfn = nwfn + 1
998 4 : CALL compose(string, "<PP_RELWFC", counter=nwfn)
999 4 : WRITE (iw, '(T8,A)') TRIM(string)
1000 4 : CALL compose(string, "index", ival=nwfn)
1001 4 : WRITE (iw, '(T12,A)') TRIM(string)
1002 : !!CALL compose(string, "els", cval=xxxx)
1003 : !!WRITE (iw, '(T12,A)') TRIM(string)
1004 4 : CALL compose(string, "nn", ival=nwfn)
1005 4 : WRITE (iw, '(T12,A)') TRIM(string)
1006 4 : CALL compose(string, "lchi", ival=l)
1007 4 : WRITE (iw, '(T12,A)') TRIM(string)
1008 4 : IF (l == 0) THEN
1009 2 : rlp = 0.5_dp
1010 : ELSE
1011 2 : rlp = l - 0.5_dp
1012 : END IF
1013 4 : CALL compose(string, "jchi", rval=rlp)
1014 4 : WRITE (iw, '(T12,A)') TRIM(string)
1015 4 : CALL compose(string, "oc", rval=occa)
1016 4 : WRITE (iw, '(T12,A)') TRIM(string)
1017 4 : WRITE (iw, '(T12,A)') '/>'
1018 10 : IF (l /= 0) THEN
1019 2 : nwfn = nwfn + 1
1020 2 : CALL compose(string, "<PP_RELWFC", counter=nwfn)
1021 2 : WRITE (iw, '(T8,A)') TRIM(string)
1022 2 : CALL compose(string, "index", ival=nwfn)
1023 2 : WRITE (iw, '(T12,A)') TRIM(string)
1024 : !!CALL compose(string, "els", cval=xxxx)
1025 : !!WRITE (iw, '(T12,A)') TRIM(string)
1026 2 : CALL compose(string, "nn", ival=nwfn)
1027 2 : WRITE (iw, '(T12,A)') TRIM(string)
1028 2 : CALL compose(string, "lchi", ival=l)
1029 2 : WRITE (iw, '(T12,A)') TRIM(string)
1030 2 : rlp = l + 0.5_dp
1031 2 : CALL compose(string, "jchi", rval=rlp)
1032 2 : WRITE (iw, '(T12,A)') TRIM(string)
1033 2 : CALL compose(string, "oc", rval=occa)
1034 2 : WRITE (iw, '(T12,A)') TRIM(string)
1035 2 : WRITE (iw, '(T12,A)') '/>'
1036 : END IF
1037 : END DO
1038 : END DO
1039 1 : ibeta = 0
1040 7 : DO l = 0, lmat
1041 6 : IF (pot%nl(l) == 0) CYCLE
1042 12 : DO i = 1, pot%nl(l)
1043 8 : ibeta = ibeta + 1
1044 8 : CALL compose(string, "<PP_RELBETA", counter=ibeta)
1045 8 : WRITE (iw, '(T8,A)') TRIM(string)
1046 8 : CALL compose(string, "index", ival=ibeta)
1047 8 : WRITE (iw, '(T12,A)') TRIM(string)
1048 8 : CALL compose(string, "lll", ival=l)
1049 8 : WRITE (iw, '(T12,A)') TRIM(string)
1050 8 : IF (l == 0) THEN
1051 3 : rlp = 0.5_dp
1052 : ELSE
1053 5 : rlp = l - 0.5_dp
1054 : END IF
1055 8 : CALL compose(string, "jjj", rval=rlp)
1056 8 : WRITE (iw, '(T12,A)') TRIM(string)
1057 8 : WRITE (iw, '(T12,A)') '/>'
1058 14 : IF (l > 0) THEN
1059 5 : ibeta = ibeta + 1
1060 5 : CALL compose(string, "<PP_RELBETA", counter=ibeta)
1061 5 : WRITE (iw, '(T8,A)') TRIM(string)
1062 5 : CALL compose(string, "index", ival=ibeta)
1063 5 : WRITE (iw, '(T12,A)') TRIM(string)
1064 5 : CALL compose(string, "lll", ival=l)
1065 5 : WRITE (iw, '(T12,A)') TRIM(string)
1066 5 : rlp = l + 0.5_dp
1067 5 : CALL compose(string, "jjj", rval=rlp)
1068 5 : WRITE (iw, '(T12,A)') TRIM(string)
1069 5 : WRITE (iw, '(T12,A)') '/>'
1070 : END IF
1071 : END DO
1072 : END DO
1073 1 : WRITE (iw, '(T4,A)') '</PP_SPIN_ORB>'
1074 : END IF
1075 :
1076 3 : WRITE (iw, '(A)') '</UPF>'
1077 :
1078 6 : END SUBROUTINE atom_write_upf
1079 :
1080 : ! **************************************************************************************************
1081 : !> \brief Produce an XML attribute string 'tag="counter/rval/ival/cval"'
1082 : !> \param string composed string
1083 : !> \param tag attribute tag
1084 : !> \param counter counter
1085 : !> \param rval real variable
1086 : !> \param ival integer variable
1087 : !> \param cval string variable
1088 : !> \param isfinal close the current XML element if this is the last attribute
1089 : !> \par History
1090 : !> * 09.2012 created [Juerg Hutter]
1091 : ! **************************************************************************************************
1092 242 : SUBROUTINE compose(string, tag, counter, rval, ival, cval, isfinal)
1093 : CHARACTER(len=*), INTENT(OUT) :: string
1094 : CHARACTER(len=*), INTENT(IN) :: tag
1095 : INTEGER, INTENT(IN), OPTIONAL :: counter
1096 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rval
1097 : INTEGER, INTENT(IN), OPTIONAL :: ival
1098 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: cval
1099 : LOGICAL, INTENT(IN), OPTIONAL :: isfinal
1100 :
1101 : CHARACTER(LEN=default_string_length) :: str
1102 : LOGICAL :: fin
1103 :
1104 242 : IF (PRESENT(counter)) THEN
1105 71 : WRITE (str, "(I12)") counter
1106 171 : ELSEIF (PRESENT(rval)) THEN
1107 54 : WRITE (str, "(G18.8)") rval
1108 117 : ELSEIF (PRESENT(ival)) THEN
1109 114 : WRITE (str, "(I12)") ival
1110 3 : ELSEIF (PRESENT(cval)) THEN
1111 3 : WRITE (str, "(A)") TRIM(ADJUSTL(cval))
1112 : ELSE
1113 0 : WRITE (str, "(A)") ""
1114 : END IF
1115 242 : fin = .FALSE.
1116 242 : IF (PRESENT(isfinal)) fin = isfinal
1117 242 : IF (PRESENT(counter)) THEN
1118 71 : IF (fin) THEN
1119 26 : WRITE (string, "(A,A1,A,A1)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str)), '>'
1120 : ELSE
1121 45 : WRITE (string, "(A,A1,A)") TRIM(ADJUSTL(tag)), '.', TRIM(ADJUSTL(str))
1122 : END IF
1123 : ELSE
1124 171 : IF (fin) THEN
1125 0 : WRITE (string, "(A,A2,A,A2)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '>"'
1126 : ELSE
1127 171 : WRITE (string, "(A,A2,A,A1)") TRIM(ADJUSTL(tag)), '="', TRIM(ADJUSTL(str)), '"'
1128 : END IF
1129 : END IF
1130 242 : END SUBROUTINE compose
1131 :
1132 10 : END MODULE atom_energy
|
