Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !***************************************************************************************************
9 : !> \brief Interface to the SIRIUS Library
10 : !> \par History
11 : !> 07.2018 initial create
12 : !> \author JHU
13 : !***************************************************************************************************
14 : #if defined(__SIRIUS)
15 : MODULE sirius_interface
16 : USE ISO_C_BINDING, ONLY: C_DOUBLE,&
17 : C_INT,&
18 : C_LOC
19 : USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals,&
20 : gth_potential_conversion
21 : USE atom_types, ONLY: atom_gthpot_type
22 : USE atom_upf, ONLY: atom_upfpot_type
23 : USE atom_utils, ONLY: atom_local_potential
24 : USE atomic_kind_types, ONLY: atomic_kind_type,&
25 : get_atomic_kind
26 : USE cell_types, ONLY: cell_type,&
27 : real_to_scaled
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_p_file,&
32 : cp_print_key_finished_output,&
33 : cp_print_key_should_output,&
34 : cp_print_key_unit_nr
35 : USE external_potential_types, ONLY: gth_potential_type
36 : USE input_constants, ONLY: do_gapw_log
37 : USE input_cp2k_pwdft, ONLY: SIRIUS_FUNC_VDWDF,&
38 : SIRIUS_FUNC_VDWDF2,&
39 : SIRIUS_FUNC_VDWDFCX
40 : USE input_section_types, ONLY: section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_get_subs_vals2,&
43 : section_vals_type,&
44 : section_vals_val_get
45 : USE kinds, ONLY: default_string_length,&
46 : dp
47 : USE machine, ONLY: m_flush
48 : USE mathconstants, ONLY: fourpi,&
49 : gamma1
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE particle_types, ONLY: particle_type
52 : USE physcon, ONLY: massunit
53 : USE pwdft_environment_types, ONLY: pwdft_energy_type,&
54 : pwdft_env_get,&
55 : pwdft_env_set,&
56 : pwdft_environment_type
57 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
58 : create_grid_atom,&
59 : deallocate_grid_atom,&
60 : grid_atom_type
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_subsys_types, ONLY: qs_subsys_get,&
64 : qs_subsys_type
65 : USE sirius, ONLY: &
66 : SIRIUS_INTEGER_ARRAY_TYPE, SIRIUS_INTEGER_TYPE, SIRIUS_LOGICAL_ARRAY_TYPE, &
67 : SIRIUS_LOGICAL_TYPE, SIRIUS_NUMBER_ARRAY_TYPE, SIRIUS_NUMBER_TYPE, &
68 : SIRIUS_STRING_ARRAY_TYPE, SIRIUS_STRING_TYPE, sirius_add_atom, sirius_add_atom_type, &
69 : sirius_add_atom_type_radial_function, sirius_add_xc_functional, sirius_context_handler, &
70 : sirius_create_context, sirius_create_ground_state, sirius_create_kset_from_grid, &
71 : sirius_finalize, sirius_find_ground_state, sirius_get_band_energies, &
72 : sirius_get_band_occupancies, sirius_get_energy, sirius_get_forces, &
73 : sirius_get_kpoint_properties, sirius_get_num_kpoints, sirius_get_parameters, &
74 : sirius_get_stress_tensor, sirius_ground_state_handler, sirius_import_parameters, &
75 : sirius_initialize, sirius_initialize_context, sirius_is_initialized, &
76 : sirius_kpoint_set_handler, sirius_option_get_info, sirius_option_get_section_length, &
77 : sirius_option_set, sirius_set_atom_position, sirius_set_atom_type_dion, &
78 : sirius_set_atom_type_hubbard, sirius_set_atom_type_radial_grid, &
79 : sirius_set_lattice_vectors, sirius_set_mpi_grid_dims, sirius_update_ground_state
80 : #include "./base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 :
84 : PRIVATE
85 :
86 : ! Global parameters
87 :
88 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sirius_interface'
89 :
90 : ! Public subroutines
91 :
92 : PUBLIC :: cp_sirius_create_env, &
93 : cp_sirius_energy_force, &
94 : cp_sirius_finalize, &
95 : cp_sirius_init, &
96 : cp_sirius_is_initialized, &
97 : cp_sirius_update_context
98 :
99 : CONTAINS
100 :
101 : !***************************************************************************************************
102 : !> \brief Initialize the Sirius library
103 : !> \par Creation (07.2018, JHU)
104 : !> \author JHU
105 : ! **************************************************************************************************
106 20 : SUBROUTINE cp_sirius_init()
107 20 : CALL sirius_initialize(.FALSE.)
108 20 : END SUBROUTINE cp_sirius_init
109 :
110 : !***************************************************************************************************
111 : !> \brief Check the initialisation status of the Sirius library
112 : !> \return Return the initialisation status of the Sirius library as boolean
113 : !> \par Creation (03.12.2025, MK)
114 : !> \author Matthias Krack
115 : ! **************************************************************************************************
116 9842 : LOGICAL FUNCTION cp_sirius_is_initialized()
117 9842 : CALL sirius_is_initialized(cp_sirius_is_initialized)
118 9842 : END FUNCTION cp_sirius_is_initialized
119 :
120 : !***************************************************************************************************
121 : !> \brief Finalize the Sirius library
122 : !> \par Creation (07.2018, JHU)
123 : !> \author JHU
124 : ! **************************************************************************************************
125 20 : SUBROUTINE cp_sirius_finalize()
126 20 : CALL sirius_finalize(.FALSE., .FALSE., .FALSE.)
127 20 : END SUBROUTINE cp_sirius_finalize
128 :
129 : !***************************************************************************************************
130 : !> \brief ...
131 : !> \param pwdft_env ...
132 : !> \param
133 : !> \par History
134 : !> 07.2018 Create the Sirius environment
135 : !> \author JHU
136 : ! **************************************************************************************************
137 20 : SUBROUTINE cp_sirius_create_env(pwdft_env)
138 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
139 : #if defined(__SIRIUS)
140 :
141 : CHARACTER(len=2) :: element_symbol
142 : CHARACTER(len=default_string_length) :: label
143 : INTEGER :: i, ii, jj, iatom, ibeta, ifun, ikind, iwf, j, l, &
144 : n, ns, natom, nbeta, nbs, nkind, nmesh, &
145 : num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu, output_unit
146 20 : INTEGER, DIMENSION(:), POINTER :: mpi_grid_dims
147 : INTEGER(KIND=C_INT), DIMENSION(3) :: k_grid, k_shift
148 20 : INTEGER, DIMENSION(:), POINTER :: kk
149 : LOGICAL :: up, use_ref_cell
150 : LOGICAL(4) :: use_so, use_symmetry, dft_plus_u_atom
151 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: fun
152 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: dion
153 : REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v1, v2
154 : REAL(KIND=dp) :: al, angle1, angle2, cval, focc, &
155 : magnetization, mass, pf, rl, zeff, alpha_u, beta_u, &
156 : J0_u, J_u, U_u, occ_u, u_minus_J, vnlp, vnlm
157 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, ef, fe, locpot, rc, rp
158 : REAL(KIND=dp), DIMENSION(3) :: vr, vs, j_t
159 20 : REAL(KIND=dp), DIMENSION(:), POINTER :: density
160 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: wavefunction, wfninfo
161 : TYPE(atom_gthpot_type), POINTER :: gth_atompot
162 : TYPE(atom_upfpot_type), POINTER :: upf_pot
163 20 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
164 : TYPE(atomic_kind_type), POINTER :: atomic_kind
165 : TYPE(cell_type), POINTER :: my_cell
166 : TYPE(mp_para_env_type), POINTER :: para_env
167 : TYPE(grid_atom_type), POINTER :: atom_grid
168 : TYPE(gth_potential_type), POINTER :: gth_potential
169 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
170 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
171 : TYPE(qs_subsys_type), POINTER :: qs_subsys
172 : TYPE(section_vals_type), POINTER :: pwdft_section, pwdft_sub_section, &
173 : xc_fun, xc_section
174 : TYPE(sirius_context_handler) :: sctx
175 : TYPE(sirius_ground_state_handler) :: gs_handler
176 : TYPE(sirius_kpoint_set_handler) :: ks_handler
177 :
178 0 : CPASSERT(ASSOCIATED(pwdft_env))
179 :
180 20 : output_unit = cp_logger_get_default_io_unit()
181 : ! create context of simulation
182 20 : CALL pwdft_env_get(pwdft_env, para_env=para_env)
183 20 : sirius_mpi_comm = para_env%get_handle()
184 20 : CALL sirius_create_context(sirius_mpi_comm, sctx)
185 :
186 : ! the "fun" starts.
187 :
188 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_section, xc_input=xc_section)
189 :
190 : CALL section_vals_val_get(pwdft_section, "ignore_convergence_failure", &
191 20 : l_val=pwdft_env%ignore_convergence_failure)
192 : ! cp2k should *have* a function that return all xc_functionals. Doing
193 : ! manually is prone to errors
194 :
195 20 : IF (ASSOCIATED(xc_section)) THEN
196 20 : ifun = 0
197 38 : DO
198 58 : ifun = ifun + 1
199 58 : xc_fun => section_vals_get_subs_vals2(xc_section, i_section=ifun)
200 58 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
201 : ! Here, we do not have to check whether the functional name starts with XC_
202 : ! because we only allow the shorter form w/o XC_
203 58 : CALL sirius_add_xc_functional(sctx, "XC_"//TRIM(xc_fun%section%name))
204 : END DO
205 : END IF
206 :
207 : ! import control section
208 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "control")
209 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
210 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "control")
211 20 : CALL section_vals_val_get(pwdft_sub_section, "mpi_grid_dims", i_vals=mpi_grid_dims)
212 : END IF
213 :
214 : ! import parameters section
215 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "parameters")
216 :
217 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
218 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "parameters")
219 20 : CALL section_vals_val_get(pwdft_sub_section, "ngridk", i_vals=kk)
220 20 : k_grid(1) = kk(1)
221 20 : k_grid(2) = kk(2)
222 20 : k_grid(3) = kk(3)
223 :
224 20 : CALL section_vals_val_get(pwdft_sub_section, "shiftk", i_vals=kk)
225 20 : k_shift(1) = kk(1)
226 20 : k_shift(2) = kk(2)
227 20 : k_shift(3) = kk(3)
228 20 : CALL section_vals_val_get(pwdft_sub_section, "num_mag_dims", i_val=num_mag_dims)
229 20 : CALL section_vals_val_get(pwdft_sub_section, "use_symmetry", l_val=use_symmetry)
230 20 : CALL section_vals_val_get(pwdft_sub_section, "so_correction", l_val=use_so)
231 :
232 : ! now check if van der walls corrections are needed
233 : vdw_func = -1
234 : #ifdef __LIBVDWXC
235 20 : CALL section_vals_val_get(pwdft_sub_section, "vdw_functional", i_val=vdw_func)
236 0 : SELECT CASE (vdw_func)
237 : CASE (SIRIUS_FUNC_VDWDF)
238 0 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF")
239 : CASE (SIRIUS_FUNC_VDWDF2)
240 0 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF2")
241 : CASE (SIRIUS_FUNC_VDWDFCX)
242 20 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDFCX")
243 : CASE default
244 : END SELECT
245 : #endif
246 :
247 : END IF
248 :
249 : ! import mixer section
250 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "mixer")
251 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
252 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "mixer")
253 : END IF
254 :
255 : ! import settings section
256 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "settings")
257 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
258 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "settings")
259 : END IF
260 :
261 : ! import solver section
262 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "iterative_solver")
263 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
264 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "iterative_solver")
265 : END IF
266 :
267 : #if defined(__SIRIUS_DFTD3)
268 : ! import dftd3 section
269 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd3")
270 : IF (ASSOCIATED(pwdft_sub_section)) THEN
271 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd3")
272 : END IF
273 : #endif
274 : #if defined(__SIRIUS_DFTD4)
275 : ! import dftd4 section
276 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd4")
277 : IF (ASSOCIATED(pwdft_sub_section)) THEN
278 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd4")
279 : END IF
280 : #endif
281 :
282 : #if defined(__SIRIUS_NLCG)
283 : ! import nlcg section
284 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "nlcg")
285 : IF (ASSOCIATED(pwdft_sub_section)) THEN
286 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "nlcg")
287 : END IF
288 : #endif
289 :
290 : #if defined(__SIRIUS_VCSQNM)
291 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "vcsqnm")
292 : IF (ASSOCIATED(pwdft_sub_section)) THEN
293 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "vcsqnm")
294 : END IF
295 : #endif
296 :
297 : !CALL sirius_dump_runtime_setup(sctx, "runtime.json")
298 20 : CALL sirius_import_parameters(sctx, '{}')
299 :
300 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
301 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
302 20 : CALL qs_subsys_get(qs_subsys, cell=my_cell, use_ref_cell=use_ref_cell)
303 80 : a1(:) = my_cell%hmat(:, 1)
304 80 : a2(:) = my_cell%hmat(:, 2)
305 80 : a3(:) = my_cell%hmat(:, 3)
306 20 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
307 :
308 20 : IF (use_ref_cell) THEN
309 0 : CPWARN("SIRIUS| The specified CELL_REF will be ignored for PW_DFT calculations")
310 : END IF
311 :
312 : ! set up the atomic type definitions
313 : CALL qs_subsys_get(qs_subsys, &
314 : atomic_kind_set=atomic_kind_set, &
315 : qs_kind_set=qs_kind_set, &
316 20 : particle_set=particle_set)
317 20 : nkind = SIZE(atomic_kind_set)
318 46 : DO ikind = 1, nkind
319 : CALL get_atomic_kind(atomic_kind_set(ikind), &
320 26 : name=label, element_symbol=element_symbol, mass=mass)
321 26 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
322 26 : NULLIFY (upf_pot, gth_potential)
323 26 : CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)
324 :
325 26 : IF (ASSOCIATED(upf_pot)) THEN
326 : CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
327 : symbol=element_symbol, &
328 20 : mass=REAL(mass/massunit, KIND=C_DOUBLE))
329 :
330 6 : ELSEIF (ASSOCIATED(gth_potential)) THEN
331 : !
332 6 : NULLIFY (atom_grid)
333 6 : CALL allocate_grid_atom(atom_grid)
334 6 : nmesh = 929
335 6 : atom_grid%nr = nmesh
336 6 : CALL create_grid_atom(atom_grid, nmesh, 1, 1, 0, do_gapw_log)
337 24 : ALLOCATE (rp(nmesh), fun(nmesh))
338 6 : IF (atom_grid%rad(1) < atom_grid%rad(nmesh)) THEN
339 : up = .TRUE.
340 : ELSE
341 : up = .FALSE.
342 : END IF
343 : IF (up) THEN
344 0 : rp(1:nmesh) = atom_grid%rad(1:nmesh)
345 : ELSE
346 5580 : DO i = 1, nmesh
347 5580 : rp(i) = atom_grid%rad(nmesh - i + 1)
348 : END DO
349 : END IF
350 : ! add new atom type
351 : CALL sirius_add_atom_type(sctx, label, &
352 : zn=NINT(zeff + 0.001d0), &
353 : symbol=element_symbol, &
354 : mass=REAL(mass/massunit, KIND=C_DOUBLE), &
355 6 : spin_orbit=gth_potential%soc)
356 : !
357 2916 : ALLOCATE (gth_atompot)
358 6 : CALL gth_potential_conversion(gth_potential, gth_atompot)
359 : ! set radial grid
360 5580 : fun(1:nmesh) = rp(1:nmesh)
361 6 : CALL sirius_set_atom_type_radial_grid(sctx, label, nmesh, fun(1))
362 : ! set beta-projectors
363 : ! GTH SOC uses the same projectors, SIRIUS can use the same or different projectors for l+1/2, l-1/2 (l > 0 l+1/2 l < 0 l-/2 )
364 24 : ALLOCATE (ef(nmesh), beta(nmesh))
365 6 : ibeta = 0
366 30 : DO l = 0, 3
367 24 : IF (gth_atompot%nl(l) == 0) CYCLE
368 14 : rl = gth_atompot%rcnl(l)
369 : ! we need to multiply by r so that data transferred to sirius are r \beta(r) not beta(r)
370 13020 : ef(1:nmesh) = EXP(-0.5_dp*rp(1:nmesh)*rp(1:nmesh)/(rl*rl))
371 50 : DO i = 1, gth_atompot%nl(l)
372 30 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
373 30 : j = l + 2*i - 1
374 30 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
375 27900 : beta(:) = pf*rp**(l + 2*i - 2)*ef
376 30 : ibeta = ibeta + 1
377 27900 : fun(1:nmesh) = beta(1:nmesh)*rp(1:nmesh)
378 : CALL sirius_add_atom_type_radial_function(sctx, label, &
379 30 : "beta", fun(1), nmesh, l=l)
380 : ! we double the number of beta projectors for SO and l>0
381 54 : IF (gth_atompot%soc .AND. l /= 0) THEN
382 : CALL sirius_add_atom_type_radial_function(sctx, label, &
383 10 : "beta", fun(1), nmesh, l=-l)
384 : END IF
385 : END DO
386 : END DO
387 6 : DEALLOCATE (ef, beta)
388 6 : nbeta = ibeta
389 :
390 : ! nonlocal PP matrix elements
391 6 : IF (gth_atompot%soc) THEN
392 2 : nbs = 2*nbeta - gth_atompot%nl(0)
393 8 : ALLOCATE (dion(nbs, nbs))
394 : ELSE
395 16 : ALLOCATE (dion(nbeta, nbeta))
396 : END IF
397 458 : dion = 0.0_dp
398 6 : IF (gth_atompot%soc) THEN
399 2 : ns = gth_atompot%nl(0)
400 2 : IF (ns /= 0) THEN
401 26 : dion(1:ns, 1:ns) = gth_atompot%hnl(1:ns, 1:ns, 0)
402 : END IF
403 8 : DO l = 1, 3
404 6 : IF (gth_atompot%nl(l) == 0) CYCLE
405 16 : DO i = 1, gth_atompot%nl(l)
406 14 : ii = ns + 2*SUM(gth_atompot%nl(1:l - 1))
407 10 : ii = ii + 2*(i - 1) + 1
408 42 : DO j = 1, gth_atompot%nl(l)
409 34 : jj = ns + 2*SUM(gth_atompot%nl(1:l - 1))
410 26 : jj = jj + 2*(j - 1) + 1
411 26 : vnlp = gth_atompot%hnl(i, j, l) + 0.5_dp*l*gth_atompot%knl(i, j, l)
412 26 : vnlm = gth_atompot%hnl(i, j, l) - 0.5_dp*(l + 1)*gth_atompot%knl(i, j, l)
413 26 : dion(ii, jj) = vnlp
414 36 : dion(ii + 1, jj + 1) = vnlm
415 : END DO
416 : END DO
417 : END DO
418 2 : CALL sirius_set_atom_type_dion(sctx, label, nbs, dion(1, 1))
419 : ELSE
420 20 : DO l = 0, 3
421 16 : IF (gth_atompot%nl(l) == 0) CYCLE
422 14 : ibeta = SUM(gth_atompot%nl(0:l - 1)) + 1
423 8 : i = ibeta + gth_atompot%nl(l) - 1
424 52 : dion(ibeta:i, ibeta:i) = gth_atompot%hnl(1:gth_atompot%nl(l), 1:gth_atompot%nl(l), l)
425 : END DO
426 4 : CALL sirius_set_atom_type_dion(sctx, label, nbeta, dion(1, 1))
427 : END IF
428 :
429 6 : DEALLOCATE (dion)
430 :
431 : ! set non-linear core correction
432 6 : IF (gth_atompot%nlcc) THEN
433 0 : ALLOCATE (corden(nmesh), fe(nmesh), rc(nmesh))
434 0 : corden(:) = 0.0_dp
435 0 : n = gth_atompot%nexp_nlcc
436 0 : DO i = 1, n
437 0 : al = gth_atompot%alpha_nlcc(i)
438 0 : rc(:) = rp(:)/al
439 0 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
440 0 : DO j = 1, gth_atompot%nct_nlcc(i)
441 0 : cval = gth_atompot%cval_nlcc(j, i)
442 0 : corden(:) = corden(:) + fe(:)*rc(:)**(2*j - 2)*cval
443 : END DO
444 : END DO
445 0 : fun(1:nmesh) = corden(1:nmesh)*rp(1:nmesh)
446 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_core", &
447 0 : fun(1), nmesh)
448 0 : DEALLOCATE (corden, fe, rc)
449 : END IF
450 :
451 : ! local potential
452 18 : ALLOCATE (locpot(nmesh))
453 5580 : locpot(:) = 0.0_dp
454 6 : CALL atom_local_potential(locpot, gth_atompot, rp)
455 5580 : fun(1:nmesh) = locpot(1:nmesh)
456 : CALL sirius_add_atom_type_radial_function(sctx, label, "vloc", &
457 6 : fun(1), nmesh)
458 6 : DEALLOCATE (locpot)
459 : !
460 6 : NULLIFY (density, wavefunction, wfninfo)
461 : CALL calculate_atomic_orbitals(atomic_kind_set(ikind), qs_kind_set(ikind), &
462 : density=density, wavefunction=wavefunction, &
463 6 : wfninfo=wfninfo, agrid=atom_grid)
464 :
465 : ! set the atomic radial functions
466 22 : DO iwf = 1, SIZE(wavefunction, 2)
467 16 : focc = wfninfo(1, iwf)
468 16 : l = NINT(wfninfo(2, iwf))
469 : ! we can not easily get the principal quantum number
470 16 : nu = -1
471 16 : IF (up) THEN
472 0 : fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
473 : ELSE
474 14880 : DO i = 1, nmesh
475 14880 : fun(i) = wavefunction(nmesh - i + 1, iwf)*rp(i)
476 : END DO
477 : END IF
478 : CALL sirius_add_atom_type_radial_function(sctx, &
479 : label, "ps_atomic_wf", &
480 22 : fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
481 : END DO
482 :
483 : ! set total charge density of a free atom (to compute initial rho(r))
484 6 : IF (up) THEN
485 0 : fun(1:nmesh) = fourpi*density(1:nmesh)*atom_grid%rad(1:nmesh)**2
486 : ELSE
487 5580 : DO i = 1, nmesh
488 5580 : fun(i) = fourpi*density(nmesh - i + 1)*atom_grid%rad(nmesh - i + 1)**2
489 : END DO
490 : END IF
491 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
492 6 : fun(1), nmesh)
493 :
494 6 : IF (ASSOCIATED(density)) DEALLOCATE (density)
495 6 : IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
496 6 : IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
497 :
498 6 : CALL deallocate_grid_atom(atom_grid)
499 6 : DEALLOCATE (rp, fun)
500 6 : DEALLOCATE (gth_atompot)
501 : !
502 : ELSE
503 : CALL cp_abort(__LOCATION__, &
504 0 : "CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
505 : END IF
506 :
507 : CALL get_qs_kind(qs_kind_set(ikind), &
508 : dft_plus_u_atom=dft_plus_u_atom, &
509 : l_of_dft_plus_u=lu, &
510 : n_of_dft_plus_u=nu, &
511 : u_minus_j_target=u_minus_j, &
512 : U_of_dft_plus_u=U_u, &
513 : J_of_dft_plus_u=J_u, &
514 : alpha_of_dft_plus_u=alpha_u, &
515 : beta_of_dft_plus_u=beta_u, &
516 : J0_of_dft_plus_u=J0_u, &
517 26 : occupation_of_dft_plus_u=occ_u)
518 :
519 72 : IF (dft_plus_u_atom) THEN
520 0 : IF (nu < 1) THEN
521 0 : CPABORT("CP2K/SIRIUS (hubbard): principal quantum number not specified")
522 : END IF
523 :
524 0 : IF (lu < 0) THEN
525 0 : CPABORT("CP2K/SIRIUS (hubbard): l can not be negative.")
526 : END IF
527 :
528 0 : IF (occ_u < 0.0) THEN
529 0 : CPABORT("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
530 : END IF
531 :
532 0 : j_t(:) = 0.0
533 0 : IF (ABS(u_minus_j) < 1e-8) THEN
534 0 : j_t(1) = J_u
535 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
536 0 : occ_u, U_u, j_t, alpha_u, beta_u, J0_u)
537 : ELSE
538 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
539 0 : occ_u, u_minus_j, j_t, alpha_u, beta_u, J0_u)
540 : END IF
541 : END IF
542 :
543 : END DO
544 :
545 : ! add atoms to the unit cell
546 : ! WARNING: sirius accepts only fractional coordinates;
547 20 : natom = SIZE(particle_set)
548 58 : DO iatom = 1, natom
549 152 : vr(1:3) = particle_set(iatom)%r(1:3)
550 38 : CALL real_to_scaled(vs, vr, my_cell)
551 38 : atomic_kind => particle_set(iatom)%atomic_kind
552 38 : ikind = atomic_kind%kind_number
553 38 : CALL get_atomic_kind(atomic_kind, name=label)
554 38 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, magnetization=magnetization)
555 : ! angle of magnetization might come from input Atom x y z mx my mz
556 : ! or as an angle?
557 : ! Answer : SIRIUS only accept the magnetization as mx, my, mz
558 38 : IF (num_mag_dims == 3) THEN
559 4 : angle1 = 0.0_dp
560 4 : angle2 = 0.0_dp
561 4 : v1(1) = magnetization*SIN(angle1)*COS(angle2)
562 4 : v1(2) = magnetization*SIN(angle1)*SIN(angle2)
563 4 : v1(3) = magnetization*COS(angle1)
564 : ELSE
565 34 : v1 = 0._dp
566 34 : v1(3) = magnetization
567 : END IF
568 38 : v2(1:3) = vs(1:3)
569 96 : CALL sirius_add_atom(sctx, label, v2(1), v1(1))
570 : END DO
571 :
572 20 : CALL sirius_set_mpi_grid_dims(sctx, 2, mpi_grid_dims)
573 :
574 : ! initialize global variables/indices/arrays/etc. of the simulation
575 20 : CALL sirius_initialize_context(sctx)
576 :
577 : ! strictly speaking the parameter use_symmetry is initialized at the
578 : ! beginning but it does no harm to do it that way
579 20 : IF (use_symmetry) THEN
580 16 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.TRUE., kset_handler=ks_handler)
581 : ELSE
582 4 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.FALSE., kset_handler=ks_handler)
583 : END IF
584 : ! create ground-state class
585 20 : CALL sirius_create_ground_state(ks_handler, gs_handler)
586 :
587 20 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler, ks_handler=ks_handler)
588 : #endif
589 40 : END SUBROUTINE cp_sirius_create_env
590 :
591 : !***************************************************************************************************
592 : !> \brief ...
593 : !> \param pwdft_env ...
594 : !> \param
595 : !> \par History
596 : !> 07.2018 Update the Sirius environment
597 : !> \author JHU
598 : ! **************************************************************************************************
599 20 : SUBROUTINE cp_sirius_update_context(pwdft_env)
600 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
601 :
602 : INTEGER :: iatom, natom
603 : REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v2
604 : REAL(KIND=dp), DIMENSION(3) :: vr, vs
605 : TYPE(cell_type), POINTER :: my_cell
606 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
607 : TYPE(qs_subsys_type), POINTER :: qs_subsys
608 : TYPE(sirius_context_handler) :: sctx
609 : TYPE(sirius_ground_state_handler) :: gs_handler
610 :
611 0 : CPASSERT(ASSOCIATED(pwdft_env))
612 20 : CALL pwdft_env_get(pwdft_env, sctx=sctx, gs_handler=gs_handler)
613 :
614 : ! get current positions and lattice vectors
615 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
616 :
617 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
618 20 : CALL qs_subsys_get(qs_subsys, cell=my_cell)
619 80 : a1(:) = my_cell%hmat(:, 1)
620 80 : a2(:) = my_cell%hmat(:, 2)
621 80 : a3(:) = my_cell%hmat(:, 3)
622 20 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
623 :
624 : ! new atomic positions
625 20 : CALL qs_subsys_get(qs_subsys, particle_set=particle_set)
626 20 : natom = SIZE(particle_set)
627 58 : DO iatom = 1, natom
628 152 : vr(1:3) = particle_set(iatom)%r(1:3)
629 38 : CALL real_to_scaled(vs, vr, my_cell)
630 38 : v2(1:3) = vs(1:3)
631 58 : CALL sirius_set_atom_position(sctx, iatom, v2(1))
632 : END DO
633 :
634 : ! update ground-state class
635 20 : CALL sirius_update_ground_state(gs_handler)
636 :
637 20 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler)
638 :
639 20 : END SUBROUTINE cp_sirius_update_context
640 :
641 : ! **************************************************************************************************
642 : !> \brief ...
643 : !> \param sctx ...
644 : !> \param section ...
645 : !> \param section_name ...
646 : ! **************************************************************************************************
647 100 : SUBROUTINE cp_sirius_fill_in_section(sctx, section, section_name)
648 : TYPE(sirius_context_handler), INTENT(INOUT) :: sctx
649 : TYPE(section_vals_type), POINTER :: section
650 : CHARACTER(*), INTENT(in) :: section_name
651 :
652 : CHARACTER(len=256), TARGET :: option_name
653 : CHARACTER(len=4096) :: description, usage
654 100 : CHARACTER(len=80), DIMENSION(:), POINTER :: tmp
655 : CHARACTER(len=80), TARGET :: str
656 : INTEGER :: ctype, elem, ic, j
657 100 : INTEGER, DIMENSION(:), POINTER :: ivals
658 : INTEGER, TARGET :: enum_length, ival, length, &
659 : num_possible_values, number_of_options
660 : LOGICAL :: explicit
661 100 : LOGICAL, DIMENSION(:), POINTER :: lvals
662 : LOGICAL, TARGET :: found, lval
663 100 : REAL(kind=dp), DIMENSION(:), POINTER :: rvals
664 : REAL(kind=dp), TARGET :: rval
665 :
666 100 : NULLIFY (rvals)
667 100 : NULLIFY (ivals)
668 100 : CALL sirius_option_get_section_length(section_name, number_of_options)
669 :
670 2220 : DO elem = 1, number_of_options
671 2120 : option_name = ''
672 : CALL sirius_option_get_info(section_name, &
673 : elem, &
674 : option_name, &
675 : 256, &
676 : ctype, &
677 : num_possible_values, &
678 : enum_length, &
679 : description, &
680 : 4096, &
681 : usage, &
682 2120 : 4096)
683 :
684 : ! ignore the keyword parametes for sections dftd3 and dftd4.
685 2120 : IF (((section_name == 'dftd3') .OR. (section_name == 'dftd4')) .AND. (option_name == 'parameters')) THEN
686 : CYCLE
687 : END IF
688 :
689 4340 : IF ((option_name /= 'memory_usage') .AND. (option_name /= 'xc_functionals') .AND. (option_name /= 'vk')) THEN
690 2080 : CALL section_vals_val_get(section, option_name, explicit=found)
691 2080 : IF (found) THEN
692 160 : SELECT CASE (ctype)
693 : CASE (SIRIUS_INTEGER_TYPE)
694 160 : CALL section_vals_val_get(section, option_name, i_val=ival)
695 160 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ival))
696 : CASE (SIRIUS_NUMBER_TYPE)
697 124 : CALL section_vals_val_get(section, option_name, r_val=rval)
698 124 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rval))
699 : CASE (SIRIUS_LOGICAL_TYPE)
700 30 : CALL section_vals_val_get(section, option_name, l_val=lval)
701 30 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lval))
702 : CASE (SIRIUS_STRING_TYPE) ! string nightmare
703 124 : str = ''
704 124 : CALL section_vals_val_get(section, option_name, explicit=explicit, c_val=str)
705 124 : str = TRIM(ADJUSTL(str))
706 10044 : DO j = 1, LEN(str)
707 9920 : ic = ICHAR(str(j:j))
708 10044 : IF (ic >= 65 .AND. ic < 90) str(j:j) = CHAR(ic + 32)
709 : END DO
710 :
711 124 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), max_length=LEN_TRIM(str))
712 : CASE (SIRIUS_INTEGER_ARRAY_TYPE)
713 20 : CALL section_vals_val_get(section, option_name, i_vals=ivals)
714 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ivals(1)), &
715 20 : max_length=num_possible_values)
716 : CASE (SIRIUS_NUMBER_ARRAY_TYPE)
717 0 : CALL section_vals_val_get(section, option_name, r_vals=rvals)
718 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rvals(1)), &
719 0 : max_length=num_possible_values)
720 : CASE (SIRIUS_LOGICAL_ARRAY_TYPE)
721 0 : CALL section_vals_val_get(section, option_name, l_vals=lvals)
722 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lvals(1)), &
723 0 : max_length=num_possible_values)
724 : CASE (SIRIUS_STRING_ARRAY_TYPE)
725 0 : CALL section_vals_val_get(section, option_name, explicit=explicit, n_rep_val=length)
726 458 : DO j = 1, length
727 0 : str = ''
728 0 : CALL section_vals_val_get(section, option_name, i_rep_val=j, explicit=explicit, c_vals=tmp)
729 0 : str = TRIM(ADJUSTL(tmp(j)))
730 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), &
731 0 : max_length=LEN_TRIM(str), append=.TRUE.)
732 : END DO
733 : CASE DEFAULT
734 : END SELECT
735 : END IF
736 : END IF
737 : END DO
738 100 : END SUBROUTINE cp_sirius_fill_in_section
739 :
740 : !***************************************************************************************************
741 : !> \brief ...
742 : !> \param pwdft_env ...
743 : !> \param calculate_forces ...
744 : !> \param calculate_stress_tensor ...
745 : !> \param
746 : !> \par History
747 : !> 07.2018 start the Sirius library
748 : !> \author JHU
749 : ! **************************************************************************************************
750 20 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
751 : TYPE(pwdft_environment_type), INTENT(INOUT), &
752 : POINTER :: pwdft_env
753 : LOGICAL, INTENT(IN) :: calculate_forces, calculate_stress_tensor
754 :
755 : INTEGER :: iw, n1, n2
756 : LOGICAL :: do_print, gs_converged
757 : REAL(KIND=C_DOUBLE) :: etotal
758 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: cforces
759 : REAL(KIND=C_DOUBLE), DIMENSION(3, 3) :: cstress
760 : REAL(KIND=dp), DIMENSION(3, 3) :: stress
761 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
762 : TYPE(cp_logger_type), POINTER :: logger
763 : TYPE(pwdft_energy_type), POINTER :: energy
764 : TYPE(section_vals_type), POINTER :: print_section, pwdft_input
765 : TYPE(sirius_ground_state_handler) :: gs_handler
766 :
767 0 : CPASSERT(ASSOCIATED(pwdft_env))
768 :
769 20 : NULLIFY (logger)
770 20 : logger => cp_get_default_logger()
771 20 : iw = cp_logger_get_default_io_unit(logger)
772 :
773 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
774 20 : CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
775 :
776 20 : IF (gs_converged) THEN
777 20 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
778 : ELSE
779 0 : IF (pwdft_env%ignore_convergence_failure) THEN
780 0 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
781 : ELSE
782 0 : CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
783 : END IF
784 : END IF
785 20 : IF (iw > 0) CALL m_flush(iw)
786 :
787 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
788 : etotal = 0.0_C_DOUBLE
789 :
790 20 : CALL sirius_get_energy(gs_handler, 'band-gap', etotal)
791 20 : energy%band_gap = etotal
792 :
793 : etotal = 0.0_C_DOUBLE
794 20 : CALL sirius_get_energy(gs_handler, 'total', etotal)
795 20 : energy%etotal = etotal
796 :
797 : ! extract entropy (TS returned by sirius is always negative, sign
798 : ! convention in QE)
799 : etotal = 0.0_C_DOUBLE
800 20 : CALL sirius_get_energy(gs_handler, 'demet', etotal)
801 20 : energy%entropy = -etotal
802 :
803 20 : IF (calculate_forces) THEN
804 14 : CALL pwdft_env_get(pwdft_env=pwdft_env, forces=forces)
805 14 : n1 = SIZE(forces, 1)
806 14 : n2 = SIZE(forces, 2)
807 :
808 56 : ALLOCATE (cforces(n2, n1))
809 142 : cforces = 0.0_C_DOUBLE
810 14 : CALL sirius_get_forces(gs_handler, 'total', cforces)
811 : ! Sirius computes the forces but cp2k use the gradient everywhere
812 : ! so a minus sign is needed.
813 : ! note also that sirius and cp2k store the forces transpose to each other
814 : ! sirius : forces(coordinates, atoms)
815 : ! cp2k : forces(atoms, coordinates)
816 152 : forces = -TRANSPOSE(cforces(:, :))
817 14 : DEALLOCATE (cforces)
818 : END IF
819 :
820 20 : IF (calculate_stress_tensor) THEN
821 0 : cstress = 0.0_C_DOUBLE
822 0 : CALL sirius_get_stress_tensor(gs_handler, 'total', cstress)
823 0 : stress(1:3, 1:3) = cstress(1:3, 1:3)
824 0 : CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
825 : END IF
826 :
827 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_input)
828 20 : print_section => section_vals_get_subs_vals(pwdft_input, "PRINT")
829 20 : CALL section_vals_get(print_section, explicit=do_print)
830 20 : IF (do_print) THEN
831 2 : CALL cp_sirius_print_results(pwdft_env, print_section)
832 : END IF
833 40 : END SUBROUTINE cp_sirius_energy_force
834 :
835 : !***************************************************************************************************
836 : !> \brief ...
837 : !> \param pwdft_env ...
838 : !> \param print_section ...
839 : !> \param
840 : !> \par History
841 : !> 12.2019 init
842 : !> \author JHU
843 : ! **************************************************************************************************
844 2 : SUBROUTINE cp_sirius_print_results(pwdft_env, print_section)
845 : TYPE(pwdft_environment_type), INTENT(INOUT), &
846 : POINTER :: pwdft_env
847 : TYPE(section_vals_type), POINTER :: print_section
848 :
849 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
850 : INTEGER :: i, ik, iounit, ispn, iterstep, iv, iw, &
851 : nbands, nhist, nkpts, nspins
852 : INTEGER(KIND=C_INT) :: cint
853 : LOGICAL :: append, dos, ionode
854 : REAL(KIND=C_DOUBLE) :: creal
855 2 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: slist
856 : REAL(KIND=dp) :: de, e_fermi(2), emax, emin, eval
857 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkpt
858 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
859 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: energies, occupations
860 : TYPE(cp_logger_type), POINTER :: logger
861 : TYPE(sirius_context_handler) :: sctx
862 : TYPE(sirius_ground_state_handler) :: gs_handler
863 : TYPE(sirius_kpoint_set_handler) :: ks_handler
864 :
865 2 : NULLIFY (logger)
866 4 : logger => cp_get_default_logger()
867 : ionode = logger%para_env%is_source()
868 2 : iounit = cp_logger_get_default_io_unit(logger)
869 :
870 : ! Density of States
871 2 : dos = BTEST(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
872 2 : IF (dos) THEN
873 2 : CALL pwdft_env_get(pwdft_env, ks_handler=ks_handler)
874 2 : CALL pwdft_env_get(pwdft_env, gs_handler=gs_handler)
875 2 : CALL pwdft_env_get(pwdft_env, sctx=sctx)
876 :
877 2 : CALL section_vals_val_get(print_section, "DOS%DELTA_E", r_val=de)
878 2 : CALL section_vals_val_get(print_section, "DOS%APPEND", l_val=append)
879 :
880 2 : CALL sirius_get_num_kpoints(ks_handler, cint)
881 2 : nkpts = cint
882 2 : CALL sirius_get_parameters(sctx, num_bands=cint)
883 2 : nbands = cint
884 2 : CALL sirius_get_parameters(sctx, num_spins=cint)
885 2 : nspins = cint
886 2 : e_fermi(:) = 0.0_dp
887 10 : ALLOCATE (energies(nbands, nspins, nkpts))
888 442 : energies = 0.0_dp
889 8 : ALLOCATE (occupations(nbands, nspins, nkpts))
890 442 : occupations = 0.0_dp
891 6 : ALLOCATE (wkpt(nkpts))
892 6 : ALLOCATE (slist(nbands))
893 10 : DO ik = 1, nkpts
894 8 : CALL sirius_get_kpoint_properties(ks_handler, ik, creal)
895 10 : wkpt(ik) = creal
896 : END DO
897 10 : DO ik = 1, nkpts
898 26 : DO ispn = 1, nspins
899 16 : CALL sirius_get_band_energies(ks_handler, ik, ispn, slist)
900 432 : energies(1:nbands, ispn, ik) = slist(1:nbands)
901 16 : CALL sirius_get_band_occupancies(ks_handler, ik, ispn, slist)
902 440 : occupations(1:nbands, ispn, ik) = slist(1:nbands)
903 : END DO
904 : END DO
905 442 : emin = MINVAL(energies)
906 442 : emax = MAXVAL(energies)
907 2 : nhist = NINT((emax - emin)/de) + 1
908 16 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
909 16110 : hist = 0.0_dp
910 16110 : occval = 0.0_dp
911 16110 : ehist = 0.0_dp
912 :
913 10 : DO ik = 1, nkpts
914 26 : DO ispn = 1, nspins
915 440 : DO i = 1, nbands
916 416 : eval = energies(i, ispn, ik) - emin
917 416 : iv = NINT(eval/de) + 1
918 416 : CPASSERT((iv > 0) .AND. (iv <= nhist))
919 416 : hist(iv, ispn) = hist(iv, ispn) + wkpt(ik)
920 432 : occval(iv, ispn) = occval(iv, ispn) + wkpt(ik)*occupations(i, ispn, ik)
921 : END DO
922 : END DO
923 : END DO
924 16110 : hist = hist/REAL(nbands, KIND=dp)
925 8054 : DO i = 1, nhist
926 24158 : ehist(i, 1:nspins) = emin + (i - 1)*de
927 : END DO
928 :
929 2 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
930 2 : my_act = "WRITE"
931 2 : IF (append .AND. iterstep > 1) THEN
932 0 : my_pos = "APPEND"
933 : ELSE
934 2 : my_pos = "REWIND"
935 : END IF
936 :
937 : iw = cp_print_key_unit_nr(logger, print_section, "DOS", &
938 : extension=".dos", file_position=my_pos, file_action=my_act, &
939 2 : file_form="FORMATTED")
940 2 : IF (iw > 0) THEN
941 1 : IF (nspins == 2) THEN
942 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
943 1 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
944 1 : WRITE (UNIT=iw, FMT="(T2,A, A)") " Energy[a.u.] Alpha_Density Occupation", &
945 2 : " Beta_Density Occupation"
946 : ELSE
947 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
948 0 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
949 0 : WRITE (UNIT=iw, FMT="(T2,A)") " Energy[a.u.] Density Occupation"
950 : END IF
951 4027 : DO i = 1, nhist
952 4026 : eval = emin + (i - 1)*de
953 4027 : IF (nspins == 2) THEN
954 4026 : WRITE (UNIT=iw, FMT="(F15.8,4F15.4)") eval, hist(i, 1), occval(i, 1), &
955 8052 : hist(i, 2), occval(i, 2)
956 : ELSE
957 0 : WRITE (UNIT=iw, FMT="(F15.8,2F15.4)") eval, hist(i, 1), occval(i, 1)
958 : END IF
959 : END DO
960 : END IF
961 2 : CALL cp_print_key_finished_output(iw, logger, print_section, "DOS")
962 :
963 2 : DEALLOCATE (energies, occupations, wkpt, slist)
964 4 : DEALLOCATE (hist, occval, ehist)
965 : END IF
966 4 : END SUBROUTINE cp_sirius_print_results
967 :
968 : END MODULE sirius_interface
969 :
970 : #else
971 :
972 : !***************************************************************************************************
973 : !> \brief Empty implementation in case SIRIUS is not compiled in.
974 : !***************************************************************************************************
975 : MODULE sirius_interface
976 : USE pwdft_environment_types, ONLY: pwdft_environment_type
977 : #include "./base/base_uses.f90"
978 :
979 : IMPLICIT NONE
980 :
981 : PRIVATE
982 :
983 : ! Public subroutines
984 :
985 : PUBLIC :: cp_sirius_create_env, &
986 : cp_sirius_energy_force, &
987 : cp_sirius_finalize, &
988 : cp_sirius_init, &
989 : cp_sirius_is_initialized, &
990 : cp_sirius_update_context
991 :
992 : CONTAINS
993 :
994 : ! **************************************************************************************************
995 : !> \brief Empty implementation in case SIRIUS is not compiled in.
996 : ! **************************************************************************************************
997 : SUBROUTINE cp_sirius_init()
998 : END SUBROUTINE cp_sirius_init
999 :
1000 : ! **************************************************************************************************
1001 : !> \brief Return always .FALSE. because the Sirius library is not compiled in.
1002 : !> \return Return the initialisation status of the Sirius library as boolean
1003 : ! **************************************************************************************************
1004 : LOGICAL FUNCTION cp_sirius_is_initialized()
1005 : cp_sirius_is_initialized = .FALSE.
1006 : END FUNCTION cp_sirius_is_initialized
1007 :
1008 : ! **************************************************************************************************
1009 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1010 : ! **************************************************************************************************
1011 : SUBROUTINE cp_sirius_finalize()
1012 : END SUBROUTINE cp_sirius_finalize
1013 :
1014 : ! **************************************************************************************************
1015 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1016 : !> \param pwdft_env ...
1017 : ! **************************************************************************************************
1018 : SUBROUTINE cp_sirius_create_env(pwdft_env)
1019 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1020 :
1021 : MARK_USED(pwdft_env)
1022 : CPABORT("Sirius library is missing")
1023 : END SUBROUTINE cp_sirius_create_env
1024 :
1025 : ! **************************************************************************************************
1026 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1027 : !> \param pwdft_env ...
1028 : !> \param calculate_forces ...
1029 : !> \param calculate_stress ...
1030 : ! **************************************************************************************************
1031 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
1032 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1033 : LOGICAL :: calculate_forces, calculate_stress
1034 :
1035 : MARK_USED(pwdft_env)
1036 : MARK_USED(calculate_forces)
1037 : MARK_USED(calculate_stress)
1038 : CPABORT("Sirius library is missing")
1039 : END SUBROUTINE cp_sirius_energy_force
1040 :
1041 : ! **************************************************************************************************
1042 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1043 : !> \param pwdft_env ...
1044 : ! **************************************************************************************************
1045 : SUBROUTINE cp_sirius_update_context(pwdft_env)
1046 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1047 :
1048 : MARK_USED(pwdft_env)
1049 : CPABORT("Sirius library is missing")
1050 : END SUBROUTINE cp_sirius_update_context
1051 :
1052 : END MODULE sirius_interface
1053 :
1054 : #endif
|
