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 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 0 : CALL ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr)
158 0 : ec_env%mo_occ => mo_ref
159 0 : ec_env%cpmos => cpmos
160 : END IF
161 :
162 12 : IF (calculate_forces) THEN
163 : ! check for orbital rotations
164 4 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
165 8 : DO ispin = 1, nspins
166 : CALL align_vectors(ec_env%cpmos(ispin), ec_env%mo_occ(ispin), mo_occ(ispin), &
167 8 : matrix_s(1)%matrix, unit_nr)
168 : END DO
169 : ! set up matrices for response
170 4 : CALL ec_ext_setup(qs_env, ec_env, .TRUE., unit_nr)
171 : ! orthogonality force
172 : CALL matrix_r_forces(qs_env, ec_env%cpmos, ec_env%mo_occ, &
173 4 : ec_env%matrix_w(1, 1)%matrix, unit_nr)
174 : ELSE
175 8 : CALL ec_ext_setup(qs_env, ec_env, .FALSE., unit_nr)
176 : END IF
177 :
178 12 : CALL cp_fm_release(mo_occ)
179 :
180 12 : CALL timestop(handle)
181 :
182 12 : END SUBROUTINE ec_ext_energy
183 :
184 : ! **************************************************************************************************
185 :
186 : ! **************************************************************************************************
187 : !> \brief ...
188 : !> \param qs_env ...
189 : !> \param ec_env ...
190 : !> \param mo_ref ...
191 : !> \param cpmos ...
192 : !> \param calculate_forces ...
193 : !> \param unit_nr ...
194 : ! **************************************************************************************************
195 0 : SUBROUTINE ec_ext_interface(qs_env, ec_env, mo_ref, cpmos, calculate_forces, unit_nr)
196 : TYPE(qs_environment_type), POINTER :: qs_env
197 : TYPE(energy_correction_type), POINTER :: ec_env
198 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_ref, cpmos
199 : LOGICAL, INTENT(IN) :: calculate_forces
200 : INTEGER, INTENT(IN) :: unit_nr
201 :
202 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_interface'
203 :
204 : INTEGER :: handle, ispin, nao, nocc(2), nspins
205 : TYPE(dft_control_type), POINTER :: dft_control
206 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
207 : TYPE(mp_para_env_type), POINTER :: para_env
208 :
209 0 : CALL timeset(routineN, handle)
210 :
211 : MARK_USED(qs_env)
212 : MARK_USED(ec_env)
213 : MARK_USED(mo_ref)
214 : MARK_USED(cpmos)
215 : MARK_USED(calculate_forces)
216 : MARK_USED(unit_nr)
217 :
218 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
219 :
220 0 : nspins = dft_control%nspins
221 0 : CALL get_qs_env(qs_env, mos=mos)
222 0 : nocc = 0
223 0 : DO ispin = 1, nspins
224 0 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, homo=nocc(ispin))
225 : END DO
226 :
227 0 : IF (unit_nr > 0) THEN
228 0 : WRITE (unit_nr, '(T2,A)') "Read EXTERNAL Response from file: "//TRIM(ec_env%exresp_fn)
229 : END IF
230 :
231 0 : CALL timestop(handle)
232 :
233 0 : END SUBROUTINE ec_ext_interface
234 :
235 : ! **************************************************************************************************
236 :
237 : ! **************************************************************************************************
238 : !> \brief ...
239 : !> \param qs_env ...
240 : !> \param ec_env ...
241 : !> \param calculate_forces ...
242 : !> \param unit_nr ...
243 : ! **************************************************************************************************
244 12 : SUBROUTINE ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr)
245 : TYPE(qs_environment_type), POINTER :: qs_env
246 : TYPE(energy_correction_type), POINTER :: ec_env
247 : LOGICAL, INTENT(IN) :: calculate_forces
248 : INTEGER, INTENT(IN) :: unit_nr
249 :
250 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_debug'
251 :
252 : CHARACTER(LEN=default_string_length) :: headline
253 : INTEGER :: handle, ispin, nocc, nspins
254 : TYPE(cp_fm_type), POINTER :: mo_coeff
255 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s
256 : TYPE(dft_control_type), POINTER :: dft_control
257 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
258 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
259 12 : POINTER :: sab_orb
260 : TYPE(qs_rho_type), POINTER :: rho
261 :
262 12 : CALL timeset(routineN, handle)
263 :
264 : CALL get_qs_env(qs_env=qs_env, &
265 : dft_control=dft_control, &
266 : sab_orb=sab_orb, &
267 : rho=rho, &
268 : matrix_s_kp=matrix_s, &
269 12 : matrix_h_kp=matrix_h)
270 :
271 12 : nspins = dft_control%nspins
272 12 : CALL get_qs_env(qs_env, mos=mos)
273 24 : DO ispin = 1, nspins
274 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
275 24 : CALL cp_fm_to_fm(mo_coeff, ec_env%mo_occ(ispin), nocc)
276 : END DO
277 :
278 : ! Core Hamiltonian matrix
279 12 : IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
280 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
281 12 : headline = "CORE HAMILTONIAN MATRIX"
282 12 : ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
283 : CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=TRIM(headline), &
284 12 : template=matrix_h(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
285 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
286 12 : CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
287 :
288 : ! Get density matrix of reference calculation
289 12 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
290 : ! Use Core energy as model energy
291 12 : CALL calculate_ptrace(ec_env%matrix_h, matrix_p, ec_env%ex, nspins)
292 :
293 12 : IF (calculate_forces) THEN
294 : ! force of model energy
295 4 : CALL ec_debug_force(qs_env, matrix_p, unit_nr)
296 : END IF
297 :
298 12 : CALL timestop(handle)
299 :
300 12 : END SUBROUTINE ec_ext_debug
301 :
302 : ! **************************************************************************************************
303 : !> \brief ...
304 : !> \param qs_env ...
305 : !> \param matrix_p ...
306 : !> \param unit_nr ...
307 : ! **************************************************************************************************
308 4 : SUBROUTINE ec_debug_force(qs_env, matrix_p, unit_nr)
309 : TYPE(qs_environment_type), POINTER :: qs_env
310 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
311 : INTEGER, INTENT(IN) :: unit_nr
312 :
313 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force'
314 :
315 : INTEGER :: handle, iounit, nder, nimages
316 4 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
317 : LOGICAL :: calculate_forces, debug_forces, &
318 : debug_stress, use_virial
319 : REAL(KIND=dp) :: eps_ppnl, fconv
320 : REAL(KIND=dp), DIMENSION(3) :: fodeb
321 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
322 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
323 : TYPE(cell_type), POINTER :: cell
324 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
325 : TYPE(dft_control_type), POINTER :: dft_control
326 : TYPE(mp_para_env_type), POINTER :: para_env
327 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
328 4 : POINTER :: sab_orb, sac_ppl, sap_ppnl
329 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
330 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
331 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
332 : TYPE(qs_ks_env_type), POINTER :: ks_env
333 : TYPE(virial_type), POINTER :: virial
334 :
335 4 : CALL timeset(routineN, handle)
336 :
337 4 : debug_forces = .TRUE.
338 4 : debug_stress = .TRUE.
339 4 : iounit = unit_nr
340 :
341 4 : calculate_forces = .TRUE.
342 :
343 : ! no k-points possible
344 4 : NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
345 : CALL get_qs_env(qs_env=qs_env, &
346 : cell=cell, &
347 : dft_control=dft_control, &
348 : force=force, &
349 : ks_env=ks_env, &
350 : para_env=para_env, &
351 4 : virial=virial)
352 4 : nimages = dft_control%nimages
353 4 : IF (nimages /= 1) THEN
354 0 : CPABORT("K-points not implemented")
355 : END IF
356 :
357 : ! check for virial
358 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
359 :
360 4 : fconv = 1.0E-9_dp*pascal/cell%deth
361 4 : IF (debug_stress .AND. use_virial) THEN
362 0 : sttot = virial%pv_virial
363 : END IF
364 :
365 : ! check for GAPW/GAPW_XC
366 4 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
367 0 : CPABORT("GAPW not implemented")
368 : END IF
369 :
370 : ! get neighbor lists
371 4 : NULLIFY (sab_orb, sac_ppl, sap_ppnl)
372 : CALL get_qs_env(qs_env=qs_env, &
373 4 : sab_orb=sab_orb, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
374 :
375 : ! initialize src matrix
376 4 : NULLIFY (scrm)
377 4 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
378 4 : ALLOCATE (scrm(1, 1)%matrix)
379 4 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix)
380 4 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
381 :
382 4 : nder = 1
383 4 : IF (SIZE(matrix_p, 1) == 2) THEN
384 : CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
385 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
386 : END IF
387 :
388 : ! kinetic energy
389 16 : IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
390 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic
391 : CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
392 : matrix_name="KINETIC ENERGY MATRIX", &
393 : basis_type="ORB", &
394 : sab_nl=sab_orb, calculate_forces=.TRUE., &
395 4 : matrixkp_p=matrix_p)
396 : IF (debug_forces) THEN
397 16 : fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
398 4 : CALL para_env%sum(fodeb)
399 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb
400 : END IF
401 4 : IF (debug_stress .AND. use_virial) THEN
402 0 : stdeb = fconv*(virial%pv_ekinetic - stdeb)
403 0 : CALL para_env%sum(stdeb)
404 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
405 0 : 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
406 : END IF
407 4 : IF (SIZE(matrix_p, 1) == 2) THEN
408 : CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
409 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
410 : END IF
411 :
412 : ! compute the ppl contribution to the core hamiltonian
413 4 : NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
414 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
415 4 : atomic_kind_set=atomic_kind_set)
416 :
417 4 : IF (ASSOCIATED(sac_ppl)) THEN
418 16 : IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
419 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
420 : CALL build_core_ppl(scrm, matrix_p, force, &
421 : virial, calculate_forces, use_virial, nder, &
422 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
423 4 : nimages, cell_to_index, "ORB")
424 : IF (calculate_forces .AND. debug_forces) THEN
425 16 : fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
426 4 : CALL para_env%sum(fodeb)
427 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb
428 : END IF
429 4 : IF (debug_stress .AND. use_virial) THEN
430 0 : stdeb = fconv*(virial%pv_ppl - stdeb)
431 0 : CALL para_env%sum(stdeb)
432 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
433 0 : 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb)
434 : END IF
435 : END IF
436 :
437 : ! compute the ppnl contribution to the core hamiltonian ***
438 4 : eps_ppnl = dft_control%qs_control%eps_ppnl
439 4 : IF (ASSOCIATED(sap_ppnl)) THEN
440 8 : IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
441 2 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
442 : CALL build_core_ppnl(scrm, matrix_p, force, &
443 : virial, calculate_forces, use_virial, nder, &
444 : qs_kind_set, atomic_kind_set, particle_set, &
445 : sab_orb, sap_ppnl, eps_ppnl, &
446 2 : nimages, cell_to_index, "ORB")
447 : IF (calculate_forces .AND. debug_forces) THEN
448 8 : fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
449 2 : CALL para_env%sum(fodeb)
450 2 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb
451 : END IF
452 2 : IF (debug_stress .AND. use_virial) THEN
453 0 : stdeb = fconv*(virial%pv_ppnl - stdeb)
454 0 : CALL para_env%sum(stdeb)
455 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
456 0 : 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb)
457 : END IF
458 : END IF
459 :
460 4 : IF (debug_stress .AND. use_virial) THEN
461 0 : stdeb = fconv*(virial%pv_virial - sttot)
462 0 : CALL para_env%sum(stdeb)
463 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
464 0 : 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
465 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") ' '
466 : END IF
467 :
468 : ! delete scr matrix
469 4 : CALL dbcsr_deallocate_matrix_set(scrm)
470 :
471 4 : CALL timestop(handle)
472 :
473 4 : END SUBROUTINE ec_debug_force
474 :
475 : ! **************************************************************************************************
476 :
477 : ! **************************************************************************************************
478 : !> \brief ...
479 : !> \param qs_env ...
480 : !> \param ec_env ...
481 : !> \param calc_forces ...
482 : !> \param unit_nr ...
483 : ! **************************************************************************************************
484 12 : SUBROUTINE ec_ext_setup(qs_env, ec_env, calc_forces, unit_nr)
485 : TYPE(qs_environment_type), POINTER :: qs_env
486 : TYPE(energy_correction_type), POINTER :: ec_env
487 : LOGICAL, INTENT(IN) :: calc_forces
488 : INTEGER, INTENT(IN) :: unit_nr
489 :
490 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_setup'
491 :
492 : CHARACTER(LEN=default_string_length) :: headline
493 : INTEGER :: handle, ispin, nao, nocc, nspins
494 : REAL(KIND=dp) :: a_max, c_max
495 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
496 : TYPE(cp_fm_type) :: emat, ksmo, smo
497 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks, matrix_p, matrix_s
498 : TYPE(dft_control_type), POINTER :: dft_control
499 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
500 12 : POINTER :: sab_orb
501 : TYPE(qs_rho_type), POINTER :: rho
502 :
503 12 : CALL timeset(routineN, handle)
504 :
505 : CALL get_qs_env(qs_env=qs_env, &
506 : dft_control=dft_control, &
507 : sab_orb=sab_orb, &
508 : rho=rho, &
509 : matrix_s_kp=matrix_s, &
510 : matrix_ks_kp=matrix_ks, &
511 12 : matrix_h_kp=matrix_h)
512 :
513 12 : nspins = dft_control%nspins
514 :
515 : ! KS Hamiltonian matrix
516 12 : IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
517 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
518 12 : headline = "HAMILTONIAN MATRIX"
519 24 : DO ispin = 1, nspins
520 12 : ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
521 : CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), &
522 12 : template=matrix_ks(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
523 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, sab_orb)
524 24 : CALL dbcsr_copy(ec_env%matrix_ks(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix)
525 : END DO
526 :
527 : ! Overlap matrix
528 12 : IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
529 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_s, 1, 1)
530 12 : headline = "OVERLAP MATRIX"
531 12 : ALLOCATE (ec_env%matrix_s(1, 1)%matrix)
532 : CALL dbcsr_create(ec_env%matrix_s(1, 1)%matrix, name=TRIM(headline), &
533 12 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
534 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_s(1, 1)%matrix, sab_orb)
535 12 : CALL dbcsr_copy(ec_env%matrix_s(1, 1)%matrix, matrix_s(1, 1)%matrix)
536 :
537 : ! density matrix
538 : ! Get density matrix of reference calculation
539 12 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
540 12 : IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
541 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
542 12 : headline = "DENSITY MATRIX"
543 24 : DO ispin = 1, nspins
544 12 : ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
545 : CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), &
546 12 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
547 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, sab_orb)
548 24 : CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
549 : END DO
550 :
551 12 : IF (calc_forces) THEN
552 : ! energy weighted density matrix
553 : ! for security, we recalculate W here (stored in qs_env)
554 4 : IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
555 4 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
556 4 : headline = "ENERGY WEIGHTED DENSITY MATRIX"
557 8 : DO ispin = 1, nspins
558 4 : ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
559 : CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), &
560 4 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
561 4 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, sab_orb)
562 8 : CALL dbcsr_set(ec_env%matrix_w(ispin, 1)%matrix, 0.0_dp)
563 : END DO
564 :
565 : ! hz matrix
566 4 : IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
567 4 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
568 4 : headline = "Hz MATRIX"
569 8 : DO ispin = 1, nspins
570 4 : ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
571 : CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, name=TRIM(headline), &
572 4 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
573 4 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_hz(ispin)%matrix, sab_orb)
574 8 : CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
575 : END DO
576 :
577 : ! Test for consistency of orbitals and KS matrix
578 8 : DO ispin = 1, nspins
579 4 : mat_struct => ec_env%mo_occ(ispin)%matrix_struct
580 4 : CALL cp_fm_create(ksmo, mat_struct)
581 4 : CALL cp_fm_get_info(ksmo, nrow_global=nao, ncol_global=nocc)
582 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%mo_occ(ispin), &
583 4 : ksmo, nocc, alpha=1.0_dp, beta=0.0_dp)
584 4 : CALL cp_fm_create(smo, mat_struct)
585 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_s(1, 1)%matrix, ec_env%mo_occ(ispin), &
586 4 : smo, nocc, alpha=1.0_dp, beta=0.0_dp)
587 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
588 4 : template_fmstruct=mat_struct)
589 4 : CALL cp_fm_create(emat, fm_struct)
590 4 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, ec_env%mo_occ(ispin), ksmo, 0.0_dp, emat)
591 4 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, -1.0_dp, smo, emat, 1.0_dp, ksmo)
592 4 : CALL cp_fm_maxabsval(ksmo, a_max)
593 4 : CALL cp_fm_struct_release(fm_struct)
594 4 : CALL cp_fm_release(smo)
595 4 : CALL cp_fm_release(ksmo)
596 4 : CALL cp_fm_release(emat)
597 4 : CALL cp_fm_maxabsval(ec_env%mo_occ(ispin), c_max)
598 24 : IF (unit_nr > 0) THEN
599 2 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO coeficients", ispin, c_max
600 2 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO gradients", ispin, a_max
601 : END IF
602 : END DO
603 : END IF
604 :
605 12 : CALL timestop(handle)
606 :
607 12 : END SUBROUTINE ec_ext_setup
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param cpmos ...
612 : !> \param mo_ref ...
613 : !> \param mo_occ ...
614 : !> \param matrix_s ...
615 : !> \param unit_nr ...
616 : ! **************************************************************************************************
617 4 : SUBROUTINE align_vectors(cpmos, mo_ref, mo_occ, matrix_s, unit_nr)
618 : TYPE(cp_fm_type), INTENT(IN) :: cpmos, mo_ref, mo_occ
619 : TYPE(dbcsr_type), POINTER :: matrix_s
620 : INTEGER, INTENT(IN) :: unit_nr
621 :
622 : CHARACTER(LEN=*), PARAMETER :: routineN = 'align_vectors'
623 :
624 : INTEGER :: handle, i, nao, nocc, scg
625 : REAL(KIND=dp) :: a_max
626 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag
627 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
628 : TYPE(cp_fm_type) :: emat, smo
629 :
630 4 : CALL timeset(routineN, handle)
631 :
632 4 : mat_struct => cpmos%matrix_struct
633 4 : CALL cp_fm_create(smo, mat_struct)
634 4 : CALL cp_fm_get_info(smo, nrow_global=nao, ncol_global=nocc)
635 4 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_occ, smo, nocc, alpha=1.0_dp, beta=0.0_dp)
636 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
637 4 : template_fmstruct=mat_struct)
638 4 : CALL cp_fm_create(emat, fm_struct)
639 4 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, mo_ref, smo, 0.0_dp, emat)
640 4 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, cpmos, emat, 0.0_dp, smo)
641 4 : CALL cp_fm_to_fm(smo, cpmos)
642 4 : CALL cp_fm_to_fm(mo_occ, mo_ref)
643 : !
644 12 : ALLOCATE (diag(nocc))
645 4 : CALL cp_fm_get_diag(emat, diag)
646 14 : a_max = nocc - SUM(diag)
647 4 : scg = 0
648 14 : DO i = 1, nocc
649 14 : IF (ABS(diag(i) + 1.0_dp) < 0.001) scg = scg + 1
650 : END DO
651 4 : IF (unit_nr > 0) THEN
652 2 : WRITE (unit_nr, "(T3,A,T61,F20.8)") "External:: Orbital rotation index", a_max
653 2 : WRITE (unit_nr, "(T3,A,T71,I10)") "External:: Number of orbital phase changes", scg
654 : END IF
655 :
656 4 : DEALLOCATE (diag)
657 4 : CALL cp_fm_struct_release(fm_struct)
658 4 : CALL cp_fm_release(smo)
659 4 : CALL cp_fm_release(emat)
660 :
661 4 : CALL timestop(handle)
662 :
663 8 : END SUBROUTINE align_vectors
664 :
665 : ! **************************************************************************************************
666 : !> \brief ...
667 : !> \param qs_env ...
668 : !> \param matrix_w ...
669 : !> \param unit_nr ...
670 : ! **************************************************************************************************
671 0 : SUBROUTINE matrix_w_forces(qs_env, matrix_w, unit_nr)
672 : TYPE(qs_environment_type), POINTER :: qs_env
673 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w
674 : INTEGER, INTENT(IN) :: unit_nr
675 :
676 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_w_forces'
677 :
678 : INTEGER :: handle, iounit, nder, nimages
679 : LOGICAL :: debug_forces, debug_stress, use_virial
680 : REAL(KIND=dp) :: fconv
681 : REAL(KIND=dp), DIMENSION(3) :: fodeb
682 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
683 : TYPE(cell_type), POINTER :: cell
684 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
685 : TYPE(dft_control_type), POINTER :: dft_control
686 : TYPE(mp_para_env_type), POINTER :: para_env
687 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
688 0 : POINTER :: sab_orb
689 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
690 : TYPE(qs_ks_env_type), POINTER :: ks_env
691 : TYPE(virial_type), POINTER :: virial
692 :
693 0 : CALL timeset(routineN, handle)
694 :
695 0 : debug_forces = .TRUE.
696 0 : debug_stress = .TRUE.
697 :
698 0 : iounit = unit_nr
699 :
700 : ! no k-points possible
701 : CALL get_qs_env(qs_env=qs_env, &
702 : cell=cell, &
703 : dft_control=dft_control, &
704 : force=force, &
705 : ks_env=ks_env, &
706 : sab_orb=sab_orb, &
707 : para_env=para_env, &
708 0 : virial=virial)
709 0 : nimages = dft_control%nimages
710 0 : IF (nimages /= 1) THEN
711 0 : CPABORT("K-points not implemented")
712 : END IF
713 :
714 : ! check for virial
715 0 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
716 :
717 0 : fconv = 1.0E-9_dp*pascal/cell%deth
718 0 : IF (debug_stress .AND. use_virial) THEN
719 0 : sttot = virial%pv_virial
720 : END IF
721 :
722 : ! initialize src matrix
723 0 : NULLIFY (scrm)
724 0 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
725 0 : ALLOCATE (scrm(1, 1)%matrix)
726 0 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_w(1, 1)%matrix)
727 0 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
728 :
729 0 : nder = 1
730 0 : IF (SIZE(matrix_w, 1) == 2) THEN
731 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
732 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
733 : END IF
734 :
735 : ! Overlap and kinetic energy matrices
736 0 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
737 0 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
738 : CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
739 : matrix_name="OVERLAP MATRIX", &
740 : basis_type_a="ORB", &
741 : basis_type_b="ORB", &
742 : sab_nl=sab_orb, calculate_forces=.TRUE., &
743 0 : matrixkp_p=matrix_w)
744 :
745 : IF (debug_forces) THEN
746 0 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
747 0 : CALL para_env%sum(fodeb)
748 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
749 : END IF
750 0 : IF (debug_stress .AND. use_virial) THEN
751 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
752 0 : CALL para_env%sum(stdeb)
753 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
754 0 : 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
755 : END IF
756 0 : IF (SIZE(matrix_w, 1) == 2) THEN
757 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
758 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
759 : END IF
760 :
761 : ! delete scrm matrix
762 0 : CALL dbcsr_deallocate_matrix_set(scrm)
763 :
764 0 : CALL timestop(handle)
765 :
766 0 : END SUBROUTINE matrix_w_forces
767 :
768 : ! **************************************************************************************************
769 : !> \brief ...
770 : !> \param qs_env ...
771 : !> \param cpmos ...
772 : !> \param mo_occ ...
773 : !> \param matrix_r ...
774 : !> \param unit_nr ...
775 : ! **************************************************************************************************
776 4 : SUBROUTINE matrix_r_forces(qs_env, cpmos, mo_occ, matrix_r, unit_nr)
777 : TYPE(qs_environment_type), POINTER :: qs_env
778 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ
779 : TYPE(dbcsr_type), POINTER :: matrix_r
780 : INTEGER, INTENT(IN) :: unit_nr
781 :
782 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_r_forces'
783 :
784 : INTEGER :: handle, iounit, ispin, nao, nocc, nspins
785 : LOGICAL :: debug_forces, debug_stress, use_virial
786 : REAL(KIND=dp) :: fconv, focc
787 : REAL(KIND=dp), DIMENSION(3) :: fodeb
788 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
789 : TYPE(cell_type), POINTER :: cell
790 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
791 : TYPE(cp_fm_type) :: chcmat, rcvec
792 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: scrm
793 : TYPE(dft_control_type), POINTER :: dft_control
794 : TYPE(mp_para_env_type), POINTER :: para_env
795 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
796 4 : POINTER :: sab_orb
797 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
798 : TYPE(qs_ks_env_type), POINTER :: ks_env
799 : TYPE(virial_type), POINTER :: virial
800 :
801 4 : CALL timeset(routineN, handle)
802 :
803 4 : debug_forces = .TRUE.
804 4 : debug_stress = .TRUE.
805 4 : iounit = unit_nr
806 :
807 4 : nspins = SIZE(mo_occ)
808 4 : focc = 1.0_dp
809 4 : IF (nspins == 1) focc = 2.0_dp
810 4 : focc = 0.25_dp*focc
811 :
812 4 : CALL dbcsr_set(matrix_r, 0.0_dp)
813 8 : DO ispin = 1, nspins
814 4 : CALL cp_fm_get_info(cpmos(ispin), matrix_struct=fm_struct, nrow_global=nao, ncol_global=nocc)
815 4 : CALL cp_fm_create(rcvec, fm_struct)
816 4 : CALL cp_fm_struct_create(mat_struct, nrow_global=nocc, ncol_global=nocc, template_fmstruct=fm_struct)
817 4 : CALL cp_fm_create(chcmat, mat_struct)
818 4 : CALL parallel_gemm("T", "N", nocc, nocc, nao, 1.0_dp, mo_occ(ispin), cpmos(ispin), 0.0_dp, chcmat)
819 4 : CALL parallel_gemm("N", "N", nao, nocc, nocc, 1.0_dp, mo_occ(ispin), chcmat, 0.0_dp, rcvec)
820 4 : CALL cp_dbcsr_plus_fm_fm_t(matrix_r, matrix_v=rcvec, matrix_g=mo_occ(ispin), ncol=nocc, alpha=focc)
821 4 : CALL cp_fm_struct_release(mat_struct)
822 4 : CALL cp_fm_release(rcvec)
823 16 : CALL cp_fm_release(chcmat)
824 : END DO
825 :
826 : CALL get_qs_env(qs_env=qs_env, &
827 : cell=cell, &
828 : dft_control=dft_control, &
829 : force=force, &
830 : ks_env=ks_env, &
831 : sab_orb=sab_orb, &
832 : para_env=para_env, &
833 4 : virial=virial)
834 : ! check for virial
835 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
836 4 : fconv = 1.0E-9_dp*pascal/cell%deth
837 :
838 16 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
839 4 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
840 4 : NULLIFY (scrm)
841 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
842 : matrix_name="OVERLAP MATRIX", &
843 : basis_type_a="ORB", basis_type_b="ORB", &
844 : sab_nl=sab_orb, calculate_forces=.TRUE., &
845 4 : matrix_p=matrix_r)
846 : IF (debug_forces) THEN
847 16 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
848 4 : CALL para_env%sum(fodeb)
849 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
850 : END IF
851 4 : IF (debug_stress .AND. use_virial) THEN
852 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
853 0 : CALL para_env%sum(stdeb)
854 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
855 0 : 'STRESS| Wz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
856 : END IF
857 4 : CALL dbcsr_deallocate_matrix_set(scrm)
858 :
859 4 : CALL timestop(handle)
860 :
861 4 : END SUBROUTINE matrix_r_forces
862 :
863 : END MODULE ec_external
|
