1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Energy correction environment setup and handling
8!> \par History
9!>       2019.09 created
10!> \author JGH
11! **************************************************************************************************
12MODULE ec_environment
13   USE atomic_kind_types,               ONLY: atomic_kind_type
14   USE basis_set_container_types,       ONLY: add_basis_set_to_container,&
15                                              remove_basis_from_container
16   USE basis_set_types,                 ONLY: copy_gto_basis_set,&
17                                              create_primitive_basis_set,&
18                                              gto_basis_set_type
19   USE cp_control_types,                ONLY: dft_control_type
20   USE cp_para_types,                   ONLY: cp_para_env_type
21   USE ec_env_types,                    ONLY: energy_correction_type
22   USE input_constants,                 ONLY: ec_functional_harris,&
23                                              xc_vdw_fun_nonloc,&
24                                              xc_vdw_fun_pairpot
25   USE input_cp2k_check,                ONLY: xc_functionals_expand
26   USE input_section_types,             ONLY: section_vals_get,&
27                                              section_vals_get_subs_vals,&
28                                              section_vals_type,&
29                                              section_vals_val_get
30   USE kinds,                           ONLY: dp
31   USE orbital_pointers,                ONLY: init_orbital_pointers
32   USE qs_dispersion_nonloc,            ONLY: qs_dispersion_nonloc_init
33   USE qs_dispersion_pairpot,           ONLY: qs_dispersion_pairpot_init
34   USE qs_dispersion_types,             ONLY: qs_dispersion_type
35   USE qs_dispersion_utils,             ONLY: qs_dispersion_env_set
36   USE qs_environment_types,            ONLY: get_qs_env,&
37                                              qs_environment_type
38   USE qs_interactions,                 ONLY: init_interaction_radii_orb_basis
39   USE qs_kind_types,                   ONLY: get_qs_kind,&
40                                              get_qs_kind_set,&
41                                              qs_kind_type
42   USE string_utilities,                ONLY: uppercase
43#include "./base/base_uses.f90"
44
45   IMPLICIT NONE
46
47   PRIVATE
48
49   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_environment'
50
51   PUBLIC :: ec_env_create
52
53CONTAINS
54
55! **************************************************************************************************
56!> \brief Allocates and intitializes ec_env
57!> \param qs_env ...
58!> \param ec_env the object to create
59!> \param dft_section ...
60!> \par History
61!>       2019.09 created
62!> \author JGH
63! **************************************************************************************************
64   SUBROUTINE ec_env_create(qs_env, ec_env, dft_section)
65      TYPE(qs_environment_type), POINTER                 :: qs_env
66      TYPE(energy_correction_type), POINTER              :: ec_env
67      TYPE(section_vals_type), POINTER                   :: dft_section
68
69      CHARACTER(len=*), PARAMETER :: routineN = 'ec_env_create', routineP = moduleN//':'//routineN
70
71      CPASSERT(.NOT. ASSOCIATED(ec_env))
72      ALLOCATE (ec_env)
73      CALL init_ec_env(qs_env, ec_env, dft_section)
74
75   END SUBROUTINE ec_env_create
76
77! **************************************************************************************************
78!> \brief Initializes ec_env
79!> \param qs_env ...
80!> \param ec_env ...
81!> \param dft_section ...
82!> \par History
83!>       2019.09 created
84!> \author JGH
85! **************************************************************************************************
86   SUBROUTINE init_ec_env(qs_env, ec_env, dft_section)
87      TYPE(qs_environment_type), POINTER                 :: qs_env
88      TYPE(energy_correction_type), POINTER              :: ec_env
89      TYPE(section_vals_type), OPTIONAL, POINTER         :: dft_section
90
91      CHARACTER(LEN=*), PARAMETER :: routineN = 'init_ec_env', routineP = moduleN//':'//routineN
92
93      INTEGER                                            :: ikind, maxlgto, nkind
94      LOGICAL                                            :: explicit
95      REAL(KIND=dp)                                      :: eps_pgf_orb
96      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
97      TYPE(cp_para_env_type), POINTER                    :: para_env
98      TYPE(dft_control_type), POINTER                    :: dft_control
99      TYPE(gto_basis_set_type), POINTER                  :: basis_set, harris_basis
100      TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
101      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
102      TYPE(qs_kind_type), POINTER                        :: qs_kind
103      TYPE(section_vals_type), POINTER                   :: nl_section, pp_section, section1, &
104                                                            section2, xc_section
105
106      NULLIFY (atomic_kind_set, dispersion_env, para_env)
107      NULLIFY (ec_env%sab_orb, ec_env%sac_ppl, ec_env%sap_ppnl)
108      NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s)
109      NULLIFY (ec_env%matrix_t, ec_env%matrix_p, ec_env%matrix_w)
110      NULLIFY (ec_env%task_list)
111      NULLIFY (ec_env%force)
112      NULLIFY (ec_env%mao_coef)
113      NULLIFY (ec_env%dispersion_env)
114      NULLIFY (ec_env%xc_section)
115      NULLIFY (ec_env%cpmos)
116      NULLIFY (ec_env%matrix_hz)
117      NULLIFY (ec_env%p_env)
118      NULLIFY (ec_env%vh_rspace)
119      NULLIFY (ec_env%vxc_rspace)
120      NULLIFY (ec_env%vtau_rspace)
121      ec_env%should_update = .TRUE.
122      ec_env%mao = .FALSE.
123
124      IF (qs_env%energy_correction) THEN
125         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%ALGORITHM", &
126                                   i_val=ec_env%ks_solver)
127         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%ENERGY_FUNCTIONAL", &
128                                   i_val=ec_env%energy_functional)
129         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%FACTORIZATION", &
130                                   i_val=ec_env%factorization)
131         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%EPS_DEFAULT", &
132                                   r_val=ec_env%eps_default)
133         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%HARRIS_BASIS", &
134                                   c_val=ec_env%basis)
135         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO", &
136                                   l_val=ec_env%mao)
137         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO_MAX_ITER", &
138                                   i_val=ec_env%mao_max_iter)
139         CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO_EPS_GRAD", &
140                                   r_val=ec_env%mao_eps_grad)
141
142         ! set basis
143         CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
144         CALL uppercase(ec_env%basis)
145         SELECT CASE (ec_env%basis)
146         CASE ("ORBITAL")
147            DO ikind = 1, nkind
148               qs_kind => qs_kind_set(ikind)
149               CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
150               IF (ASSOCIATED(basis_set)) THEN
151                  NULLIFY (harris_basis)
152                  CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
153                  IF (ASSOCIATED(harris_basis)) THEN
154                     CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
155                  END IF
156                  NULLIFY (harris_basis)
157                  CALL copy_gto_basis_set(basis_set, harris_basis)
158                  CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
159               END IF
160            END DO
161         CASE ("PRIMITIVE")
162            DO ikind = 1, nkind
163               qs_kind => qs_kind_set(ikind)
164               CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
165               IF (ASSOCIATED(basis_set)) THEN
166                  NULLIFY (harris_basis)
167                  CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
168                  IF (ASSOCIATED(harris_basis)) THEN
169                     CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
170                  END IF
171                  NULLIFY (harris_basis)
172                  CALL create_primitive_basis_set(basis_set, harris_basis)
173                  CALL get_qs_env(qs_env, dft_control=dft_control)
174                  eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
175                  CALL init_interaction_radii_orb_basis(harris_basis, eps_pgf_orb)
176                  harris_basis%kind_radius = basis_set%kind_radius
177                  CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
178               END IF
179            END DO
180         CASE ("HARRIS")
181            DO ikind = 1, nkind
182               qs_kind => qs_kind_set(ikind)
183               NULLIFY (harris_basis)
184               CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
185               IF (.NOT. ASSOCIATED(harris_basis)) THEN
186                  CPWARN("Harris Basis not defined for all types of atoms.")
187               END IF
188            END DO
189         CASE DEFAULT
190            CPABORT("Unknown basis set for energy correction (Harris functional)")
191         END SELECT
192         !
193         CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS")
194         CALL init_orbital_pointers(maxlgto + 1)
195         ! set functional
196         SELECT CASE (ec_env%energy_functional)
197         CASE (ec_functional_harris)
198            ec_env%ec_name = "Harris"
199         CASE DEFAULT
200            CPABORT("unknown energy correction")
201         END SELECT
202         ! select the XC section
203         NULLIFY (xc_section)
204         xc_section => section_vals_get_subs_vals(dft_section, "XC")
205         section1 => section_vals_get_subs_vals(dft_section, "ENERGY_CORRECTION%XC")
206         section2 => section_vals_get_subs_vals(dft_section, "ENERGY_CORRECTION%XC%XC_FUNCTIONAL")
207         CALL section_vals_get(section2, explicit=explicit)
208         IF (explicit) THEN
209            CALL xc_functionals_expand(section2, section1)
210            ec_env%xc_section => section1
211         ELSE
212            ec_env%xc_section => xc_section
213         END IF
214         ! dispersion
215         ALLOCATE (dispersion_env)
216         NULLIFY (xc_section)
217         xc_section => ec_env%xc_section
218         CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env)
219         CALL qs_dispersion_env_set(dispersion_env, xc_section)
220         IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
221            NULLIFY (pp_section)
222            pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
223            CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
224         ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
225            NULLIFY (nl_section)
226            nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
227            CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
228         END IF
229         ec_env%dispersion_env => dispersion_env
230      END IF
231
232   END SUBROUTINE init_ec_env
233
234END MODULE ec_environment
235