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 Fractional occupation number weighted density (Grimme and Hansen).
10 : !> \author Falk William Seidel (initial implementation)
11 : ! **************************************************************************************************
12 : MODULE qs_fod
13 : USE bibliography, ONLY: Grimme2015FOD,&
14 : cite_reference
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 : dbcsr_p_type,&
18 : dbcsr_release,&
19 : dbcsr_type
20 : USE cp_log_handling, ONLY: cp_logger_type
21 : USE cp_output_handling, ONLY: cp_p_file,&
22 : cp_print_key_finished_output,&
23 : cp_print_key_should_output,&
24 : cp_print_key_unit_nr
25 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
26 : USE input_constants, ONLY: smear_fermi_dirac
27 : USE input_section_types, ONLY: section_get_ivals,&
28 : section_vals_get_subs_vals,&
29 : section_vals_type,&
30 : section_vals_val_get
31 : USE kinds, ONLY: dp
32 : USE particle_list_types, ONLY: particle_list_type
33 : USE pw_env_types, ONLY: pw_env_get,&
34 : pw_env_type
35 : USE pw_pool_types, ONLY: pw_pool_type
36 : USE pw_types, ONLY: pw_c1d_gs_type,&
37 : pw_r3d_rs_type
38 : USE qs_collocate_density, ONLY: calculate_rho_elec
39 : USE qs_density_matrices, ONLY: calculate_density_matrix
40 : USE qs_environment_types, ONLY: get_qs_env,&
41 : qs_environment_type
42 : USE qs_kind_types, ONLY: get_qs_kind,&
43 : qs_kind_type
44 : USE qs_ks_types, ONLY: qs_ks_env_type
45 : USE qs_mo_types, ONLY: deallocate_mo_set,&
46 : duplicate_mo_set,&
47 : mo_set_type
48 : USE qs_rho_types, ONLY: qs_rho_get,&
49 : qs_rho_type
50 : USE qs_subsys_types, ONLY: qs_subsys_get,&
51 : qs_subsys_type
52 : USE scf_control_types, ONLY: scf_control_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fod'
58 : PUBLIC :: qs_scf_post_fod, qs_fod_validate, fod_weights
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief Hole weights below the chemical potential, particle weights above it.
64 : !> \param occupation Fermi-Dirac occupations (0 to maxocc).
65 : !> \param eigenvalues Orbital energies in the same units as mu.
66 : !> \param mu Chemical potential for this spin channel.
67 : !> \param maxocc Maximum occupation, two for RKS and one for UKS.
68 : !> \return Nonnegative FOD weights; both branches agree at half occupation.
69 : ! **************************************************************************************************
70 16 : PURE FUNCTION fod_weights(occupation, eigenvalues, mu, maxocc) RESULT(weights)
71 : REAL(dp), DIMENSION(:), INTENT(IN) :: occupation, eigenvalues
72 : REAL(dp), INTENT(IN) :: mu, maxocc
73 : REAL(dp), DIMENSION(SIZE(occupation)) :: weights
74 :
75 104 : WHERE (eigenvalues <= mu)
76 : weights = maxocc - occupation
77 : ELSE WHERE
78 : weights = occupation
79 : END WHERE
80 16 : END FUNCTION fod_weights
81 :
82 : ! **************************************************************************************************
83 : !> \brief Reject unsupported FOD settings before starting an expensive SCF calculation.
84 : !> \param input FORCE_EVAL input section.
85 : !> \param logger Output logger.
86 : !> \param qs_env Quickstep environment.
87 : ! **************************************************************************************************
88 24795 : SUBROUTINE qs_fod_validate(input, logger, qs_env)
89 : TYPE(section_vals_type), POINTER :: input
90 : TYPE(cp_logger_type), POINTER :: logger
91 : TYPE(qs_environment_type), POINTER :: qs_env
92 :
93 : LOGICAL :: do_kpoints
94 : TYPE(dft_control_type), POINTER :: dft_control
95 : TYPE(scf_control_type), POINTER :: scf_control
96 : TYPE(section_vals_type), POINTER :: fod_section
97 :
98 24795 : fod_section => section_vals_get_subs_vals(input, "DFT%PRINT%FOD")
99 24795 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, fod_section, ""), cp_p_file)) RETURN
100 8 : CALL get_qs_env(qs_env, dft_control=dft_control, scf_control=scf_control, do_kpoints=do_kpoints)
101 8 : IF (do_kpoints) CPABORT("FOD currently supports Gamma-point calculations only.")
102 8 : IF (.NOT. dft_control%qs_control%gpw) CPABORT("FOD currently requires the GPW method.")
103 8 : IF (qs_env%run_rtp) CPABORT("FOD is not supported for real-time propagation.")
104 8 : IF (.NOT. scf_control%smear%do_smear .OR. scf_control%smear%method /= smear_fermi_dirac) THEN
105 0 : CPABORT("FOD requires SCF%SMEAR with METHOD FERMI_DIRAC.")
106 : END IF
107 8 : IF (scf_control%smear%electronic_temperature <= 0.0_dp) THEN
108 0 : CPABORT("FOD requires a positive electronic temperature.")
109 : END IF
110 : END SUBROUTINE qs_fod_validate
111 :
112 : ! **************************************************************************************************
113 : !> \brief Print the orbital FOD sum and, optionally, its real-space density.
114 : !> \param input FORCE_EVAL input section.
115 : !> \param logger Output logger.
116 : !> \param qs_env Quickstep environment; the SCF orbitals and density are not modified.
117 : !> \param output_unit Main output unit, nonpositive on non-writing ranks.
118 : ! **************************************************************************************************
119 12381 : SUBROUTINE qs_scf_post_fod(input, logger, qs_env, output_unit)
120 : TYPE(section_vals_type), POINTER :: input
121 : TYPE(cp_logger_type), POINTER :: logger
122 : TYPE(qs_environment_type), POINTER :: qs_env
123 : INTEGER, INTENT(IN) :: output_unit
124 :
125 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_fod'
126 :
127 : INTEGER :: handle, iatom, ikind, ispin, unit_nr
128 : LOGICAL :: write_cube
129 : REAL(dp) :: grid_integral, nfod, spin_integral
130 12381 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: weights, zcharge
131 12381 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
132 : TYPE(dbcsr_type), POINTER :: matrix_fod
133 : TYPE(mo_set_type) :: mo_fod
134 12381 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
135 : TYPE(particle_list_type), POINTER :: particles
136 : TYPE(pw_c1d_gs_type) :: tmp_g
137 : TYPE(pw_env_type), POINTER :: pw_env
138 : TYPE(pw_pool_type), POINTER :: pool
139 : TYPE(pw_r3d_rs_type) :: fod_r, tmp_r
140 12381 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
141 : TYPE(qs_ks_env_type), POINTER :: ks_env
142 : TYPE(qs_rho_type), POINTER :: rho
143 : TYPE(qs_subsys_type), POINTER :: subsys
144 : TYPE(scf_control_type), POINTER :: scf_control
145 : TYPE(section_vals_type), POINTER :: fod_section
146 :
147 24762 : fod_section => section_vals_get_subs_vals(input, "DFT%PRINT%FOD")
148 12381 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, fod_section, ""), cp_p_file)) RETURN
149 4 : CALL timeset(routineN, handle)
150 4 : CALL qs_fod_validate(input, logger, qs_env)
151 4 : CALL get_qs_env(qs_env, mos=mos, scf_control=scf_control)
152 4 : CPASSERT(ASSOCIATED(mos))
153 4 : CALL cite_reference(Grimme2015FOD)
154 4 : CALL section_vals_val_get(fod_section, "CUBE", l_val=write_cube)
155 :
156 4 : NULLIFY (rho, rho_ao, ks_env, pw_env, pool, subsys, particles, matrix_fod)
157 4 : IF (write_cube) THEN
158 2 : CALL get_qs_env(qs_env, rho=rho, ks_env=ks_env, pw_env=pw_env, subsys=subsys, qs_kind_set=qs_kind_set)
159 2 : CALL qs_rho_get(rho, rho_ao=rho_ao)
160 2 : CALL qs_subsys_get(subsys, particles=particles)
161 6 : ALLOCATE (zcharge(particles%n_els))
162 6 : DO iatom = 1, particles%n_els
163 4 : ikind = particles%els(iatom)%atomic_kind%kind_number
164 6 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zcharge(iatom))
165 : END DO
166 2 : CALL pw_env_get(pw_env, auxbas_pw_pool=pool)
167 2 : CALL pool%create_pw(fod_r)
168 2 : CALL pool%create_pw(tmp_r)
169 2 : CALL pool%create_pw(tmp_g)
170 95267 : fod_r%array = 0.0_dp
171 : END IF
172 :
173 4 : nfod = 0.0_dp
174 4 : grid_integral = 0.0_dp
175 10 : DO ispin = 1, SIZE(mos)
176 6 : CPASSERT(mos(ispin)%nmo > 0)
177 6 : IF (mos(ispin)%occupation_numbers(mos(ispin)%nmo) > &
178 : mos(ispin)%maxocc*scf_control%smear%eps_fermi_dirac) THEN
179 0 : CPWARN("FOD: the highest available orbital is occupied; increase ADDED_MOS and check convergence.")
180 : END IF
181 18 : ALLOCATE (weights(mos(ispin)%nmo))
182 : weights(:) = fod_weights(mos(ispin)%occupation_numbers, mos(ispin)%eigenvalues, &
183 6 : mos(ispin)%mu, mos(ispin)%maxocc)
184 54 : nfod = nfod + SUM(weights)
185 6 : IF (write_cube) THEN
186 2 : CALL duplicate_mo_set(mo_fod, mos(ispin))
187 18 : mo_fod%occupation_numbers = weights
188 2 : mo_fod%uniform_occupation = .FALSE.
189 : ! Include every supplied orbital, including the fractional virtual tail.
190 2 : mo_fod%homo = mo_fod%nmo
191 2 : ALLOCATE (matrix_fod)
192 2 : CALL dbcsr_copy(matrix_fod, rho_ao(ispin)%matrix)
193 2 : CALL calculate_density_matrix(mo_fod, matrix_fod)
194 : CALL calculate_rho_elec(matrix_p=matrix_fod, rho=tmp_r, rho_gspace=tmp_g, &
195 2 : total_rho=spin_integral, ks_env=ks_env)
196 95267 : fod_r%array = fod_r%array + tmp_r%array
197 : ! calculate_rho_elec returns electronic charge, not electron count.
198 2 : grid_integral = grid_integral - spin_integral
199 2 : CALL dbcsr_release(matrix_fod)
200 2 : DEALLOCATE (matrix_fod)
201 2 : CALL deallocate_mo_set(mo_fod)
202 : END IF
203 10 : DEALLOCATE (weights)
204 : END DO
205 :
206 4 : IF (output_unit > 0) THEN
207 2 : WRITE (output_unit, '(T2,A,T61,F20.10)') "FOD| N_FOD (orbital sum)", nfod
208 2 : IF (write_cube) WRITE (output_unit, '(T2,A,T61,F20.10)') "FOD| Grid integral", grid_integral
209 : END IF
210 4 : IF (write_cube) THEN
211 : unit_nr = cp_print_key_unit_nr(logger, fod_section, "", extension=".cube", &
212 2 : middle_name="FOD", file_position="REWIND", log_filename=.FALSE.)
213 : ! Collective call: non-writing ranks must participate even with unit_nr=-1.
214 : CALL cp_pw_to_cube(fod_r, unit_nr, "FRACTIONAL OCCUPATION DENSITY [e/bohr^3]", &
215 2 : particles=particles, zeff=zcharge, stride=section_get_ivals(fod_section, "STRIDE"))
216 2 : CALL cp_print_key_finished_output(unit_nr, logger, fod_section, "")
217 2 : CALL pool%give_back_pw(fod_r)
218 2 : CALL pool%give_back_pw(tmp_r)
219 2 : CALL pool%give_back_pw(tmp_g)
220 2 : DEALLOCATE (zcharge)
221 : END IF
222 4 : CALL timestop(handle)
223 12385 : END SUBROUTINE qs_scf_post_fod
224 : END MODULE qs_fod
|
