Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 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 44 : 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 : !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 0 : DO ispin = 1, nspins
209 0 : CALL cp_fm_to_fm(efg_psi0(ispin, ijdir), h1_psi0(ispin))
210 : END DO
211 : !
212 : !
213 0 : linres_control%lr_triplet = .FALSE.
214 0 : linres_control%do_kernel = .FALSE.
215 0 : linres_control%converged = .FALSE.
216 0 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
217 : !
218 : !
219 : ! copy the response
220 0 : DO ispin = 1, nspins
221 0 : CALL cp_fm_to_fm(psi1(ispin), psi1_efg(ispin, ijdir))
222 0 : fro = cp_fm_frobenius_norm(psi1(ispin))
223 0 : 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 0 : IF (output_unit > 0) THEN
240 0 : 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 44 : IF (do_pso) THEN
250 136 : DO idir = 1, 3
251 216 : DO ispin = 1, nspins
252 216 : CALL cp_fm_set_all(psi1_pso(ispin, idir), 0.0_dp)
253 : END DO
254 102 : IF (output_unit > 0) THEN
255 51 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator pso_"//ACHAR(idir + 119)
256 : END IF
257 : !
258 : !Initial guess for psi1
259 216 : DO ispin = 1, nspins
260 216 : 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 216 : DO ispin = 1, nspins
267 216 : CALL cp_fm_to_fm(pso_psi0(ispin, idir), h1_psi0(ispin))
268 : END DO
269 : !
270 : !
271 102 : linres_control%lr_triplet = .FALSE. ! we do singlet response
272 102 : linres_control%do_kernel = .FALSE. ! we do uncoupled response
273 102 : linres_control%converged = .FALSE.
274 102 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
275 : !
276 : !
277 : ! copy the response
278 216 : DO ispin = 1, nspins
279 114 : CALL cp_fm_to_fm(psi1(ispin), psi1_pso(ispin, idir))
280 114 : fro = cp_fm_frobenius_norm(psi1(ispin))
281 216 : 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 238 : IF (output_unit > 0) THEN
298 51 : 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 44 : IF (do_fc) THEN
307 0 : DO ispin = 1, nspins
308 0 : CALL cp_fm_set_all(psi1_fc(ispin), 0.0_dp)
309 : END DO
310 0 : IF (output_unit > 0) THEN
311 0 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator fc"
312 : END IF
313 : !
314 : !Initial guess for psi1
315 0 : DO ispin = 1, nspins
316 0 : 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 0 : DO ispin = 1, nspins
323 0 : CALL cp_fm_to_fm(fc_psi0(ispin), h1_psi0(ispin))
324 : END DO
325 : !
326 : !
327 0 : linres_control%lr_triplet = .TRUE. ! we do triplet response
328 0 : linres_control%do_kernel = .TRUE. ! we do coupled response
329 0 : linres_control%converged = .FALSE.
330 0 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
331 : !
332 : !
333 : ! copy the response
334 0 : DO ispin = 1, nspins
335 0 : CALL cp_fm_to_fm(psi1(ispin), psi1_fc(ispin))
336 0 : fro = cp_fm_frobenius_norm(psi1(ispin))
337 0 : 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 0 : IF (output_unit > 0) THEN
354 0 : 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 44 : IF (do_dso) THEN
364 8 : DO idir = 1, 3
365 12 : DO ispin = 1, nspins
366 12 : CALL cp_fm_set_all(psi1_dso(ispin, idir), 0.0_dp)
367 : END DO
368 6 : IF (output_unit > 0) THEN
369 3 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator r_"//ACHAR(idir + 119)
370 : END IF
371 : !
372 : !Initial guess for psi1
373 12 : DO ispin = 1, nspins
374 12 : 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 12 : DO ispin = 1, nspins
381 12 : CALL cp_fm_to_fm(dso_psi0(ispin, idir), h1_psi0(ispin))
382 : END DO
383 : !
384 : !
385 6 : linres_control%lr_triplet = .FALSE. ! we do singlet response
386 6 : linres_control%do_kernel = .TRUE. ! we do uncoupled response
387 6 : linres_control%converged = .FALSE.
388 6 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
389 : !
390 : !
391 : ! copy the response
392 12 : DO ispin = 1, nspins
393 6 : CALL cp_fm_to_fm(psi1(ispin), psi1_dso(ispin, idir))
394 6 : fro = cp_fm_frobenius_norm(psi1(ispin))
395 12 : 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 14 : IF (output_unit > 0) THEN
412 3 : 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 44 : IF (output_unit > 0) THEN
424 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| response: CheckSum =', chk
425 : END IF
426 : !
427 : !
428 : ! clean up
429 44 : CALL cp_fm_release(psi1)
430 44 : CALL cp_fm_release(h1_psi0)
431 44 : DEALLOCATE (psi0_order)
432 : !
433 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
434 44 : & "PRINT%PROGRAM_RUN_INFO")
435 : !
436 44 : CALL timestop(handle)
437 : !
438 88 : END SUBROUTINE issc_response
439 :
440 : ! **************************************************************************************************
441 : !> \brief ...
442 : !> \param issc_env ...
443 : !> \param qs_env ...
444 : !> \param iatom ...
445 : ! **************************************************************************************************
446 44 : 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 44 : REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc
462 : TYPE(cell_type), POINTER :: cell
463 44 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc
464 44 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_dso, psi1_efg, psi1_pso
465 : TYPE(cp_fm_type), POINTER :: mo_coeff
466 44 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, &
467 44 : matrix_pso
468 : TYPE(dft_control_type), POINTER :: dft_control
469 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
470 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
471 : TYPE(section_vals_type), POINTER :: issc_section
472 :
473 44 : CALL timeset(routineN, handle)
474 :
475 44 : NULLIFY (cell, dft_control, particle_set, issc, psi1_fc, psi1_efg, psi1_pso)
476 44 : 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 44 : mos=mos)
483 :
484 44 : gapw = dft_control%qs_control%gapw
485 44 : natom = SIZE(particle_set, 1)
486 44 : 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 44 : do_dso=do_dso)
503 :
504 44 : g = e_mass/(2.0_dp*p_mass)
505 44 : facfc = hertz*g**2*a_fine**4
506 44 : facpso = hertz*g**2*a_fine**4
507 44 : facsd = hertz*g**2*a_fine**4
508 44 : facdso = hertz*g**2*a_fine**4
509 :
510 : !
511 : !
512 : issc_section => section_vals_get_subs_vals(qs_env%input, &
513 44 : & "PROPERTIES%LINRES%SPINSPIN")
514 : !
515 : ! Initialize
516 92 : DO ispin = 1, nspins
517 48 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, maxocc=maxocc)
518 48 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
519 :
520 292 : DO jatom = 1, natom
521 800 : r_i = particle_set(iatom)%r
522 800 : r_j = particle_set(jatom)%r
523 800 : 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 200 : IF (do_fc .AND. iatom .NE. jatom) THEN
533 : !
534 : ! build the integral for the jatom
535 0 : CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp)
536 0 : 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 0 : & alpha=1.0_dp)
540 :
541 0 : CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
542 0 : WRITE (*, *) ' jatom', jatom, 'tr(P*fc)=', buf
543 :
544 0 : CALL cp_fm_trace(fc_psi0(ispin), psi1_fc(ispin), buf)
545 0 : issc_fc = 2.0_dp*2.0_dp*maxocc*facfc*buf
546 0 : issc(1, 1, iatom, jatom, 1) = issc(1, 1, iatom, jatom, 1) + issc_fc
547 0 : issc(2, 2, iatom, jatom, 1) = issc(2, 2, iatom, jatom, 1) + issc_fc
548 0 : issc(3, 3, iatom, jatom, 1) = issc(3, 3, iatom, jatom, 1) + issc_fc
549 : END IF
550 : !
551 : ! SD term
552 : !
553 200 : IF (do_sd .AND. iatom .NE. jatom) THEN
554 : !
555 : ! build the integral for the jatom
556 0 : CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp)
557 0 : CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp)
558 0 : CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp)
559 0 : CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp)
560 0 : CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp)
561 0 : CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp)
562 0 : CALL build_efg_matrix(qs_env, matrix_efg, r_j)
563 0 : 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 0 : & alpha=1.0_dp, beta=0.0_dp)
567 0 : CALL cp_fm_trace(fc_psi0(ispin), mo_coeff, buf)
568 0 : WRITE (*, *) ' jatom', jatom, ixyz, 'tr(P*efg)=', buf
569 0 : DO jxyz = 1, 6
570 0 : CALL cp_fm_trace(fc_psi0(ispin), psi1_efg(ispin, jxyz), buf)
571 0 : 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 200 : IF (do_pso .AND. iatom .NE. jatom) THEN
581 : !
582 : ! build the integral for the jatom
583 128 : CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp)
584 128 : CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp)
585 128 : CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp)
586 128 : CALL build_pso_matrix(qs_env, matrix_pso, r_j)
587 512 : 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 384 : & alpha=1.0_dp, beta=0.0_dp)
591 1664 : DO jxyz = 1, 3
592 1152 : CALL cp_fm_trace(fc_psi0(ispin), psi1_pso(ispin, jxyz), buf)
593 1152 : issc_pso = -2.0_dp*maxocc*facpso*buf
594 1536 : 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 248 : 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 8 : 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 6 : & alpha=1.0_dp, beta=0.0_dp)
628 26 : DO jxyz = 1, 3
629 18 : CALL cp_fm_trace(psi1_dso(ispin, jxyz), fc_psi0(ispin), buf)
630 : ! we save the polarizability for a checksum later on !
631 18 : issc_dso = 2.0_dp*maxocc*buf
632 24 : 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 44 : CALL timestop(handle)
644 : !
645 44 : END SUBROUTINE issc_issc
646 :
647 : ! **************************************************************************************************
648 : !> \brief ...
649 : !> \param issc_env ...
650 : !> \param qs_env ...
651 : ! **************************************************************************************************
652 12 : 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 12 : 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 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
670 : TYPE(section_vals_type), POINTER :: issc_section
671 :
672 12 : NULLIFY (logger, particle_set, atom_kind_i, atom_kind_j, dft_control)
673 :
674 24 : logger => cp_get_default_logger()
675 12 : output_unit = cp_logger_get_default_io_unit(logger)
676 :
677 : issc_section => section_vals_get_subs_vals(qs_env%input, &
678 12 : "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 12 : do_dso=do_dso)
686 : !
687 : CALL get_qs_env(qs_env=qs_env, &
688 : dft_control=dft_control, &
689 12 : particle_set=particle_set)
690 :
691 12 : natom = SIZE(particle_set, 1)
692 12 : gapw = dft_control%qs_control%gapw
693 :
694 : !
695 12 : IF (output_unit > 0) THEN
696 6 : WRITE (output_unit, '(T2,A,E14.6)') 'ISSC| CheckSum K =', &
697 42 : SQRT(DDOT(SIZE(issc), issc, 1, issc, 1))
698 : END IF
699 : !
700 12 : 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 12 : extension=".data", middle_name="K", log_filename=.FALSE.)
705 :
706 12 : IF (unit_atoms > 0) THEN
707 6 : WRITE (unit_atoms, *)
708 6 : WRITE (unit_atoms, *)
709 6 : WRITE (title, '(A)') "Indirect spin-spin coupling matrix"
710 6 : WRITE (unit_atoms, '(T2,A)') title
711 28 : DO iatom = 1, natom
712 22 : atom_kind_i => particle_set(iatom)%atomic_kind
713 22 : CALL get_atomic_kind(atom_kind_i, name=name_i, element_symbol=element_symbol_i)
714 124 : DO jatom = 1, natom
715 96 : atom_kind_j => particle_set(jatom)%atomic_kind
716 96 : CALL get_atomic_kind(atom_kind_j, name=name_j, element_symbol=element_symbol_j)
717 : !
718 96 : IF (iatom .EQ. jatom .AND. .NOT. do_dso) CYCLE
719 : !
720 : !
721 : ! FC
722 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 1)
723 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
724 75 : CALL diamat_all(issc_tmp, eig)
725 75 : issc_iso_fc = (eig(1) + eig(2) + eig(3))/3.0_dp
726 : !
727 : ! SD
728 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 2)
729 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
730 75 : CALL diamat_all(issc_tmp, eig)
731 75 : issc_iso_sd = (eig(1) + eig(2) + eig(3))/3.0_dp
732 : !
733 : ! PSO
734 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 3)
735 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
736 75 : CALL diamat_all(issc_tmp, eig)
737 75 : issc_iso_pso = (eig(1) + eig(2) + eig(3))/3.0_dp
738 : !
739 : ! DSO
740 975 : issc_tmp(:, :) = issc(:, :, iatom, jatom, 4)
741 1875 : issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :)))
742 75 : CALL diamat_all(issc_tmp, eig)
743 75 : issc_iso_dso = (eig(1) + eig(2) + eig(3))/3.0_dp
744 : !
745 : ! TOT
746 75 : issc_iso_tot = issc_iso_fc + issc_iso_sd + issc_iso_dso + issc_iso_pso
747 : !
748 : !
749 75 : WRITE (unit_atoms, *)
750 75 : WRITE (unit_atoms, '(T2,2(A,I5,A,2X,A2))') 'Indirect spin-spin coupling between ', &
751 75 : iatom, TRIM(name_i), element_symbol_i, ' and ', &
752 150 : jatom, TRIM(name_j), element_symbol_j
753 : !
754 75 : IF (do_fc) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic FC contribution = ', issc_iso_fc, ' Hz'
755 75 : IF (do_sd) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic SD contribution = ', issc_iso_sd, ' Hz'
756 75 : 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 75 : IF (do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' !!! POLARIZABILITY (for the moment) = ', issc_iso_dso, ' Hz'
759 97 : 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 12 : & "PRINT%K_MATRIX")
765 : END IF
766 : !
767 : !
768 12 : END SUBROUTINE issc_print
769 :
770 : ! **************************************************************************************************
771 : !> \brief Initialize the issc environment
772 : !> \param issc_env ...
773 : !> \param qs_env ...
774 : ! **************************************************************************************************
775 12 : 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 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
786 12 : 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 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
795 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
796 12 : POINTER :: sab_orb
797 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
798 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
799 : TYPE(section_vals_type), POINTER :: issc_section, lr_section
800 :
801 : !
802 :
803 12 : CALL timeset(routineN, handle)
804 :
805 12 : NULLIFY (linres_control)
806 12 : NULLIFY (logger, issc_section)
807 12 : NULLIFY (tmp_fm_struct)
808 12 : NULLIFY (particle_set, qs_kind_set)
809 12 : NULLIFY (sab_orb)
810 :
811 12 : logger => cp_get_default_logger()
812 12 : 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 12 : extension=".linresLog")
816 :
817 12 : CALL issc_env_cleanup(issc_env)
818 :
819 12 : IF (output_unit > 0) THEN
820 6 : WRITE (output_unit, "(/,T20,A,/)") "*** Start indirect spin-spin coupling Calculation ***"
821 6 : 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 12 : & "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 12 : dbcsr_dist=dbcsr_dist)
837 : !
838 : !
839 12 : gapw = dft_control%qs_control%gapw
840 12 : nspins = dft_control%nspins
841 12 : natom = SIZE(particle_set, 1)
842 : !
843 : ! check that the psi0 are localized and you have all the centers
844 12 : IF (.NOT. linres_control%localized_psi0) &
845 : CALL cp_warn(__LOCATION__, 'To get indirect spin-spin coupling parameters within '// &
846 0 : 'PBC you need to localize zero order orbitals')
847 : !
848 : !
849 : ! read terms need to be calculated
850 : ! FC
851 12 : CALL section_vals_val_get(issc_section, "DO_FC", l_val=issc_env%do_fc)
852 : ! SD
853 12 : CALL section_vals_val_get(issc_section, "DO_SD", l_val=issc_env%do_sd)
854 : ! PSO
855 12 : CALL section_vals_val_get(issc_section, "DO_PSO", l_val=issc_env%do_pso)
856 : ! DSO
857 12 : 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 12 : CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", n_rep_val=n_rep)
862 : !
863 : !
864 12 : NULLIFY (issc_env%issc_on_atom_list)
865 12 : n = 0
866 16 : DO ir = 1, n_rep
867 4 : NULLIFY (list)
868 4 : CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", i_rep_val=ir, i_vals=list)
869 16 : IF (ASSOCIATED(list)) THEN
870 4 : CALL reallocate(issc_env%issc_on_atom_list, 1, n + SIZE(list))
871 14 : DO ini = 1, SIZE(list)
872 14 : issc_env%issc_on_atom_list(ini + n) = list(ini)
873 : END DO
874 4 : n = n + SIZE(list)
875 : END IF
876 : END DO
877 : !
878 12 : IF (.NOT. ASSOCIATED(issc_env%issc_on_atom_list)) THEN
879 30 : ALLOCATE (issc_env%issc_on_atom_list(natom), STAT=istat)
880 10 : CPASSERT(istat .EQ. 0)
881 44 : DO iatom = 1, natom
882 44 : issc_env%issc_on_atom_list(iatom) = iatom
883 : END DO
884 : END IF
885 12 : 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 84 : STAT=istat)
892 12 : CPASSERT(istat == 0)
893 10220 : issc_env%issc(:, :, :, :, :) = 0.0_dp
894 10220 : 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 796 : STAT=istat)
901 12 : CPASSERT(istat == 0)
902 26 : DO ispin = 1, nspins
903 : !mo_coeff => current_env%psi0_order(ispin)%matrix
904 14 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
905 14 : CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao)
906 :
907 14 : NULLIFY (tmp_fm_struct)
908 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
909 : ncol_global=m, &
910 14 : context=mo_coeff%matrix_struct%context)
911 98 : DO idir = 1, 6
912 84 : CALL cp_fm_create(issc_env%psi1_efg(ispin, idir), tmp_fm_struct)
913 98 : CALL cp_fm_create(issc_env%efg_psi0(ispin, idir), tmp_fm_struct)
914 : END DO
915 56 : DO idir = 1, 3
916 42 : CALL cp_fm_create(issc_env%psi1_pso(ispin, idir), tmp_fm_struct)
917 42 : CALL cp_fm_create(issc_env%pso_psi0(ispin, idir), tmp_fm_struct)
918 42 : CALL cp_fm_create(issc_env%psi1_dso(ispin, idir), tmp_fm_struct)
919 56 : CALL cp_fm_create(issc_env%dso_psi0(ispin, idir), tmp_fm_struct)
920 : END DO
921 14 : CALL cp_fm_create(issc_env%psi1_fc(ispin), tmp_fm_struct)
922 14 : CALL cp_fm_create(issc_env%fc_psi0(ispin), tmp_fm_struct)
923 40 : CALL cp_fm_struct_release(tmp_fm_struct)
924 : END DO
925 : !
926 : ! prepare for allocation
927 36 : ALLOCATE (first_sgf(natom))
928 24 : ALLOCATE (last_sgf(natom))
929 : CALL get_particle_set(particle_set, qs_kind_set, &
930 : first_sgf=first_sgf, &
931 12 : last_sgf=last_sgf)
932 24 : ALLOCATE (row_blk_sizes(natom))
933 12 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
934 12 : DEALLOCATE (first_sgf)
935 12 : DEALLOCATE (last_sgf)
936 :
937 : !
938 : ! efg, pso and fc operators
939 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_efg, 6)
940 12 : 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 12 : nze=0, mutable_work=.TRUE.)
946 12 : 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 12 : 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 12 : 'efg xy')
953 : CALL dbcsr_copy(issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(1)%matrix, &
954 12 : 'efg xz')
955 : CALL dbcsr_copy(issc_env%matrix_efg(4)%matrix, issc_env%matrix_efg(1)%matrix, &
956 12 : 'efg (3yy-rr)/3')
957 : CALL dbcsr_copy(issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(1)%matrix, &
958 12 : 'efg yz')
959 : CALL dbcsr_copy(issc_env%matrix_efg(6)%matrix, issc_env%matrix_efg(1)%matrix, &
960 12 : 'efg (3zz-rr)/3')
961 :
962 12 : CALL dbcsr_set(issc_env%matrix_efg(1)%matrix, 0.0_dp)
963 12 : CALL dbcsr_set(issc_env%matrix_efg(2)%matrix, 0.0_dp)
964 12 : CALL dbcsr_set(issc_env%matrix_efg(3)%matrix, 0.0_dp)
965 12 : CALL dbcsr_set(issc_env%matrix_efg(4)%matrix, 0.0_dp)
966 12 : CALL dbcsr_set(issc_env%matrix_efg(5)%matrix, 0.0_dp)
967 12 : CALL dbcsr_set(issc_env%matrix_efg(6)%matrix, 0.0_dp)
968 : !
969 : ! PSO
970 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_pso, 3)
971 12 : 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 12 : nze=0, mutable_work=.TRUE.)
977 12 : CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_pso(1)%matrix, sab_orb)
978 :
979 12 : 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 12 : 'pso y')
982 : CALL dbcsr_copy(issc_env%matrix_pso(3)%matrix, issc_env%matrix_pso(1)%matrix, &
983 12 : 'pso z')
984 12 : CALL dbcsr_set(issc_env%matrix_pso(1)%matrix, 0.0_dp)
985 12 : CALL dbcsr_set(issc_env%matrix_pso(2)%matrix, 0.0_dp)
986 12 : CALL dbcsr_set(issc_env%matrix_pso(3)%matrix, 0.0_dp)
987 : !
988 : ! DSO
989 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_dso, 3)
990 12 : 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 12 : 'dso x')
993 : CALL dbcsr_copy(issc_env%matrix_dso(2)%matrix, issc_env%matrix_efg(1)%matrix, &
994 12 : 'dso y')
995 : CALL dbcsr_copy(issc_env%matrix_dso(3)%matrix, issc_env%matrix_efg(1)%matrix, &
996 12 : 'dso z')
997 12 : CALL dbcsr_set(issc_env%matrix_dso(1)%matrix, 0.0_dp)
998 12 : CALL dbcsr_set(issc_env%matrix_dso(2)%matrix, 0.0_dp)
999 12 : CALL dbcsr_set(issc_env%matrix_dso(3)%matrix, 0.0_dp)
1000 : !
1001 : ! FC
1002 12 : CALL dbcsr_allocate_matrix_set(issc_env%matrix_fc, 1)
1003 12 : ALLOCATE (issc_env%matrix_fc(1)%matrix)
1004 : CALL dbcsr_copy(issc_env%matrix_fc(1)%matrix, issc_env%matrix_efg(1)%matrix, &
1005 12 : 'fc')
1006 12 : CALL dbcsr_set(issc_env%matrix_fc(1)%matrix, 0.0_dp)
1007 :
1008 12 : DEALLOCATE (row_blk_sizes)
1009 : !
1010 : ! Conversion factors
1011 12 : IF (output_unit > 0) THEN
1012 : WRITE (output_unit, "(T2,A,T60,I4,A)")&
1013 6 : & "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 12 : & "PRINT%PROGRAM_RUN_INFO")
1018 :
1019 12 : CALL timestop(handle)
1020 :
1021 24 : END SUBROUTINE issc_env_init
1022 :
1023 : ! **************************************************************************************************
1024 : !> \brief Deallocate the issc environment
1025 : !> \param issc_env ...
1026 : !> \par History
1027 : ! **************************************************************************************************
1028 24 : SUBROUTINE issc_env_cleanup(issc_env)
1029 :
1030 : TYPE(issc_env_type), INTENT(INOUT) :: issc_env
1031 :
1032 24 : IF (ASSOCIATED(issc_env%issc_on_atom_list)) THEN
1033 12 : DEALLOCATE (issc_env%issc_on_atom_list)
1034 : END IF
1035 24 : IF (ASSOCIATED(issc_env%issc)) THEN
1036 12 : DEALLOCATE (issc_env%issc)
1037 : END IF
1038 24 : IF (ASSOCIATED(issc_env%issc_loc)) THEN
1039 12 : DEALLOCATE (issc_env%issc_loc)
1040 : END IF
1041 : !
1042 : !efg_psi0
1043 24 : CALL cp_fm_release(issc_env%efg_psi0)
1044 : !
1045 : !pso_psi0
1046 24 : CALL cp_fm_release(issc_env%pso_psi0)
1047 : !
1048 : !dso_psi0
1049 24 : CALL cp_fm_release(issc_env%dso_psi0)
1050 : !
1051 : !fc_psi0
1052 24 : CALL cp_fm_release(issc_env%fc_psi0)
1053 : !
1054 : !psi1_efg
1055 24 : CALL cp_fm_release(issc_env%psi1_efg)
1056 : !
1057 : !psi1_pso
1058 24 : CALL cp_fm_release(issc_env%psi1_pso)
1059 : !
1060 : !psi1_dso
1061 24 : CALL cp_fm_release(issc_env%psi1_dso)
1062 : !
1063 : !psi1_fc
1064 24 : CALL cp_fm_release(issc_env%psi1_fc)
1065 : !
1066 : !matrix_efg
1067 24 : IF (ASSOCIATED(issc_env%matrix_efg)) THEN
1068 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_efg)
1069 : END IF
1070 : !
1071 : !matrix_pso
1072 24 : IF (ASSOCIATED(issc_env%matrix_pso)) THEN
1073 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_pso)
1074 : END IF
1075 : !
1076 : !matrix_dso
1077 24 : IF (ASSOCIATED(issc_env%matrix_dso)) THEN
1078 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_dso)
1079 : END IF
1080 : !
1081 : !matrix_fc
1082 24 : IF (ASSOCIATED(issc_env%matrix_fc)) THEN
1083 12 : CALL dbcsr_deallocate_matrix_set(issc_env%matrix_fc)
1084 : END IF
1085 :
1086 24 : END SUBROUTINE issc_env_cleanup
1087 :
1088 : END MODULE qs_linres_issc_utils
|
