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