1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Contains some functions used in the context of adiabatic hybrid functionals
8!> \par History
9!>      01.2008 created [Manuel Guidon]
10!> \author Manuel Guidon
11! **************************************************************************************************
12MODULE xc_adiabatic_methods
13
14   USE input_constants,                 ONLY: do_potential_coulomb
15   USE kinds,                           ONLY: dp
16   USE mathconstants,                   ONLY: oorootpi
17   USE qs_energy_types,                 ONLY: qs_energy_type
18   USE qs_environment_types,            ONLY: get_qs_env,&
19                                              qs_environment_type
20   USE qs_mo_types,                     ONLY: get_mo_set,&
21                                              mo_set_p_type
22#include "./base/base_uses.f90"
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC :: rescale_MCY3_pade
28
29CONTAINS
30
31! **************************************************************************************************
32!> \brief - Calculates rescaling factors for XC potentials and energy expression
33!> \param qs_env ...
34!> \param hf_energy Array of size 2 containing the two HF energies (Ex^{HF} and Ex^{HF,LR}
35!> \param energy QS energy type
36!> \param adiabatic_lambda , adiabatic_omega: Parameters for adiabatic connection
37!> \param adiabatic_omega ...
38!> \param scale_dEx1 scaling coefficient for xc-potentials to be calculated
39!> \param scale_ddW0 scaling coefficient for xc-potentials to be calculated
40!> \param scale_dDFA scaling coefficient for xc-potentials to be calculated
41!> \param scale_dEx2 scaling coefficient for xc-potentials to be calculated
42!> \param total_energy_xc will contain the full xc energy
43!> \par History
44!>      09.2007 created [Manuel Guidon]
45!> \author Manuel Guidon
46!> \note
47!>      - Model for adiabatic connection:
48!>
49!>         W_lambda = a + b*lambda/(1+c*lambda)
50!>         Exc = a + b*(c-log(1+c)/c^2)
51!>         a = Ex^{HF}
52!>         b = -c1*2*omega/sqrt(PI)*nelectron
53!>         c = -1/lambda - b/(Ex^{HF}-W_lambda^{BLYP} + c2*W_lambda^{B88,LR}-c3*W_lambda^{HF,LR}
54! **************************************************************************************************
55   SUBROUTINE rescale_MCY3_pade(qs_env, hf_energy, energy, adiabatic_lambda, &
56                                adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, &
57                                scale_dEx2, total_energy_xc)
58
59      TYPE(qs_environment_type), POINTER                 :: qs_env
60      REAL(dp), INTENT(INOUT)                            :: hf_energy(*)
61      TYPE(qs_energy_type), POINTER                      :: energy
62      REAL(dp), INTENT(IN)                               :: adiabatic_lambda, adiabatic_omega
63      REAL(dp), INTENT(INOUT)                            :: scale_dEx1, scale_ddW0, scale_dDFA, &
64                                                            scale_dEx2, total_energy_xc
65
66      INTEGER                                            :: nelec_a, nelec_b, nelectron, nspins
67      LOGICAL                                            :: do_swap_hf
68      REAL(dp) :: a, b, c, c1, da_dDFA, da_ddW0, da_dEx1, da_dEx2, db_dDFA, db_ddW0, db_dEx1, &
69         db_dEx2, dc_dDFA, dc_ddW0, dc_dEx1, dc_dEx2, dExc_da, dExc_db, dExc_dc, dfa_energy, &
70         swap_value
71      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
72
73      do_swap_hf = .FALSE.
74      !! Assume the first HF section is the Coulomb one
75      IF (qs_env%x_data(1, 1)%potential_parameter%potential_type /= do_potential_coulomb) do_swap_hf = .TRUE.
76
77      IF (do_swap_hf) THEN
78         swap_value = hf_energy(1)
79         hf_energy(1) = hf_energy(2)
80         hf_energy(2) = swap_value
81      END IF
82
83      c1 = 0.23163_dp
84
85      CALL get_qs_env(qs_env=qs_env, mos=mos)
86      CALL get_mo_set(mo_set=mos(1)%mo_set, nelectron=nelec_a)
87      nspins = SIZE(mos)
88      IF (nspins == 2) THEN
89         CALL get_mo_set(mo_set=mos(2)%mo_set, nelectron=nelec_b)
90      ELSE
91         nelec_b = 0
92      END IF
93      nelectron = nelec_a + nelec_b
94      dfa_energy = energy%exc + energy%exc1
95      a = hf_energy(1)
96      b = -c1*2.0_dp*adiabatic_omega*oorootpi*nelectron !-0.23163_dp*2.0_dp*0.2_dp*oorootpi*nelectron
97      c = -1.0_dp/adiabatic_lambda - b/(hf_energy(1) - dfa_energy - hf_energy(2))
98
99      dExc_da = 1.0_dp
100      dExc_db = 1.0_dp/c - (LOG(ABS(1.0_dp + c))/(c*c))
101      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)
102
103      da_dEx1 = 1.0_dp
104      da_ddW0 = 0.0_dp
105      da_dDFA = 0.0_dp
106      da_dEx2 = 0.0_dp
107
108      db_dEx1 = 0.0_dp
109      db_ddW0 = 1.0_dp
110      db_dDFA = 0.0_dp
111      db_dEx2 = 0.0_dp
112
113      dc_dEx1 = b/(hf_energy(1) - dfa_energy - hf_energy(2))**2
114      dc_ddW0 = -1.0_dp/(hf_energy(1) - dfa_energy - hf_energy(2))
115      dc_dDFA = -dc_dEx1
116      dc_dEx2 = -dc_dEx1
117
118      scale_dEx1 = dExc_da*da_dEx1 + dExc_db*db_dEx1 + dExc_dc*dc_dEx1
119      scale_ddW0 = dExc_da*da_ddW0 + dExc_db*db_ddW0 + dExc_dc*dc_ddW0
120      scale_dDFA = dExc_da*da_dDFA + dExc_db*db_dDFA + dExc_dc*dc_dDFA
121      scale_dEx2 = dExc_da*da_dEx2 + dExc_db*db_dEx2 + dExc_dc*dc_dEx2
122
123      total_energy_xc = a + b/(c*c)*(c - LOG(ABS(1.0_dp + c)))
124      IF (do_swap_hf) THEN
125         swap_value = scale_dEx1
126         scale_dEx1 = scale_dEx2
127         scale_dEx2 = swap_value
128      END IF
129   END SUBROUTINE rescale_MCY3_pade
130
131END MODULE xc_adiabatic_methods
132