63 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_dos'
65 CHARACTER(LEN=20) :: fmtstr_data
66 CHARACTER(LEN=default_string_length) :: my_act, my_pos
67 INTEGER :: handle, i, iounit, ispin, iterstep, iv, &
68 iw, ndigits, nhist, nmo(2), nspins
69 LOGICAL :: append, ionode, should_output
70 REAL(kind=
dp) :: de, e1, e2, e_fermi(2), emax, emin, eval
71 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ehist, hist, occval
72 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues, occupation_numbers
78 ionode = logger%para_env%is_source()
82 IF ((.NOT. should_output))
RETURN
84 CALL timeset(routinen, handle)
85 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
87 IF (iounit > 0)
WRITE (unit=iounit, fmt=
'(/,(T3,A,T61,I10))') &
88 " Calculate DOS at iteration step ", iterstep
93 IF (append .AND. iterstep > 1)
THEN
98 ndigits = min(max(ndigits, 1), 10)
107 CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
108 eigenvalues => mo_set%eigenvalues
109 e1 = minval(eigenvalues(1:nmo(ispin)))
110 e2 = maxval(eigenvalues(1:nmo(ispin)))
115 IF (de > 0.0_dp)
THEN
116 nhist = nint((emax - emin)/de) + 1
117 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
123 occupation_numbers => mo_set%occupation_numbers
124 eigenvalues => mo_set%eigenvalues
126 eval = eigenvalues(i) - emin
127 iv = nint(eval/de) + 1
128 cpassert((iv > 0) .AND. (iv <= nhist))
129 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
130 occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
132 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=
dp)
135 ehist(i, 1:nspins) = emin + (i - 1)*de
139 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
145 occupation_numbers => mo_set%occupation_numbers
146 eigenvalues => mo_set%eigenvalues
148 ehist(i, ispin) = eigenvalues(i)
149 hist(i, ispin) = 1.0_dp
150 occval(i, ispin) = occupation_numbers(i)
152 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=
dp)
158 extension=
".dos", file_position=my_pos, file_action=my_act, &
159 file_form=
"FORMATTED")
161 IF (nspins == 2)
THEN
162 WRITE (unit=iw, fmt=
"(T2,A,I0,A,2F12.6)") &
163 "# DOS at iteration step i = ", iterstep,
", E_Fermi[a.u.] = ", e_fermi(1:2)
164 WRITE (unit=iw, fmt=
"(T2,A, A)")
"# Energy[a.u.] Alpha_Density Occupation", &
165 " Energy[a.u.] Beta_Density Occupation"
167 WRITE (unit=fmtstr_data, fmt=
"(A,I0,A)")
"(2(F15.8,2F15.", ndigits,
"))"
169 WRITE (unit=iw, fmt=
"(T2,A,I0,A,F12.6)") &
170 "# DOS at iteration step i = ", iterstep,
", E_Fermi[a.u.] = ", e_fermi(1)
171 WRITE (unit=iw, fmt=
"(T2,A)")
"# Energy[a.u.] Density Occupation"
173 WRITE (unit=fmtstr_data, fmt=
"(A,I0,A)")
"(F15.8,2F15.", ndigits,
")"
176 IF (nspins == 2)
THEN
180 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
181 e2, hist(i, 2), occval(i, 2)
185 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
190 DEALLOCATE (hist, occval, ehist)
192 CALL timestop(handle)
210 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_dos_kp'
212 CHARACTER(LEN=16) :: fmtstr_data
213 CHARACTER(LEN=default_string_length) :: my_act, my_pos
214 INTEGER :: handle, i, ik, iounit, ispin, iterstep, &
215 iv, iw, ndigits, nhist, nmo(2), &
217 INTEGER,
DIMENSION(:),
POINTER :: nkp_grid
218 LOGICAL :: append, explicit, ionode, should_output
219 REAL(kind=
dp) :: de, e1, e2, emax, emin, eval, wkp
220 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ehist, hist, occval
221 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues, occupation_numbers
229 NULLIFY (logger, kpoints)
231 ionode = logger%para_env%is_source()
235 IF ((.NOT. should_output))
RETURN
237 CALL timeset(routinen, handle)
238 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
246 IF (nkp_grid(i) < 1)
THEN
247 WRITE (unit=iounit, fmt=
'(T4,A,I3,A,I1)') &
248 "Invalid kpoint grid for DOS ", nkp_grid(i),
" in dimension ", i
259 IF (iounit > 0)
WRITE (unit=iounit, fmt=
'(/,(T3,A,T61,I10))') &
260 " Calculate DOS at iteration step ", iterstep
261 IF (iounit > 0)
WRITE (unit=iounit, fmt=
'((T3,A,3I3,A))') &
262 " Using a", kpoints%nkp_grid(:),
' '//trim(kpoints%kp_scheme)//
' grid'
268 de = max(de, 0.00001_dp)
269 IF (append .AND. iterstep > 1)
THEN
274 ndigits = min(max(ndigits, 1), 10)
276 CALL get_qs_env(qs_env, dft_control=dft_control)
277 nspins = dft_control%nspins
278 para_env => kpoints%para_env_inter_kp
283 IF (kpoints%nkp /= 0)
THEN
284 DO ik = 1,
SIZE(kpoints%kp_env)
285 mos => kpoints%kp_env(ik)%kpoint_env%mos
286 cpassert(
ASSOCIATED(mos))
288 mo_set => mos(1, ispin)
290 eigenvalues => mo_set%eigenvalues
291 e1 = minval(eigenvalues(1:nmo_kp))
292 e2 = maxval(eigenvalues(1:nmo_kp))
295 nmo(ispin) = max(nmo(ispin), nmo_kp)
299 CALL para_env%min(emin)
300 CALL para_env%max(emax)
301 CALL para_env%max(nmo)
303 nhist = nint((emax - emin)/de) + 1
304 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
309 IF (kpoints%nkp /= 0)
THEN
310 DO ik = 1,
SIZE(kpoints%kp_env)
311 mos => kpoints%kp_env(ik)%kpoint_env%mos
312 wkp = kpoints%kp_env(ik)%kpoint_env%wkp
314 mo_set => mos(1, ispin)
315 occupation_numbers => mo_set%occupation_numbers
316 eigenvalues => mo_set%eigenvalues
318 eval = eigenvalues(i) - emin
319 iv = nint(eval/de) + 1
320 cpassert((iv > 0) .AND. (iv <= nhist))
321 hist(iv, ispin) = hist(iv, ispin) + wkp
322 occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
327 CALL para_env%sum(hist)
328 CALL para_env%sum(occval)
330 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=
dp)
333 ehist(i, 1:nspins) = emin + (i - 1)*de
338 extension=
".dos", file_position=my_pos, file_action=my_act, &
339 file_form=
"FORMATTED")
341 IF (nspins == 2)
THEN
342 WRITE (unit=iw, fmt=
"(T2,A,I0)")
"# DOS at iteration step i = ", iterstep
343 WRITE (unit=iw, fmt=
"(T2,A,A)")
"# Energy[a.u.] Alpha_Density Occupation", &
344 " Beta_Density Occupation"
346 WRITE (unit=fmtstr_data, fmt=
"(A,I0,A)")
"(F15.8,4F15.", ndigits,
")"
348 WRITE (unit=iw, fmt=
"(T2,A,I0)")
"# DOS at iteration step i = ", iterstep
349 WRITE (unit=iw, fmt=
"(T2,A)")
"# Energy[a.u.] Density Occupation"
351 WRITE (unit=fmtstr_data, fmt=
"(A,I0,A)")
"(F15.8,2F15.", ndigits,
")"
354 eval = emin + (i - 1)*de
355 IF (nspins == 2)
THEN
357 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
358 hist(i, 2), occval(i, 2)
361 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
366 DEALLOCATE (hist, occval)
373 CALL timestop(handle)
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.