Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief
10 : !>
11 : !>
12 : !> \par History
13 : !> refactoring 03-2011 [MI]
14 : !> \author MI
15 : ! **************************************************************************************************
16 : MODULE xc_adiabatic_utils
17 :
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
20 : USE hfx_communication, ONLY: scale_and_add_fock_to_ks_matrix
21 : USE hfx_derivatives, ONLY: derivatives_four_center
22 : USE hfx_types, ONLY: hfx_type
23 : USE input_constants, ONLY: do_adiabatic_hybrid_mcy3,&
24 : do_adiabatic_model_pade
25 : USE input_section_types, ONLY: section_vals_get,&
26 : section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: dp
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE pw_types, ONLY: pw_r3d_rs_type
32 : USE qs_energy_types, ONLY: qs_energy_type
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_ks_types, ONLY: qs_ks_env_type
36 : USE qs_rho_types, ONLY: qs_rho_get,&
37 : qs_rho_type
38 : USE qs_vxc, ONLY: qs_vxc_create
39 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
40 : USE xc_adiabatic_methods, ONLY: rescale_MCY3_pade
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : ! *** Public subroutines ***
48 : PUBLIC :: rescale_xc_potential
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_adiabatic_utils'
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief
56 : !>
57 : !> \param qs_env ...
58 : !> \param ks_matrix ...
59 : !> \param rho ...
60 : !> \param energy ...
61 : !> \param v_rspace_new ...
62 : !> \param v_tau_rspace ...
63 : !> \param hf_energy ...
64 : !> \param just_energy ...
65 : !> \param calculate_forces ...
66 : !> \param use_virial ...
67 : ! **************************************************************************************************
68 504 : SUBROUTINE rescale_xc_potential(qs_env, ks_matrix, rho, energy, v_rspace_new, v_tau_rspace, &
69 56 : hf_energy, just_energy, calculate_forces, use_virial)
70 :
71 : TYPE(qs_environment_type), POINTER :: qs_env
72 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
73 : TYPE(qs_rho_type), POINTER :: rho
74 : TYPE(qs_energy_type), POINTER :: energy
75 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
76 : REAL(dp), DIMENSION(:) :: hf_energy
77 : LOGICAL, INTENT(in) :: just_energy, calculate_forces, use_virial
78 :
79 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rescale_xc_potential'
80 :
81 : INTEGER :: adiabatic_functional, adiabatic_model, &
82 : handle, n_rep_hf, nimages
83 : LOGICAL :: do_adiabatic_rescaling, do_hfx, gapw, &
84 : gapw_xc
85 : REAL(dp) :: adiabatic_lambda, adiabatic_omega, &
86 : scale_dDFA, scale_ddW0, scale_dEx1, &
87 : scale_dEx2, total_energy_xc
88 56 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_resp
89 56 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
90 : TYPE(dft_control_type), POINTER :: dft_control
91 56 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
92 : TYPE(mp_para_env_type), POINTER :: para_env
93 : TYPE(qs_ks_env_type), POINTER :: ks_env
94 : TYPE(qs_rho_type), POINTER :: rho_xc
95 : TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
96 : hfx_sections, input, xc_section
97 :
98 56 : CALL timeset(routineN, handle)
99 56 : NULLIFY (para_env, dft_control, adiabatic_rescaling_section, hfx_sections, &
100 56 : input, xc_section, rho_xc, ks_env, rho_ao, rho_ao_resp, x_data)
101 :
102 : CALL get_qs_env(qs_env, &
103 : dft_control=dft_control, &
104 : para_env=para_env, &
105 : input=input, &
106 : rho_xc=rho_xc, &
107 : ks_env=ks_env, &
108 56 : x_data=x_data)
109 :
110 56 : IF (x_data(1, 1)%do_hfx_ri) CPABORT("RI-HFX not compatible with this kinf of functionals")
111 56 : nimages = dft_control%nimages
112 56 : CPASSERT(nimages == 1)
113 :
114 56 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
115 :
116 56 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
117 56 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
118 56 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
119 56 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
120 56 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
121 :
122 56 : gapw = dft_control%qs_control%gapw
123 56 : gapw_xc = dft_control%qs_control%gapw_xc
124 :
125 : CALL section_vals_val_get(adiabatic_rescaling_section, "FUNCTIONAL_TYPE", &
126 56 : i_val=adiabatic_functional)
127 : CALL section_vals_val_get(adiabatic_rescaling_section, "FUNCTIONAL_MODEL", &
128 56 : i_val=adiabatic_model)
129 : CALL section_vals_val_get(adiabatic_rescaling_section, "LAMBDA", &
130 56 : r_val=adiabatic_lambda)
131 : CALL section_vals_val_get(adiabatic_rescaling_section, "OMEGA", &
132 56 : r_val=adiabatic_omega)
133 56 : SELECT CASE (adiabatic_functional)
134 : CASE (do_adiabatic_hybrid_mcy3)
135 56 : SELECT CASE (adiabatic_model)
136 : CASE (do_adiabatic_model_pade)
137 56 : IF (n_rep_hf /= 2) &
138 : CALL cp_abort(__LOCATION__, &
139 : "For this kind of adiababatic hybrid functional 2 HF sections have to be provided. "// &
140 0 : "Please check your input file.")
141 : CALL rescale_MCY3_pade(qs_env, hf_energy, energy, adiabatic_lambda, &
142 : adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, &
143 56 : scale_dEx2, total_energy_xc)
144 :
145 : !! Scale and add Fock matrix to KS matrix
146 56 : IF (do_hfx) THEN
147 : CALL scale_and_add_fock_to_ks_matrix(para_env, qs_env, ks_matrix, 1, &
148 56 : scale_dEx1)
149 : CALL scale_and_add_fock_to_ks_matrix(para_env, qs_env, ks_matrix, 2, &
150 56 : scale_dEx2)
151 : END IF
152 56 : IF (calculate_forces) THEN
153 0 : CPASSERT(.NOT. use_virial)
154 : !! we also have to scale the forces!!!!
155 : CALL derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_sections, &
156 : para_env, 1, use_virial, &
157 0 : adiabatic_rescale_factor=scale_dEx1)
158 : CALL derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_sections, &
159 : para_env, 2, use_virial, &
160 0 : adiabatic_rescale_factor=scale_dEx2)
161 : END IF
162 :
163 : ! Calculate vxc and rescale it
164 56 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
165 56 : IF (gapw_xc) THEN
166 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, &
167 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
168 0 : just_energy=just_energy, adiabatic_rescale_factor=scale_dDFA)
169 : ELSE
170 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
171 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
172 56 : just_energy=just_energy, adiabatic_rescale_factor=scale_dDFA)
173 : END IF
174 :
175 : ! Calculate vxc and rescale it
176 56 : IF (gapw .OR. gapw_xc) THEN
177 56 : CALL calculate_vxc_atom(qs_env, just_energy, energy%exc1, adiabatic_rescale_factor=scale_dDFA)
178 : END IF
179 : !! Hack for the total energy expression
180 56 : energy%ex = 0.0_dp
181 56 : energy%exc1 = 0.0_dp
182 112 : energy%exc = total_energy_xc
183 :
184 : END SELECT
185 : END SELECT
186 56 : CALL timestop(handle)
187 :
188 56 : END SUBROUTINE rescale_xc_potential
189 :
190 : END MODULE xc_adiabatic_utils
191 :
|
