Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Chemical shift calculation by dfpt
10 : !> Initialization of the issc_env, creation of the special neighbor lists
11 : !> Perturbation Hamiltonians by application of the p and rxp oprtators to psi0
12 : !> Write output
13 : !> Deallocate everything
14 : !> \note
15 : !> The psi0 should be localized
16 : !> the Sebastiani method works within the assumption that the orbitals are
17 : !> completely contained in the simulation box
18 : ! **************************************************************************************************
19 : MODULE qs_linres_issc_utils
20 : USE atomic_kind_types, ONLY: atomic_kind_type,&
21 : get_atomic_kind
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
26 : dbcsr_copy,&
27 : dbcsr_create,&
28 : dbcsr_distribution_type,&
29 : dbcsr_p_type,&
30 : dbcsr_set,&
31 : dbcsr_type_antisymmetric,&
32 : dbcsr_type_symmetric
33 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
34 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
35 : dbcsr_allocate_matrix_set,&
36 : dbcsr_deallocate_matrix_set
37 : USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,&
38 : cp_fm_trace
39 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
40 : cp_fm_struct_release,&
41 : cp_fm_struct_type
42 : USE cp_fm_types, ONLY: cp_fm_create,&
43 : cp_fm_get_info,&
44 : cp_fm_release,&
45 : cp_fm_set_all,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_get_default_io_unit,&
50 : cp_logger_type
51 : USE cp_output_handling, ONLY: cp_p_file,&
52 : cp_print_key_finished_output,&
53 : cp_print_key_should_output,&
54 : cp_print_key_unit_nr
55 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
56 : section_vals_type,&
57 : section_vals_val_get
58 : USE kinds, ONLY: default_string_length,&
59 : dp
60 : USE mathlib, ONLY: diamat_all
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE particle_methods, ONLY: get_particle_set
64 : USE particle_types, ONLY: particle_type
65 : USE physcon, ONLY: a_fine,&
66 : e_mass,&
67 : hertz,&
68 : p_mass
69 : USE qs_elec_field, ONLY: build_efg_matrix
70 : USE qs_environment_types, ONLY: get_qs_env,&
71 : qs_environment_type
72 : USE qs_fermi_contact, ONLY: build_fermi_contact_matrix
73 : USE qs_kind_types, ONLY: qs_kind_type
74 : USE qs_linres_methods, ONLY: linres_solver
75 : USE qs_linres_types, ONLY: get_issc_env,&
76 : issc_env_type,&
77 : linres_control_type
78 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
79 : USE qs_mo_types, ONLY: get_mo_set,&
80 : mo_set_type
81 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
82 : USE qs_p_env_types, ONLY: qs_p_env_type
83 : USE qs_spin_orbit, ONLY: build_pso_matrix
84 : #include "./base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 : PUBLIC :: issc_env_cleanup, issc_env_init, issc_response, issc_issc, issc_print
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_issc_utils'
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief Initialize the issc environment
97 : !> \param issc_env ...
98 : !> \param p_env ...
99 : !> \param qs_env ...
100 : ! **************************************************************************************************
101 44 : SUBROUTINE issc_response(issc_env, p_env, qs_env)
102 : !
103 : TYPE(issc_env_type) :: issc_env
104 : TYPE(qs_p_env_type) :: p_env
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 :
107 : CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_response'
108 :
109 : INTEGER :: handle, idir, ijdir, ispin, jdir, nao, &
110 : nmo, nspins, output_unit
111 : LOGICAL :: do_dso, do_fc, do_pso, do_sd, should_stop
112 : REAL(dp) :: chk, fro
113 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
114 44 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
115 44 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc
116 44 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dso_psi0, efg_psi0, psi1_dso, psi1_efg, &
117 44 : psi1_pso, pso_psi0
118 : TYPE(cp_fm_type), POINTER :: mo_coeff
119 : TYPE(cp_logger_type), POINTER :: logger
120 : TYPE(dft_control_type), POINTER :: dft_control
121 : TYPE(linres_control_type), POINTER :: linres_control
122 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
123 : TYPE(mp_para_env_type), POINTER :: para_env
124 : TYPE(qs_matrix_pools_type), POINTER :: mpools
125 : TYPE(section_vals_type), POINTER :: issc_section, lr_section
126 :
127 44 : CALL timeset(routineN, handle)
128 : !
129 44 : NULLIFY (dft_control, linres_control, lr_section, issc_section)
130 44 : NULLIFY (logger, mpools, mo_coeff, para_env)
131 44 : NULLIFY (tmp_fm_struct, psi1_fc, psi1_efg, psi1_pso, pso_psi0, fc_psi0, efg_psi0)
132 :
133 44 : logger => cp_get_default_logger()
134 44 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
135 : issc_section => section_vals_get_subs_vals(qs_env%input, &
136 44 : "PROPERTIES%LINRES%SPINSPIN")
137 :
138 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
139 44 : extension=".linresLog")
140 44 : IF (output_unit > 0) THEN
141 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
142 22 : "*** Self consistent optimization of the response wavefunctions ***"
143 : END IF
144 :
145 : CALL get_qs_env(qs_env=qs_env, &
146 : dft_control=dft_control, &
147 : mpools=mpools, &
148 : linres_control=linres_control, &
149 : mos=mos, &
150 44 : para_env=para_env)
151 :
152 44 : nspins = dft_control%nspins
153 :
154 : CALL get_issc_env(issc_env=issc_env, &
155 : !list_cubes=list_cubes, &
156 : psi1_efg=psi1_efg, &
157 : psi1_pso=psi1_pso, &
158 : psi1_dso=psi1_dso, &
159 : psi1_fc=psi1_fc, &
160 : efg_psi0=efg_psi0, &
161 : pso_psi0=pso_psi0, &
162 : dso_psi0=dso_psi0, &
163 : fc_psi0=fc_psi0, &
164 : do_fc=do_fc, &
165 : do_sd=do_sd, &
166 : do_pso=do_pso, &
167 44 : do_dso=do_dso)
168 : !
169 : ! allocate the vectors
170 180 : ALLOCATE (psi0_order(nspins))
171 228 : ALLOCATE (psi1(nspins), h1_psi0(nspins))
172 92 : DO ispin = 1, nspins
173 48 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
174 48 : psi0_order(ispin) = mo_coeff
175 48 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
176 48 : NULLIFY (tmp_fm_struct)
177 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
178 : ncol_global=nmo, &
179 48 : context=mo_coeff%matrix_struct%context)
180 48 : CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
181 48 : CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
182 140 : CALL cp_fm_struct_release(tmp_fm_struct)
183 : END DO
184 44 : chk = 0.0_dp
185 : should_stop = .FALSE.
186 : !
187 : ! operator efg
188 44 : IF (do_sd) THEN
189 : ijdir = 0
190 0 : DO idir = 1, 3
191 0 : DO jdir = idir, 3
192 0 : ijdir = ijdir + 1
193 0 : DO ispin = 1, nspins
194 0 : CALL cp_fm_set_all(psi1_efg(ispin, ijdir), 0.0_dp)
195 : END DO
196 0 : IF (output_unit > 0) THEN
197 0 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator efg_"//ACHAR(idir + 119)//ACHAR(jdir + 119)
198 : END IF
199 : !
200 : !Initial guess for psi1
201 0 : DO ispin = 1, nspins
202 0 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
203 : END DO
204 : !
205 0 : DO ispin = 1, nspins
206 0 : CALL cp_fm_to_fm(efg_psi0(ispin, ijdir), h1_psi0(ispin))
207 : END DO
208 : !
209 : !
210 0 : linres_control%lr_triplet = .FALSE.
211 0 : linres_control%do_kernel = .FALSE.
212 0 : linres_control%converged = .FALSE.
213 0 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
214 : !
215 : !
216 : ! copy the response
217 0 : DO ispin = 1, nspins
218 0 : CALL cp_fm_to_fm(psi1(ispin), psi1_efg(ispin, ijdir))
219 0 : fro = cp_fm_frobenius_norm(psi1(ispin))
220 0 : chk = chk + fro
221 : END DO
222 : !
223 : !
224 0 : IF (output_unit > 0) THEN
225 0 : WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
226 : END IF
227 : !
228 : ! Write the result in the restart file
229 : END DO ! jdir
230 : END DO ! idir
231 : END IF
232 : !
233 : ! operator pso
234 44 : IF (do_pso) THEN
235 136 : DO idir = 1, 3
236 216 : DO ispin = 1, nspins
237 216 : CALL cp_fm_set_all(psi1_pso(ispin, idir), 0.0_dp)
238 : END DO
239 102 : IF (output_unit > 0) THEN
240 51 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator pso_"//ACHAR(idir + 119)
241 : END IF
242 : !
243 : !Initial guess for psi1
244 216 : DO ispin = 1, nspins
245 216 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
246 : !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin))
247 : !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
248 : END DO
249 : !
250 : !DO scf cycle to optimize psi1
251 216 : DO ispin = 1, nspins
252 216 : CALL cp_fm_to_fm(pso_psi0(ispin, idir), h1_psi0(ispin))
253 : END DO
254 : !
255 : !
256 102 : linres_control%lr_triplet = .FALSE. ! we do singlet response
257 102 : linres_control%do_kernel = .FALSE. ! we do uncoupled response
258 102 : linres_control%converged = .FALSE.
259 102 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
260 : !
261 : !
262 : ! copy the response
263 216 : DO ispin = 1, nspins
264 114 : CALL cp_fm_to_fm(psi1(ispin), psi1_pso(ispin, idir))
265 114 : fro = cp_fm_frobenius_norm(psi1(ispin))
266 216 : chk = chk + fro
267 : END DO
268 : !
269 : !
270 238 : IF (output_unit > 0) THEN
271 51 : WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
272 : END IF
273 : !
274 : ! Write the result in the restart file
275 : END DO ! idir
276 : END IF
277 : !
278 : ! operator fc
279 44 : IF (do_fc) THEN
280 0 : DO ispin = 1, nspins
281 0 : CALL cp_fm_set_all(psi1_fc(ispin), 0.0_dp)
282 : END DO
283 0 : IF (output_unit > 0) THEN
284 0 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator fc"
285 : END IF
286 : !
287 : !Initial guess for psi1
288 0 : DO ispin = 1, nspins
289 0 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
290 : END DO
291 : !
292 0 : DO ispin = 1, nspins
293 0 : CALL cp_fm_to_fm(fc_psi0(ispin), h1_psi0(ispin))
294 : END DO
295 : !
296 : !
297 0 : linres_control%lr_triplet = .TRUE. ! we do triplet response
298 0 : linres_control%do_kernel = .TRUE. ! we do coupled response
299 0 : linres_control%converged = .FALSE.
300 0 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
301 : !
302 : !
303 : ! copy the response
304 0 : DO ispin = 1, nspins
305 0 : CALL cp_fm_to_fm(psi1(ispin), psi1_fc(ispin))
306 0 : fro = cp_fm_frobenius_norm(psi1(ispin))
307 0 : chk = chk + fro
308 : END DO
309 : !
310 : !
311 0 : IF (output_unit > 0) THEN
312 0 : WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
313 : END IF
314 : !
315 : ! Write the result in the restart file
316 : END IF
317 :
318 : !>>>> debugging only
319 : !
320 : ! here we have the operator r and compute the polarizability for debugging the kernel only
321 44 : IF (do_dso) THEN
322 8 : DO idir = 1, 3
323 12 : DO ispin = 1, nspins
324 12 : CALL cp_fm_set_all(psi1_dso(ispin, idir), 0.0_dp)
325 : END DO
326 6 : IF (output_unit > 0) THEN
327 3 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator r_"//ACHAR(idir + 119)
328 : END IF
329 : !
330 : !Initial guess for psi1
331 12 : DO ispin = 1, nspins
332 12 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
333 : !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin))
334 : !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
335 : END DO
336 : !
337 : !DO scf cycle to optimize psi1
338 12 : DO ispin = 1, nspins
339 12 : CALL cp_fm_to_fm(dso_psi0(ispin, idir), h1_psi0(ispin))
340 : END DO
341 : !
342 : !
343 6 : linres_control%lr_triplet = .FALSE. ! we do singlet response
344 6 : linres_control%do_kernel = .TRUE. ! we do uncoupled response
345 6 : linres_control%converged = .FALSE.
346 6 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
347 : !
348 : !
349 : ! copy the response
350 12 : DO ispin = 1, nspins
351 6 : CALL cp_fm_to_fm(psi1(ispin), psi1_dso(ispin, idir))
352 6 : fro = cp_fm_frobenius_norm(psi1(ispin))
353 12 : chk = chk + fro
354 : END DO
355 14 : IF (output_unit > 0) THEN
356 3 : WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
357 : END IF
358 : !
359 : ! Write the result in the restart file
360 : END DO ! idir
361 : END IF
362 : !<<<< debugging only
363 :
364 : !
365 : !
366 : ! print the checksum
367 44 : IF (output_unit > 0) THEN
368 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| response: CheckSum =', chk
369 : END IF
370 : !
371 : !
372 : ! clean up
373 44 : CALL cp_fm_release(psi1)
374 44 : CALL cp_fm_release(h1_psi0)
375 44 : DEALLOCATE (psi0_order)
376 : !
377 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
378 44 : & "PRINT%PROGRAM_RUN_INFO")
379 : !
380 44 : CALL timestop(handle)
381 : !
382 88 : END SUBROUTINE issc_response
383 :
384 : ! **************************************************************************************************
385 : !> \brief ...
386 : !> \param issc_env ...
387 : !> \param qs_env ...
388 : !> \param iatom ...
389 : ! **************************************************************************************************
390 44 : SUBROUTINE issc_issc(issc_env, qs_env, iatom)
391 :
392 : TYPE(issc_env_type) :: issc_env
393 : TYPE(qs_environment_type), POINTER :: qs_env
394 : INTEGER, INTENT(IN) :: iatom
395 :
396 : CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_issc'
397 :
398 : INTEGER :: handle, ispin, ixyz, jatom, jxyz, natom, &
399 : nmo, nspins
400 : LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw
401 : REAL(dp) :: buf, facdso, facfc, facpso, facsd, g, &
402 : issc_dso, issc_fc, issc_pso, issc_sd, &
403 : maxocc
404 : REAL(dp), DIMENSION(3) :: r_i, r_j
405 44 : REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc
406 : TYPE(cell_type), POINTER :: cell
407 44 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc
408 44 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_dso, psi1_efg, psi1_pso
409 : TYPE(cp_fm_type), POINTER :: mo_coeff
410 44 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, &
411 44 : matrix_pso
412 : TYPE(dft_control_type), POINTER :: dft_control
413 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
414 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
415 : TYPE(section_vals_type), POINTER :: issc_section
416 :
417 44 : CALL timeset(routineN, handle)
418 :
419 44 : NULLIFY (cell, dft_control, particle_set, issc, psi1_fc, psi1_efg, psi1_pso)
420 44 : NULLIFY (matrix_efg, matrix_fc, matrix_pso, mos, mo_coeff, fc_psi0)
421 :
422 : CALL get_qs_env(qs_env=qs_env, &
423 : cell=cell, &
424 : dft_control=dft_control, &
425 : particle_set=particle_set, &
426 44 : mos=mos)
427 :
428 44 : gapw = dft_control%qs_control%gapw
429 44 : natom = SIZE(particle_set, 1)
430 44 : nspins = dft_control%nspins
431 :
432 : CALL get_issc_env(issc_env=issc_env, &
433 : matrix_efg=matrix_efg, &
434 : matrix_pso=matrix_pso, &
435 : matrix_fc=matrix_fc, &
436 : matrix_dso=matrix_dso, &
437 : psi1_fc=psi1_fc, &
438 : psi1_efg=psi1_efg, &
439 : psi1_pso=psi1_pso, &
440 : psi1_dso=psi1_dso, &
441 : fc_psi0=fc_psi0, &
442 : issc=issc, &
443 : do_fc=do_fc, &
444 : do_sd=do_sd, &
445 : do_pso=do_pso, &
446 44 : do_dso=do_dso)
447 :
448 44 : g = e_mass/(2.0_dp*p_mass)
449 44 : facfc = hertz*g**2*a_fine**4
450 44 : facpso = hertz*g**2*a_fine**4
451 44 : facsd = hertz*g**2*a_fine**4
452 44 : facdso = hertz*g**2*a_fine**4
453 :
454 : !
455 : !
456 : issc_section => section_vals_get_subs_vals(qs_env%input, &
457 44 : & "PROPERTIES%LINRES%SPINSPIN")
458 : !
459 : ! Initialize
460 92 : DO ispin = 1, nspins
461 48 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, maxocc=maxocc)
462 48 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
463 :
464 292 : DO jatom = 1, natom
465 800 : r_i = particle_set(iatom)%r
466 800 : r_j = particle_set(jatom)%r
467 800 : r_j = pbc(r_i, r_j, cell) + r_i
468 : !
469 : !
470 : !
471 : !write(*,*) 'iatom =',iatom,' r_i=',r_i
472 : !write(*,*) 'jatom =',jatom,' r_j=',r_j
473 : !
474 : ! FC term
475 : !
476 200 : IF (do_fc .AND. iatom /= jatom) THEN
477 : !
478 : ! build the integral for the jatom
479 0 : CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp)
480 0 : CALL build_fermi_contact_matrix(qs_env, matrix_fc, r_j)
481 : CALL cp_dbcsr_sm_fm_multiply(matrix_fc(1)%matrix, mo_coeff, &
482 : fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
483 0 : & alpha=1.0_dp)
484 :
485 0 : CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
486 0 : WRITE (*, *) ' jatom', jatom, 'tr(P*fc)=', buf
487 :
488 0 : CALL cp_fm_trace(fc_psi0(ispin), psi1_fc(ispin), buf)
489 0 : issc_fc = 2.0_dp*2.0_dp*maxocc*facfc*buf
490 0 : issc(1, 1, iatom, jatom, 1) = issc(1, 1, iatom, jatom, 1) + issc_fc
491 0 : issc(2, 2, iatom, jatom, 1) = issc(2, 2, iatom, jatom, 1) + issc_fc
492 0 : issc(3, 3, iatom, jatom, 1) = issc(3, 3, iatom, jatom, 1) + issc_fc
493 : END IF
494 : !
495 : ! SD term
496 : !
497 200 : IF (do_sd .AND. iatom /= jatom) THEN
498 : !
499 : ! build the integral for the jatom
500 0 : CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp)
501 0 : CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp)
502 0 : CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp)
503 0 : CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp)
504 0 : CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp)
505 0 : CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp)
506 0 : CALL build_efg_matrix(qs_env, matrix_efg, r_j)
507 0 : DO ixyz = 1, 6
508 : CALL cp_dbcsr_sm_fm_multiply(matrix_efg(ixyz)%matrix, mo_coeff, &
509 : fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
510 0 : & alpha=1.0_dp, beta=0.0_dp)
511 0 : CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
512 0 : WRITE (*, *) ' jatom', jatom, ixyz, 'tr(P*efg)=', buf
513 0 : DO jxyz = 1, 6
514 0 : CALL cp_fm_trace(fc_psi0(ispin), psi1_efg(ispin, jxyz), buf)
515 0 : issc_sd = 2.0_dp*maxocc*facsd*buf
516 : !issc(ixyz,jxyz,iatom,jatom) = issc_sd
517 : !write(*,*) 'pso_',ixyz,jxyz,' iatom',iatom,'jatom',jatom,issc_pso
518 : END DO
519 : END DO
520 : END IF
521 : !
522 : ! PSO term
523 : !
524 200 : IF (do_pso .AND. iatom /= jatom) THEN
525 : !
526 : ! build the integral for the jatom
527 128 : CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp)
528 128 : CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp)
529 128 : CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp)
530 128 : CALL build_pso_matrix(qs_env, matrix_pso, r_j)
531 512 : DO ixyz = 1, 3
532 : CALL cp_dbcsr_sm_fm_multiply(matrix_pso(ixyz)%matrix, mo_coeff, &
533 : fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
534 384 : & alpha=1.0_dp, beta=0.0_dp)
535 1664 : DO jxyz = 1, 3
536 1152 : CALL cp_fm_trace(fc_psi0(ispin), psi1_pso(ispin, jxyz), buf)
537 1152 : issc_pso = -2.0_dp*maxocc*facpso*buf
538 1536 : issc(ixyz, jxyz, iatom, jatom, 3) = issc(ixyz, jxyz, iatom, jatom, 3) + issc_pso
539 : END DO
540 : END DO
541 : END IF
542 : !
543 : ! DSO term
544 : !
545 : !>>>>> for debugging we compute here the polarizability and NOT the DSO term!
546 248 : IF (do_dso .AND. iatom == natom .AND. jatom == natom) THEN
547 8 : DO ixyz = 1, 3
548 : CALL cp_dbcsr_sm_fm_multiply(matrix_dso(ixyz)%matrix, mo_coeff, &
549 : fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
550 6 : & alpha=1.0_dp, beta=0.0_dp)
551 26 : DO jxyz = 1, 3
552 18 : CALL cp_fm_trace(psi1_dso(ispin, jxyz), fc_psi0(ispin), buf)
553 : ! we save the polarizability for a checksum later on !
554 18 : issc_dso = 2.0_dp*maxocc*buf
555 24 : issc(ixyz, jxyz, iatom, jatom, 4) = issc(ixyz, jxyz, iatom, jatom, 4) + issc_dso
556 : END DO
557 : END DO
558 :
559 : END IF
560 : !
561 : END DO ! jatom
562 : END DO ! ispin
563 : !
564 : !
565 : ! Finalize
566 44 : CALL timestop(handle)
567 : !
568 44 : END SUBROUTINE issc_issc
569 :
570 : ! **************************************************************************************************
571 : !> \brief ...
572 : !> \param issc_env ...
573 : !> \param qs_env ...
574 : ! **************************************************************************************************
575 12 : SUBROUTINE issc_print(issc_env, qs_env)
576 : TYPE(issc_env_type) :: issc_env
577 : TYPE(qs_environment_type), POINTER :: qs_env
578 :
579 : CHARACTER(LEN=2) :: element_symbol_i, element_symbol_j
580 : CHARACTER(LEN=default_string_length) :: name_i, name_j, title
581 : INTEGER :: iatom, jatom, natom, output_unit, &
582 : unit_atoms
583 : LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw
584 : REAL(dp) :: eig(3), issc_iso_dso, issc_iso_fc, &
585 : issc_iso_pso, issc_iso_sd, &
586 : issc_iso_tot, issc_tmp(3, 3)
587 12 : REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc
588 : REAL(dp), EXTERNAL :: DDOT
589 : TYPE(atomic_kind_type), POINTER :: atom_kind_i, atom_kind_j
590 : TYPE(cp_logger_type), POINTER :: logger
591 : TYPE(dft_control_type), POINTER :: dft_control
592 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
593 : TYPE(section_vals_type), POINTER :: issc_section
594 :
595 12 : NULLIFY (logger, particle_set, atom_kind_i, atom_kind_j, dft_control)
596 :
597 24 : logger => cp_get_default_logger()
598 12 : output_unit = cp_logger_get_default_io_unit(logger)
599 :
600 : issc_section => section_vals_get_subs_vals(qs_env%input, &
601 12 : "PROPERTIES%LINRES%SPINSPIN")
602 :
603 : CALL get_issc_env(issc_env=issc_env, &
604 : issc=issc, &
605 : do_fc=do_fc, &
606 : do_sd=do_sd, &
607 : do_pso=do_pso, &
608 12 : do_dso=do_dso)
609 : !
610 : CALL get_qs_env(qs_env=qs_env, &
611 : dft_control=dft_control, &
612 12 : particle_set=particle_set)
613 :
614 12 : natom = SIZE(particle_set, 1)
615 12 : gapw = dft_control%qs_control%gapw
616 :
617 : !
618 12 : IF (output_unit > 0) THEN
619 6 : WRITE (output_unit, '(T2,A,E14.6)') 'ISSC| CheckSum K =', &
620 42 : SQRT(DDOT(SIZE(issc), issc, 1, issc, 1))
621 : END IF
622 : !
623 12 : IF (BTEST(cp_print_key_should_output(logger%iter_info, issc_section, &
624 : "PRINT%K_MATRIX"), cp_p_file)) THEN
625 :
626 : unit_atoms = cp_print_key_unit_nr(logger, issc_section, "PRINT%K_MATRIX", &
627 12 : extension=".data", middle_name="K", log_filename=.FALSE.)
628 :
629 12 : IF (unit_atoms > 0) THEN
630 6 : WRITE (unit_atoms, *)
631 6 : WRITE (unit_atoms, *)
632 6 : WRITE (title, '(A)') "Indirect spin-spin coupling matrix"
633 6 : WRITE (unit_atoms, '(T2,A)') title
634 28 : DO iatom = 1, natom
635 22 : atom_kind_i => particle_set(iatom)%atomic_kind
636 22 : CALL get_atomic_kind(atom_kind_i, name=name_i, element_symbol=element_symbol_i)
637 124 : DO jatom = 1, natom
638 96 : atom_kind_j => particle_set(jatom)%atomic_kind
639 96 : CALL get_atomic_kind(atom_kind_j, name=name_j, element_symbol=element_symbol_j)
640 : !
641 96 : IF (iatom == jatom .AND. .NOT. do_dso) CYCLE
642 : !
643 : !
644 : ! FC
645 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 1)
646 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
647 75 : CALL diamat_all(issc_tmp, eig)
648 75 : issc_iso_fc = (eig(1) + eig(2) + eig(3))/3.0_dp
649 : !
650 : ! SD
651 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 2)
652 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
653 75 : CALL diamat_all(issc_tmp, eig)
654 75 : issc_iso_sd = (eig(1) + eig(2) + eig(3))/3.0_dp
655 : !
656 : ! PSO
657 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 3)
658 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
659 75 : CALL diamat_all(issc_tmp, eig)
660 75 : issc_iso_pso = (eig(1) + eig(2) + eig(3))/3.0_dp
661 : !
662 : ! DSO
663 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 4)
664 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
665 75 : CALL diamat_all(issc_tmp, eig)
666 75 : issc_iso_dso = (eig(1) + eig(2) + eig(3))/3.0_dp
667 : !
668 : ! TOT
669 75 : issc_iso_tot = issc_iso_fc + issc_iso_sd + issc_iso_dso + issc_iso_pso
670 : !
671 : !
672 75 : WRITE (unit_atoms, *)
673 75 : WRITE (unit_atoms, '(T2,2(A,I5,A,2X,A2))') 'Indirect spin-spin coupling between ', &
674 75 : iatom, TRIM(name_i), element_symbol_i, ' and ', &
675 150 : jatom, TRIM(name_j), element_symbol_j
676 : !
677 75 : IF (do_fc) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic FC contribution = ', issc_iso_fc, ' Hz'
678 75 : IF (do_sd) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic SD contribution = ', issc_iso_sd, ' Hz'
679 75 : IF (do_pso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic PSO contribution = ', issc_iso_pso, ' Hz'
680 : !IF(do_dso) WRITE(unit_atoms,'(T1,A,f12.4,A)') ' Isotropic DSO contribution = ',issc_iso_dso,' Hz'
681 75 : IF (do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' !!! POLARIZABILITY (for the moment) = ', issc_iso_dso, ' Hz'
682 97 : IF (.NOT. do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic coupling = ', issc_iso_tot, ' Hz'
683 : END DO
684 : END DO
685 : END IF
686 : CALL cp_print_key_finished_output(unit_atoms, logger, issc_section,&
687 12 : & "PRINT%K_MATRIX")
688 : END IF
689 : !
690 : !
691 12 : END SUBROUTINE issc_print
692 :
693 : ! **************************************************************************************************
694 : !> \brief Initialize the issc environment
695 : !> \param issc_env ...
696 : !> \param qs_env ...
697 : ! **************************************************************************************************
698 12 : SUBROUTINE issc_env_init(issc_env, qs_env)
699 : !
700 : TYPE(issc_env_type) :: issc_env
701 : TYPE(qs_environment_type), POINTER :: qs_env
702 :
703 : CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_env_init'
704 :
705 : INTEGER :: handle, iatom, idir, ini, ir, ispin, m, &
706 : n, n_rep, nao, natom, nspins, &
707 : output_unit
708 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
709 12 : INTEGER, DIMENSION(:), POINTER :: list, row_blk_sizes
710 : LOGICAL :: gapw
711 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
712 : TYPE(cp_fm_type), POINTER :: mo_coeff
713 : TYPE(cp_logger_type), POINTER :: logger
714 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
715 : TYPE(dft_control_type), POINTER :: dft_control
716 : TYPE(linres_control_type), POINTER :: linres_control
717 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
718 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
719 12 : POINTER :: sab_orb
720 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
721 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
722 : TYPE(section_vals_type), POINTER :: issc_section, lr_section
723 :
724 : !
725 :
726 12 : CALL timeset(routineN, handle)
727 :
728 12 : NULLIFY (linres_control)
729 12 : NULLIFY (logger, issc_section)
730 12 : NULLIFY (tmp_fm_struct)
731 12 : NULLIFY (particle_set, qs_kind_set)
732 12 : NULLIFY (sab_orb)
733 :
734 12 : logger => cp_get_default_logger()
735 12 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
736 :
737 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
738 12 : extension=".linresLog")
739 :
740 12 : CALL issc_env_cleanup(issc_env)
741 :
742 12 : IF (output_unit > 0) THEN
743 6 : WRITE (output_unit, "(/,T20,A,/)") "*** Start indirect spin-spin coupling Calculation ***"
744 6 : WRITE (output_unit, "(T10,A,/)") "Inizialization of the ISSC environment"
745 : END IF
746 :
747 : issc_section => section_vals_get_subs_vals(qs_env%input, &
748 12 : & "PROPERTIES%LINRES%SPINSPIN")
749 : !CALL section_vals_val_get(nmr_section,"INTERPOLATE_SHIFT",l_val=nmr_env%interpolate_shift)
750 : !CALL section_vals_val_get(nmr_section,"SHIFT_GAPW_RADIUS",r_val=nmr_env%shift_gapw_radius)
751 :
752 : CALL get_qs_env(qs_env=qs_env, &
753 : dft_control=dft_control, &
754 : linres_control=linres_control, &
755 : mos=mos, &
756 : sab_orb=sab_orb, &
757 : particle_set=particle_set, &
758 : qs_kind_set=qs_kind_set, &
759 12 : dbcsr_dist=dbcsr_dist)
760 : !
761 : !
762 12 : gapw = dft_control%qs_control%gapw
763 12 : nspins = dft_control%nspins
764 12 : natom = SIZE(particle_set, 1)
765 : !
766 : ! check that the psi0 are localized and you have all the centers
767 12 : IF (.NOT. linres_control%localized_psi0) THEN
768 : CALL cp_warn(__LOCATION__, 'To get indirect spin-spin coupling parameters within '// &
769 0 : 'PBC you need to localize zero order orbitals')
770 : END IF
771 : !
772 : !
773 : ! read terms need to be calculated
774 : ! FC
775 12 : CALL section_vals_val_get(issc_section, "DO_FC", l_val=issc_env%do_fc)
776 : ! SD
777 12 : CALL section_vals_val_get(issc_section, "DO_SD", l_val=issc_env%do_sd)
778 : ! PSO
779 12 : CALL section_vals_val_get(issc_section, "DO_PSO", l_val=issc_env%do_pso)
780 : ! DSO
781 12 : CALL section_vals_val_get(issc_section, "DO_DSO", l_val=issc_env%do_dso)
782 : !
783 : !
784 : ! read the list of atoms on which the issc need to be calculated
785 12 : CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", n_rep_val=n_rep)
786 : !
787 : !
788 12 : NULLIFY (issc_env%issc_on_atom_list)
789 12 : n = 0
790 16 : DO ir = 1, n_rep
791 4 : NULLIFY (list)
792 4 : CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", i_rep_val=ir, i_vals=list)
793 16 : IF (ASSOCIATED(list)) THEN
794 4 : CALL reallocate(issc_env%issc_on_atom_list, 1, n + SIZE(list))
795 14 : DO ini = 1, SIZE(list)
796 14 : issc_env%issc_on_atom_list(ini + n) = list(ini)
797 : END DO
798 4 : n = n + SIZE(list)
799 : END IF
800 : END DO
801 : !
802 12 : IF (.NOT. ASSOCIATED(issc_env%issc_on_atom_list)) THEN
803 30 : ALLOCATE (issc_env%issc_on_atom_list(natom))
804 44 : DO iatom = 1, natom
805 44 : issc_env%issc_on_atom_list(iatom) = iatom
806 : END DO
807 : END IF
808 12 : issc_env%issc_natms = SIZE(issc_env%issc_on_atom_list)
809 : !
810 : !
811 : ! Initialize the issc tensor
812 : ALLOCATE (issc_env%issc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4), &
813 84 : issc_env%issc_loc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4))
814 10220 : issc_env%issc(:, :, :, :, :) = 0.0_dp
815 10220 : issc_env%issc_loc(:, :, :, :, :) = 0.0_dp
816 : !
817 : ! allocation
818 : ALLOCATE (issc_env%efg_psi0(nspins, 6), issc_env%pso_psi0(nspins, 3), issc_env%fc_psi0(nspins), &
819 : issc_env%psi1_efg(nspins, 6), issc_env%psi1_pso(nspins, 3), issc_env%psi1_fc(nspins), &
820 796 : issc_env%dso_psi0(nspins, 3), issc_env%psi1_dso(nspins, 3))
821 26 : DO ispin = 1, nspins
822 14 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
823 14 : CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
824 :
825 14 : NULLIFY (tmp_fm_struct)
826 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
827 : ncol_global=m, &
828 14 : context=mo_coeff%matrix_struct%context)
829 98 : DO idir = 1, 6
830 84 : CALL cp_fm_create(issc_env%psi1_efg(ispin, idir), tmp_fm_struct)
831 98 : CALL cp_fm_create(issc_env%efg_psi0(ispin, idir), tmp_fm_struct)
832 : END DO
833 56 : DO idir = 1, 3
834 42 : CALL cp_fm_create(issc_env%psi1_pso(ispin, idir), tmp_fm_struct)
835 42 : CALL cp_fm_create(issc_env%pso_psi0(ispin, idir), tmp_fm_struct)
836 42 : CALL cp_fm_create(issc_env%psi1_dso(ispin, idir), tmp_fm_struct)
837 56 : CALL cp_fm_create(issc_env%dso_psi0(ispin, idir), tmp_fm_struct)
838 : END DO
839 14 : CALL cp_fm_create(issc_env%psi1_fc(ispin), tmp_fm_struct)
840 14 : CALL cp_fm_create(issc_env%fc_psi0(ispin), tmp_fm_struct)
841 40 : CALL cp_fm_struct_release(tmp_fm_struct)
842 : END DO
843 : !
844 : ! prepare for allocation
845 36 : ALLOCATE (first_sgf(natom))
846 24 : ALLOCATE (last_sgf(natom))
847 : CALL get_particle_set(particle_set, qs_kind_set, &
848 : first_sgf=first_sgf, &
849 12 : last_sgf=last_sgf)
850 24 : ALLOCATE (row_blk_sizes(natom))
851 12 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
852 12 : DEALLOCATE (first_sgf)
853 12 : DEALLOCATE (last_sgf)
854 :
855 : !
856 : ! efg, pso and fc operators
857 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_efg, 6)
858 12 : ALLOCATE (issc_env%matrix_efg(1)%matrix)
859 : CALL dbcsr_create(matrix=issc_env%matrix_efg(1)%matrix, &
860 : name="efg (3xx-rr)/3", &
861 : dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
862 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
863 12 : mutable_work=.TRUE.)
864 12 : CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_efg(1)%matrix, sab_orb)
865 :
866 : ALLOCATE (issc_env%matrix_efg(2)%matrix, &
867 : issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(4)%matrix, &
868 12 : issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(6)%matrix)
869 : CALL dbcsr_copy(issc_env%matrix_efg(2)%matrix, issc_env%matrix_efg(1)%matrix, &
870 12 : 'efg xy')
871 : CALL dbcsr_copy(issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(1)%matrix, &
872 12 : 'efg xz')
873 : CALL dbcsr_copy(issc_env%matrix_efg(4)%matrix, issc_env%matrix_efg(1)%matrix, &
874 12 : 'efg (3yy-rr)/3')
875 : CALL dbcsr_copy(issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(1)%matrix, &
876 12 : 'efg yz')
877 : CALL dbcsr_copy(issc_env%matrix_efg(6)%matrix, issc_env%matrix_efg(1)%matrix, &
878 12 : 'efg (3zz-rr)/3')
879 :
880 12 : CALL dbcsr_set(issc_env%matrix_efg(1)%matrix, 0.0_dp)
881 12 : CALL dbcsr_set(issc_env%matrix_efg(2)%matrix, 0.0_dp)
882 12 : CALL dbcsr_set(issc_env%matrix_efg(3)%matrix, 0.0_dp)
883 12 : CALL dbcsr_set(issc_env%matrix_efg(4)%matrix, 0.0_dp)
884 12 : CALL dbcsr_set(issc_env%matrix_efg(5)%matrix, 0.0_dp)
885 12 : CALL dbcsr_set(issc_env%matrix_efg(6)%matrix, 0.0_dp)
886 : !
887 : ! PSO
888 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_pso, 3)
889 12 : ALLOCATE (issc_env%matrix_pso(1)%matrix)
890 : CALL dbcsr_create(matrix=issc_env%matrix_pso(1)%matrix, &
891 : name="pso x", &
892 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
893 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
894 12 : mutable_work=.TRUE.)
895 12 : CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_pso(1)%matrix, sab_orb)
896 :
897 12 : ALLOCATE (issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(3)%matrix)
898 : CALL dbcsr_copy(issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(1)%matrix, &
899 12 : 'pso y')
900 : CALL dbcsr_copy(issc_env%matrix_pso(3)%matrix, issc_env%matrix_pso(1)%matrix, &
901 12 : 'pso z')
902 12 : CALL dbcsr_set(issc_env%matrix_pso(1)%matrix, 0.0_dp)
903 12 : CALL dbcsr_set(issc_env%matrix_pso(2)%matrix, 0.0_dp)
904 12 : CALL dbcsr_set(issc_env%matrix_pso(3)%matrix, 0.0_dp)
905 : !
906 : ! DSO
907 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_dso, 3)
908 12 : ALLOCATE (issc_env%matrix_dso(1)%matrix, issc_env%matrix_dso(2)%matrix, issc_env%matrix_dso(3)%matrix)
909 : CALL dbcsr_copy(issc_env%matrix_dso(1)%matrix, issc_env%matrix_efg(1)%matrix, &
910 12 : 'dso x')
911 : CALL dbcsr_copy(issc_env%matrix_dso(2)%matrix, issc_env%matrix_efg(1)%matrix, &
912 12 : 'dso y')
913 : CALL dbcsr_copy(issc_env%matrix_dso(3)%matrix, issc_env%matrix_efg(1)%matrix, &
914 12 : 'dso z')
915 12 : CALL dbcsr_set(issc_env%matrix_dso(1)%matrix, 0.0_dp)
916 12 : CALL dbcsr_set(issc_env%matrix_dso(2)%matrix, 0.0_dp)
917 12 : CALL dbcsr_set(issc_env%matrix_dso(3)%matrix, 0.0_dp)
918 : !
919 : ! FC
920 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_fc, 1)
921 12 : ALLOCATE (issc_env%matrix_fc(1)%matrix)
922 : CALL dbcsr_copy(issc_env%matrix_fc(1)%matrix, issc_env%matrix_efg(1)%matrix, &
923 12 : 'fc')
924 12 : CALL dbcsr_set(issc_env%matrix_fc(1)%matrix, 0.0_dp)
925 :
926 12 : DEALLOCATE (row_blk_sizes)
927 : !
928 : ! Conversion factors
929 12 : IF (output_unit > 0) THEN
930 : WRITE (output_unit, "(T2,A,T60,I4,A)")&
931 6 : & "ISSC| spin-spin coupling computed for ", issc_env%issc_natms, ' atoms'
932 : END IF
933 :
934 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
935 12 : & "PRINT%PROGRAM_RUN_INFO")
936 :
937 12 : CALL timestop(handle)
938 :
939 24 : END SUBROUTINE issc_env_init
940 :
941 : ! **************************************************************************************************
942 : !> \brief Deallocate the issc environment
943 : !> \param issc_env ...
944 : !> \par History
945 : ! **************************************************************************************************
946 24 : SUBROUTINE issc_env_cleanup(issc_env)
947 :
948 : TYPE(issc_env_type), INTENT(INOUT) :: issc_env
949 :
950 24 : IF (ASSOCIATED(issc_env%issc_on_atom_list)) THEN
951 12 : DEALLOCATE (issc_env%issc_on_atom_list)
952 : END IF
953 24 : IF (ASSOCIATED(issc_env%issc)) THEN
954 12 : DEALLOCATE (issc_env%issc)
955 : END IF
956 24 : IF (ASSOCIATED(issc_env%issc_loc)) THEN
957 12 : DEALLOCATE (issc_env%issc_loc)
958 : END IF
959 : !
960 : !efg_psi0
961 24 : CALL cp_fm_release(issc_env%efg_psi0)
962 : !
963 : !pso_psi0
964 24 : CALL cp_fm_release(issc_env%pso_psi0)
965 : !
966 : !dso_psi0
967 24 : CALL cp_fm_release(issc_env%dso_psi0)
968 : !
969 : !fc_psi0
970 24 : CALL cp_fm_release(issc_env%fc_psi0)
971 : !
972 : !psi1_efg
973 24 : CALL cp_fm_release(issc_env%psi1_efg)
974 : !
975 : !psi1_pso
976 24 : CALL cp_fm_release(issc_env%psi1_pso)
977 : !
978 : !psi1_dso
979 24 : CALL cp_fm_release(issc_env%psi1_dso)
980 : !
981 : !psi1_fc
982 24 : CALL cp_fm_release(issc_env%psi1_fc)
983 : !
984 : !matrix_efg
985 24 : IF (ASSOCIATED(issc_env%matrix_efg)) THEN
986 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_efg)
987 : END IF
988 : !
989 : !matrix_pso
990 24 : IF (ASSOCIATED(issc_env%matrix_pso)) THEN
991 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_pso)
992 : END IF
993 : !
994 : !matrix_dso
995 24 : IF (ASSOCIATED(issc_env%matrix_dso)) THEN
996 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_dso)
997 : END IF
998 : !
999 : !matrix_fc
1000 24 : IF (ASSOCIATED(issc_env%matrix_fc)) THEN
1001 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_fc)
1002 : END IF
1003 :
1004 24 : END SUBROUTINE issc_env_cleanup
1005 :
1006 : END MODULE qs_linres_issc_utils
|
