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