(git:1f285aa)
xc_adiabatic_methods.F
Go to the documentation of this file.
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 Contains some functions used in the context of adiabatic hybrid functionals
10 !> \par History
11 !> 01.2008 created [Manuel Guidon]
12 !> \author Manuel Guidon
13 ! **************************************************************************************************
15 
17  USE kinds, ONLY: dp
18  USE mathconstants, ONLY: oorootpi
19  USE qs_energy_types, ONLY: qs_energy_type
20  USE qs_environment_types, ONLY: get_qs_env,&
21  qs_environment_type
22  USE qs_mo_types, ONLY: get_mo_set,&
23  mo_set_type
24 #include "./base/base_uses.f90"
25 
26  IMPLICIT NONE
27  PRIVATE
28 
29  PUBLIC :: rescale_mcy3_pade
30 
31 CONTAINS
32 
33 ! **************************************************************************************************
34 !> \brief - Calculates rescaling factors for XC potentials and energy expression
35 !> \param qs_env ...
36 !> \param hf_energy Array of size 2 containing the two HF energies (Ex^{HF} and Ex^{HF,LR}
37 !> \param energy QS energy type
38 !> \param adiabatic_lambda , adiabatic_omega: Parameters for adiabatic connection
39 !> \param adiabatic_omega ...
40 !> \param scale_dEx1 scaling coefficient for xc-potentials to be calculated
41 !> \param scale_ddW0 scaling coefficient for xc-potentials to be calculated
42 !> \param scale_dDFA scaling coefficient for xc-potentials to be calculated
43 !> \param scale_dEx2 scaling coefficient for xc-potentials to be calculated
44 !> \param total_energy_xc will contain the full xc energy
45 !> \par History
46 !> 09.2007 created [Manuel Guidon]
47 !> \author Manuel Guidon
48 !> \note
49 !> - Model for adiabatic connection:
50 !>
51 !> W_lambda = a + b*lambda/(1+c*lambda)
52 !> Exc = a + b*(c-log(1+c)/c^2)
53 !> a = Ex^{HF}
54 !> b = -c1*2*omega/sqrt(PI)*nelectron
55 !> c = -1/lambda - b/(Ex^{HF}-W_lambda^{BLYP} + c2*W_lambda^{B88,LR}-c3*W_lambda^{HF,LR}
56 ! **************************************************************************************************
57  SUBROUTINE rescale_mcy3_pade(qs_env, hf_energy, energy, adiabatic_lambda, &
58  adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, &
59  scale_dEx2, total_energy_xc)
60 
61  TYPE(qs_environment_type), POINTER :: qs_env
62  REAL(dp), INTENT(INOUT) :: hf_energy(*)
63  TYPE(qs_energy_type), POINTER :: energy
64  REAL(dp), INTENT(IN) :: adiabatic_lambda, adiabatic_omega
65  REAL(dp), INTENT(INOUT) :: scale_dex1, scale_ddw0, scale_ddfa, &
66  scale_dex2, total_energy_xc
67 
68  INTEGER :: nelec_a, nelec_b, nelectron, nspins
69  LOGICAL :: do_swap_hf
70  REAL(dp) :: a, b, c, c1, da_ddfa, da_ddw0, da_dex1, da_dex2, db_ddfa, db_ddw0, db_dex1, &
71  db_dex2, dc_ddfa, dc_ddw0, dc_dex1, dc_dex2, dexc_da, dexc_db, dexc_dc, dfa_energy, &
72  swap_value
73  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
74 
75  do_swap_hf = .false.
76  !! Assume the first HF section is the Coulomb one
77  IF (qs_env%x_data(1, 1)%potential_parameter%potential_type /= do_potential_coulomb) do_swap_hf = .true.
78 
79  IF (do_swap_hf) THEN
80  swap_value = hf_energy(1)
81  hf_energy(1) = hf_energy(2)
82  hf_energy(2) = swap_value
83  END IF
84 
85  c1 = 0.23163_dp
86 
87  CALL get_qs_env(qs_env=qs_env, mos=mos)
88  CALL get_mo_set(mo_set=mos(1), nelectron=nelec_a)
89  nspins = SIZE(mos)
90  IF (nspins == 2) THEN
91  CALL get_mo_set(mo_set=mos(2), nelectron=nelec_b)
92  ELSE
93  nelec_b = 0
94  END IF
95  nelectron = nelec_a + nelec_b
96  dfa_energy = energy%exc + energy%exc1
97  a = hf_energy(1)
98  b = -c1*2.0_dp*adiabatic_omega*oorootpi*nelectron !-0.23163_dp*2.0_dp*0.2_dp*oorootpi*nelectron
99  c = -1.0_dp/adiabatic_lambda - b/(hf_energy(1) - dfa_energy - hf_energy(2))
100 
101  dexc_da = 1.0_dp
102  dexc_db = 1.0_dp/c - (log(abs(1.0_dp + c))/(c*c))
103  dexc_dc = -b/(c*c*c*(1.0_dp + c))*(2.0_dp*c + c*c - 2.0_dp*log(abs(1.0_dp + c)) - 2.0_dp*log(abs(1.0_dp + c))*c)
104 
105  da_dex1 = 1.0_dp
106  da_ddw0 = 0.0_dp
107  da_ddfa = 0.0_dp
108  da_dex2 = 0.0_dp
109 
110  db_dex1 = 0.0_dp
111  db_ddw0 = 1.0_dp
112  db_ddfa = 0.0_dp
113  db_dex2 = 0.0_dp
114 
115  dc_dex1 = b/(hf_energy(1) - dfa_energy - hf_energy(2))**2
116  dc_ddw0 = -1.0_dp/(hf_energy(1) - dfa_energy - hf_energy(2))
117  dc_ddfa = -dc_dex1
118  dc_dex2 = -dc_dex1
119 
120  scale_dex1 = dexc_da*da_dex1 + dexc_db*db_dex1 + dexc_dc*dc_dex1
121  scale_ddw0 = dexc_da*da_ddw0 + dexc_db*db_ddw0 + dexc_dc*dc_ddw0
122  scale_ddfa = dexc_da*da_ddfa + dexc_db*db_ddfa + dexc_dc*dc_ddfa
123  scale_dex2 = dexc_da*da_dex2 + dexc_db*db_dex2 + dexc_dc*dc_dex2
124 
125  total_energy_xc = a + b/(c*c)*(c - log(abs(1.0_dp + c)))
126  IF (do_swap_hf) THEN
127  swap_value = scale_dex1
128  scale_dex1 = scale_dex2
129  scale_dex2 = swap_value
130  END IF
131  END SUBROUTINE rescale_mcy3_pade
132 
133 END MODULE xc_adiabatic_methods
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_potential_coulomb
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public oorootpi
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
Contains some functions used in the context of adiabatic hybrid functionals.
subroutine, public rescale_mcy3_pade(qs_env, hf_energy, energy, adiabatic_lambda, adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, scale_dEx2, total_energy_xc)
Calculates rescaling factors for XC potentials and energy expression