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 4900 : SUBROUTINE create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
88 :
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc, vxc_tau
91 : TYPE(qs_rho_type), INTENT(IN), POINTER :: rho
92 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), &
93 : POINTER :: rho1_r
94 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), &
95 : POINTER :: rho1_g
96 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), &
97 : POINTER :: tau1_r
98 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
99 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
100 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
101 : OPTIONAL :: virial_xc
102 :
103 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kernel'
104 :
105 : INTEGER :: handle
106 : TYPE(pw_env_type), POINTER :: pw_env
107 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
108 4900 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
109 : TYPE(pw_r3d_rs_type), POINTER :: weights
110 : TYPE(xc_derivative_set_type) :: deriv_set
111 : TYPE(xc_rho_set_type) :: rho_set
112 :
113 2450 : CALL timeset(routineN, handle)
114 :
115 2450 : NULLIFY (auxbas_pw_pool, pw_env, rho_r, vxc, vxc_tau)
116 :
117 2450 : CALL get_qs_env(qs_env, pw_env=pw_env, xcint_weights=weights)
118 2450 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
119 : ! Get density
120 2450 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
121 :
122 : CALL xc_prep_2nd_deriv(deriv_set=deriv_set, & ! containing potentials
123 : rho_set=rho_set, & ! density at which derivs are calculated
124 : rho_r=rho_r, & ! place where derivative is evaluated
125 : tau_r=tau_r, &
126 : pw_pool=auxbas_pw_pool, & ! pool for grids
127 : weights=weights, &
128 2450 : xc_section=xc_section)
129 :
130 : ! folding of second deriv with density in rho1_set
131 : CALL xc_calc_2nd_deriv(v_xc=vxc, & ! XC-potential, rho-dependent
132 : v_xc_tau=vxc_tau, & ! XC-potential, tau-dependent
133 : deriv_set=deriv_set, & ! deriv of xc-potential
134 : rho_set=rho_set, & ! density at which deriv are calculated
135 : rho1_r=rho1_r, & ! density with which to fold
136 : rho1_g=rho1_g, & ! density with which to fold
137 : tau1_r=tau1_r, &
138 : pw_pool=auxbas_pw_pool, & ! pool for grids
139 : weights=weights, &
140 : xc_section=xc_section, &
141 : gapw=.FALSE., &
142 : compute_virial=compute_virial, &
143 2450 : virial_xc=virial_xc)
144 :
145 : ! Release second deriv stuff
146 2450 : CALL xc_dset_release(deriv_set)
147 2450 : CALL xc_rho_set_release(rho_set=rho_set, pw_pool=auxbas_pw_pool)
148 :
149 2450 : CALL timestop(handle)
150 :
151 53900 : END SUBROUTINE create_kernel
152 :
153 : ! **************************************************************************************************
154 : !> \brief Allocate and initiate molecular orbitals environment
155 : !>
156 : !> \param qs_env ...
157 : !> \param matrix_s Used as template
158 : !> \param
159 : !>
160 : !> \par History
161 : !> 2020.10 created [Fabian Belleflamme]
162 : !> \author Fabian Belleflamme
163 : ! **************************************************************************************************
164 10 : SUBROUTINE ec_mos_init(qs_env, matrix_s)
165 : TYPE(qs_environment_type), POINTER :: qs_env
166 : TYPE(dbcsr_type) :: matrix_s
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_mos_init'
169 :
170 : INTEGER :: handle, ispin, multiplicity, n_ao, &
171 : nelectron, nmo, nspins
172 : INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
173 : REAL(dp) :: maxocc
174 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
175 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
176 : TYPE(cp_fm_type), POINTER :: mo_coeff
177 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
178 : TYPE(dft_control_type), POINTER :: dft_control
179 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
180 : TYPE(mp_para_env_type), POINTER :: para_env
181 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
182 : TYPE(qs_matrix_pools_type), POINTER :: my_mpools
183 :
184 10 : CALL timeset(routineN, handle)
185 :
186 10 : NULLIFY (blacs_env, dft_control, mo_coeff, mo_coeff_b, mos, my_mpools, qs_kind_set)
187 :
188 : CALL get_qs_env(qs_env=qs_env, &
189 : dft_control=dft_control, &
190 : blacs_env=blacs_env, &
191 : qs_kind_set=qs_kind_set, &
192 : nelectron_spin=nelectron_spin, &
193 10 : para_env=para_env)
194 10 : nspins = dft_control%nspins
195 :
196 : ! Start setup
197 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)
198 :
199 : ! the total number of electrons
200 10 : nelectron = nelectron - dft_control%charge
201 10 : multiplicity = dft_control%multiplicity
202 :
203 : ! setting maxocc and n_mo
204 10 : IF (dft_control%nspins == 1) THEN
205 10 : maxocc = 2.0_dp
206 10 : nelectron_spin(1) = nelectron
207 10 : nelectron_spin(2) = 0
208 10 : IF (MODULO(nelectron, 2) == 0) THEN
209 10 : n_mo(1) = nelectron/2
210 : ELSE
211 0 : n_mo(1) = INT(nelectron/2._dp) + 1
212 : END IF
213 10 : n_mo(2) = 0
214 : ELSE
215 0 : maxocc = 1.0_dp
216 :
217 : ! The simplist spin distribution is written here. Special cases will
218 : ! need additional user input
219 0 : IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN
220 0 : CPABORT("LSD: try to use a different multiplicity")
221 : END IF
222 :
223 0 : nelectron_spin(1) = (nelectron + multiplicity - 1)/2
224 0 : nelectron_spin(2) = (nelectron - multiplicity + 1)/2
225 :
226 0 : IF (nelectron_spin(2) < 0) THEN
227 0 : CPABORT("LSD: too few electrons for this multiplicity")
228 : END IF
229 :
230 0 : n_mo(1) = nelectron_spin(1)
231 0 : n_mo(2) = nelectron_spin(2)
232 :
233 : END IF
234 :
235 : ! Allocate MO set
236 40 : ALLOCATE (mos(nspins))
237 20 : DO ispin = 1, nspins
238 : CALL allocate_mo_set(mo_set=mos(ispin), &
239 : nao=n_ao, &
240 : nmo=n_mo(ispin), &
241 : nelectron=nelectron_spin(ispin), &
242 : n_el_f=REAL(nelectron_spin(ispin), dp), &
243 : maxocc=maxocc, &
244 20 : flexible_electron_count=dft_control%relax_multiplicity)
245 : END DO
246 :
247 10 : CALL set_qs_env(qs_env, mos=mos)
248 :
249 : ! finish initialization of the MOs
250 10 : NULLIFY (mo_coeff, mo_coeff_b)
251 20 : DO ispin = 1, SIZE(mos)
252 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
253 10 : nmo=nmo, nao=n_ao)
254 :
255 10 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
256 : CALL cp_fm_struct_create(fm_struct, nrow_global=n_ao, &
257 : ncol_global=nmo, para_env=para_env, &
258 10 : context=blacs_env)
259 :
260 : CALL init_mo_set(mos(ispin), &
261 : fm_struct=fm_struct, &
262 10 : name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
263 10 : CALL cp_fm_struct_release(fm_struct)
264 : END IF
265 :
266 30 : IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
267 10 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=nmo)
268 10 : CALL dbcsr_init_p(mos(ispin)%mo_coeff_b)
269 : CALL cp_dbcsr_m_by_n_from_row_template(mos(ispin)%mo_coeff_b, &
270 : template=matrix_s, &
271 : n=nmo, &
272 10 : sym=dbcsr_type_no_symmetry)
273 : END IF
274 : END DO
275 :
276 10 : CALL mpools_release(mpools=my_mpools)
277 :
278 10 : CALL timestop(handle)
279 :
280 20 : END SUBROUTINE ec_mos_init
281 :
282 : END MODULE ec_methods
|
