Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Test routines for HFX caclulations using PW
10 : !>
11 : !>
12 : !> \par History
13 : !> refactoring 03-2011 [MI]
14 : !> Made GAPW compatible 12.2019 (A. Bussy)
15 : !> Refactored from hfx_admm_utils (JGH)
16 : !> \author MI
17 : ! **************************************************************************************************
18 : MODULE hfx_pw_methods
19 : USE atomic_kind_types, ONLY: atomic_kind_type
20 : USE cell_types, ONLY: cell_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_type
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
24 : USE cp_fm_types, ONLY: cp_fm_get_info,&
25 : cp_fm_type
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_type
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_unit_nr
30 : USE input_constants, ONLY: do_potential_coulomb,&
31 : do_potential_short,&
32 : do_potential_truncated
33 : USE input_section_types, ONLY: section_vals_get,&
34 : section_vals_get_subs_vals,&
35 : section_vals_type,&
36 : section_vals_val_get
37 : USE kinds, ONLY: dp
38 : USE mathconstants, ONLY: fourpi
39 : USE particle_types, ONLY: particle_type
40 : USE pw_env_types, ONLY: pw_env_type
41 : USE pw_grid_types, ONLY: pw_grid_type
42 : USE pw_methods, ONLY: pw_copy,&
43 : pw_multiply,&
44 : pw_transfer,&
45 : pw_zero
46 : USE pw_poisson_methods, ONLY: pw_poisson_solve
47 : USE pw_poisson_types, ONLY: pw_poisson_type
48 : USE pw_pool_types, ONLY: pw_pool_type
49 : USE pw_types, ONLY: pw_c1d_gs_type,&
50 : pw_r3d_rs_type
51 : USE qs_collocate_density, ONLY: calculate_wavefunction
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_kind_types, ONLY: qs_kind_type
55 : USE qs_mo_types, ONLY: get_mo_set,&
56 : mo_set_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : ! *** Public subroutines ***
64 : PUBLIC :: pw_hfx
65 :
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_pw_methods'
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief computes the Hartree-Fock energy brute force in a pw basis
72 : !> \param qs_env ...
73 : !> \param ehfx ...
74 : !> \param hfx_section ...
75 : !> \param poisson_env ...
76 : !> \param auxbas_pw_pool ...
77 : !> \param irep ...
78 : !> \par History
79 : !> 12.2007 created [Joost VandeVondele]
80 : !> \note
81 : !> only computes the HFX energy, no derivatives as yet
82 : ! **************************************************************************************************
83 25282 : SUBROUTINE pw_hfx(qs_env, ehfx, hfx_section, poisson_env, auxbas_pw_pool, irep)
84 : TYPE(qs_environment_type), POINTER :: qs_env
85 : REAL(KIND=dp), INTENT(IN) :: ehfx
86 : TYPE(section_vals_type), POINTER :: hfx_section
87 : TYPE(pw_poisson_type), POINTER :: poisson_env
88 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
89 : INTEGER :: irep
90 :
91 : CHARACTER(*), PARAMETER :: routineN = 'pw_hfx'
92 :
93 : INTEGER :: blocksize, handle, ig, iloc, iorb, &
94 : iorb_block, ispin, iw, jloc, jorb, &
95 : jorb_block, norb, potential_type
96 : LOGICAL :: do_kpoints, do_pw_hfx, explicit
97 : REAL(KIND=dp) :: exchange_energy, fraction, g2, g3d, gg, &
98 : omega, pair_energy, rcut, scaling
99 25282 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
100 : TYPE(cell_type), POINTER :: cell
101 : TYPE(cp_fm_type), POINTER :: mo_coeff
102 : TYPE(cp_logger_type), POINTER :: logger
103 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
104 : TYPE(dft_control_type), POINTER :: dft_control
105 25282 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
106 25282 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 : TYPE(pw_c1d_gs_type) :: greenfn, pot_g, rho_g
108 : TYPE(pw_env_type), POINTER :: pw_env
109 : TYPE(pw_grid_type), POINTER :: grid
110 : TYPE(pw_r3d_rs_type) :: rho_r
111 25282 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_i, rho_j
112 25282 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
113 : TYPE(section_vals_type), POINTER :: ip_section
114 :
115 25282 : CALL timeset(routineN, handle)
116 25282 : logger => cp_get_default_logger()
117 :
118 25282 : CALL section_vals_val_get(hfx_section, "PW_HFX", l_val=do_pw_hfx, i_rep_section=irep)
119 :
120 25282 : IF (do_pw_hfx) THEN
121 32 : CALL section_vals_val_get(hfx_section, "FRACTION", r_val=fraction)
122 32 : CALL section_vals_val_get(hfx_section, "PW_HFX_BLOCKSIZE", i_val=blocksize)
123 :
124 : CALL get_qs_env(qs_env, mos=mo_array, pw_env=pw_env, dft_control=dft_control, &
125 : cell=cell, particle_set=particle_set, do_kpoints=do_kpoints, &
126 32 : atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
127 32 : IF (do_kpoints) CPABORT("PW HFX not implemented with K-points")
128 :
129 : ! limit the blocksize by the number of orbitals
130 32 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
131 32 : CALL cp_fm_get_info(mo_coeff, ncol_global=norb)
132 32 : blocksize = MAX(1, MIN(blocksize, norb))
133 :
134 32 : CALL auxbas_pw_pool%create_pw(rho_r)
135 32 : CALL auxbas_pw_pool%create_pw(rho_g)
136 32 : CALL auxbas_pw_pool%create_pw(pot_g)
137 :
138 218 : ALLOCATE (rho_i(blocksize))
139 186 : ALLOCATE (rho_j(blocksize))
140 :
141 32 : CALL auxbas_pw_pool%create_pw(greenfn)
142 32 : ip_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL")
143 32 : CALL section_vals_get(ip_section, explicit=explicit)
144 32 : potential_type = do_potential_coulomb
145 32 : IF (explicit) THEN
146 22 : CALL section_vals_val_get(ip_section, "POTENTIAL_TYPE", i_val=potential_type)
147 : END IF
148 32 : IF (potential_type == do_potential_coulomb) THEN
149 12 : CALL pw_copy(poisson_env%green_fft%influence_fn, greenfn)
150 20 : ELSEIF (potential_type == do_potential_truncated) THEN
151 10 : CALL section_vals_val_get(ip_section, "CUTOFF_RADIUS", r_val=rcut)
152 10 : grid => poisson_env%green_fft%influence_fn%pw_grid
153 1866245 : DO ig = grid%first_gne0, grid%ngpts_cut_local
154 1866235 : g2 = grid%gsq(ig)
155 1866235 : gg = SQRT(g2)
156 1866235 : g3d = fourpi/g2
157 1866245 : greenfn%array(ig) = g3d*(1.0_dp - COS(rcut*gg))
158 : END DO
159 10 : IF (grid%have_g0) &
160 5 : greenfn%array(1) = 0.5_dp*fourpi*rcut*rcut
161 10 : ELSEIF (potential_type == do_potential_short) THEN
162 10 : CALL section_vals_val_get(ip_section, "OMEGA", r_val=omega)
163 10 : IF (omega > 0.0_dp) omega = 0.25_dp/(omega*omega)
164 10 : grid => poisson_env%green_fft%influence_fn%pw_grid
165 1866245 : DO ig = grid%first_gne0, grid%ngpts_cut_local
166 1866235 : g2 = grid%gsq(ig)
167 1866235 : gg = -omega*g2
168 1866235 : g3d = fourpi/g2
169 1866245 : greenfn%array(ig) = g3d*(1.0_dp - EXP(gg))
170 : END DO
171 10 : IF (grid%have_g0) greenfn%array(1) = 0.0_dp
172 : ELSE
173 0 : CPWARN("PW_SCF: Potential type not supported, calculation uses Coulomb potential.")
174 : END IF
175 :
176 154 : DO iorb_block = 1, blocksize
177 122 : CALL rho_i(iorb_block)%create(rho_r%pw_grid)
178 154 : CALL rho_j(iorb_block)%create(rho_r%pw_grid)
179 : END DO
180 :
181 32 : exchange_energy = 0.0_dp
182 :
183 64 : DO ispin = 1, SIZE(mo_array)
184 32 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b)
185 :
186 32 : IF (mo_array(ispin)%use_mo_coeff_b) THEN !fm->dbcsr
187 32 : CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff) !fm->dbcsr
188 : END IF !fm->dbcsr
189 :
190 32 : CALL cp_fm_get_info(mo_coeff, ncol_global=norb)
191 :
192 96 : DO iorb_block = 1, norb, blocksize
193 :
194 154 : DO iorb = iorb_block, MIN(iorb_block + blocksize - 1, norb)
195 :
196 122 : iloc = iorb - iorb_block + 1
197 : CALL calculate_wavefunction(mo_coeff, iorb, rho_i(iloc), rho_g, &
198 : atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
199 154 : pw_env)
200 :
201 : END DO
202 :
203 96 : DO jorb_block = iorb_block, norb, blocksize
204 :
205 154 : DO jorb = jorb_block, MIN(jorb_block + blocksize - 1, norb)
206 :
207 122 : jloc = jorb - jorb_block + 1
208 : CALL calculate_wavefunction(mo_coeff, jorb, rho_j(jloc), rho_g, &
209 : atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
210 154 : pw_env)
211 :
212 : END DO
213 :
214 186 : DO iorb = iorb_block, MIN(iorb_block + blocksize - 1, norb)
215 122 : iloc = iorb - iorb_block + 1
216 636 : DO jorb = jorb_block, MIN(jorb_block + blocksize - 1, norb)
217 482 : jloc = jorb - jorb_block + 1
218 482 : IF (jorb < iorb) CYCLE
219 :
220 : ! compute the pair density
221 302 : CALL pw_zero(rho_r)
222 302 : CALL pw_multiply(rho_r, rho_i(iloc), rho_j(jloc))
223 :
224 : ! go the g-space and compute hartree energy
225 302 : CALL pw_transfer(rho_r, rho_g)
226 : CALL pw_poisson_solve(poisson_env, rho_g, pair_energy, pot_g, &
227 302 : greenfn=greenfn)
228 :
229 : ! sum up to the full energy
230 302 : scaling = fraction
231 302 : IF (SIZE(mo_array) == 1) scaling = scaling*2.0_dp
232 302 : IF (iorb /= jorb) scaling = scaling*2.0_dp
233 :
234 906 : exchange_energy = exchange_energy - scaling*pair_energy
235 :
236 : END DO
237 : END DO
238 :
239 : END DO
240 : END DO
241 : END DO
242 :
243 154 : DO iorb_block = 1, blocksize
244 122 : CALL rho_i(iorb_block)%release()
245 154 : CALL rho_j(iorb_block)%release()
246 : END DO
247 :
248 32 : CALL auxbas_pw_pool%give_back_pw(rho_r)
249 32 : CALL auxbas_pw_pool%give_back_pw(rho_g)
250 32 : CALL auxbas_pw_pool%give_back_pw(pot_g)
251 32 : CALL auxbas_pw_pool%give_back_pw(greenfn)
252 :
253 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
254 32 : extension=".scfLog")
255 :
256 32 : IF (iw > 0) THEN
257 : WRITE (UNIT=iw, FMT="((T3,A,T61,F20.10))") &
258 16 : "HF_PW_HFX| PW exchange energy:", exchange_energy
259 : WRITE (UNIT=iw, FMT="((T3,A,T61,F20.10),/)") &
260 16 : "HF_PW_HFX| Gaussian exchange energy:", ehfx
261 : END IF
262 :
263 128 : CALL cp_print_key_finished_output(iw, logger, hfx_section, "HF_INFO")
264 :
265 : END IF
266 :
267 25282 : CALL timestop(handle)
268 :
269 50564 : END SUBROUTINE pw_hfx
270 :
271 : END MODULE hfx_pw_methods
|
