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 2 : SUBROUTINE cp_sirius_init()
107 2 : CALL sirius_initialize(.FALSE.)
108 2 : 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 1496 : LOGICAL FUNCTION cp_sirius_is_initialized()
117 1496 : CALL sirius_is_initialized(cp_sirius_is_initialized)
118 1496 : 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 2 : SUBROUTINE cp_sirius_finalize()
126 2 : CALL sirius_finalize(call_mpi_fin=.FALSE., call_device_reset=.FALSE.)
127 2 : 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 24 : 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 24 : INTEGER, DIMENSION(:), POINTER :: mpi_grid_dims
147 : INTEGER(KIND=C_INT), DIMENSION(3) :: k_grid, k_shift
148 24 : INTEGER, DIMENSION(:), POINTER :: kk
149 : LOGICAL :: up, use_ref_cell
150 : LOGICAL(4) :: use_so, use_symmetry, dft_plus_u_atom
151 24 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: fun
152 24 : 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 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, ef, fe, locpot, rc, rp
158 : REAL(KIND=dp), DIMENSION(3) :: vr, vs, j_t
159 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: density
160 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: wavefunction, wfninfo
161 : TYPE(atom_gthpot_type), POINTER :: gth_atompot
162 : TYPE(atom_upfpot_type), POINTER :: upf_pot
163 24 : 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 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
170 24 : 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 24 : output_unit = cp_logger_get_default_io_unit()
181 : ! create context of simulation
182 24 : CALL pwdft_env_get(pwdft_env, para_env=para_env)
183 24 : sirius_mpi_comm = para_env%get_handle()
184 24 : CALL sirius_create_context(sirius_mpi_comm, sctx)
185 :
186 : ! the "fun" starts.
187 :
188 24 : 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 24 : 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 24 : IF (ASSOCIATED(xc_section)) THEN
196 24 : ifun = 0
197 46 : DO
198 70 : ifun = ifun + 1
199 70 : xc_fun => section_vals_get_subs_vals2(xc_section, i_section=ifun)
200 70 : 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 70 : CALL sirius_add_xc_functional(sctx, "XC_"//TRIM(xc_fun%section%name))
204 : END DO
205 : END IF
206 :
207 : ! import control section
208 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "control")
209 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
210 24 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "control")
211 24 : 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 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "parameters")
216 :
217 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
218 24 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "parameters")
219 24 : CALL section_vals_val_get(pwdft_sub_section, "ngridk", i_vals=kk)
220 24 : k_grid(1) = kk(1)
221 24 : k_grid(2) = kk(2)
222 24 : k_grid(3) = kk(3)
223 :
224 24 : CALL section_vals_val_get(pwdft_sub_section, "shiftk", i_vals=kk)
225 24 : k_shift(1) = kk(1)
226 24 : k_shift(2) = kk(2)
227 24 : k_shift(3) = kk(3)
228 24 : CALL section_vals_val_get(pwdft_sub_section, "num_mag_dims", i_val=num_mag_dims)
229 24 : CALL section_vals_val_get(pwdft_sub_section, "use_symmetry", l_val=use_symmetry)
230 24 : 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 24 : 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 24 : 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 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "mixer")
251 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
252 24 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "mixer")
253 : END IF
254 :
255 : ! import settings section
256 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "settings")
257 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
258 24 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "settings")
259 : END IF
260 :
261 : ! import solver section
262 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "iterative_solver")
263 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
264 24 : 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 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd3")
270 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
271 24 : 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 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd4")
277 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
278 24 : 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 24 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "vcsqnm")
292 24 : IF (ASSOCIATED(pwdft_sub_section)) THEN
293 24 : 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 24 : CALL sirius_import_parameters(sctx, '{}')
299 :
300 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
301 24 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
302 24 : CALL qs_subsys_get(qs_subsys, cell=my_cell, use_ref_cell=use_ref_cell)
303 96 : a1(:) = my_cell%hmat(:, 1)
304 96 : a2(:) = my_cell%hmat(:, 2)
305 96 : a3(:) = my_cell%hmat(:, 3)
306 24 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
307 :
308 24 : 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 24 : particle_set=particle_set)
317 24 : nkind = SIZE(atomic_kind_set)
318 54 : DO ikind = 1, nkind
319 : CALL get_atomic_kind(atomic_kind_set(ikind), &
320 30 : name=label, element_symbol=element_symbol, mass=mass)
321 30 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
322 30 : NULLIFY (upf_pot, gth_potential)
323 30 : CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)
324 :
325 30 : IF (ASSOCIATED(upf_pot)) THEN
326 : CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
327 : symbol=element_symbol, &
328 22 : mass=REAL(mass/massunit, KIND=C_DOUBLE))
329 :
330 8 : ELSE IF (ASSOCIATED(gth_potential)) THEN
331 : !
332 8 : NULLIFY (atom_grid)
333 8 : CALL allocate_grid_atom(atom_grid)
334 8 : nmesh = 929
335 8 : atom_grid%nr = nmesh
336 8 : CALL create_grid_atom(atom_grid, nmesh, 1, 1, 0, do_gapw_log)
337 32 : ALLOCATE (rp(nmesh), fun(nmesh))
338 8 : 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 7440 : DO i = 1, nmesh
347 7440 : 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 8 : spin_orbit=gth_potential%soc)
356 : !
357 3888 : ALLOCATE (gth_atompot)
358 8 : CALL gth_potential_conversion(gth_potential, gth_atompot)
359 : ! set radial grid
360 7440 : fun(1:nmesh) = rp(1:nmesh)
361 8 : 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 32 : ALLOCATE (ef(nmesh), beta(nmesh))
365 8 : ibeta = 0
366 40 : DO l = 0, 3
367 32 : IF (gth_atompot%nl(l) == 0) CYCLE
368 18 : 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 16740 : ef(1:nmesh) = EXP(-0.5_dp*rp(1:nmesh)*rp(1:nmesh)/(rl*rl))
371 60 : DO i = 1, gth_atompot%nl(l)
372 34 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
373 34 : j = l + 2*i - 1
374 34 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
375 31620 : beta(:) = pf*rp**(l + 2*i - 2)*ef
376 34 : ibeta = ibeta + 1
377 31620 : fun(1:nmesh) = beta(1:nmesh)*rp(1:nmesh)
378 : CALL sirius_add_atom_type_radial_function(sctx, label, &
379 34 : "beta", fun(1), nmesh, l=l)
380 : ! we double the number of beta projectors for SO and l>0
381 66 : 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 8 : DEALLOCATE (ef, beta)
388 8 : nbeta = ibeta
389 :
390 : ! nonlocal PP matrix elements
391 8 : IF (gth_atompot%soc) THEN
392 2 : nbs = 2*nbeta - gth_atompot%nl(0)
393 8 : ALLOCATE (dion(nbs, nbs))
394 : ELSE
395 24 : ALLOCATE (dion(nbeta, nbeta))
396 : END IF
397 8 : dion = 0.0_dp
398 8 : 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 30 : DO l = 0, 3
421 24 : IF (gth_atompot%nl(l) == 0) CYCLE
422 20 : ibeta = SUM(gth_atompot%nl(0:l - 1)) + 1
423 12 : i = ibeta + gth_atompot%nl(l) - 1
424 66 : dion(ibeta:i, ibeta:i) = gth_atompot%hnl(1:gth_atompot%nl(l), 1:gth_atompot%nl(l), l)
425 : END DO
426 6 : CALL sirius_set_atom_type_dion(sctx, label, nbeta, dion(1, 1))
427 : END IF
428 :
429 8 : DEALLOCATE (dion)
430 :
431 : ! set non-linear core correction
432 8 : IF (gth_atompot%nlcc) THEN
433 10 : ALLOCATE (corden(nmesh), fe(nmesh), rc(nmesh))
434 2 : corden(:) = 0.0_dp
435 2 : n = gth_atompot%nexp_nlcc
436 4 : DO i = 1, n
437 2 : al = gth_atompot%alpha_nlcc(i)
438 1860 : rc(:) = rp(:)/al
439 1860 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
440 6 : DO j = 1, gth_atompot%nct_nlcc(i)
441 2 : cval = gth_atompot%cval_nlcc(j, i)
442 1862 : corden(:) = corden(:) + fe(:)*rc(:)**(2*j - 2)*cval
443 : END DO
444 : END DO
445 : ! SIRIUS expects the volume density rho_core(r), not r*rho_core(r).
446 1860 : fun(1:nmesh) = corden(1:nmesh)
447 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_core", &
448 2 : fun(1), nmesh)
449 2 : DEALLOCATE (corden, fe, rc)
450 : END IF
451 :
452 : ! local potential
453 24 : ALLOCATE (locpot(nmesh))
454 8 : locpot(:) = 0.0_dp
455 8 : CALL atom_local_potential(locpot, gth_atompot, rp)
456 7440 : fun(1:nmesh) = locpot(1:nmesh)
457 : CALL sirius_add_atom_type_radial_function(sctx, label, "vloc", &
458 8 : fun(1), nmesh)
459 8 : DEALLOCATE (locpot)
460 : !
461 8 : NULLIFY (density, wavefunction, wfninfo)
462 : CALL calculate_atomic_orbitals(atomic_kind_set(ikind), qs_kind_set(ikind), &
463 : density=density, wavefunction=wavefunction, &
464 8 : wfninfo=wfninfo, agrid=atom_grid)
465 :
466 : ! set the atomic radial functions
467 28 : DO iwf = 1, SIZE(wavefunction, 2)
468 20 : focc = wfninfo(1, iwf)
469 20 : l = NINT(wfninfo(2, iwf))
470 : ! we can not easily get the principal quantum number
471 20 : nu = -1
472 20 : IF (up) THEN
473 0 : fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
474 : ELSE
475 18600 : DO i = 1, nmesh
476 18600 : fun(i) = wavefunction(nmesh - i + 1, iwf)*rp(i)
477 : END DO
478 : END IF
479 : CALL sirius_add_atom_type_radial_function(sctx, &
480 : label, "ps_atomic_wf", &
481 28 : fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
482 : END DO
483 :
484 : ! set total charge density of a free atom (to compute initial rho(r))
485 8 : IF (up) THEN
486 0 : fun(1:nmesh) = fourpi*density(1:nmesh)*atom_grid%rad(1:nmesh)**2
487 : ELSE
488 7440 : DO i = 1, nmesh
489 7440 : fun(i) = fourpi*density(nmesh - i + 1)*atom_grid%rad(nmesh - i + 1)**2
490 : END DO
491 : END IF
492 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
493 8 : fun(1), nmesh)
494 :
495 8 : IF (ASSOCIATED(density)) DEALLOCATE (density)
496 8 : IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
497 8 : IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
498 :
499 8 : CALL deallocate_grid_atom(atom_grid)
500 8 : DEALLOCATE (rp, fun)
501 8 : DEALLOCATE (gth_atompot)
502 : !
503 : ELSE
504 : CALL cp_abort(__LOCATION__, &
505 0 : "CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
506 : END IF
507 :
508 : CALL get_qs_kind(qs_kind_set(ikind), &
509 : dft_plus_u_atom=dft_plus_u_atom, &
510 : l_of_dft_plus_u=lu, &
511 : n_of_dft_plus_u=nu, &
512 : u_minus_j_target=u_minus_j, &
513 : U_of_dft_plus_u=U_u, &
514 : J_of_dft_plus_u=J_u, &
515 : alpha_of_dft_plus_u=alpha_u, &
516 : beta_of_dft_plus_u=beta_u, &
517 : J0_of_dft_plus_u=J0_u, &
518 30 : occupation_of_dft_plus_u=occ_u)
519 :
520 84 : IF (dft_plus_u_atom) THEN
521 0 : IF (nu < 1) THEN
522 0 : CPABORT("CP2K/SIRIUS (hubbard): principal quantum number not specified")
523 : END IF
524 :
525 0 : IF (lu < 0) THEN
526 0 : CPABORT("CP2K/SIRIUS (hubbard): l can not be negative.")
527 : END IF
528 :
529 0 : IF (occ_u < 0.0) THEN
530 0 : CPABORT("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
531 : END IF
532 :
533 0 : j_t(:) = 0.0
534 0 : IF (ABS(u_minus_j) < 1e-8) THEN
535 0 : j_t(1) = J_u
536 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
537 0 : occ_u, U_u, j_t, alpha_u, beta_u, J0_u)
538 : ELSE
539 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
540 0 : occ_u, u_minus_j, j_t, alpha_u, beta_u, J0_u)
541 : END IF
542 : END IF
543 :
544 : END DO
545 :
546 : ! add atoms to the unit cell
547 : ! WARNING: sirius accepts only fractional coordinates;
548 24 : natom = SIZE(particle_set)
549 70 : DO iatom = 1, natom
550 184 : vr(1:3) = particle_set(iatom)%r(1:3)
551 46 : CALL real_to_scaled(vs, vr, my_cell)
552 46 : atomic_kind => particle_set(iatom)%atomic_kind
553 46 : ikind = atomic_kind%kind_number
554 46 : CALL get_atomic_kind(atomic_kind, name=label)
555 46 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, magnetization=magnetization)
556 : ! angle of magnetization might come from input Atom x y z mx my mz
557 : ! or as an angle?
558 : ! Answer : SIRIUS only accept the magnetization as mx, my, mz
559 46 : IF (num_mag_dims == 3) THEN
560 4 : angle1 = 0.0_dp
561 4 : angle2 = 0.0_dp
562 4 : v1(1) = magnetization*SIN(angle1)*COS(angle2)
563 4 : v1(2) = magnetization*SIN(angle1)*SIN(angle2)
564 4 : v1(3) = magnetization*COS(angle1)
565 : ELSE
566 42 : v1 = 0._dp
567 42 : v1(3) = magnetization
568 : END IF
569 46 : v2(1:3) = vs(1:3)
570 116 : CALL sirius_add_atom(sctx, label, v2(1), v1(1))
571 : END DO
572 :
573 24 : CALL sirius_set_mpi_grid_dims(sctx, 2, mpi_grid_dims)
574 :
575 : ! initialize global variables/indices/arrays/etc. of the simulation
576 24 : CALL sirius_initialize_context(sctx)
577 :
578 : ! strictly speaking the parameter use_symmetry is initialized at the
579 : ! beginning but it does no harm to do it that way
580 24 : IF (use_symmetry) THEN
581 20 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.TRUE., kset_handler=ks_handler)
582 : ELSE
583 4 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.FALSE., kset_handler=ks_handler)
584 : END IF
585 : ! create ground-state class
586 24 : CALL sirius_create_ground_state(ks_handler, gs_handler)
587 :
588 24 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler, ks_handler=ks_handler)
589 : #endif
590 48 : END SUBROUTINE cp_sirius_create_env
591 :
592 : !***************************************************************************************************
593 : !> \brief ...
594 : !> \param pwdft_env ...
595 : !> \param
596 : !> \par History
597 : !> 07.2018 Update the Sirius environment
598 : !> \author JHU
599 : ! **************************************************************************************************
600 24 : SUBROUTINE cp_sirius_update_context(pwdft_env)
601 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
602 :
603 : INTEGER :: iatom, natom
604 : REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v2
605 : REAL(KIND=dp), DIMENSION(3) :: vr, vs
606 : TYPE(cell_type), POINTER :: my_cell
607 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
608 : TYPE(qs_subsys_type), POINTER :: qs_subsys
609 : TYPE(sirius_context_handler) :: sctx
610 : TYPE(sirius_ground_state_handler) :: gs_handler
611 :
612 0 : CPASSERT(ASSOCIATED(pwdft_env))
613 24 : CALL pwdft_env_get(pwdft_env, sctx=sctx, gs_handler=gs_handler)
614 :
615 : ! get current positions and lattice vectors
616 24 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
617 :
618 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
619 24 : CALL qs_subsys_get(qs_subsys, cell=my_cell)
620 96 : a1(:) = my_cell%hmat(:, 1)
621 96 : a2(:) = my_cell%hmat(:, 2)
622 96 : a3(:) = my_cell%hmat(:, 3)
623 24 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
624 :
625 : ! new atomic positions
626 24 : CALL qs_subsys_get(qs_subsys, particle_set=particle_set)
627 24 : natom = SIZE(particle_set)
628 70 : DO iatom = 1, natom
629 184 : vr(1:3) = particle_set(iatom)%r(1:3)
630 46 : CALL real_to_scaled(vs, vr, my_cell)
631 46 : v2(1:3) = vs(1:3)
632 70 : CALL sirius_set_atom_position(sctx, iatom, v2(1))
633 : END DO
634 :
635 : ! update ground-state class
636 24 : CALL sirius_update_ground_state(gs_handler)
637 :
638 24 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler)
639 :
640 24 : END SUBROUTINE cp_sirius_update_context
641 :
642 : ! **************************************************************************************************
643 : !> \brief ...
644 : !> \param sctx ...
645 : !> \param section ...
646 : !> \param section_name ...
647 : ! **************************************************************************************************
648 192 : SUBROUTINE cp_sirius_fill_in_section(sctx, section, section_name)
649 : TYPE(sirius_context_handler), INTENT(INOUT) :: sctx
650 : TYPE(section_vals_type), POINTER :: section
651 : CHARACTER(*), INTENT(in) :: section_name
652 :
653 : CHARACTER(len=256), TARGET :: option_name
654 : CHARACTER(len=4096) :: description, usage
655 192 : CHARACTER(len=80), DIMENSION(:), POINTER :: tmp
656 : CHARACTER(len=80), TARGET :: str
657 : INTEGER :: ctype, elem, ic, j
658 192 : INTEGER, DIMENSION(:), POINTER :: ivals
659 : INTEGER, TARGET :: enum_length, ival, length, &
660 : num_possible_values, number_of_options
661 : LOGICAL :: explicit
662 192 : LOGICAL, DIMENSION(:), POINTER :: lvals
663 : LOGICAL, TARGET :: found, lval
664 192 : REAL(kind=dp), DIMENSION(:), POINTER :: rvals
665 : REAL(kind=dp), TARGET :: rval
666 :
667 192 : NULLIFY (rvals)
668 192 : NULLIFY (ivals)
669 192 : CALL sirius_option_get_section_length(section_name, number_of_options)
670 :
671 3192 : DO elem = 1, number_of_options
672 3000 : option_name = ''
673 : CALL sirius_option_get_info(section_name, &
674 : elem, &
675 : option_name, &
676 : 256, &
677 : ctype, &
678 : num_possible_values, &
679 : enum_length, &
680 : description, &
681 : 4096, &
682 : usage, &
683 3000 : 4096)
684 :
685 : ! ignore the keyword parametes for sections dftd3 and dftd4.
686 3000 : IF (((section_name == 'dftd3') .OR. (section_name == 'dftd4')) .AND. (option_name == 'parameters')) THEN
687 : CYCLE
688 : END IF
689 :
690 6144 : IF ((option_name /= 'memory_usage') .AND. (option_name /= 'xc_functionals') .AND. (option_name /= 'vk')) THEN
691 2904 : CALL section_vals_val_get(section, option_name, explicit=found)
692 2904 : IF (found) THEN
693 180 : SELECT CASE (ctype)
694 : CASE (SIRIUS_INTEGER_TYPE)
695 180 : CALL section_vals_val_get(section, option_name, i_val=ival)
696 180 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ival))
697 : CASE (SIRIUS_NUMBER_TYPE)
698 152 : CALL section_vals_val_get(section, option_name, r_val=rval)
699 152 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rval))
700 : CASE (SIRIUS_LOGICAL_TYPE)
701 30 : CALL section_vals_val_get(section, option_name, l_val=lval)
702 30 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lval))
703 : CASE (SIRIUS_STRING_TYPE) ! string nightmare
704 148 : str = ''
705 148 : CALL section_vals_val_get(section, option_name, explicit=explicit, c_val=str)
706 148 : str = TRIM(ADJUSTL(str))
707 11988 : DO j = 1, LEN(str)
708 11840 : ic = ICHAR(str(j:j))
709 11988 : IF (ic >= 65 .AND. ic < 90) str(j:j) = CHAR(ic + 32)
710 : END DO
711 :
712 148 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), max_length=LEN_TRIM(str))
713 : CASE (SIRIUS_INTEGER_ARRAY_TYPE)
714 24 : CALL section_vals_val_get(section, option_name, i_vals=ivals)
715 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ivals(1)), &
716 24 : max_length=num_possible_values)
717 : CASE (SIRIUS_NUMBER_ARRAY_TYPE)
718 0 : CALL section_vals_val_get(section, option_name, r_vals=rvals)
719 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rvals(1)), &
720 0 : max_length=num_possible_values)
721 : CASE (SIRIUS_LOGICAL_ARRAY_TYPE)
722 0 : CALL section_vals_val_get(section, option_name, l_vals=lvals)
723 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lvals(1)), &
724 0 : max_length=num_possible_values)
725 : CASE (SIRIUS_STRING_ARRAY_TYPE)
726 0 : CALL section_vals_val_get(section, option_name, explicit=explicit, n_rep_val=length)
727 534 : DO j = 1, length
728 0 : str = ''
729 0 : CALL section_vals_val_get(section, option_name, i_rep_val=j, explicit=explicit, c_vals=tmp)
730 0 : str = TRIM(ADJUSTL(tmp(j)))
731 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), &
732 0 : max_length=LEN_TRIM(str), append=.TRUE.)
733 : END DO
734 : CASE DEFAULT
735 : END SELECT
736 : END IF
737 : END IF
738 : END DO
739 192 : END SUBROUTINE cp_sirius_fill_in_section
740 :
741 : !***************************************************************************************************
742 : !> \brief ...
743 : !> \param pwdft_env ...
744 : !> \param calculate_forces ...
745 : !> \param calculate_stress_tensor ...
746 : !> \param
747 : !> \par History
748 : !> 07.2018 start the Sirius library
749 : !> \author JHU
750 : ! **************************************************************************************************
751 24 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
752 : TYPE(pwdft_environment_type), INTENT(INOUT), &
753 : POINTER :: pwdft_env
754 : LOGICAL, INTENT(IN) :: calculate_forces, calculate_stress_tensor
755 :
756 : INTEGER :: iw, n1, n2
757 : LOGICAL :: do_print, gs_converged
758 : REAL(KIND=C_DOUBLE) :: etotal
759 24 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: cforces
760 : REAL(KIND=C_DOUBLE), DIMENSION(3, 3) :: cstress
761 : REAL(KIND=dp), DIMENSION(3, 3) :: stress
762 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
763 : TYPE(cp_logger_type), POINTER :: logger
764 : TYPE(pwdft_energy_type), POINTER :: energy
765 : TYPE(section_vals_type), POINTER :: print_section, pwdft_input
766 : TYPE(sirius_ground_state_handler) :: gs_handler
767 :
768 0 : CPASSERT(ASSOCIATED(pwdft_env))
769 :
770 24 : NULLIFY (logger)
771 24 : logger => cp_get_default_logger()
772 24 : iw = cp_logger_get_default_io_unit(logger)
773 :
774 24 : CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
775 24 : CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
776 :
777 24 : IF (gs_converged) THEN
778 24 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
779 : ELSE
780 0 : IF (pwdft_env%ignore_convergence_failure) THEN
781 0 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
782 : ELSE
783 0 : CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
784 : END IF
785 : END IF
786 24 : IF (iw > 0) CALL m_flush(iw)
787 :
788 24 : CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
789 : etotal = 0.0_C_DOUBLE
790 :
791 24 : CALL sirius_get_energy(gs_handler, 'band-gap', etotal)
792 24 : energy%band_gap = etotal
793 :
794 : etotal = 0.0_C_DOUBLE
795 24 : CALL sirius_get_energy(gs_handler, 'total', etotal)
796 24 : energy%etotal = etotal
797 :
798 : ! extract entropy (TS returned by sirius is always negative, sign
799 : ! convention in QE)
800 : etotal = 0.0_C_DOUBLE
801 24 : CALL sirius_get_energy(gs_handler, 'demet', etotal)
802 24 : energy%entropy = -etotal
803 :
804 24 : IF (calculate_forces) THEN
805 18 : CALL pwdft_env_get(pwdft_env=pwdft_env, forces=forces)
806 18 : n1 = SIZE(forces, 1)
807 18 : n2 = SIZE(forces, 2)
808 :
809 72 : ALLOCATE (cforces(n2, n1))
810 18 : cforces = 0.0_C_DOUBLE
811 18 : CALL sirius_get_forces(gs_handler, 'total', cforces)
812 : ! Sirius computes the forces but cp2k use the gradient everywhere
813 : ! so a minus sign is needed.
814 : ! note also that sirius and cp2k store the forces transpose to each other
815 : ! sirius : forces(coordinates, atoms)
816 : ! cp2k : forces(atoms, coordinates)
817 192 : forces = -TRANSPOSE(cforces(:, :))
818 18 : DEALLOCATE (cforces)
819 : END IF
820 :
821 24 : IF (calculate_stress_tensor) THEN
822 0 : cstress = 0.0_C_DOUBLE
823 0 : CALL sirius_get_stress_tensor(gs_handler, 'total', cstress)
824 0 : stress(1:3, 1:3) = cstress(1:3, 1:3)
825 0 : CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
826 : END IF
827 :
828 24 : CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_input)
829 24 : print_section => section_vals_get_subs_vals(pwdft_input, "PRINT")
830 24 : CALL section_vals_get(print_section, explicit=do_print)
831 24 : IF (do_print) THEN
832 2 : CALL cp_sirius_print_results(pwdft_env, print_section)
833 : END IF
834 48 : END SUBROUTINE cp_sirius_energy_force
835 :
836 : !***************************************************************************************************
837 : !> \brief ...
838 : !> \param pwdft_env ...
839 : !> \param print_section ...
840 : !> \param
841 : !> \par History
842 : !> 12.2019 init
843 : !> \author JHU
844 : ! **************************************************************************************************
845 2 : SUBROUTINE cp_sirius_print_results(pwdft_env, print_section)
846 : TYPE(pwdft_environment_type), INTENT(INOUT), &
847 : POINTER :: pwdft_env
848 : TYPE(section_vals_type), POINTER :: print_section
849 :
850 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
851 : INTEGER :: i, ik, iounit, ispn, iterstep, iv, iw, &
852 : nbands, nhist, nkpts, nspins
853 : INTEGER(KIND=C_INT) :: cint
854 : LOGICAL :: append, dos, ionode
855 : REAL(KIND=C_DOUBLE) :: creal
856 2 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: slist
857 : REAL(KIND=dp) :: de, e_fermi(2), emax, emin, eval
858 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkpt
859 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
860 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: energies, occupations
861 : TYPE(cp_logger_type), POINTER :: logger
862 : TYPE(sirius_context_handler) :: sctx
863 : TYPE(sirius_ground_state_handler) :: gs_handler
864 : TYPE(sirius_kpoint_set_handler) :: ks_handler
865 :
866 2 : NULLIFY (logger)
867 4 : logger => cp_get_default_logger()
868 : ionode = logger%para_env%is_source()
869 2 : iounit = cp_logger_get_default_io_unit(logger)
870 :
871 : ! Density of States
872 2 : dos = BTEST(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
873 2 : IF (dos) THEN
874 2 : CALL pwdft_env_get(pwdft_env, ks_handler=ks_handler)
875 2 : CALL pwdft_env_get(pwdft_env, gs_handler=gs_handler)
876 2 : CALL pwdft_env_get(pwdft_env, sctx=sctx)
877 :
878 2 : CALL section_vals_val_get(print_section, "DOS%DELTA_E", r_val=de)
879 2 : CALL section_vals_val_get(print_section, "DOS%APPEND", l_val=append)
880 :
881 2 : CALL sirius_get_num_kpoints(ks_handler, cint)
882 2 : nkpts = cint
883 2 : CALL sirius_get_parameters(sctx, num_bands=cint)
884 2 : nbands = cint
885 2 : CALL sirius_get_parameters(sctx, num_spins=cint)
886 2 : nspins = cint
887 2 : e_fermi(:) = 0.0_dp
888 10 : ALLOCATE (energies(nbands, nspins, nkpts))
889 2 : energies = 0.0_dp
890 8 : ALLOCATE (occupations(nbands, nspins, nkpts))
891 2 : occupations = 0.0_dp
892 6 : ALLOCATE (wkpt(nkpts))
893 6 : ALLOCATE (slist(nbands))
894 10 : DO ik = 1, nkpts
895 8 : CALL sirius_get_kpoint_properties(ks_handler, ik, creal)
896 10 : wkpt(ik) = creal
897 : END DO
898 10 : DO ik = 1, nkpts
899 26 : DO ispn = 1, nspins
900 16 : CALL sirius_get_band_energies(ks_handler, ik, ispn, slist)
901 432 : energies(1:nbands, ispn, ik) = slist(1:nbands)
902 16 : CALL sirius_get_band_occupancies(ks_handler, ik, ispn, slist)
903 440 : occupations(1:nbands, ispn, ik) = slist(1:nbands)
904 : END DO
905 : END DO
906 442 : emin = MINVAL(energies)
907 442 : emax = MAXVAL(energies)
908 2 : nhist = NINT((emax - emin)/de) + 1
909 16 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
910 2 : hist = 0.0_dp
911 2 : occval = 0.0_dp
912 : ehist = 0.0_dp
913 :
914 10 : DO ik = 1, nkpts
915 26 : DO ispn = 1, nspins
916 440 : DO i = 1, nbands
917 416 : eval = energies(i, ispn, ik) - emin
918 416 : iv = NINT(eval/de) + 1
919 416 : CPASSERT((iv > 0) .AND. (iv <= nhist))
920 416 : hist(iv, ispn) = hist(iv, ispn) + wkpt(ik)
921 432 : occval(iv, ispn) = occval(iv, ispn) + wkpt(ik)*occupations(i, ispn, ik)
922 : END DO
923 : END DO
924 : END DO
925 16110 : hist = hist/REAL(nbands, KIND=dp)
926 8054 : DO i = 1, nhist
927 24158 : ehist(i, 1:nspins) = emin + (i - 1)*de
928 : END DO
929 :
930 2 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
931 2 : my_act = "WRITE"
932 2 : IF (append .AND. iterstep > 1) THEN
933 0 : my_pos = "APPEND"
934 : ELSE
935 2 : my_pos = "REWIND"
936 : END IF
937 :
938 : iw = cp_print_key_unit_nr(logger, print_section, "DOS", &
939 : extension=".dos", file_position=my_pos, file_action=my_act, &
940 2 : file_form="FORMATTED")
941 2 : IF (iw > 0) THEN
942 1 : IF (nspins == 2) THEN
943 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
944 1 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
945 1 : WRITE (UNIT=iw, FMT="(T2,A, A)") " Energy[a.u.] Alpha_Density Occupation", &
946 2 : " Beta_Density Occupation"
947 : ELSE
948 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
949 0 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
950 0 : WRITE (UNIT=iw, FMT="(T2,A)") " Energy[a.u.] Density Occupation"
951 : END IF
952 4027 : DO i = 1, nhist
953 4026 : eval = emin + (i - 1)*de
954 4027 : IF (nspins == 2) THEN
955 4026 : WRITE (UNIT=iw, FMT="(F15.8,4F15.4)") eval, hist(i, 1), occval(i, 1), &
956 8052 : hist(i, 2), occval(i, 2)
957 : ELSE
958 0 : WRITE (UNIT=iw, FMT="(F15.8,2F15.4)") eval, hist(i, 1), occval(i, 1)
959 : END IF
960 : END DO
961 : END IF
962 2 : CALL cp_print_key_finished_output(iw, logger, print_section, "DOS")
963 :
964 2 : DEALLOCATE (energies, occupations, wkpt, slist)
965 4 : DEALLOCATE (hist, occval, ehist)
966 : END IF
967 4 : END SUBROUTINE cp_sirius_print_results
968 :
969 : END MODULE sirius_interface
970 :
971 : #else
972 :
973 : !***************************************************************************************************
974 : !> \brief Empty implementation in case SIRIUS is not compiled in.
975 : !***************************************************************************************************
976 : MODULE sirius_interface
977 : USE pwdft_environment_types, ONLY: pwdft_environment_type
978 : #include "./base/base_uses.f90"
979 :
980 : IMPLICIT NONE
981 :
982 : PRIVATE
983 :
984 : ! Public subroutines
985 :
986 : PUBLIC :: cp_sirius_create_env, &
987 : cp_sirius_energy_force, &
988 : cp_sirius_finalize, &
989 : cp_sirius_init, &
990 : cp_sirius_is_initialized, &
991 : cp_sirius_update_context
992 :
993 : CONTAINS
994 :
995 : ! **************************************************************************************************
996 : !> \brief Empty implementation in case SIRIUS is not compiled in.
997 : ! **************************************************************************************************
998 : SUBROUTINE cp_sirius_init()
999 : END SUBROUTINE cp_sirius_init
1000 :
1001 : ! **************************************************************************************************
1002 : !> \brief Return always .FALSE. because the Sirius library is not compiled in.
1003 : !> \return Return the initialisation status of the Sirius library as boolean
1004 : ! **************************************************************************************************
1005 : LOGICAL FUNCTION cp_sirius_is_initialized()
1006 : cp_sirius_is_initialized = .FALSE.
1007 : END FUNCTION cp_sirius_is_initialized
1008 :
1009 : ! **************************************************************************************************
1010 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1011 : ! **************************************************************************************************
1012 : SUBROUTINE cp_sirius_finalize()
1013 : END SUBROUTINE cp_sirius_finalize
1014 :
1015 : ! **************************************************************************************************
1016 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1017 : !> \param pwdft_env ...
1018 : ! **************************************************************************************************
1019 : SUBROUTINE cp_sirius_create_env(pwdft_env)
1020 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1021 :
1022 : MARK_USED(pwdft_env)
1023 : CPABORT("Sirius library is missing")
1024 : END SUBROUTINE cp_sirius_create_env
1025 :
1026 : ! **************************************************************************************************
1027 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1028 : !> \param pwdft_env ...
1029 : !> \param calculate_forces ...
1030 : !> \param calculate_stress ...
1031 : ! **************************************************************************************************
1032 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
1033 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1034 : LOGICAL :: calculate_forces, calculate_stress
1035 :
1036 : MARK_USED(pwdft_env)
1037 : MARK_USED(calculate_forces)
1038 : MARK_USED(calculate_stress)
1039 : CPABORT("Sirius library is missing")
1040 : END SUBROUTINE cp_sirius_energy_force
1041 :
1042 : ! **************************************************************************************************
1043 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1044 : !> \param pwdft_env ...
1045 : ! **************************************************************************************************
1046 : SUBROUTINE cp_sirius_update_context(pwdft_env)
1047 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1048 :
1049 : MARK_USED(pwdft_env)
1050 : CPABORT("Sirius library is missing")
1051 : END SUBROUTINE cp_sirius_update_context
1052 :
1053 : END MODULE sirius_interface
1054 :
1055 : #endif
|
