(git:34ef472)
qs_linres_issc_utils.F
Go to the documentation of this file.
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 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
26  USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
27  dbcsr_allocate_matrix_set,&
28  dbcsr_deallocate_matrix_set
29  USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,&
30  cp_fm_trace
31  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32  cp_fm_struct_release,&
33  cp_fm_struct_type
34  USE cp_fm_types, ONLY: cp_fm_create,&
35  cp_fm_get_info,&
36  cp_fm_release,&
37  cp_fm_set_all,&
38  cp_fm_to_fm,&
39  cp_fm_type
40  USE cp_log_handling, ONLY: cp_get_default_logger,&
41  cp_logger_get_default_io_unit,&
42  cp_logger_type
43  USE cp_output_handling, ONLY: cp_p_file,&
44  cp_print_key_finished_output,&
45  cp_print_key_should_output,&
46  cp_print_key_unit_nr
47  USE dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
48  dbcsr_copy,&
49  dbcsr_create,&
50  dbcsr_distribution_type,&
51  dbcsr_p_type,&
52  dbcsr_set,&
53  dbcsr_type_antisymmetric,&
54  dbcsr_type_symmetric
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  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
115  TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc
116  TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dso_psi0, efg_psi0, psi1_dso, psi1_efg, &
117  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  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  CALL timeset(routinen, handle)
128  !
129  NULLIFY (dft_control, linres_control, lr_section, issc_section)
130  NULLIFY (logger, mpools, mo_coeff, para_env)
131  NULLIFY (tmp_fm_struct, psi1_fc, psi1_efg, psi1_pso, pso_psi0, fc_psi0, efg_psi0)
132 
133  logger => cp_get_default_logger()
134  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
135  issc_section => section_vals_get_subs_vals(qs_env%input, &
136  "PROPERTIES%LINRES%SPINSPIN")
137 
138  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
139  extension=".linresLog")
140  IF (output_unit > 0) THEN
141  WRITE (unit=output_unit, fmt="(T10,A,/)") &
142  "*** 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  para_env=para_env)
151 
152  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  do_dso=do_dso)
168  !
169  ! allocate the vectors
170  ALLOCATE (psi0_order(nspins))
171  ALLOCATE (psi1(nspins), h1_psi0(nspins))
172  DO ispin = 1, nspins
173  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
174  psi0_order(ispin) = mo_coeff
175  CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
176  NULLIFY (tmp_fm_struct)
177  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
178  ncol_global=nmo, &
179  context=mo_coeff%matrix_struct%context)
180  CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
181  CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
182  CALL cp_fm_struct_release(tmp_fm_struct)
183  END DO
184  chk = 0.0_dp
185  should_stop = .false.
186  !
187  ! operator efg
188  IF (do_sd) THEN
189  ijdir = 0
190  DO idir = 1, 3
191  DO jdir = idir, 3
192  ijdir = ijdir + 1
193  DO ispin = 1, nspins
194  CALL cp_fm_set_all(psi1_efg(ispin, ijdir), 0.0_dp)
195  END DO
196  IF (output_unit > 0) THEN
197  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  DO ispin = 1, nspins
202  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
203  !CALL cp_fm_to_fm(p_psi0(ispin,ijdir)%matrix, psi1(ispin))
204  !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
205  END DO
206  !
207  !DO scf cycle to optimize psi1
208  DO ispin = 1, nspins
209  CALL cp_fm_to_fm(efg_psi0(ispin, ijdir), h1_psi0(ispin))
210  END DO
211  !
212  !
213  linres_control%lr_triplet = .false.
214  linres_control%do_kernel = .false.
215  linres_control%converged = .false.
216  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
217  !
218  !
219  ! copy the response
220  DO ispin = 1, nspins
221  CALL cp_fm_to_fm(psi1(ispin), psi1_efg(ispin, ijdir))
222  fro = cp_fm_frobenius_norm(psi1(ispin))
223  chk = chk + fro
224  END DO
225  !
226  ! print response functions
227  !IF(BTEST(cp_print_key_should_output(logger%iter_info,issc_section,&
228  ! & "PRINT%RESPONSE_FUNCTION_CUBES"),cp_p_file)) THEN
229  ! ncubes = SIZE(list_cubes,1)
230  ! print_key => section_vals_get_subs_vals(issc_section,"PRINT%RESPONSE_FUNCTION_CUBES")
231  ! DO ispin = 1,nspins
232  ! CALL qs_print_cubes(qs_env,psi1(ispin),ncubes,list_cubes,&
233  ! centers_set(ispin)%array,print_key,'psi1_efg',&
234  ! idir=ijdir,ispin=ispin)
235  ! ENDDO ! ispin
236  !ENDIF ! print response functions
237  !
238  !
239  IF (output_unit > 0) THEN
240  WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
241  END IF
242  !
243  ! Write the result in the restart file
244  END DO ! jdir
245  END DO ! idir
246  END IF
247  !
248  ! operator pso
249  IF (do_pso) THEN
250  DO idir = 1, 3
251  DO ispin = 1, nspins
252  CALL cp_fm_set_all(psi1_pso(ispin, idir), 0.0_dp)
253  END DO
254  IF (output_unit > 0) THEN
255  WRITE (output_unit, "(T10,A)") "Response to the perturbation operator pso_"//achar(idir + 119)
256  END IF
257  !
258  !Initial guess for psi1
259  DO ispin = 1, nspins
260  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
261  !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin))
262  !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
263  END DO
264  !
265  !DO scf cycle to optimize psi1
266  DO ispin = 1, nspins
267  CALL cp_fm_to_fm(pso_psi0(ispin, idir), h1_psi0(ispin))
268  END DO
269  !
270  !
271  linres_control%lr_triplet = .false. ! we do singlet response
272  linres_control%do_kernel = .false. ! we do uncoupled response
273  linres_control%converged = .false.
274  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
275  !
276  !
277  ! copy the response
278  DO ispin = 1, nspins
279  CALL cp_fm_to_fm(psi1(ispin), psi1_pso(ispin, idir))
280  fro = cp_fm_frobenius_norm(psi1(ispin))
281  chk = chk + fro
282  END DO
283  !
284  ! print response functions
285  !IF(BTEST(cp_print_key_should_output(logger%iter_info,issc_section,&
286  ! & "PRINT%RESPONSE_FUNCTION_CUBES"),cp_p_file)) THEN
287  ! ncubes = SIZE(list_cubes,1)
288  ! print_key => section_vals_get_subs_vals(issc_section,"PRINT%RESPONSE_FUNCTION_CUBES")
289  ! DO ispin = 1,nspins
290  ! CALL qs_print_cubes(qs_env,psi1(ispin),ncubes,list_cubes,&
291  ! centers_set(ispin)%array,print_key,'psi1_pso',&
292  ! idir=idir,ispin=ispin)
293  ! ENDDO ! ispin
294  !ENDIF ! print response functions
295  !
296  !
297  IF (output_unit > 0) THEN
298  WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
299  END IF
300  !
301  ! Write the result in the restart file
302  END DO ! idir
303  END IF
304  !
305  ! operator fc
306  IF (do_fc) THEN
307  DO ispin = 1, nspins
308  CALL cp_fm_set_all(psi1_fc(ispin), 0.0_dp)
309  END DO
310  IF (output_unit > 0) THEN
311  WRITE (output_unit, "(T10,A)") "Response to the perturbation operator fc"
312  END IF
313  !
314  !Initial guess for psi1
315  DO ispin = 1, nspins
316  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
317  !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin))
318  !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
319  END DO
320  !
321  !DO scf cycle to optimize psi1
322  DO ispin = 1, nspins
323  CALL cp_fm_to_fm(fc_psi0(ispin), h1_psi0(ispin))
324  END DO
325  !
326  !
327  linres_control%lr_triplet = .true. ! we do triplet response
328  linres_control%do_kernel = .true. ! we do coupled response
329  linres_control%converged = .false.
330  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
331  !
332  !
333  ! copy the response
334  DO ispin = 1, nspins
335  CALL cp_fm_to_fm(psi1(ispin), psi1_fc(ispin))
336  fro = cp_fm_frobenius_norm(psi1(ispin))
337  chk = chk + fro
338  END DO
339  !
340  ! print response functions
341  !IF(BTEST(cp_print_key_should_output(logger%iter_info,issc_section,&
342  ! & "PRINT%RESPONSE_FUNCTION_CUBES"),cp_p_file)) THEN
343  ! ncubes = SIZE(list_cubes,1)
344  ! print_key => section_vals_get_subs_vals(issc_section,"PRINT%RESPONSE_FUNCTION_CUBES")
345  ! DO ispin = 1,nspins
346  ! CALL qs_print_cubes(qs_env,psi1(ispin),ncubes,list_cubes,&
347  ! centers_set(ispin)%array,print_key,'psi1_pso',&
348  ! idir=idir,ispin=ispin)
349  ! ENDDO ! ispin
350  !ENDIF ! print response functions
351  !
352  !
353  IF (output_unit > 0) THEN
354  WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
355  END IF
356  !
357  ! Write the result in the restart file
358  END IF
359 
360  !>>>> debugging only
361  !
362  ! here we have the operator r and compute the polarizability for debugging the kernel only
363  IF (do_dso) THEN
364  DO idir = 1, 3
365  DO ispin = 1, nspins
366  CALL cp_fm_set_all(psi1_dso(ispin, idir), 0.0_dp)
367  END DO
368  IF (output_unit > 0) THEN
369  WRITE (output_unit, "(T10,A)") "Response to the perturbation operator r_"//achar(idir + 119)
370  END IF
371  !
372  !Initial guess for psi1
373  DO ispin = 1, nspins
374  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
375  !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin))
376  !CALL cp_fm_scale(-1.0_dp,psi1(ispin))
377  END DO
378  !
379  !DO scf cycle to optimize psi1
380  DO ispin = 1, nspins
381  CALL cp_fm_to_fm(dso_psi0(ispin, idir), h1_psi0(ispin))
382  END DO
383  !
384  !
385  linres_control%lr_triplet = .false. ! we do singlet response
386  linres_control%do_kernel = .true. ! we do uncoupled response
387  linres_control%converged = .false.
388  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
389  !
390  !
391  ! copy the response
392  DO ispin = 1, nspins
393  CALL cp_fm_to_fm(psi1(ispin), psi1_dso(ispin, idir))
394  fro = cp_fm_frobenius_norm(psi1(ispin))
395  chk = chk + fro
396  END DO
397  !
398  ! print response functions
399  !IF(BTEST(cp_print_key_should_output(logger%iter_info,issc_section,&
400  ! & "PRINT%RESPONSE_FUNCTION_CUBES"),cp_p_file)) THEN
401  ! ncubes = SIZE(list_cubes,1)
402  ! print_key => section_vals_get_subs_vals(issc_section,"PRINT%RESPONSE_FUNCTION_CUBES")
403  ! DO ispin = 1,nspins
404  ! CALL qs_print_cubes(qs_env,psi1(ispin),ncubes,list_cubes,&
405  ! centers_set(ispin)%array,print_key,'psi1_pso',&
406  ! idir=idir,ispin=ispin)
407  ! ENDDO ! ispin
408  !ENDIF ! print response functions
409  !
410  !
411  IF (output_unit > 0) THEN
412  WRITE (output_unit, "(T10,A)") "Write the resulting psi1 in restart file... not implemented yet"
413  END IF
414  !
415  ! Write the result in the restart file
416  END DO ! idir
417  END IF
418  !<<<< debugging only
419 
420  !
421  !
422  ! print the checksum
423  IF (output_unit > 0) THEN
424  WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| response: CheckSum =', chk
425  END IF
426  !
427  !
428  ! clean up
429  CALL cp_fm_release(psi1)
430  CALL cp_fm_release(h1_psi0)
431  DEALLOCATE (psi0_order)
432  !
433  CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
434  & "PRINT%PROGRAM_RUN_INFO")
435  !
436  CALL timestop(handle)
437  !
438  END SUBROUTINE issc_response
439 
440 ! **************************************************************************************************
441 !> \brief ...
442 !> \param issc_env ...
443 !> \param qs_env ...
444 !> \param iatom ...
445 ! **************************************************************************************************
446  SUBROUTINE issc_issc(issc_env, qs_env, iatom)
447 
448  TYPE(issc_env_type) :: issc_env
449  TYPE(qs_environment_type), POINTER :: qs_env
450  INTEGER, INTENT(IN) :: iatom
451 
452  CHARACTER(LEN=*), PARAMETER :: routinen = 'issc_issc'
453 
454  INTEGER :: handle, ispin, ixyz, jatom, jxyz, natom, &
455  nmo, nspins
456  LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw
457  REAL(dp) :: buf, facdso, facfc, facpso, facsd, g, &
458  issc_dso, issc_fc, issc_pso, issc_sd, &
459  maxocc
460  REAL(dp), DIMENSION(3) :: r_i, r_j
461  REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc
462  TYPE(cell_type), POINTER :: cell
463  TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc
464  TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_dso, psi1_efg, psi1_pso
465  TYPE(cp_fm_type), POINTER :: mo_coeff
466  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, &
467  matrix_pso
468  TYPE(dft_control_type), POINTER :: dft_control
469  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
470  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
471  TYPE(section_vals_type), POINTER :: issc_section
472 
473  CALL timeset(routinen, handle)
474 
475  NULLIFY (cell, dft_control, particle_set, issc, psi1_fc, psi1_efg, psi1_pso)
476  NULLIFY (matrix_efg, matrix_fc, matrix_pso, mos, mo_coeff, fc_psi0)
477 
478  CALL get_qs_env(qs_env=qs_env, &
479  cell=cell, &
480  dft_control=dft_control, &
481  particle_set=particle_set, &
482  mos=mos)
483 
484  gapw = dft_control%qs_control%gapw
485  natom = SIZE(particle_set, 1)
486  nspins = dft_control%nspins
487 
488  CALL get_issc_env(issc_env=issc_env, &
489  matrix_efg=matrix_efg, &
490  matrix_pso=matrix_pso, &
491  matrix_fc=matrix_fc, &
492  matrix_dso=matrix_dso, &
493  psi1_fc=psi1_fc, &
494  psi1_efg=psi1_efg, &
495  psi1_pso=psi1_pso, &
496  psi1_dso=psi1_dso, &
497  fc_psi0=fc_psi0, &
498  issc=issc, &
499  do_fc=do_fc, &
500  do_sd=do_sd, &
501  do_pso=do_pso, &
502  do_dso=do_dso)
503 
504  g = e_mass/(2.0_dp*p_mass)
505  facfc = hertz*g**2*a_fine**4
506  facpso = hertz*g**2*a_fine**4
507  facsd = hertz*g**2*a_fine**4
508  facdso = hertz*g**2*a_fine**4
509 
510  !
511  !
512  issc_section => section_vals_get_subs_vals(qs_env%input, &
513  & "PROPERTIES%LINRES%SPINSPIN")
514  !
515  ! Initialize
516  DO ispin = 1, nspins
517  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, maxocc=maxocc)
518  CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
519 
520  DO jatom = 1, natom
521  r_i = particle_set(iatom)%r
522  r_j = particle_set(jatom)%r
523  r_j = pbc(r_i, r_j, cell) + r_i
524  !
525  !
526  !
527  !write(*,*) 'iatom =',iatom,' r_i=',r_i
528  !write(*,*) 'jatom =',jatom,' r_j=',r_j
529  !
530  ! FC term
531  !
532  IF (do_fc .AND. iatom .NE. jatom) THEN
533  !
534  ! build the integral for the jatom
535  CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp)
536  CALL build_fermi_contact_matrix(qs_env, matrix_fc, r_j)
537  CALL cp_dbcsr_sm_fm_multiply(matrix_fc(1)%matrix, mo_coeff, &
538  fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
539  & alpha=1.0_dp)
540 
541  CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
542  WRITE (*, *) ' jatom', jatom, 'tr(P*fc)=', buf
543 
544  CALL cp_fm_trace(fc_psi0(ispin), psi1_fc(ispin), buf)
545  issc_fc = 2.0_dp*2.0_dp*maxocc*facfc*buf
546  issc(1, 1, iatom, jatom, 1) = issc(1, 1, iatom, jatom, 1) + issc_fc
547  issc(2, 2, iatom, jatom, 1) = issc(2, 2, iatom, jatom, 1) + issc_fc
548  issc(3, 3, iatom, jatom, 1) = issc(3, 3, iatom, jatom, 1) + issc_fc
549  END IF
550  !
551  ! SD term
552  !
553  IF (do_sd .AND. iatom .NE. jatom) THEN
554  !
555  ! build the integral for the jatom
556  CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp)
557  CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp)
558  CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp)
559  CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp)
560  CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp)
561  CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp)
562  CALL build_efg_matrix(qs_env, matrix_efg, r_j)
563  DO ixyz = 1, 6
564  CALL cp_dbcsr_sm_fm_multiply(matrix_efg(ixyz)%matrix, mo_coeff, &
565  fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
566  & alpha=1.0_dp, beta=0.0_dp)
567  CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
568  WRITE (*, *) ' jatom', jatom, ixyz, 'tr(P*efg)=', buf
569  DO jxyz = 1, 6
570  CALL cp_fm_trace(fc_psi0(ispin), psi1_efg(ispin, jxyz), buf)
571  issc_sd = 2.0_dp*maxocc*facsd*buf
572  !issc(ixyz,jxyz,iatom,jatom) = issc_sd
573  !write(*,*) 'pso_',ixyz,jxyz,' iatom',iatom,'jatom',jatom,issc_pso
574  END DO
575  END DO
576  END IF
577  !
578  ! PSO term
579  !
580  IF (do_pso .AND. iatom .NE. jatom) THEN
581  !
582  ! build the integral for the jatom
583  CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp)
584  CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp)
585  CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp)
586  CALL build_pso_matrix(qs_env, matrix_pso, r_j)
587  DO ixyz = 1, 3
588  CALL cp_dbcsr_sm_fm_multiply(matrix_pso(ixyz)%matrix, mo_coeff, &
589  fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
590  & alpha=1.0_dp, beta=0.0_dp)
591  DO jxyz = 1, 3
592  CALL cp_fm_trace(fc_psi0(ispin), psi1_pso(ispin, jxyz), buf)
593  issc_pso = -2.0_dp*maxocc*facpso*buf
594  issc(ixyz, jxyz, iatom, jatom, 3) = issc(ixyz, jxyz, iatom, jatom, 3) + issc_pso
595  END DO
596  END DO
597  END IF
598  !
599  ! DSO term
600  !
601  !>>>>> for debugging we compute here the polarizability and NOT the DSO term!
602  IF (do_dso .AND. iatom .EQ. natom .AND. jatom .EQ. natom) THEN
603  !
604  ! build the integral for the jatom
605  !CALL dbcsr_set(matrix_dso(1)%matrix,0.0_dp)
606  !CALL dbcsr_set(matrix_dso(2)%matrix,0.0_dp)
607  !CALL dbcsr_set(matrix_dso(3)%matrix,0.0_dp)
608  !CALL dbcsr_set(matrix_dso(4)%matrix,0.0_dp)
609  !CALL dbcsr_set(matrix_dso(5)%matrix,0.0_dp)
610  !CALL dbcsr_set(matrix_dso(6)%matrix,0.0_dp)
611  !CALL build_dso_matrix(qs_env,matrix_dso,r_i,r_j)
612  !DO ixyz = 1,6
613  ! CALL cp_sm_fm_multiply(matrix_dso(ixyz)%matrix,mo_coeff,&
614  ! & fc_psi0(ispin),ncol=nmo,& ! fc_psi0 a buffer
615  ! & alpha=1.0_dp,beta=0.0_dp)
616  ! CALL cp_fm_trace(fc_psi0(ispin),mo_coeff,buf)
617  ! issc_dso = 2.0_dp * maxocc * facdso * buf
618  ! issc(ixyz,jxyz,iatom,jatom,4) = issc_dso
619  !ENDDO
620  !CALL dbcsr_set(matrix_dso(1)%matrix,0.0_dp)
621  !CALL dbcsr_set(matrix_dso(2)%matrix,0.0_dp)
622  !CALL dbcsr_set(matrix_dso(3)%matrix,0.0_dp)
623  !CALL rRc_xyz_ao(matrix_dso,qs_env,(/0.0_dp,0.0_dp,0.0_dp/),1)
624  DO ixyz = 1, 3
625  CALL cp_dbcsr_sm_fm_multiply(matrix_dso(ixyz)%matrix, mo_coeff, &
626  fc_psi0(ispin), ncol=nmo,& ! fc_psi0 a buffer
627  & alpha=1.0_dp, beta=0.0_dp)
628  DO jxyz = 1, 3
629  CALL cp_fm_trace(psi1_dso(ispin, jxyz), fc_psi0(ispin), buf)
630  ! we save the polarizability for a checksum later on !
631  issc_dso = 2.0_dp*maxocc*buf
632  issc(ixyz, jxyz, iatom, jatom, 4) = issc(ixyz, jxyz, iatom, jatom, 4) + issc_dso
633  END DO
634  END DO
635 
636  END IF
637  !
638  END DO ! jatom
639  END DO ! ispin
640  !
641  !
642  ! Finalize
643  CALL timestop(handle)
644  !
645  END SUBROUTINE issc_issc
646 
647 ! **************************************************************************************************
648 !> \brief ...
649 !> \param issc_env ...
650 !> \param qs_env ...
651 ! **************************************************************************************************
652  SUBROUTINE issc_print(issc_env, qs_env)
653  TYPE(issc_env_type) :: issc_env
654  TYPE(qs_environment_type), POINTER :: qs_env
655 
656  CHARACTER(LEN=2) :: element_symbol_i, element_symbol_j
657  CHARACTER(LEN=default_string_length) :: name_i, name_j, title
658  INTEGER :: iatom, jatom, natom, output_unit, &
659  unit_atoms
660  LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw
661  REAL(dp) :: eig(3), issc_iso_dso, issc_iso_fc, &
662  issc_iso_pso, issc_iso_sd, &
663  issc_iso_tot, issc_tmp(3, 3)
664  REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc
665  REAL(dp), EXTERNAL :: ddot
666  TYPE(atomic_kind_type), POINTER :: atom_kind_i, atom_kind_j
667  TYPE(cp_logger_type), POINTER :: logger
668  TYPE(dft_control_type), POINTER :: dft_control
669  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
670  TYPE(section_vals_type), POINTER :: issc_section
671 
672  NULLIFY (logger, particle_set, atom_kind_i, atom_kind_j, dft_control)
673 
674  logger => cp_get_default_logger()
675  output_unit = cp_logger_get_default_io_unit(logger)
676 
677  issc_section => section_vals_get_subs_vals(qs_env%input, &
678  "PROPERTIES%LINRES%SPINSPIN")
679 
680  CALL get_issc_env(issc_env=issc_env, &
681  issc=issc, &
682  do_fc=do_fc, &
683  do_sd=do_sd, &
684  do_pso=do_pso, &
685  do_dso=do_dso)
686  !
687  CALL get_qs_env(qs_env=qs_env, &
688  dft_control=dft_control, &
689  particle_set=particle_set)
690 
691  natom = SIZE(particle_set, 1)
692  gapw = dft_control%qs_control%gapw
693 
694  !
695  IF (output_unit > 0) THEN
696  WRITE (output_unit, '(T2,A,E14.6)') 'ISSC| CheckSum K =', &
697  sqrt(ddot(SIZE(issc), issc, 1, issc, 1))
698  END IF
699  !
700  IF (btest(cp_print_key_should_output(logger%iter_info, issc_section, &
701  "PRINT%K_MATRIX"), cp_p_file)) THEN
702 
703  unit_atoms = cp_print_key_unit_nr(logger, issc_section, "PRINT%K_MATRIX", &
704  extension=".data", middle_name="K", log_filename=.false.)
705 
706  IF (unit_atoms > 0) THEN
707  WRITE (unit_atoms, *)
708  WRITE (unit_atoms, *)
709  WRITE (title, '(A)') "Indirect spin-spin coupling matrix"
710  WRITE (unit_atoms, '(T2,A)') title
711  DO iatom = 1, natom
712  atom_kind_i => particle_set(iatom)%atomic_kind
713  CALL get_atomic_kind(atom_kind_i, name=name_i, element_symbol=element_symbol_i)
714  DO jatom = 1, natom
715  atom_kind_j => particle_set(jatom)%atomic_kind
716  CALL get_atomic_kind(atom_kind_j, name=name_j, element_symbol=element_symbol_j)
717  !
718  IF (iatom .EQ. jatom .AND. .NOT. do_dso) cycle
719  !
720  !
721  ! FC
722  issc_tmp(:, :) = issc(:, :, iatom, jatom, 1)
723  issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + transpose(issc_tmp(:, :)))
724  CALL diamat_all(issc_tmp, eig)
725  issc_iso_fc = (eig(1) + eig(2) + eig(3))/3.0_dp
726  !
727  ! SD
728  issc_tmp(:, :) = issc(:, :, iatom, jatom, 2)
729  issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + transpose(issc_tmp(:, :)))
730  CALL diamat_all(issc_tmp, eig)
731  issc_iso_sd = (eig(1) + eig(2) + eig(3))/3.0_dp
732  !
733  ! PSO
734  issc_tmp(:, :) = issc(:, :, iatom, jatom, 3)
735  issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + transpose(issc_tmp(:, :)))
736  CALL diamat_all(issc_tmp, eig)
737  issc_iso_pso = (eig(1) + eig(2) + eig(3))/3.0_dp
738  !
739  ! DSO
740  issc_tmp(:, :) = issc(:, :, iatom, jatom, 4)
741  issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + transpose(issc_tmp(:, :)))
742  CALL diamat_all(issc_tmp, eig)
743  issc_iso_dso = (eig(1) + eig(2) + eig(3))/3.0_dp
744  !
745  ! TOT
746  issc_iso_tot = issc_iso_fc + issc_iso_sd + issc_iso_dso + issc_iso_pso
747  !
748  !
749  WRITE (unit_atoms, *)
750  WRITE (unit_atoms, '(T2,2(A,I5,A,2X,A2))') 'Indirect spin-spin coupling between ', &
751  iatom, trim(name_i), element_symbol_i, ' and ', &
752  jatom, trim(name_j), element_symbol_j
753  !
754  IF (do_fc) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic FC contribution = ', issc_iso_fc, ' Hz'
755  IF (do_sd) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic SD contribution = ', issc_iso_sd, ' Hz'
756  IF (do_pso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic PSO contribution = ', issc_iso_pso, ' Hz'
757  !IF(do_dso) WRITE(unit_atoms,'(T1,A,f12.4,A)') ' Isotropic DSO contribution = ',issc_iso_dso,' Hz'
758  IF (do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' !!! POLARIZABILITY (for the moment) = ', issc_iso_dso, ' Hz'
759  IF (.NOT. do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic coupling = ', issc_iso_tot, ' Hz'
760  END DO
761  END DO
762  END IF
763  CALL cp_print_key_finished_output(unit_atoms, logger, issc_section,&
764  & "PRINT%K_MATRIX")
765  END IF
766  !
767  !
768  END SUBROUTINE issc_print
769 
770 ! **************************************************************************************************
771 !> \brief Initialize the issc environment
772 !> \param issc_env ...
773 !> \param qs_env ...
774 ! **************************************************************************************************
775  SUBROUTINE issc_env_init(issc_env, qs_env)
776  !
777  TYPE(issc_env_type) :: issc_env
778  TYPE(qs_environment_type), POINTER :: qs_env
779 
780  CHARACTER(LEN=*), PARAMETER :: routinen = 'issc_env_init'
781 
782  INTEGER :: handle, iatom, idir, ini, ir, ispin, &
783  istat, m, n, n_rep, nao, natom, &
784  nspins, output_unit
785  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
786  INTEGER, DIMENSION(:), POINTER :: list, row_blk_sizes
787  LOGICAL :: gapw
788  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
789  TYPE(cp_fm_type), POINTER :: mo_coeff
790  TYPE(cp_logger_type), POINTER :: logger
791  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
792  TYPE(dft_control_type), POINTER :: dft_control
793  TYPE(linres_control_type), POINTER :: linres_control
794  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
795  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
796  POINTER :: sab_orb
797  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
798  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
799  TYPE(section_vals_type), POINTER :: issc_section, lr_section
800 
801 !
802 
803  CALL timeset(routinen, handle)
804 
805  NULLIFY (linres_control)
806  NULLIFY (logger, issc_section)
807  NULLIFY (tmp_fm_struct)
808  NULLIFY (particle_set, qs_kind_set)
809  NULLIFY (sab_orb)
810 
811  logger => cp_get_default_logger()
812  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
813 
814  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
815  extension=".linresLog")
816 
817  CALL issc_env_cleanup(issc_env)
818 
819  IF (output_unit > 0) THEN
820  WRITE (output_unit, "(/,T20,A,/)") "*** Start indirect spin-spin coupling Calculation ***"
821  WRITE (output_unit, "(T10,A,/)") "Inizialization of the ISSC environment"
822  END IF
823 
824  issc_section => section_vals_get_subs_vals(qs_env%input, &
825  & "PROPERTIES%LINRES%SPINSPIN")
826  !CALL section_vals_val_get(nmr_section,"INTERPOLATE_SHIFT",l_val=nmr_env%interpolate_shift)
827  !CALL section_vals_val_get(nmr_section,"SHIFT_GAPW_RADIUS",r_val=nmr_env%shift_gapw_radius)
828 
829  CALL get_qs_env(qs_env=qs_env, &
830  dft_control=dft_control, &
831  linres_control=linres_control, &
832  mos=mos, &
833  sab_orb=sab_orb, &
834  particle_set=particle_set, &
835  qs_kind_set=qs_kind_set, &
836  dbcsr_dist=dbcsr_dist)
837  !
838  !
839  gapw = dft_control%qs_control%gapw
840  nspins = dft_control%nspins
841  natom = SIZE(particle_set, 1)
842  !
843  ! check that the psi0 are localized and you have all the centers
844  IF (.NOT. linres_control%localized_psi0) &
845  CALL cp_warn(__location__, 'To get indirect spin-spin coupling parameters within '// &
846  'PBC you need to localize zero order orbitals')
847  !
848  !
849  ! read terms need to be calculated
850  ! FC
851  CALL section_vals_val_get(issc_section, "DO_FC", l_val=issc_env%do_fc)
852  ! SD
853  CALL section_vals_val_get(issc_section, "DO_SD", l_val=issc_env%do_sd)
854  ! PSO
855  CALL section_vals_val_get(issc_section, "DO_PSO", l_val=issc_env%do_pso)
856  ! DSO
857  CALL section_vals_val_get(issc_section, "DO_DSO", l_val=issc_env%do_dso)
858  !
859  !
860  ! read the list of atoms on which the issc need to be calculated
861  CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", n_rep_val=n_rep)
862  !
863  !
864  NULLIFY (issc_env%issc_on_atom_list)
865  n = 0
866  DO ir = 1, n_rep
867  NULLIFY (list)
868  CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", i_rep_val=ir, i_vals=list)
869  IF (ASSOCIATED(list)) THEN
870  CALL reallocate(issc_env%issc_on_atom_list, 1, n + SIZE(list))
871  DO ini = 1, SIZE(list)
872  issc_env%issc_on_atom_list(ini + n) = list(ini)
873  END DO
874  n = n + SIZE(list)
875  END IF
876  END DO
877  !
878  IF (.NOT. ASSOCIATED(issc_env%issc_on_atom_list)) THEN
879  ALLOCATE (issc_env%issc_on_atom_list(natom), stat=istat)
880  cpassert(istat .EQ. 0)
881  DO iatom = 1, natom
882  issc_env%issc_on_atom_list(iatom) = iatom
883  END DO
884  END IF
885  issc_env%issc_natms = SIZE(issc_env%issc_on_atom_list)
886  !
887  !
888  ! Initialize the issc tensor
889  ALLOCATE (issc_env%issc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4), &
890  issc_env%issc_loc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4), &
891  stat=istat)
892  cpassert(istat == 0)
893  issc_env%issc(:, :, :, :, :) = 0.0_dp
894  issc_env%issc_loc(:, :, :, :, :) = 0.0_dp
895  !
896  ! allocation
897  ALLOCATE (issc_env%efg_psi0(nspins, 6), issc_env%pso_psi0(nspins, 3), issc_env%fc_psi0(nspins), &
898  issc_env%psi1_efg(nspins, 6), issc_env%psi1_pso(nspins, 3), issc_env%psi1_fc(nspins), &
899  issc_env%dso_psi0(nspins, 3), issc_env%psi1_dso(nspins, 3), &
900  stat=istat)
901  cpassert(istat == 0)
902  DO ispin = 1, nspins
903  !mo_coeff => current_env%psi0_order(ispin)%matrix
904  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
905  CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
906 
907  NULLIFY (tmp_fm_struct)
908  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
909  ncol_global=m, &
910  context=mo_coeff%matrix_struct%context)
911  DO idir = 1, 6
912  CALL cp_fm_create(issc_env%psi1_efg(ispin, idir), tmp_fm_struct)
913  CALL cp_fm_create(issc_env%efg_psi0(ispin, idir), tmp_fm_struct)
914  END DO
915  DO idir = 1, 3
916  CALL cp_fm_create(issc_env%psi1_pso(ispin, idir), tmp_fm_struct)
917  CALL cp_fm_create(issc_env%pso_psi0(ispin, idir), tmp_fm_struct)
918  CALL cp_fm_create(issc_env%psi1_dso(ispin, idir), tmp_fm_struct)
919  CALL cp_fm_create(issc_env%dso_psi0(ispin, idir), tmp_fm_struct)
920  END DO
921  CALL cp_fm_create(issc_env%psi1_fc(ispin), tmp_fm_struct)
922  CALL cp_fm_create(issc_env%fc_psi0(ispin), tmp_fm_struct)
923  CALL cp_fm_struct_release(tmp_fm_struct)
924  END DO
925  !
926  ! prepare for allocation
927  ALLOCATE (first_sgf(natom))
928  ALLOCATE (last_sgf(natom))
929  CALL get_particle_set(particle_set, qs_kind_set, &
930  first_sgf=first_sgf, &
931  last_sgf=last_sgf)
932  ALLOCATE (row_blk_sizes(natom))
933  CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
934  DEALLOCATE (first_sgf)
935  DEALLOCATE (last_sgf)
936 
937  !
938  ! efg, pso and fc operators
939  CALL dbcsr_allocate_matrix_set(issc_env%matrix_efg, 6)
940  ALLOCATE (issc_env%matrix_efg(1)%matrix)
941  CALL dbcsr_create(matrix=issc_env%matrix_efg(1)%matrix, &
942  name="efg (3xx-rr)/3", &
943  dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
944  row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
945  nze=0, mutable_work=.true.)
946  CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_efg(1)%matrix, sab_orb)
947 
948  ALLOCATE (issc_env%matrix_efg(2)%matrix, &
949  issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(4)%matrix, &
950  issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(6)%matrix)
951  CALL dbcsr_copy(issc_env%matrix_efg(2)%matrix, issc_env%matrix_efg(1)%matrix, &
952  'efg xy')
953  CALL dbcsr_copy(issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(1)%matrix, &
954  'efg xz')
955  CALL dbcsr_copy(issc_env%matrix_efg(4)%matrix, issc_env%matrix_efg(1)%matrix, &
956  'efg (3yy-rr)/3')
957  CALL dbcsr_copy(issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(1)%matrix, &
958  'efg yz')
959  CALL dbcsr_copy(issc_env%matrix_efg(6)%matrix, issc_env%matrix_efg(1)%matrix, &
960  'efg (3zz-rr)/3')
961 
962  CALL dbcsr_set(issc_env%matrix_efg(1)%matrix, 0.0_dp)
963  CALL dbcsr_set(issc_env%matrix_efg(2)%matrix, 0.0_dp)
964  CALL dbcsr_set(issc_env%matrix_efg(3)%matrix, 0.0_dp)
965  CALL dbcsr_set(issc_env%matrix_efg(4)%matrix, 0.0_dp)
966  CALL dbcsr_set(issc_env%matrix_efg(5)%matrix, 0.0_dp)
967  CALL dbcsr_set(issc_env%matrix_efg(6)%matrix, 0.0_dp)
968  !
969  ! PSO
970  CALL dbcsr_allocate_matrix_set(issc_env%matrix_pso, 3)
971  ALLOCATE (issc_env%matrix_pso(1)%matrix)
972  CALL dbcsr_create(matrix=issc_env%matrix_pso(1)%matrix, &
973  name="pso x", &
974  dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
975  row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
976  nze=0, mutable_work=.true.)
977  CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_pso(1)%matrix, sab_orb)
978 
979  ALLOCATE (issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(3)%matrix)
980  CALL dbcsr_copy(issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(1)%matrix, &
981  'pso y')
982  CALL dbcsr_copy(issc_env%matrix_pso(3)%matrix, issc_env%matrix_pso(1)%matrix, &
983  'pso z')
984  CALL dbcsr_set(issc_env%matrix_pso(1)%matrix, 0.0_dp)
985  CALL dbcsr_set(issc_env%matrix_pso(2)%matrix, 0.0_dp)
986  CALL dbcsr_set(issc_env%matrix_pso(3)%matrix, 0.0_dp)
987  !
988  ! DSO
989  CALL dbcsr_allocate_matrix_set(issc_env%matrix_dso, 3)
990  ALLOCATE (issc_env%matrix_dso(1)%matrix, issc_env%matrix_dso(2)%matrix, issc_env%matrix_dso(3)%matrix)
991  CALL dbcsr_copy(issc_env%matrix_dso(1)%matrix, issc_env%matrix_efg(1)%matrix, &
992  'dso x')
993  CALL dbcsr_copy(issc_env%matrix_dso(2)%matrix, issc_env%matrix_efg(1)%matrix, &
994  'dso y')
995  CALL dbcsr_copy(issc_env%matrix_dso(3)%matrix, issc_env%matrix_efg(1)%matrix, &
996  'dso z')
997  CALL dbcsr_set(issc_env%matrix_dso(1)%matrix, 0.0_dp)
998  CALL dbcsr_set(issc_env%matrix_dso(2)%matrix, 0.0_dp)
999  CALL dbcsr_set(issc_env%matrix_dso(3)%matrix, 0.0_dp)
1000  !
1001  ! FC
1002  CALL dbcsr_allocate_matrix_set(issc_env%matrix_fc, 1)
1003  ALLOCATE (issc_env%matrix_fc(1)%matrix)
1004  CALL dbcsr_copy(issc_env%matrix_fc(1)%matrix, issc_env%matrix_efg(1)%matrix, &
1005  'fc')
1006  CALL dbcsr_set(issc_env%matrix_fc(1)%matrix, 0.0_dp)
1007 
1008  DEALLOCATE (row_blk_sizes)
1009  !
1010  ! Conversion factors
1011  IF (output_unit > 0) THEN
1012  WRITE (output_unit, "(T2,A,T60,I4,A)")&
1013  & "ISSC| spin-spin coupling computed for ", issc_env%issc_natms, ' atoms'
1014  END IF
1015 
1016  CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
1017  & "PRINT%PROGRAM_RUN_INFO")
1018 
1019  CALL timestop(handle)
1020 
1021  END SUBROUTINE issc_env_init
1022 
1023 ! **************************************************************************************************
1024 !> \brief Deallocate the issc environment
1025 !> \param issc_env ...
1026 !> \par History
1027 ! **************************************************************************************************
1028  SUBROUTINE issc_env_cleanup(issc_env)
1029 
1030  TYPE(issc_env_type), INTENT(INOUT) :: issc_env
1031 
1032  IF (ASSOCIATED(issc_env%issc_on_atom_list)) THEN
1033  DEALLOCATE (issc_env%issc_on_atom_list)
1034  END IF
1035  IF (ASSOCIATED(issc_env%issc)) THEN
1036  DEALLOCATE (issc_env%issc)
1037  END IF
1038  IF (ASSOCIATED(issc_env%issc_loc)) THEN
1039  DEALLOCATE (issc_env%issc_loc)
1040  END IF
1041  !
1042  !efg_psi0
1043  CALL cp_fm_release(issc_env%efg_psi0)
1044  !
1045  !pso_psi0
1046  CALL cp_fm_release(issc_env%pso_psi0)
1047  !
1048  !dso_psi0
1049  CALL cp_fm_release(issc_env%dso_psi0)
1050  !
1051  !fc_psi0
1052  CALL cp_fm_release(issc_env%fc_psi0)
1053  !
1054  !psi1_efg
1055  CALL cp_fm_release(issc_env%psi1_efg)
1056  !
1057  !psi1_pso
1058  CALL cp_fm_release(issc_env%psi1_pso)
1059  !
1060  !psi1_dso
1061  CALL cp_fm_release(issc_env%psi1_dso)
1062  !
1063  !psi1_fc
1064  CALL cp_fm_release(issc_env%psi1_fc)
1065  !
1066  !matrix_efg
1067  IF (ASSOCIATED(issc_env%matrix_efg)) THEN
1068  CALL dbcsr_deallocate_matrix_set(issc_env%matrix_efg)
1069  END IF
1070  !
1071  !matrix_pso
1072  IF (ASSOCIATED(issc_env%matrix_pso)) THEN
1073  CALL dbcsr_deallocate_matrix_set(issc_env%matrix_pso)
1074  END IF
1075  !
1076  !matrix_dso
1077  IF (ASSOCIATED(issc_env%matrix_dso)) THEN
1078  CALL dbcsr_deallocate_matrix_set(issc_env%matrix_dso)
1079  END IF
1080  !
1081  !matrix_fc
1082  IF (ASSOCIATED(issc_env%matrix_fc)) THEN
1083  CALL dbcsr_deallocate_matrix_set(issc_env%matrix_fc)
1084  END IF
1085 
1086  END SUBROUTINE issc_env_cleanup
1087 
1088 END MODULE qs_linres_issc_utils
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_frobenius_norm
real(kind=dp) function, public cp_fm_frobenius_norm(matrix_a)
computes the Frobenius norm of matrix_a
Definition: cp_fm_basic_linalg.F:629
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:372
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::a_fine
real(kind=dp), parameter, public a_fine
Definition: physcon.F:119
physcon::hertz
real(kind=dp), parameter, public hertz
Definition: physcon.F:186
physcon::p_mass
real(kind=dp), parameter, public p_mass
Definition: physcon.F:112
physcon::e_mass
real(kind=dp), parameter, public e_mass
Definition: physcon.F:109
qs_elec_field
Distribution of the electric field gradient integral matrix.
Definition: qs_elec_field.F:13
qs_elec_field::build_efg_matrix
subroutine, public build_efg_matrix(qs_env, matrix_efg, rc)
Calculation of the electric field gradient matrix over Cartesian Gaussian functions.
Definition: qs_elec_field.F:75
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_fermi_contact
Distribution of the Fermi contact integral matrix.
Definition: qs_fermi_contact.F:13
qs_fermi_contact::build_fermi_contact_matrix
subroutine, public build_fermi_contact_matrix(qs_env, matrix_fc, rc)
Calculation of the Fermi contact matrix over Cartesian Gaussian functions.
Definition: qs_fermi_contact.F:75
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_linres_issc_utils
Chemical shift calculation by dfpt Initialization of the issc_env, creation of the special neighbor l...
Definition: qs_linres_issc_utils.F:19
qs_linres_issc_utils::issc_issc
subroutine, public issc_issc(issc_env, qs_env, iatom)
...
Definition: qs_linres_issc_utils.F:447
qs_linres_issc_utils::issc_response
subroutine, public issc_response(issc_env, p_env, qs_env)
Initialize the issc environment.
Definition: qs_linres_issc_utils.F:102
qs_linres_issc_utils::issc_env_cleanup
subroutine, public issc_env_cleanup(issc_env)
Deallocate the issc environment.
Definition: qs_linres_issc_utils.F:1029
qs_linres_issc_utils::issc_env_init
subroutine, public issc_env_init(issc_env, qs_env)
Initialize the issc environment.
Definition: qs_linres_issc_utils.F:776
qs_linres_issc_utils::issc_print
subroutine, public issc_print(issc_env, qs_env)
...
Definition: qs_linres_issc_utils.F:653
qs_linres_methods
localize wavefunctions linear response scf
Definition: qs_linres_methods.F:15
qs_linres_methods::linres_solver
subroutine, public linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition: qs_linres_methods.F:216
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_linres_types::get_issc_env
subroutine, public get_issc_env(issc_env, issc_on_atom_list, issc_gapw_radius, issc_loc, do_fc, do_sd, do_pso, do_dso, issc, interpolate_issc, psi1_efg, psi1_pso, psi1_dso, psi1_fc, efg_psi0, pso_psi0, dso_psi0, fc_psi0, matrix_efg, matrix_pso, matrix_dso, matrix_fc)
...
Definition: qs_linres_types.F:558
qs_matrix_pools
wrapper for the pools of matrixes
Definition: qs_matrix_pools.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_spin_orbit
Distribution of the spin orbit integral matrix.
Definition: qs_spin_orbit.F:13
qs_spin_orbit::build_pso_matrix
subroutine, public build_pso_matrix(qs_env, matrix_so, rc)
Calculation of the paramagnetic spin orbit matrix over Cartesian Gaussian functions.
Definition: qs_spin_orbit.F:74