Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines used for Harris functional
10 : !> Kohn-Sham calculation
11 : !> \par History
12 : !> 10.2020 created
13 : !> \author Fabian Belleflamme
14 : ! **************************************************************************************************
15 : MODULE ec_methods
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_init_p,&
19 : dbcsr_type,&
20 : dbcsr_type_no_symmetry
21 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_m_by_n_from_row_template
22 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 : cp_fm_struct_release,&
24 : cp_fm_struct_type
25 : USE cp_fm_types, ONLY: cp_fm_get_info,&
26 : cp_fm_type
27 : USE cp_log_handling, ONLY: cp_to_string
28 : USE input_section_types, ONLY: section_vals_type
29 : USE kinds, ONLY: dp
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE pw_env_types, ONLY: pw_env_get,&
32 : pw_env_type
33 : USE pw_pool_types, ONLY: pw_pool_type
34 : USE pw_types, ONLY: pw_c1d_gs_type,&
35 : pw_r3d_rs_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type,&
38 : set_qs_env
39 : USE qs_kind_types, ONLY: get_qs_kind_set,&
40 : qs_kind_type
41 : USE qs_matrix_pools, ONLY: mpools_release,&
42 : qs_matrix_pools_type
43 : USE qs_mo_types, ONLY: allocate_mo_set,&
44 : get_mo_set,&
45 : init_mo_set,&
46 : mo_set_type
47 : USE qs_rho_types, ONLY: qs_rho_get,&
48 : qs_rho_type
49 : USE xc, ONLY: xc_calc_2nd_deriv,&
50 : xc_prep_2nd_deriv
51 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
52 : xc_dset_release
53 : USE xc_rho_set_types, ONLY: xc_rho_set_release,&
54 : xc_rho_set_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : ! *** Global parameters ***
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_methods'
64 :
65 : PUBLIC :: create_kernel
66 : PUBLIC :: ec_mos_init
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief Creation of second derivative xc-potential
72 : !> \param qs_env ...
73 : !> \param vxc will contain the partially integrated second derivative
74 : !> taken with respect to rho, evaluated in rho and folded with rho1
75 : !> vxc is allocated here and needs to be deallocated by the caller.
76 : !> \param vxc_tau ...
77 : !> \param rho density at which derivatives were calculated
78 : !> \param rho1_r density in r-space with which to fold
79 : !> \param rho1_g density in g-space with which to fold
80 : !> \param tau1_r ...
81 : !> \param xc_section XC parameters of this potential
82 : !> \param compute_virial Enable stress-tensor calculation
83 : !> \param virial_xc Will contain GGA contribution of XC-kernel to stress-tensor
84 : !> \date 11.2019
85 : !> \author fbelle
86 : ! **************************************************************************************************
87 4920 : SUBROUTINE create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, &
88 : compute_virial, virial_xc)
89 :
90 : TYPE(qs_environment_type), POINTER :: qs_env
91 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc, vxc_tau
92 : TYPE(qs_rho_type), INTENT(IN), POINTER :: rho
93 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), &
94 : POINTER :: rho1_r
95 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), &
96 : POINTER :: rho1_g
97 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), &
98 : POINTER :: tau1_r
99 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
100 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
101 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
102 : OPTIONAL :: virial_xc
103 :
104 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kernel'
105 :
106 : INTEGER :: handle
107 : TYPE(pw_env_type), POINTER :: pw_env
108 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
109 4920 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
110 : TYPE(pw_r3d_rs_type), POINTER :: weights
111 : TYPE(xc_derivative_set_type) :: deriv_set
112 : TYPE(xc_rho_set_type) :: rho_set
113 :
114 2460 : CALL timeset(routineN, handle)
115 :
116 2460 : NULLIFY (auxbas_pw_pool, pw_env, rho_r, vxc, vxc_tau)
117 :
118 2460 : CALL get_qs_env(qs_env, pw_env=pw_env, xcint_weights=weights)
119 2460 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
120 : ! Get density
121 2460 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
122 :
123 : CALL xc_prep_2nd_deriv(deriv_set=deriv_set, & ! containing potentials
124 : rho_set=rho_set, & ! density at which derivs are calculated
125 : rho_r=rho_r, & ! place where derivative is evaluated
126 : tau_r=tau_r, &
127 : pw_pool=auxbas_pw_pool, & ! pool for grids
128 : weights=weights, &
129 2460 : xc_section=xc_section)
130 :
131 : ! folding of second deriv with density in rho1_set
132 : CALL xc_calc_2nd_deriv(v_xc=vxc, & ! XC-potential, rho-dependent
133 : v_xc_tau=vxc_tau, & ! XC-potential, tau-dependent
134 : deriv_set=deriv_set, & ! deriv of xc-potential
135 : rho_set=rho_set, & ! density at which deriv are calculated
136 : rho1_r=rho1_r, & ! density with which to fold
137 : rho1_g=rho1_g, & ! density with which to fold
138 : tau1_r=tau1_r, &
139 : pw_pool=auxbas_pw_pool, & ! pool for grids
140 : weights=weights, &
141 : xc_section=xc_section, &
142 : gapw=.FALSE., &
143 : compute_virial=compute_virial, &
144 2460 : virial_xc=virial_xc)
145 :
146 : ! Release second deriv stuff
147 2460 : CALL xc_dset_release(deriv_set)
148 2460 : CALL xc_rho_set_release(rho_set=rho_set, pw_pool=auxbas_pw_pool)
149 :
150 2460 : CALL timestop(handle)
151 :
152 54120 : END SUBROUTINE create_kernel
153 :
154 : ! **************************************************************************************************
155 : !> \brief Allocate and initiate molecular orbitals environment
156 : !>
157 : !> \param qs_env ...
158 : !> \param matrix_s Used as template
159 : !> \param
160 : !>
161 : !> \par History
162 : !> 2020.10 created [Fabian Belleflamme]
163 : !> \author Fabian Belleflamme
164 : ! **************************************************************************************************
165 10 : SUBROUTINE ec_mos_init(qs_env, matrix_s)
166 : TYPE(qs_environment_type), POINTER :: qs_env
167 : TYPE(dbcsr_type) :: matrix_s
168 :
169 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_mos_init'
170 :
171 : INTEGER :: handle, ispin, multiplicity, n_ao, &
172 : nelectron, nmo, nspins
173 : INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
174 : REAL(dp) :: maxocc
175 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
176 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
177 : TYPE(cp_fm_type), POINTER :: mo_coeff
178 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
179 : TYPE(dft_control_type), POINTER :: dft_control
180 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
181 : TYPE(mp_para_env_type), POINTER :: para_env
182 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
183 : TYPE(qs_matrix_pools_type), POINTER :: my_mpools
184 :
185 10 : CALL timeset(routineN, handle)
186 :
187 10 : NULLIFY (blacs_env, dft_control, mo_coeff, mo_coeff_b, mos, my_mpools, qs_kind_set)
188 :
189 : CALL get_qs_env(qs_env=qs_env, &
190 : dft_control=dft_control, &
191 : blacs_env=blacs_env, &
192 : qs_kind_set=qs_kind_set, &
193 : nelectron_spin=nelectron_spin, &
194 10 : para_env=para_env)
195 10 : nspins = dft_control%nspins
196 :
197 : ! Start setup
198 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)
199 :
200 : ! the total number of electrons
201 10 : nelectron = nelectron - dft_control%charge
202 10 : multiplicity = dft_control%multiplicity
203 :
204 : ! setting maxocc and n_mo
205 10 : IF (dft_control%nspins == 1) THEN
206 10 : maxocc = 2.0_dp
207 10 : nelectron_spin(1) = nelectron
208 10 : nelectron_spin(2) = 0
209 10 : IF (MODULO(nelectron, 2) == 0) THEN
210 10 : n_mo(1) = nelectron/2
211 : ELSE
212 0 : n_mo(1) = INT(nelectron/2._dp) + 1
213 : END IF
214 10 : n_mo(2) = 0
215 : ELSE
216 0 : maxocc = 1.0_dp
217 :
218 : ! The simplist spin distribution is written here. Special cases will
219 : ! need additional user input
220 0 : IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN
221 0 : CPABORT("LSD: try to use a different multiplicity")
222 : END IF
223 :
224 0 : nelectron_spin(1) = (nelectron + multiplicity - 1)/2
225 0 : nelectron_spin(2) = (nelectron - multiplicity + 1)/2
226 :
227 0 : IF (nelectron_spin(2) < 0) THEN
228 0 : CPABORT("LSD: too few electrons for this multiplicity")
229 : END IF
230 :
231 0 : n_mo(1) = nelectron_spin(1)
232 0 : n_mo(2) = nelectron_spin(2)
233 :
234 : END IF
235 :
236 : ! Allocate MO set
237 40 : ALLOCATE (mos(nspins))
238 20 : DO ispin = 1, nspins
239 : CALL allocate_mo_set(mo_set=mos(ispin), &
240 : nao=n_ao, &
241 : nmo=n_mo(ispin), &
242 : nelectron=nelectron_spin(ispin), &
243 : n_el_f=REAL(nelectron_spin(ispin), dp), &
244 : maxocc=maxocc, &
245 20 : flexible_electron_count=dft_control%relax_multiplicity)
246 : END DO
247 :
248 10 : CALL set_qs_env(qs_env, mos=mos)
249 :
250 : ! finish initialization of the MOs
251 10 : NULLIFY (mo_coeff, mo_coeff_b)
252 20 : DO ispin = 1, SIZE(mos)
253 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
254 10 : nmo=nmo, nao=n_ao)
255 :
256 10 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
257 : CALL cp_fm_struct_create(fm_struct, nrow_global=n_ao, &
258 : ncol_global=nmo, para_env=para_env, &
259 10 : context=blacs_env)
260 :
261 : CALL init_mo_set(mos(ispin), &
262 : fm_struct=fm_struct, &
263 10 : name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
264 10 : CALL cp_fm_struct_release(fm_struct)
265 : END IF
266 :
267 30 : IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
268 10 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=nmo)
269 10 : CALL dbcsr_init_p(mos(ispin)%mo_coeff_b)
270 : CALL cp_dbcsr_m_by_n_from_row_template(mos(ispin)%mo_coeff_b, &
271 : template=matrix_s, &
272 : n=nmo, &
273 10 : sym=dbcsr_type_no_symmetry)
274 : END IF
275 : END DO
276 :
277 10 : CALL mpools_release(mpools=my_mpools)
278 :
279 10 : CALL timestop(handle)
280 :
281 20 : END SUBROUTINE ec_mos_init
282 :
283 : END MODULE ec_methods
|
