Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE optimize_dmfet_potential
9 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
10 : cp_blacs_env_release,&
11 : cp_blacs_env_type
12 : USE cp_control_types, ONLY: dft_control_type
13 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr
14 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
15 : cp_fm_struct_release,&
16 : cp_fm_struct_type
17 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
18 : cp_fm_create,&
19 : cp_fm_maxabsval,&
20 : cp_fm_release,&
21 : cp_fm_set_all,&
22 : cp_fm_type
23 : USE dbcsr_api, ONLY: dbcsr_create,&
24 : dbcsr_init_p,&
25 : dbcsr_multiply,&
26 : dbcsr_p_type,&
27 : dbcsr_release,&
28 : dbcsr_trace,&
29 : dbcsr_type,&
30 : dbcsr_type_no_symmetry
31 : USE embed_types, ONLY: opt_dmfet_pot_type
32 : USE input_section_types, ONLY: section_vals_type,&
33 : section_vals_val_get
34 : USE kinds, ONLY: dp
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE parallel_gemm_api, ONLY: parallel_gemm
37 : USE qs_environment_types, ONLY: get_qs_env,&
38 : qs_environment_type
39 : USE qs_mo_types, ONLY: get_mo_set,&
40 : mo_set_type
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_dmfet_potential'
48 :
49 : PUBLIC :: build_full_dm, prepare_dmfet_opt, release_dmfet_opt, subsys_spin, check_dmfet
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief ...
55 : !> \param opt_dmfet ...
56 : !> \param opt_dmfet_section ...
57 : ! **************************************************************************************************
58 0 : SUBROUTINE read_opt_dmfet_section(opt_dmfet, opt_dmfet_section)
59 : TYPE(opt_dmfet_pot_type) :: opt_dmfet
60 : TYPE(section_vals_type), POINTER :: opt_dmfet_section
61 :
62 : ! Read keywords
63 :
64 : CALL section_vals_val_get(opt_dmfet_section, "N_ITER", &
65 0 : i_val=opt_dmfet%n_iter)
66 :
67 : CALL section_vals_val_get(opt_dmfet_section, "TRUST_RAD", &
68 0 : r_val=opt_dmfet%trust_rad)
69 :
70 : CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_MAX", &
71 0 : r_val=opt_dmfet%conv_max)
72 :
73 : CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_INT", &
74 0 : r_val=opt_dmfet%conv_int)
75 :
76 : CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_MAX", &
77 0 : r_val=opt_dmfet%conv_max_beta)
78 :
79 : CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_INT", &
80 0 : r_val=opt_dmfet%conv_int_beta)
81 :
82 : CALL section_vals_val_get(opt_dmfet_section, "READ_DMFET_POT", &
83 0 : l_val=opt_dmfet%read_dmfet_pot)
84 :
85 0 : END SUBROUTINE read_opt_dmfet_section
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param qs_env ...
90 : !> \return ...
91 : ! **************************************************************************************************
92 0 : FUNCTION subsys_spin(qs_env) RESULT(subsys_open_shell)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : LOGICAL :: subsys_open_shell
95 :
96 : TYPE(dft_control_type), POINTER :: dft_control
97 :
98 0 : NULLIFY (dft_control)
99 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
100 0 : subsys_open_shell = .FALSE.
101 0 : IF (dft_control%nspins == 2) subsys_open_shell = .TRUE.
102 :
103 0 : END FUNCTION subsys_spin
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param qs_env ...
108 : !> \param opt_dmfet ...
109 : !> \param opt_dmfet_section ...
110 : ! **************************************************************************************************
111 0 : SUBROUTINE prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section)
112 : TYPE(qs_environment_type), POINTER :: qs_env
113 : TYPE(opt_dmfet_pot_type) :: opt_dmfet
114 : TYPE(section_vals_type), POINTER :: opt_dmfet_section
115 :
116 : INTEGER :: diff_size, nao
117 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
118 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
119 : TYPE(cp_fm_type), POINTER :: mo_coeff
120 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
121 : TYPE(mp_para_env_type), POINTER :: para_env
122 :
123 : ! Read the input
124 0 : CALL read_opt_dmfet_section(opt_dmfet, opt_dmfet_section)
125 :
126 : ! Get the orbital coefficients
127 0 : CALL get_qs_env(qs_env=qs_env, mos=mos, para_env=para_env)
128 0 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, nao=nao)
129 :
130 : ! Make cp_fm matrices
131 0 : NULLIFY (blacs_env)
132 0 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
133 :
134 0 : NULLIFY (opt_dmfet%dmfet_pot, opt_dmfet%dm_1, opt_dmfet%dm_2, opt_dmfet%dm_total, opt_dmfet%dm_diff)
135 0 : DEALLOCATE (opt_dmfet%dmfet_pot, opt_dmfet%dm_subsys, opt_dmfet%dm_total, opt_dmfet%dm_diff)
136 : NULLIFY (fm_struct)
137 :
138 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
139 : nrow_global=nao, ncol_global=nao)
140 : CALL cp_fm_create(opt_dmfet%dmfet_pot, fm_struct, name="dmfet_pot")
141 : CALL cp_fm_create(opt_dmfet%dm_subsys, fm_struct, name="dm_subsys")
142 : CALL cp_fm_create(opt_dmfet%dm_total, fm_struct, name="dm_total")
143 : CALL cp_fm_create(opt_dmfet%dm_diff, fm_struct, name="dm_diff")
144 :
145 : CALL cp_fm_set_all(opt_dmfet%dmfet_pot, 0.0_dp)
146 : CALL cp_fm_set_all(opt_dmfet%dm_subsys, 0.0_dp)
147 : CALL cp_fm_set_all(opt_dmfet%dm_total, 0.0_dp)
148 : CALL cp_fm_set_all(opt_dmfet%dm_diff, 0.0_dp)
149 :
150 : ! Beta spin counterparts
151 : IF (opt_dmfet%open_shell_embed) THEN
152 : NULLIFY (opt_dmfet%dmfet_pot_beta, opt_dmfet%dm_subsys_beta, &
153 : opt_dmfet%dm_total_beta, opt_dmfet%dm_diff_beta)
154 : ALLOCATE (opt_dmfet%dmfet_pot_beta, opt_dmfet%dm_subsys_beta, &
155 : opt_dmfet%dm_total_beta, opt_dmfet%dm_diff_beta)
156 : CALL cp_fm_create(opt_dmfet%dmfet_pot_beta, fm_struct, name="dmfet_pot_beta")
157 : CALL cp_fm_create(opt_dmfet%dm_subsys_beta, fm_struct, name="dm_subsys_beta")
158 : CALL cp_fm_create(opt_dmfet%dm_total_beta, fm_struct, name="dm_total_beta")
159 : CALL cp_fm_create(opt_dmfet%dm_diff_beta, fm_struct, name="dm_diff_beta")
160 :
161 : CALL cp_fm_set_all(opt_dmfet%dmfet_pot_beta, 0.0_dp)
162 : CALL cp_fm_set_all(opt_dmfet%dm_subsys_beta, 0.0_dp)
163 : CALL cp_fm_set_all(opt_dmfet%dm_total_beta, 0.0_dp)
164 : CALL cp_fm_set_all(opt_dmfet%dm_diff_beta, 0.0_dp)
165 : END IF
166 :
167 : ! Release structure and context
168 : CALL cp_fm_struct_release(fm_struct)
169 : CALL cp_blacs_env_release(blacs_env)
170 :
171 : ! Array to store functional values
172 : ALLOCATE (opt_dmfet%w_func(opt_dmfet%n_iter))
173 : opt_dmfet%w_func = 0.0_dp
174 :
175 : ! Allocate max_diff and int_diff
176 : diff_size = 1
177 : IF (opt_dmfet%open_shell_embed) diff_size = 2
178 : ALLOCATE (opt_dmfet%max_diff(diff_size))
179 : ALLOCATE (opt_dmfet%int_diff(diff_size))
180 :
181 0 : END SUBROUTINE prepare_dmfet_opt
182 :
183 : ! **************************************************************************************************
184 : !> \brief ...
185 : !> \param opt_dmfet ...
186 : ! **************************************************************************************************
187 0 : SUBROUTINE release_dmfet_opt(opt_dmfet)
188 : TYPE(opt_dmfet_pot_type) :: opt_dmfet
189 :
190 0 : IF (ASSOCIATED(opt_dmfet%dmfet_pot)) THEN
191 0 : CALL cp_fm_release(opt_dmfet%dmfet_pot)
192 0 : DEALLOCATE (opt_dmfet%dmfet_pot)
193 : NULLIFY (opt_dmfet%dmfet_pot)
194 : END IF
195 0 : IF (ASSOCIATED(opt_dmfet%dm_subsys)) THEN
196 0 : CALL cp_fm_release(opt_dmfet%dm_subsys)
197 0 : DEALLOCATE (opt_dmfet%dm_subsys)
198 : NULLIFY (opt_dmfet%dm_subsys)
199 : END IF
200 0 : IF (ASSOCIATED(opt_dmfet%dm_total)) THEN
201 0 : CALL cp_fm_release(opt_dmfet%dm_total)
202 0 : DEALLOCATE (opt_dmfet%dm_total)
203 : NULLIFY (opt_dmfet%dm_total)
204 : END IF
205 0 : IF (ASSOCIATED(opt_dmfet%dm_diff)) THEN
206 0 : CALL cp_fm_release(opt_dmfet%dm_diff)
207 0 : DEALLOCATE (opt_dmfet%dm_diff)
208 : NULLIFY (opt_dmfet%dm_diff)
209 : END IF
210 :
211 0 : IF (opt_dmfet%open_shell_embed) THEN
212 0 : IF (ASSOCIATED(opt_dmfet%dmfet_pot_beta)) THEN
213 0 : CALL cp_fm_release(opt_dmfet%dmfet_pot_beta)
214 0 : DEALLOCATE (opt_dmfet%dmfet_pot_beta)
215 : NULLIFY (opt_dmfet%dmfet_pot_beta)
216 : END IF
217 0 : IF (ASSOCIATED(opt_dmfet%dm_subsys_beta)) THEN
218 0 : CALL cp_fm_release(opt_dmfet%dm_subsys_beta)
219 0 : DEALLOCATE (opt_dmfet%dm_subsys_beta)
220 : NULLIFY (opt_dmfet%dm_subsys_beta)
221 : END IF
222 0 : IF (ASSOCIATED(opt_dmfet%dm_total_beta)) THEN
223 0 : CALL cp_fm_release(opt_dmfet%dm_total_beta)
224 0 : DEALLOCATE (opt_dmfet%dm_total_beta)
225 : NULLIFY (opt_dmfet%dm_total_beta)
226 : END IF
227 0 : IF (ASSOCIATED(opt_dmfet%dm_diff_beta)) THEN
228 0 : CALL cp_fm_release(opt_dmfet%dm_diff_beta)
229 0 : DEALLOCATE (opt_dmfet%dm_diff_beta)
230 : NULLIFY (opt_dmfet%dm_diff_beta)
231 : END IF
232 : END IF
233 :
234 0 : DEALLOCATE (opt_dmfet%w_func)
235 0 : DEALLOCATE (opt_dmfet%max_diff)
236 0 : DEALLOCATE (opt_dmfet%int_diff)
237 0 : DEALLOCATE (opt_dmfet%all_nspins)
238 :
239 0 : END SUBROUTINE release_dmfet_opt
240 :
241 : ! **************************************************************************************************
242 : !> \brief Builds density matrices from MO coefficients in full matrix format
243 : !> \param qs_env ...
244 : !> \param dm ...
245 : !> \param open_shell Subsystem is open shell
246 : !> \param open_shell_embed Embedding is open shell
247 : !> \param dm_beta ...
248 : ! **************************************************************************************************
249 0 : SUBROUTINE build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta)
250 : TYPE(qs_environment_type), POINTER :: qs_env
251 : TYPE(cp_fm_type), INTENT(IN) :: dm
252 : LOGICAL :: open_shell, open_shell_embed
253 : TYPE(cp_fm_type), INTENT(IN) :: dm_beta
254 :
255 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_full_dm'
256 :
257 : INTEGER :: handle, homo, nao
258 : REAL(KIND=dp) :: coeff
259 : TYPE(cp_fm_type), POINTER :: mo_coeff
260 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
261 : TYPE(mp_para_env_type), POINTER :: para_env
262 :
263 0 : CALL timeset(routineN, handle)
264 :
265 0 : coeff = 2.0_dp
266 0 : IF (open_shell_embed) coeff = 1.0_dp
267 :
268 : ! Get the orbital coefficients
269 0 : CALL get_qs_env(qs_env=qs_env, mos=mos)
270 0 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, nao=nao, homo=homo)
271 :
272 : ! Build the density matrix
273 : CALL parallel_gemm(transa="N", transb="T", m=nao, n=nao, &
274 : k=homo, alpha=coeff, &
275 : matrix_a=mo_coeff, matrix_b=mo_coeff, &
276 0 : beta=0.0_dp, matrix_c=dm)
277 :
278 : ! Open shell
279 0 : IF (open_shell) THEN
280 0 : CALL get_mo_set(mo_set=mos(2), mo_coeff=mo_coeff, nao=nao, homo=homo)
281 :
282 : ! Build the density matrix
283 : CALL parallel_gemm(transa="N", transb="T", m=nao, n=nao, &
284 : k=homo, alpha=coeff, &
285 : matrix_a=mo_coeff, matrix_b=mo_coeff, &
286 0 : beta=0.0_dp, matrix_c=dm_beta)
287 : END IF
288 :
289 : ! If embedding is open shell, but subsystem is not, copy alpha-spin DM into beta-spin DM
290 0 : IF (open_shell_embed .AND. (.NOT. open_shell)) THEN
291 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
292 0 : CALL cp_fm_copy_general(dm, dm_beta, para_env)
293 :
294 : END IF
295 :
296 0 : CALL timestop(handle)
297 :
298 0 : END SUBROUTINE build_full_dm
299 :
300 : ! **************************************************************************************************
301 : !> \brief ...
302 : !> \param opt_dmfet ...
303 : !> \param qs_env ...
304 : ! **************************************************************************************************
305 0 : SUBROUTINE check_dmfet(opt_dmfet, qs_env)
306 : TYPE(opt_dmfet_pot_type) :: opt_dmfet
307 : TYPE(qs_environment_type), POINTER :: qs_env
308 :
309 : REAL(KIND=dp) :: max_diff, max_diff_beta, trace
310 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
311 : TYPE(dbcsr_type), POINTER :: diff_dbcsr, dm_s
312 :
313 0 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
314 :
315 0 : NULLIFY (diff_dbcsr)
316 0 : CALL dbcsr_init_p(diff_dbcsr)
317 : CALL dbcsr_create(matrix=diff_dbcsr, &
318 : template=matrix_s(1)%matrix, &
319 0 : matrix_type=dbcsr_type_no_symmetry)
320 :
321 0 : NULLIFY (dm_s)
322 0 : CALL dbcsr_init_p(dm_s)
323 : CALL dbcsr_create(matrix=dm_s, &
324 : template=matrix_s(1)%matrix, &
325 0 : matrix_type=dbcsr_type_no_symmetry)
326 :
327 0 : CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff, diff_dbcsr, keep_sparsity=.FALSE.)
328 :
329 : CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, &
330 0 : 0.0_dp, dm_s)
331 :
332 0 : CALL dbcsr_trace(dm_s, trace)
333 :
334 0 : IF (opt_dmfet%open_shell_embed) THEN
335 0 : CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff_beta, diff_dbcsr, keep_sparsity=.FALSE.)
336 :
337 : CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, &
338 0 : 0.0_dp, dm_s)
339 :
340 0 : CALL dbcsr_trace(dm_s, trace)
341 :
342 : END IF
343 :
344 : ! Release dbcsr objects
345 0 : CALL dbcsr_release(diff_dbcsr)
346 0 : CALL dbcsr_release(dm_s)
347 :
348 : ! Find the absolute maximum value of the DM difference
349 0 : CALL cp_fm_maxabsval(opt_dmfet%dm_diff, max_diff)
350 0 : IF (opt_dmfet%open_shell_embed) THEN
351 0 : CALL cp_fm_maxabsval(opt_dmfet%dm_diff_beta, max_diff_beta)
352 : END IF
353 :
354 0 : END SUBROUTINE check_dmfet
355 :
356 : END MODULE optimize_dmfet_potential
|
