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 : ! **************************************************************************************************
9 : !> \brief Routines for an external energy correction on top of a Kohn-Sham calculation
10 : !> \par History
11 : !> 10.2024 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ec_external
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE cell_types, ONLY: cell_type
17 : USE core_ppl, ONLY: build_core_ppl
18 : USE core_ppnl, ONLY: build_core_ppnl
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
21 : dbcsr_copy,&
22 : dbcsr_create,&
23 : dbcsr_p_type,&
24 : dbcsr_set,&
25 : dbcsr_type,&
26 : dbcsr_type_symmetric
27 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
28 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
29 : cp_dbcsr_sm_fm_multiply,&
30 : dbcsr_allocate_matrix_set,&
31 : dbcsr_deallocate_matrix_set
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_diag,&
37 : cp_fm_get_info,&
38 : cp_fm_maxabsval,&
39 : cp_fm_release,&
40 : cp_fm_set_all,&
41 : cp_fm_to_fm,&
42 : cp_fm_type
43 : USE cp_log_handling, ONLY: cp_get_default_logger,&
44 : cp_logger_get_default_unit_nr,&
45 : cp_logger_type
46 : USE ec_env_types, ONLY: energy_correction_type
47 : USE kinds, ONLY: default_string_length,&
48 : dp
49 : USE mathlib, ONLY: det_3x3
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE parallel_gemm_api, ONLY: parallel_gemm
52 : USE particle_types, ONLY: particle_type
53 : USE physcon, ONLY: pascal
54 : USE qs_core_energies, ONLY: calculate_ptrace
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_environment_type
57 : USE qs_force_types, ONLY: qs_force_type
58 : USE qs_kind_types, ONLY: qs_kind_type
59 : USE qs_kinetic, ONLY: build_kinetic_matrix
60 : USE qs_ks_types, ONLY: qs_ks_env_type
61 : USE qs_mo_types, ONLY: get_mo_set,&
62 : mo_set_type
63 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
64 : USE qs_overlap, ONLY: build_overlap_matrix
65 : USE qs_rho_types, ONLY: qs_rho_get,&
66 : qs_rho_type
67 : USE virial_methods, ONLY: one_third_sum_diag
68 : USE virial_types, ONLY: virial_type
69 : #include "./base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 :
73 : PRIVATE
74 :
75 : ! *** Global parameters ***
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_external'
78 :
79 : PUBLIC :: ec_ext_energy
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief External energy method
85 : !> \param qs_env ...
86 : !> \param ec_env ...
87 : !> \param calculate_forces ...
88 : !> \par History
89 : !> 10.2024 created
90 : !> \author JGH
91 : ! **************************************************************************************************
92 12 : SUBROUTINE ec_ext_energy(qs_env, ec_env, calculate_forces)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : TYPE(energy_correction_type), POINTER :: ec_env
95 : LOGICAL, INTENT(IN) :: calculate_forces
96 :
97 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_ext_energy'
98 :
99 : INTEGER :: handle, ispin, nocc, nspins, unit_nr
100 : REAL(KIND=dp) :: focc
101 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
102 12 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ, mo_ref
103 : TYPE(cp_fm_type), POINTER :: mo_coeff
104 : TYPE(cp_logger_type), POINTER :: logger
105 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
106 : TYPE(dft_control_type), POINTER :: dft_control
107 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
108 :
109 12 : CALL timeset(routineN, handle)
110 :
111 12 : CALL get_qs_env(qs_env, dft_control=dft_control)
112 12 : nspins = dft_control%nspins
113 :
114 12 : logger => cp_get_default_logger()
115 12 : IF (logger%para_env%is_source()) THEN
116 6 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
117 : ELSE
118 6 : unit_nr = -1
119 : END IF
120 :
121 12 : CALL get_qs_env(qs_env, mos=mos)
122 96 : ALLOCATE (cpmos(nspins), mo_ref(nspins), mo_occ(nspins))
123 24 : DO ispin = 1, nspins
124 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
125 12 : NULLIFY (fm_struct)
126 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
127 12 : template_fmstruct=mo_coeff%matrix_struct)
128 12 : CALL cp_fm_create(cpmos(ispin), fm_struct)
129 12 : CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
130 12 : CALL cp_fm_create(mo_ref(ispin), fm_struct)
131 12 : CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp)
132 12 : CALL cp_fm_create(mo_occ(ispin), fm_struct)
133 12 : CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
134 36 : CALL cp_fm_struct_release(fm_struct)
135 : END DO
136 :
137 12 : CALL cp_fm_release(ec_env%mo_occ)
138 12 : CALL cp_fm_release(ec_env%cpmos)
139 12 : IF (ec_env%debug_external) THEN
140 : !
141 12 : ec_env%mo_occ => mo_ref
142 12 : CALL ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr)
143 : !
144 12 : IF (calculate_forces) THEN
145 4 : focc = 2.0_dp
146 4 : IF (nspins == 1) focc = 4.0_dp
147 8 : DO ispin = 1, nspins
148 4 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
149 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_h(1, 1)%matrix, ec_env%mo_occ(ispin), &
150 : cpmos(ispin), nocc, &
151 8 : alpha=focc, beta=0.0_dp)
152 : END DO
153 : END IF
154 12 : ec_env%cpmos => cpmos
155 : ELSE
156 : ! get external information
157 : !deb CALL ec_ext_interface(qs_env, ec_env%ex, mo_ref, cpmos, calculate_forces, unit_nr)
158 0 : ec_env%mo_occ => mo_ref
159 0 : ec_env%cpmos => cpmos
160 0 : CPABORT("EC EXT NYA")
161 : END IF
162 :
163 12 : IF (calculate_forces) THEN
164 : ! check for orbital rotations
165 4 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
166 8 : DO ispin = 1, nspins
167 : CALL align_vectors(ec_env%cpmos(ispin), ec_env%mo_occ(ispin), mo_occ(ispin), &
168 8 : matrix_s(1)%matrix, unit_nr)
169 : END DO
170 : ! set up matrices for response
171 4 : CALL ec_ext_setup(qs_env, ec_env, .TRUE., unit_nr)
172 : ! orthogonality force
173 : CALL matrix_r_forces(qs_env, ec_env%cpmos, ec_env%mo_occ, &
174 4 : ec_env%matrix_w(1, 1)%matrix, unit_nr)
175 : ELSE
176 8 : CALL ec_ext_setup(qs_env, ec_env, .FALSE., unit_nr)
177 : END IF
178 :
179 12 : CALL cp_fm_release(mo_occ)
180 :
181 12 : CALL timestop(handle)
182 :
183 12 : END SUBROUTINE ec_ext_energy
184 :
185 : ! **************************************************************************************************
186 :
187 : ! **************************************************************************************************
188 : !> \brief ...
189 : !> \param qs_env ...
190 : !> \param ec_env ...
191 : !> \param calculate_forces ...
192 : !> \param unit_nr ...
193 : ! **************************************************************************************************
194 12 : SUBROUTINE ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr)
195 : TYPE(qs_environment_type), POINTER :: qs_env
196 : TYPE(energy_correction_type), POINTER :: ec_env
197 : LOGICAL, INTENT(IN) :: calculate_forces
198 : INTEGER, INTENT(IN) :: unit_nr
199 :
200 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_debug'
201 :
202 : CHARACTER(LEN=default_string_length) :: headline
203 : INTEGER :: handle, ispin, nocc, nspins
204 : TYPE(cp_fm_type), POINTER :: mo_coeff
205 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s
206 : TYPE(dft_control_type), POINTER :: dft_control
207 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
208 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
209 12 : POINTER :: sab_orb
210 : TYPE(qs_rho_type), POINTER :: rho
211 :
212 12 : CALL timeset(routineN, handle)
213 :
214 : CALL get_qs_env(qs_env=qs_env, &
215 : dft_control=dft_control, &
216 : sab_orb=sab_orb, &
217 : rho=rho, &
218 : matrix_s_kp=matrix_s, &
219 12 : matrix_h_kp=matrix_h)
220 :
221 12 : nspins = dft_control%nspins
222 12 : CALL get_qs_env(qs_env, mos=mos)
223 24 : DO ispin = 1, nspins
224 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
225 24 : CALL cp_fm_to_fm(mo_coeff, ec_env%mo_occ(ispin), nocc)
226 : END DO
227 :
228 : ! Core Hamiltonian matrix
229 12 : IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
230 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
231 12 : headline = "CORE HAMILTONIAN MATRIX"
232 12 : ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
233 : CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=TRIM(headline), &
234 12 : template=matrix_h(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
235 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
236 12 : CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
237 :
238 : ! Get density matrix of reference calculation
239 12 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
240 : ! Use Core energy as model energy
241 12 : CALL calculate_ptrace(ec_env%matrix_h, matrix_p, ec_env%ex, nspins)
242 :
243 12 : IF (calculate_forces) THEN
244 : ! force of model energy
245 4 : CALL ec_debug_force(qs_env, matrix_p, unit_nr)
246 : END IF
247 :
248 12 : CALL timestop(handle)
249 :
250 12 : END SUBROUTINE ec_ext_debug
251 :
252 : ! **************************************************************************************************
253 : !> \brief ...
254 : !> \param qs_env ...
255 : !> \param matrix_p ...
256 : !> \param unit_nr ...
257 : ! **************************************************************************************************
258 4 : SUBROUTINE ec_debug_force(qs_env, matrix_p, unit_nr)
259 : TYPE(qs_environment_type), POINTER :: qs_env
260 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
261 : INTEGER, INTENT(IN) :: unit_nr
262 :
263 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force'
264 :
265 : INTEGER :: handle, iounit, nder, nimages
266 4 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
267 : LOGICAL :: calculate_forces, debug_forces, &
268 : debug_stress, use_virial
269 : REAL(KIND=dp) :: eps_ppnl, fconv
270 : REAL(KIND=dp), DIMENSION(3) :: fodeb
271 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
272 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
273 : TYPE(cell_type), POINTER :: cell
274 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
275 : TYPE(dft_control_type), POINTER :: dft_control
276 : TYPE(mp_para_env_type), POINTER :: para_env
277 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
278 4 : POINTER :: sab_orb, sac_ppl, sap_ppnl
279 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
280 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
281 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
282 : TYPE(qs_ks_env_type), POINTER :: ks_env
283 : TYPE(virial_type), POINTER :: virial
284 :
285 4 : CALL timeset(routineN, handle)
286 :
287 4 : debug_forces = .TRUE.
288 4 : debug_stress = .TRUE.
289 4 : iounit = unit_nr
290 :
291 4 : calculate_forces = .TRUE.
292 :
293 : ! no k-points possible
294 4 : NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
295 : CALL get_qs_env(qs_env=qs_env, &
296 : cell=cell, &
297 : dft_control=dft_control, &
298 : force=force, &
299 : ks_env=ks_env, &
300 : para_env=para_env, &
301 4 : virial=virial)
302 4 : nimages = dft_control%nimages
303 4 : IF (nimages /= 1) THEN
304 0 : CPABORT("K-points not implemented")
305 : END IF
306 :
307 : ! check for virial
308 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
309 :
310 4 : fconv = 1.0E-9_dp*pascal/cell%deth
311 4 : IF (debug_stress .AND. use_virial) THEN
312 0 : sttot = virial%pv_virial
313 : END IF
314 :
315 : ! check for GAPW/GAPW_XC
316 4 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
317 0 : CPABORT("GAPW not implemented")
318 : END IF
319 :
320 : ! get neighbor lists
321 4 : NULLIFY (sab_orb, sac_ppl, sap_ppnl)
322 : CALL get_qs_env(qs_env=qs_env, &
323 4 : sab_orb=sab_orb, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
324 :
325 : ! initialize src matrix
326 4 : NULLIFY (scrm)
327 4 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
328 4 : ALLOCATE (scrm(1, 1)%matrix)
329 4 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix)
330 4 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
331 :
332 4 : nder = 1
333 4 : IF (SIZE(matrix_p, 1) == 2) THEN
334 : CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
335 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
336 : END IF
337 :
338 : ! kinetic energy
339 16 : IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
340 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic
341 : CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
342 : matrix_name="KINETIC ENERGY MATRIX", &
343 : basis_type="ORB", &
344 : sab_nl=sab_orb, calculate_forces=.TRUE., &
345 4 : matrixkp_p=matrix_p)
346 : IF (debug_forces) THEN
347 16 : fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
348 4 : CALL para_env%sum(fodeb)
349 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb
350 : END IF
351 4 : IF (debug_stress .AND. use_virial) THEN
352 0 : stdeb = fconv*(virial%pv_ekinetic - stdeb)
353 0 : CALL para_env%sum(stdeb)
354 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
355 0 : 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
356 : END IF
357 4 : IF (SIZE(matrix_p, 1) == 2) THEN
358 : CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
359 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
360 : END IF
361 :
362 : ! compute the ppl contribution to the core hamiltonian
363 4 : NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
364 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
365 4 : atomic_kind_set=atomic_kind_set)
366 :
367 4 : IF (ASSOCIATED(sac_ppl)) THEN
368 16 : IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
369 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
370 : CALL build_core_ppl(scrm, matrix_p, force, &
371 : virial, calculate_forces, use_virial, nder, &
372 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
373 4 : nimages, cell_to_index, "ORB")
374 : IF (calculate_forces .AND. debug_forces) THEN
375 16 : fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
376 4 : CALL para_env%sum(fodeb)
377 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb
378 : END IF
379 4 : IF (debug_stress .AND. use_virial) THEN
380 0 : stdeb = fconv*(virial%pv_ppl - stdeb)
381 0 : CALL para_env%sum(stdeb)
382 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
383 0 : 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb)
384 : END IF
385 : END IF
386 :
387 : ! compute the ppnl contribution to the core hamiltonian ***
388 4 : eps_ppnl = dft_control%qs_control%eps_ppnl
389 4 : IF (ASSOCIATED(sap_ppnl)) THEN
390 8 : IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
391 2 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
392 : CALL build_core_ppnl(scrm, matrix_p, force, &
393 : virial, calculate_forces, use_virial, nder, &
394 : qs_kind_set, atomic_kind_set, particle_set, &
395 : sab_orb, sap_ppnl, eps_ppnl, &
396 2 : nimages, cell_to_index, "ORB")
397 : IF (calculate_forces .AND. debug_forces) THEN
398 8 : fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
399 2 : CALL para_env%sum(fodeb)
400 2 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb
401 : END IF
402 2 : IF (debug_stress .AND. use_virial) THEN
403 0 : stdeb = fconv*(virial%pv_ppnl - stdeb)
404 0 : CALL para_env%sum(stdeb)
405 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
406 0 : 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb)
407 : END IF
408 : END IF
409 :
410 4 : IF (debug_stress .AND. use_virial) THEN
411 0 : stdeb = fconv*(virial%pv_virial - sttot)
412 0 : CALL para_env%sum(stdeb)
413 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
414 0 : 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
415 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") ' '
416 : END IF
417 :
418 : ! delete scr matrix
419 4 : CALL dbcsr_deallocate_matrix_set(scrm)
420 :
421 4 : CALL timestop(handle)
422 :
423 4 : END SUBROUTINE ec_debug_force
424 :
425 : ! **************************************************************************************************
426 :
427 : ! **************************************************************************************************
428 : !> \brief ...
429 : !> \param qs_env ...
430 : !> \param ec_env ...
431 : !> \param calc_forces ...
432 : !> \param unit_nr ...
433 : ! **************************************************************************************************
434 12 : SUBROUTINE ec_ext_setup(qs_env, ec_env, calc_forces, unit_nr)
435 : TYPE(qs_environment_type), POINTER :: qs_env
436 : TYPE(energy_correction_type), POINTER :: ec_env
437 : LOGICAL, INTENT(IN) :: calc_forces
438 : INTEGER, INTENT(IN) :: unit_nr
439 :
440 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_setup'
441 :
442 : CHARACTER(LEN=default_string_length) :: headline
443 : INTEGER :: handle, ispin, nao, nocc, nspins
444 : REAL(KIND=dp) :: a_max, c_max
445 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
446 : TYPE(cp_fm_type) :: emat, ksmo, smo
447 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks, matrix_p, matrix_s
448 : TYPE(dft_control_type), POINTER :: dft_control
449 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
450 12 : POINTER :: sab_orb
451 : TYPE(qs_rho_type), POINTER :: rho
452 :
453 12 : CALL timeset(routineN, handle)
454 :
455 : CALL get_qs_env(qs_env=qs_env, &
456 : dft_control=dft_control, &
457 : sab_orb=sab_orb, &
458 : rho=rho, &
459 : matrix_s_kp=matrix_s, &
460 : matrix_ks_kp=matrix_ks, &
461 12 : matrix_h_kp=matrix_h)
462 :
463 12 : nspins = dft_control%nspins
464 :
465 : ! KS Hamiltonian matrix
466 12 : IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
467 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
468 12 : headline = "HAMILTONIAN MATRIX"
469 24 : DO ispin = 1, nspins
470 12 : ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
471 : CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), &
472 12 : template=matrix_ks(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
473 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, sab_orb)
474 24 : CALL dbcsr_copy(ec_env%matrix_ks(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix)
475 : END DO
476 :
477 : ! Overlap matrix
478 12 : IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
479 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_s, 1, 1)
480 12 : headline = "OVERLAP MATRIX"
481 12 : ALLOCATE (ec_env%matrix_s(1, 1)%matrix)
482 : CALL dbcsr_create(ec_env%matrix_s(1, 1)%matrix, name=TRIM(headline), &
483 12 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
484 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_s(1, 1)%matrix, sab_orb)
485 12 : CALL dbcsr_copy(ec_env%matrix_s(1, 1)%matrix, matrix_s(1, 1)%matrix)
486 :
487 : ! density matrix
488 : ! Get density matrix of reference calculation
489 12 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
490 12 : IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
491 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
492 12 : headline = "DENSITY MATRIX"
493 24 : DO ispin = 1, nspins
494 12 : ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
495 : CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), &
496 12 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
497 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, sab_orb)
498 24 : CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
499 : END DO
500 :
501 12 : IF (calc_forces) THEN
502 : ! energy weighted density matrix
503 : ! for security, we recalculate W here (stored in qs_env)
504 4 : IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
505 4 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
506 4 : headline = "ENERGY WEIGHTED DENSITY MATRIX"
507 8 : DO ispin = 1, nspins
508 4 : ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
509 : CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), &
510 4 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
511 4 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, sab_orb)
512 8 : CALL dbcsr_set(ec_env%matrix_w(ispin, 1)%matrix, 0.0_dp)
513 : END DO
514 :
515 : ! hz matrix
516 4 : IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
517 4 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
518 4 : headline = "Hz MATRIX"
519 8 : DO ispin = 1, nspins
520 4 : ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
521 : CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, name=TRIM(headline), &
522 4 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
523 4 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_hz(ispin)%matrix, sab_orb)
524 8 : CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
525 : END DO
526 :
527 : ! Test for consistency of orbitals and KS matrix
528 8 : DO ispin = 1, nspins
529 4 : mat_struct => ec_env%mo_occ(ispin)%matrix_struct
530 4 : CALL cp_fm_create(ksmo, mat_struct)
531 4 : CALL cp_fm_get_info(ksmo, nrow_global=nao, ncol_global=nocc)
532 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%mo_occ(ispin), &
533 4 : ksmo, nocc, alpha=1.0_dp, beta=0.0_dp)
534 4 : CALL cp_fm_create(smo, mat_struct)
535 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_s(1, 1)%matrix, ec_env%mo_occ(ispin), &
536 4 : smo, nocc, alpha=1.0_dp, beta=0.0_dp)
537 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
538 4 : template_fmstruct=mat_struct)
539 4 : CALL cp_fm_create(emat, fm_struct)
540 4 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, ec_env%mo_occ(ispin), ksmo, 0.0_dp, emat)
541 4 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, -1.0_dp, smo, emat, 1.0_dp, ksmo)
542 4 : CALL cp_fm_maxabsval(ksmo, a_max)
543 4 : CALL cp_fm_struct_release(fm_struct)
544 4 : CALL cp_fm_release(smo)
545 4 : CALL cp_fm_release(ksmo)
546 4 : CALL cp_fm_release(emat)
547 4 : CALL cp_fm_maxabsval(ec_env%mo_occ(ispin), c_max)
548 24 : IF (unit_nr > 0) THEN
549 2 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO coeficients", ispin, c_max
550 2 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO gradients", ispin, a_max
551 : END IF
552 : END DO
553 : END IF
554 :
555 12 : CALL timestop(handle)
556 :
557 12 : END SUBROUTINE ec_ext_setup
558 :
559 : ! **************************************************************************************************
560 : !> \brief ...
561 : !> \param cpmos ...
562 : !> \param mo_ref ...
563 : !> \param mo_occ ...
564 : !> \param matrix_s ...
565 : !> \param unit_nr ...
566 : ! **************************************************************************************************
567 4 : SUBROUTINE align_vectors(cpmos, mo_ref, mo_occ, matrix_s, unit_nr)
568 : TYPE(cp_fm_type), INTENT(IN) :: cpmos, mo_ref, mo_occ
569 : TYPE(dbcsr_type), POINTER :: matrix_s
570 : INTEGER, INTENT(IN) :: unit_nr
571 :
572 : CHARACTER(LEN=*), PARAMETER :: routineN = 'align_vectors'
573 :
574 : INTEGER :: handle, i, nao, nocc, scg
575 : REAL(KIND=dp) :: a_max
576 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag
577 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
578 : TYPE(cp_fm_type) :: emat, smo
579 :
580 4 : CALL timeset(routineN, handle)
581 :
582 4 : mat_struct => cpmos%matrix_struct
583 4 : CALL cp_fm_create(smo, mat_struct)
584 4 : CALL cp_fm_get_info(smo, nrow_global=nao, ncol_global=nocc)
585 4 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_occ, smo, nocc, alpha=1.0_dp, beta=0.0_dp)
586 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
587 4 : template_fmstruct=mat_struct)
588 4 : CALL cp_fm_create(emat, fm_struct)
589 4 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, mo_ref, smo, 0.0_dp, emat)
590 4 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, cpmos, emat, 0.0_dp, smo)
591 4 : CALL cp_fm_to_fm(smo, cpmos)
592 4 : CALL cp_fm_to_fm(mo_occ, mo_ref)
593 : !
594 12 : ALLOCATE (diag(nocc))
595 4 : CALL cp_fm_get_diag(emat, diag)
596 14 : a_max = nocc - SUM(diag)
597 4 : scg = 0
598 14 : DO i = 1, nocc
599 14 : IF (ABS(diag(i) + 1.0_dp) < 0.001) scg = scg + 1
600 : END DO
601 4 : IF (unit_nr > 0) THEN
602 2 : WRITE (unit_nr, "(T3,A,T61,F20.8)") "External:: Orbital rotation index", a_max
603 2 : WRITE (unit_nr, "(T3,A,T71,I10)") "External:: Number of orbital phase changes", scg
604 : END IF
605 :
606 4 : DEALLOCATE (diag)
607 4 : CALL cp_fm_struct_release(fm_struct)
608 4 : CALL cp_fm_release(smo)
609 4 : CALL cp_fm_release(emat)
610 :
611 4 : CALL timestop(handle)
612 :
613 8 : END SUBROUTINE align_vectors
614 :
615 : ! **************************************************************************************************
616 : !> \brief ...
617 : !> \param qs_env ...
618 : !> \param matrix_w ...
619 : !> \param unit_nr ...
620 : ! **************************************************************************************************
621 0 : SUBROUTINE matrix_w_forces(qs_env, matrix_w, unit_nr)
622 : TYPE(qs_environment_type), POINTER :: qs_env
623 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w
624 : INTEGER, INTENT(IN) :: unit_nr
625 :
626 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_w_forces'
627 :
628 : INTEGER :: handle, iounit, nder, nimages
629 : LOGICAL :: debug_forces, debug_stress, use_virial
630 : REAL(KIND=dp) :: fconv
631 : REAL(KIND=dp), DIMENSION(3) :: fodeb
632 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
633 : TYPE(cell_type), POINTER :: cell
634 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
635 : TYPE(dft_control_type), POINTER :: dft_control
636 : TYPE(mp_para_env_type), POINTER :: para_env
637 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
638 0 : POINTER :: sab_orb
639 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
640 : TYPE(qs_ks_env_type), POINTER :: ks_env
641 : TYPE(virial_type), POINTER :: virial
642 :
643 0 : CALL timeset(routineN, handle)
644 :
645 0 : debug_forces = .TRUE.
646 0 : debug_stress = .TRUE.
647 :
648 0 : iounit = unit_nr
649 :
650 : ! no k-points possible
651 : CALL get_qs_env(qs_env=qs_env, &
652 : cell=cell, &
653 : dft_control=dft_control, &
654 : force=force, &
655 : ks_env=ks_env, &
656 : sab_orb=sab_orb, &
657 : para_env=para_env, &
658 0 : virial=virial)
659 0 : nimages = dft_control%nimages
660 0 : IF (nimages /= 1) THEN
661 0 : CPABORT("K-points not implemented")
662 : END IF
663 :
664 : ! check for virial
665 0 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
666 :
667 0 : fconv = 1.0E-9_dp*pascal/cell%deth
668 0 : IF (debug_stress .AND. use_virial) THEN
669 0 : sttot = virial%pv_virial
670 : END IF
671 :
672 : ! initialize src matrix
673 0 : NULLIFY (scrm)
674 0 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
675 0 : ALLOCATE (scrm(1, 1)%matrix)
676 0 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_w(1, 1)%matrix)
677 0 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
678 :
679 0 : nder = 1
680 0 : IF (SIZE(matrix_w, 1) == 2) THEN
681 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
682 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
683 : END IF
684 :
685 : ! Overlap and kinetic energy matrices
686 0 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
687 0 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
688 : CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
689 : matrix_name="OVERLAP MATRIX", &
690 : basis_type_a="ORB", &
691 : basis_type_b="ORB", &
692 : sab_nl=sab_orb, calculate_forces=.TRUE., &
693 0 : matrixkp_p=matrix_w)
694 :
695 : IF (debug_forces) THEN
696 0 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
697 0 : CALL para_env%sum(fodeb)
698 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
699 : END IF
700 0 : IF (debug_stress .AND. use_virial) THEN
701 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
702 0 : CALL para_env%sum(stdeb)
703 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
704 0 : 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
705 : END IF
706 0 : IF (SIZE(matrix_w, 1) == 2) THEN
707 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
708 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
709 : END IF
710 :
711 : ! delete scrm matrix
712 0 : CALL dbcsr_deallocate_matrix_set(scrm)
713 :
714 0 : CALL timestop(handle)
715 :
716 0 : END SUBROUTINE matrix_w_forces
717 :
718 : ! **************************************************************************************************
719 : !> \brief ...
720 : !> \param qs_env ...
721 : !> \param cpmos ...
722 : !> \param mo_occ ...
723 : !> \param matrix_r ...
724 : !> \param unit_nr ...
725 : ! **************************************************************************************************
726 4 : SUBROUTINE matrix_r_forces(qs_env, cpmos, mo_occ, matrix_r, unit_nr)
727 : TYPE(qs_environment_type), POINTER :: qs_env
728 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ
729 : TYPE(dbcsr_type), POINTER :: matrix_r
730 : INTEGER, INTENT(IN) :: unit_nr
731 :
732 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_r_forces'
733 :
734 : INTEGER :: handle, iounit, ispin, nao, nocc, nspins
735 : LOGICAL :: debug_forces, debug_stress, use_virial
736 : REAL(KIND=dp) :: fconv, focc
737 : REAL(KIND=dp), DIMENSION(3) :: fodeb
738 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
739 : TYPE(cell_type), POINTER :: cell
740 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
741 : TYPE(cp_fm_type) :: chcmat, rcvec
742 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: scrm
743 : TYPE(dft_control_type), POINTER :: dft_control
744 : TYPE(mp_para_env_type), POINTER :: para_env
745 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
746 4 : POINTER :: sab_orb
747 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
748 : TYPE(qs_ks_env_type), POINTER :: ks_env
749 : TYPE(virial_type), POINTER :: virial
750 :
751 4 : CALL timeset(routineN, handle)
752 :
753 4 : debug_forces = .TRUE.
754 4 : debug_stress = .TRUE.
755 4 : iounit = unit_nr
756 :
757 4 : nspins = SIZE(mo_occ)
758 4 : focc = 1.0_dp
759 4 : IF (nspins == 1) focc = 2.0_dp
760 4 : focc = 0.25_dp*focc
761 :
762 4 : CALL dbcsr_set(matrix_r, 0.0_dp)
763 8 : DO ispin = 1, nspins
764 4 : CALL cp_fm_get_info(cpmos(ispin), matrix_struct=fm_struct, nrow_global=nao, ncol_global=nocc)
765 4 : CALL cp_fm_create(rcvec, fm_struct)
766 4 : CALL cp_fm_struct_create(mat_struct, nrow_global=nocc, ncol_global=nocc, template_fmstruct=fm_struct)
767 4 : CALL cp_fm_create(chcmat, mat_struct)
768 4 : CALL parallel_gemm("T", "N", nocc, nocc, nao, 1.0_dp, mo_occ(ispin), cpmos(ispin), 0.0_dp, chcmat)
769 4 : CALL parallel_gemm("N", "N", nao, nocc, nocc, 1.0_dp, mo_occ(ispin), chcmat, 0.0_dp, rcvec)
770 4 : CALL cp_dbcsr_plus_fm_fm_t(matrix_r, matrix_v=rcvec, matrix_g=mo_occ(ispin), ncol=nocc, alpha=focc)
771 4 : CALL cp_fm_struct_release(mat_struct)
772 4 : CALL cp_fm_release(rcvec)
773 16 : CALL cp_fm_release(chcmat)
774 : END DO
775 :
776 : CALL get_qs_env(qs_env=qs_env, &
777 : cell=cell, &
778 : dft_control=dft_control, &
779 : force=force, &
780 : ks_env=ks_env, &
781 : sab_orb=sab_orb, &
782 : para_env=para_env, &
783 4 : virial=virial)
784 : ! check for virial
785 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
786 4 : fconv = 1.0E-9_dp*pascal/cell%deth
787 :
788 16 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
789 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
790 4 : NULLIFY (scrm)
791 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
792 : matrix_name="OVERLAP MATRIX", &
793 : basis_type_a="ORB", basis_type_b="ORB", &
794 : sab_nl=sab_orb, calculate_forces=.TRUE., &
795 4 : matrix_p=matrix_r)
796 : IF (debug_forces) THEN
797 16 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
798 4 : CALL para_env%sum(fodeb)
799 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
800 : END IF
801 4 : IF (debug_stress .AND. use_virial) THEN
802 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
803 0 : CALL para_env%sum(stdeb)
804 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
805 0 : 'STRESS| Wz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
806 : END IF
807 4 : CALL dbcsr_deallocate_matrix_set(scrm)
808 :
809 4 : CALL timestop(handle)
810 :
811 4 : END SUBROUTINE matrix_r_forces
812 :
813 : END MODULE ec_external
|
