1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Chemical shift calculation by dfpt 8!> Initialization of the issc_env, creation of the special neighbor lists 9!> Perturbation Hamiltonians by application of the p and rxp oprtators to psi0 10!> Write output 11!> Deallocate everything 12!> \note 13!> The psi0 should be localized 14!> the Sebastiani method works within the assumption that the orbitals are 15!> completely contained in the simulation box 16! ************************************************************************************************** 17MODULE qs_linres_issc_utils 18 USE atomic_kind_types, ONLY: atomic_kind_type,& 19 get_atomic_kind 20 USE cell_types, ONLY: cell_type,& 21 pbc 22 USE cp_control_types, ONLY: dft_control_type 23 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 24 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& 25 dbcsr_allocate_matrix_set,& 26 dbcsr_deallocate_matrix_set 27 USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,& 28 cp_fm_trace 29 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 30 cp_fm_struct_release,& 31 cp_fm_struct_type 32 USE cp_fm_types, ONLY: cp_fm_create,& 33 cp_fm_get_info,& 34 cp_fm_p_type,& 35 cp_fm_release,& 36 cp_fm_set_all,& 37 cp_fm_to_fm,& 38 cp_fm_type 39 USE cp_log_handling, ONLY: cp_get_default_logger,& 40 cp_logger_get_default_io_unit,& 41 cp_logger_type 42 USE cp_output_handling, ONLY: cp_p_file,& 43 cp_print_key_finished_output,& 44 cp_print_key_should_output,& 45 cp_print_key_unit_nr 46 USE cp_para_types, ONLY: cp_para_env_type 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 particle_methods, ONLY: get_particle_set 63 USE particle_types, ONLY: particle_type 64 USE physcon, ONLY: a_fine,& 65 e_mass,& 66 hertz,& 67 p_mass 68 USE qs_elec_field, ONLY: build_efg_matrix 69 USE qs_environment_types, ONLY: get_qs_env,& 70 qs_environment_type 71 USE qs_fermi_contact, ONLY: build_fermi_contact_matrix 72 USE qs_kind_types, ONLY: qs_kind_type 73 USE qs_linres_methods, ONLY: linres_solver 74 USE qs_linres_types, ONLY: get_issc_env,& 75 issc_env_create,& 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_p_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 93CONTAINS 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), POINTER :: p_env 105 TYPE(qs_environment_type), POINTER :: qs_env 106 107 CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_response', routineP = moduleN//':'//routineN 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_p_type), DIMENSION(:), POINTER :: fc_psi0, h1_psi0, psi0_order, psi1, & 114 psi1_fc 115 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: dso_psi0, efg_psi0, psi1_dso, psi1_efg, & 116 psi1_pso, pso_psi0 117 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 118 TYPE(cp_fm_type), POINTER :: mo_coeff 119 TYPE(cp_logger_type), POINTER :: logger 120 TYPE(cp_para_env_type), POINTER :: para_env 121 TYPE(dft_control_type), POINTER :: dft_control 122 TYPE(linres_control_type), POINTER :: linres_control 123 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 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, psi1, h1_psi0, mo_coeff, para_env) 131 NULLIFY (tmp_fm_struct, psi1_fc, psi1_efg, psi1_pso, pso_psi0, fc_psi0, efg_psi0, psi0_order) 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 ENDIF 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_set, mo_coeff=mo_coeff) 174 psi0_order(ispin)%matrix => mo_coeff 175 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao) 176 NULLIFY (tmp_fm_struct, psi1(ispin)%matrix, h1_psi0(ispin)%matrix) 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)%matrix, tmp_fm_struct) 181 CALL cp_fm_create(h1_psi0(ispin)%matrix, tmp_fm_struct) 182 CALL cp_fm_struct_release(tmp_fm_struct) 183 ENDDO 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)%matrix, 0.0_dp) 195 ENDDO 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 ENDIF 199 ! 200 !Initial guess for psi1 201 DO ispin = 1, nspins 202 CALL cp_fm_set_all(psi1(ispin)%matrix, 0.0_dp) 203 !CALL cp_fm_to_fm(p_psi0(ispin,ijdir)%matrix, psi1(ispin)%matrix) 204 !CALL cp_fm_scale(-1.0_dp,psi1(ispin)%matrix) 205 ENDDO 206 ! 207 !DO scf cycle to optimize psi1 208 DO ispin = 1, nspins 209 CALL cp_fm_to_fm(efg_psi0(ispin, ijdir)%matrix, h1_psi0(ispin)%matrix) 210 ENDDO 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)%matrix, psi1_efg(ispin, ijdir)%matrix) 222 CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) 223 chk = chk + fro 224 ENDDO 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)%matrix,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 ENDIF 242 ! 243 ! Write the result in the restart file 244 ENDDO ! jdir 245 ENDDO ! idir 246 ENDIF 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)%matrix, 0.0_dp) 253 ENDDO 254 IF (output_unit > 0) THEN 255 WRITE (output_unit, "(T10,A)") "Response to the perturbation operator pso_"//ACHAR(idir + 119) 256 ENDIF 257 ! 258 !Initial guess for psi1 259 DO ispin = 1, nspins 260 CALL cp_fm_set_all(psi1(ispin)%matrix, 0.0_dp) 261 !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin)%matrix) 262 !CALL cp_fm_scale(-1.0_dp,psi1(ispin)%matrix) 263 ENDDO 264 ! 265 !DO scf cycle to optimize psi1 266 DO ispin = 1, nspins 267 CALL cp_fm_to_fm(pso_psi0(ispin, idir)%matrix, h1_psi0(ispin)%matrix) 268 ENDDO 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)%matrix, psi1_pso(ispin, idir)%matrix) 280 CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) 281 chk = chk + fro 282 ENDDO 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)%matrix,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 ENDIF 300 ! 301 ! Write the result in the restart file 302 ENDDO ! idir 303 ENDIF 304 ! 305 ! operator fc 306 IF (do_fc) THEN 307 DO ispin = 1, nspins 308 CALL cp_fm_set_all(psi1_fc(ispin)%matrix, 0.0_dp) 309 ENDDO 310 IF (output_unit > 0) THEN 311 WRITE (output_unit, "(T10,A)") "Response to the perturbation operator fc" 312 ENDIF 313 ! 314 !Initial guess for psi1 315 DO ispin = 1, nspins 316 CALL cp_fm_set_all(psi1(ispin)%matrix, 0.0_dp) 317 !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin)%matrix) 318 !CALL cp_fm_scale(-1.0_dp,psi1(ispin)%matrix) 319 ENDDO 320 ! 321 !DO scf cycle to optimize psi1 322 DO ispin = 1, nspins 323 CALL cp_fm_to_fm(fc_psi0(ispin)%matrix, h1_psi0(ispin)%matrix) 324 ENDDO 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)%matrix, psi1_fc(ispin)%matrix) 336 CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) 337 chk = chk + fro 338 ENDDO 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)%matrix,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 ENDIF 356 ! 357 ! Write the result in the restart file 358 ENDIF 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)%matrix, 0.0_dp) 367 ENDDO 368 IF (output_unit > 0) THEN 369 WRITE (output_unit, "(T10,A)") "Response to the perturbation operator r_"//ACHAR(idir + 119) 370 ENDIF 371 ! 372 !Initial guess for psi1 373 DO ispin = 1, nspins 374 CALL cp_fm_set_all(psi1(ispin)%matrix, 0.0_dp) 375 !CALL cp_fm_to_fm(rxp_psi0(ispin,idir)%matrix, psi1(ispin)%matrix) 376 !CALL cp_fm_scale(-1.0_dp,psi1(ispin)%matrix) 377 ENDDO 378 ! 379 !DO scf cycle to optimize psi1 380 DO ispin = 1, nspins 381 CALL cp_fm_to_fm(dso_psi0(ispin, idir)%matrix, h1_psi0(ispin)%matrix) 382 ENDDO 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)%matrix, psi1_dso(ispin, idir)%matrix) 394 CALL cp_fm_frobenius_norm(psi1(ispin)%matrix, fro) 395 chk = chk + fro 396 ENDDO 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)%matrix,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 ENDIF 414 ! 415 ! Write the result in the restart file 416 ENDDO ! idir 417 ENDIF 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 ENDIF 426 ! 427 ! 428 ! clean up 429 DO ispin = 1, nspins 430 CALL cp_fm_release(psi1(ispin)%matrix) 431 CALL cp_fm_release(h1_psi0(ispin)%matrix) 432 ENDDO 433 DEALLOCATE (psi1, h1_psi0, psi0_order) 434 ! 435 CALL cp_print_key_finished_output(output_unit, logger, lr_section,& 436 & "PRINT%PROGRAM_RUN_INFO") 437 ! 438 CALL timestop(handle) 439 ! 440 END SUBROUTINE issc_response 441 442! ************************************************************************************************** 443!> \brief ... 444!> \param issc_env ... 445!> \param qs_env ... 446!> \param iatom ... 447! ************************************************************************************************** 448 SUBROUTINE issc_issc(issc_env, qs_env, iatom) 449 450 TYPE(issc_env_type) :: issc_env 451 TYPE(qs_environment_type), POINTER :: qs_env 452 INTEGER, INTENT(IN) :: iatom 453 454 CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_issc', routineP = moduleN//':'//routineN 455 456 INTEGER :: handle, ispin, ixyz, jatom, jxyz, natom, & 457 nmo, nspins 458 LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw 459 REAL(dp) :: buf, facdso, facfc, facpso, facsd, g, & 460 issc_dso, issc_fc, issc_pso, issc_sd, & 461 maxocc 462 REAL(dp), DIMENSION(3) :: r_i, r_j 463 REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc 464 TYPE(cell_type), POINTER :: cell 465 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fc_psi0, psi1_fc 466 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: psi1_dso, psi1_efg, psi1_pso 467 TYPE(cp_fm_type), POINTER :: mo_coeff 468 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, & 469 matrix_pso 470 TYPE(dft_control_type), POINTER :: dft_control 471 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 472 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 473 TYPE(section_vals_type), POINTER :: issc_section 474 475 CALL timeset(routineN, handle) 476 477 NULLIFY (cell, dft_control, particle_set, issc, psi1_fc, psi1_efg, psi1_pso) 478 NULLIFY (matrix_efg, matrix_fc, matrix_pso, mos, mo_coeff, fc_psi0) 479 480 CALL get_qs_env(qs_env=qs_env, & 481 cell=cell, & 482 dft_control=dft_control, & 483 particle_set=particle_set, & 484 mos=mos) 485 486 gapw = dft_control%qs_control%gapw 487 natom = SIZE(particle_set, 1) 488 nspins = dft_control%nspins 489 490 CALL get_issc_env(issc_env=issc_env, & 491 matrix_efg=matrix_efg, & 492 matrix_pso=matrix_pso, & 493 matrix_fc=matrix_fc, & 494 matrix_dso=matrix_dso, & 495 psi1_fc=psi1_fc, & 496 psi1_efg=psi1_efg, & 497 psi1_pso=psi1_pso, & 498 psi1_dso=psi1_dso, & 499 fc_psi0=fc_psi0, & 500 issc=issc, & 501 do_fc=do_fc, & 502 do_sd=do_sd, & 503 do_pso=do_pso, & 504 do_dso=do_dso) 505 506 g = e_mass/(2.0_dp*p_mass) 507 facfc = hertz*g**2*a_fine**4 508 facpso = hertz*g**2*a_fine**4 509 facsd = hertz*g**2*a_fine**4 510 facdso = hertz*g**2*a_fine**4 511 512 ! 513 ! 514 issc_section => section_vals_get_subs_vals(qs_env%input, & 515 & "PROPERTIES%LINRES%SPINSPIN") 516 ! 517 ! Initialize 518 DO ispin = 1, nspins 519 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, maxocc=maxocc) 520 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo) 521 522 DO jatom = 1, natom 523 r_i = particle_set(iatom)%r 524 r_j = particle_set(jatom)%r 525 r_j = pbc(r_i, r_j, cell) + r_i 526 ! 527 ! 528 ! 529 !write(*,*) 'iatom =',iatom,' r_i=',r_i 530 !write(*,*) 'jatom =',jatom,' r_j=',r_j 531 ! 532 ! FC term 533 ! 534 IF (do_fc .AND. iatom .NE. jatom) THEN 535 ! 536 ! build the integral for the jatom 537 CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp) 538 CALL build_fermi_contact_matrix(qs_env, matrix_fc, r_j) 539 CALL cp_dbcsr_sm_fm_multiply(matrix_fc(1)%matrix, mo_coeff, & 540 fc_psi0(ispin)%matrix, ncol=nmo,& ! fc_psi0 a buffer 541 & alpha=1.0_dp) 542 543 CALL cp_fm_trace(fc_psi0(ispin)%matrix, mo_coeff, buf) 544 WRITE (*, *) ' jatom', jatom, 'tr(P*fc)=', buf 545 546 CALL cp_fm_trace(fc_psi0(ispin)%matrix, psi1_fc(ispin)%matrix, buf) 547 issc_fc = 2.0_dp*2.0_dp*maxocc*facfc*buf 548 issc(1, 1, iatom, jatom, 1) = issc(1, 1, iatom, jatom, 1) + issc_fc 549 issc(2, 2, iatom, jatom, 1) = issc(2, 2, iatom, jatom, 1) + issc_fc 550 issc(3, 3, iatom, jatom, 1) = issc(3, 3, iatom, jatom, 1) + issc_fc 551 ENDIF 552 ! 553 ! SD term 554 ! 555 IF (do_sd .AND. iatom .NE. jatom) THEN 556 ! 557 ! build the integral for the jatom 558 CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp) 559 CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp) 560 CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp) 561 CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp) 562 CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp) 563 CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp) 564 CALL build_efg_matrix(qs_env, matrix_efg, r_j) 565 DO ixyz = 1, 6 566 CALL cp_dbcsr_sm_fm_multiply(matrix_efg(ixyz)%matrix, mo_coeff, & 567 fc_psi0(ispin)%matrix, ncol=nmo,& ! fc_psi0 a buffer 568 & alpha=1.0_dp, beta=0.0_dp) 569 CALL cp_fm_trace(fc_psi0(ispin)%matrix, mo_coeff, buf) 570 WRITE (*, *) ' jatom', jatom, ixyz, 'tr(P*efg)=', buf 571 DO jxyz = 1, 6 572 CALL cp_fm_trace(fc_psi0(ispin)%matrix, psi1_efg(ispin, jxyz)%matrix, buf) 573 issc_sd = 2.0_dp*maxocc*facsd*buf 574 !issc(ixyz,jxyz,iatom,jatom) = issc_sd 575 !write(*,*) 'pso_',ixyz,jxyz,' iatom',iatom,'jatom',jatom,issc_pso 576 ENDDO 577 ENDDO 578 ENDIF 579 ! 580 ! PSO term 581 ! 582 IF (do_pso .AND. iatom .NE. jatom) THEN 583 ! 584 ! build the integral for the jatom 585 CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp) 586 CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp) 587 CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp) 588 CALL build_pso_matrix(qs_env, matrix_pso, r_j) 589 DO ixyz = 1, 3 590 CALL cp_dbcsr_sm_fm_multiply(matrix_pso(ixyz)%matrix, mo_coeff, & 591 fc_psi0(ispin)%matrix, ncol=nmo,& ! fc_psi0 a buffer 592 & alpha=1.0_dp, beta=0.0_dp) 593 DO jxyz = 1, 3 594 CALL cp_fm_trace(fc_psi0(ispin)%matrix, psi1_pso(ispin, jxyz)%matrix, buf) 595 issc_pso = -2.0_dp*maxocc*facpso*buf 596 issc(ixyz, jxyz, iatom, jatom, 3) = issc(ixyz, jxyz, iatom, jatom, 3) + issc_pso 597 ENDDO 598 ENDDO 599 ENDIF 600 ! 601 ! DSO term 602 ! 603 !>>>>> for debugging we compute here the polarizability and NOT the DSO term! 604 IF (do_dso .AND. iatom .EQ. natom .AND. jatom .EQ. natom) THEN 605 ! 606 ! build the integral for the jatom 607 !CALL dbcsr_set(matrix_dso(1)%matrix,0.0_dp) 608 !CALL dbcsr_set(matrix_dso(2)%matrix,0.0_dp) 609 !CALL dbcsr_set(matrix_dso(3)%matrix,0.0_dp) 610 !CALL dbcsr_set(matrix_dso(4)%matrix,0.0_dp) 611 !CALL dbcsr_set(matrix_dso(5)%matrix,0.0_dp) 612 !CALL dbcsr_set(matrix_dso(6)%matrix,0.0_dp) 613 !CALL build_dso_matrix(qs_env,matrix_dso,r_i,r_j) 614 !DO ixyz = 1,6 615 ! CALL cp_sm_fm_multiply(matrix_dso(ixyz)%matrix,mo_coeff,& 616 ! & fc_psi0(ispin)%matrix,ncol=nmo,& ! fc_psi0 a buffer 617 ! & alpha=1.0_dp,beta=0.0_dp) 618 ! CALL cp_fm_trace(fc_psi0(ispin)%matrix,mo_coeff,buf) 619 ! issc_dso = 2.0_dp * maxocc * facdso * buf 620 ! issc(ixyz,jxyz,iatom,jatom,4) = issc_dso 621 !ENDDO 622 !CALL dbcsr_set(matrix_dso(1)%matrix,0.0_dp) 623 !CALL dbcsr_set(matrix_dso(2)%matrix,0.0_dp) 624 !CALL dbcsr_set(matrix_dso(3)%matrix,0.0_dp) 625 !CALL rRc_xyz_ao(matrix_dso,qs_env,(/0.0_dp,0.0_dp,0.0_dp/),1) 626 DO ixyz = 1, 3 627 CALL cp_dbcsr_sm_fm_multiply(matrix_dso(ixyz)%matrix, mo_coeff, & 628 fc_psi0(ispin)%matrix, ncol=nmo,& ! fc_psi0 a buffer 629 & alpha=1.0_dp, beta=0.0_dp) 630 DO jxyz = 1, 3 631 CALL cp_fm_trace(psi1_dso(ispin, jxyz)%matrix, fc_psi0(ispin)%matrix, buf) 632 ! we save the polarizability for a checksum later on ! 633 issc_dso = 2.0_dp*maxocc*buf 634 !WRITE(*,*) ixyz,jxyz,'tr(P_r*r)=',2.0_dp * maxocc * buf 635 issc(ixyz, jxyz, iatom, jatom, 4) = issc(ixyz, jxyz, iatom, jatom, 4) + issc_dso 636 ENDDO 637 ENDDO 638 639 ENDIF 640 ! 641 ENDDO ! jatom 642 ENDDO ! ispin 643 ! 644 ! 645 ! Finalize 646 CALL timestop(handle) 647 ! 648 END SUBROUTINE issc_issc 649 650! ************************************************************************************************** 651!> \brief ... 652!> \param issc_env ... 653!> \param qs_env ... 654! ************************************************************************************************** 655 SUBROUTINE issc_print(issc_env, qs_env) 656 TYPE(issc_env_type) :: issc_env 657 TYPE(qs_environment_type), POINTER :: qs_env 658 659 CHARACTER(len=*), PARAMETER :: routineN = 'issc_print', routineP = moduleN//':'//routineN 660 661 CHARACTER(LEN=2) :: element_symbol_i, element_symbol_j 662 CHARACTER(LEN=default_string_length) :: name_i, name_j, title 663 INTEGER :: iatom, jatom, natom, output_unit, & 664 unit_atoms 665 LOGICAL :: do_dso, do_fc, do_pso, do_sd, gapw 666 REAL(dp) :: eig(3), issc_iso_dso, issc_iso_fc, & 667 issc_iso_pso, issc_iso_sd, & 668 issc_iso_tot, issc_tmp(3, 3) 669 REAL(dp), DIMENSION(:, :, :, :, :), POINTER :: issc 670 REAL(dp), EXTERNAL :: DDOT 671 TYPE(atomic_kind_type), POINTER :: atom_kind_i, atom_kind_j 672 TYPE(cp_logger_type), POINTER :: logger 673 TYPE(dft_control_type), POINTER :: dft_control 674 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 675 TYPE(section_vals_type), POINTER :: issc_section 676 677 NULLIFY (logger, particle_set, atom_kind_i, atom_kind_j, dft_control) 678 679 logger => cp_get_default_logger() 680 output_unit = cp_logger_get_default_io_unit(logger) 681 682 issc_section => section_vals_get_subs_vals(qs_env%input, & 683 "PROPERTIES%LINRES%SPINSPIN") 684 685 CALL get_issc_env(issc_env=issc_env, & 686 issc=issc, & 687 do_fc=do_fc, & 688 do_sd=do_sd, & 689 do_pso=do_pso, & 690 do_dso=do_dso) 691 ! 692 CALL get_qs_env(qs_env=qs_env, & 693 dft_control=dft_control, & 694 particle_set=particle_set) 695 696 natom = SIZE(particle_set, 1) 697 gapw = dft_control%qs_control%gapw 698 699 ! 700 IF (output_unit > 0) THEN 701 WRITE (output_unit, '(T2,A,E14.6)') 'ISSC| CheckSum K =', & 702 SQRT(DDOT(SIZE(issc), issc, 1, issc, 1)) 703 ENDIF 704 ! 705 IF (BTEST(cp_print_key_should_output(logger%iter_info, issc_section, & 706 "PRINT%K_MATRIX"), cp_p_file)) THEN 707 708 unit_atoms = cp_print_key_unit_nr(logger, issc_section, "PRINT%K_MATRIX", & 709 extension=".data", middle_name="K", log_filename=.FALSE.) 710 711 IF (unit_atoms > 0) THEN 712 WRITE (unit_atoms, *) 713 WRITE (unit_atoms, *) 714 WRITE (title, '(A)') "Indirect spin-spin coupling matrix" 715 WRITE (unit_atoms, '(T2,A)') title 716 DO iatom = 1, natom 717 atom_kind_i => particle_set(iatom)%atomic_kind 718 CALL get_atomic_kind(atom_kind_i, name=name_i, element_symbol=element_symbol_i) 719 DO jatom = 1, natom 720 atom_kind_j => particle_set(jatom)%atomic_kind 721 CALL get_atomic_kind(atom_kind_j, name=name_j, element_symbol=element_symbol_j) 722 ! 723 IF (iatom .EQ. jatom .AND. .NOT. do_dso) CYCLE 724 ! 725 ! 726 ! FC 727 issc_tmp(:, :) = issc(:, :, iatom, jatom, 1) 728 issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :))) 729 CALL diamat_all(issc_tmp, eig) 730 issc_iso_fc = (eig(1) + eig(2) + eig(3))/3.0_dp 731 ! 732 ! SD 733 issc_tmp(:, :) = issc(:, :, iatom, jatom, 2) 734 issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :))) 735 CALL diamat_all(issc_tmp, eig) 736 issc_iso_sd = (eig(1) + eig(2) + eig(3))/3.0_dp 737 ! 738 ! PSO 739 issc_tmp(:, :) = issc(:, :, iatom, jatom, 3) 740 issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :))) 741 CALL diamat_all(issc_tmp, eig) 742 issc_iso_pso = (eig(1) + eig(2) + eig(3))/3.0_dp 743 ! 744 ! DSO 745 issc_tmp(:, :) = issc(:, :, iatom, jatom, 4) 746 issc_tmp(:, :) = 0.5_dp*(issc_tmp(:, :) + TRANSPOSE(issc_tmp(:, :))) 747 CALL diamat_all(issc_tmp, eig) 748 issc_iso_dso = (eig(1) + eig(2) + eig(3))/3.0_dp 749 ! 750 ! TOT 751 issc_iso_tot = issc_iso_fc + issc_iso_sd + issc_iso_dso + issc_iso_pso 752 ! 753 ! 754 WRITE (unit_atoms, *) 755 WRITE (unit_atoms, '(T2,2(A,I5,A,2X,A2))') 'Indirect spin-spin coupling between ', & 756 iatom, TRIM(name_i), element_symbol_i, ' and ', & 757 jatom, TRIM(name_j), element_symbol_j 758 ! 759 IF (do_fc) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic FC contribution = ', issc_iso_fc, ' Hz' 760 IF (do_sd) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic SD contribution = ', issc_iso_sd, ' Hz' 761 IF (do_pso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic PSO contribution = ', issc_iso_pso, ' Hz' 762 !IF(do_dso) WRITE(unit_atoms,'(T1,A,f12.4,A)') ' Isotropic DSO contribution = ',issc_iso_dso,' Hz' 763 IF (do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' !!! POLARIZABILITY (for the moment) = ', issc_iso_dso, ' Hz' 764 IF (.NOT. do_dso) WRITE (unit_atoms, '(T1,A,f12.4,A)') ' Isotropic coupling = ', issc_iso_tot, ' Hz' 765 ENDDO 766 ENDDO 767 ENDIF 768 CALL cp_print_key_finished_output(unit_atoms, logger, issc_section,& 769 & "PRINT%K_MATRIX") 770 ENDIF 771 ! 772 ! 773 END SUBROUTINE issc_print 774 775! ************************************************************************************************** 776!> \brief Initialize the issc environment 777!> \param issc_env ... 778!> \param qs_env ... 779! ************************************************************************************************** 780 SUBROUTINE issc_env_init(issc_env, qs_env) 781 ! 782 TYPE(issc_env_type) :: issc_env 783 TYPE(qs_environment_type), POINTER :: qs_env 784 785 CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_env_init', routineP = moduleN//':'//routineN 786 787 INTEGER :: handle, iatom, idir, ini, ir, ispin, & 788 istat, m, n, n_rep, nao, natom, & 789 nspins, output_unit 790 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 791 INTEGER, DIMENSION(:), POINTER :: list, row_blk_sizes 792 LOGICAL :: gapw 793 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 794 TYPE(cp_fm_type), POINTER :: mo_coeff 795 TYPE(cp_logger_type), POINTER :: logger 796 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 797 TYPE(dft_control_type), POINTER :: dft_control 798 TYPE(linres_control_type), POINTER :: linres_control 799 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 800 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 801 POINTER :: sab_orb 802 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 803 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 804 TYPE(section_vals_type), POINTER :: issc_section, lr_section 805 806! 807 808 CALL timeset(routineN, handle) 809 810 NULLIFY (linres_control) 811 NULLIFY (logger, issc_section) 812 NULLIFY (tmp_fm_struct) 813 NULLIFY (particle_set, qs_kind_set) 814 NULLIFY (sab_orb) 815 816 logger => cp_get_default_logger() 817 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 818 819 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 820 extension=".linresLog") 821 822 IF (issc_env%ref_count /= 0) THEN 823 CALL issc_env_cleanup(issc_env) 824 ENDIF 825 826 IF (output_unit > 0) THEN 827 WRITE (output_unit, "(/,T20,A,/)") "*** Start indirect spin-spin coupling Calculation ***" 828 WRITE (output_unit, "(T10,A,/)") "Inizialization of the ISSC environment" 829 ENDIF 830 831 CALL issc_env_create(issc_env) 832 ! 833 issc_section => section_vals_get_subs_vals(qs_env%input, & 834 & "PROPERTIES%LINRES%SPINSPIN") 835 !CALL section_vals_val_get(nmr_section,"INTERPOLATE_SHIFT",l_val=nmr_env%interpolate_shift) 836 !CALL section_vals_val_get(nmr_section,"SHIFT_GAPW_RADIUS",r_val=nmr_env%shift_gapw_radius) 837 838 CALL get_qs_env(qs_env=qs_env, & 839 dft_control=dft_control, & 840 linres_control=linres_control, & 841 mos=mos, & 842 sab_orb=sab_orb, & 843 particle_set=particle_set, & 844 qs_kind_set=qs_kind_set, & 845 dbcsr_dist=dbcsr_dist) 846 ! 847 ! 848 gapw = dft_control%qs_control%gapw 849 nspins = dft_control%nspins 850 natom = SIZE(particle_set, 1) 851 ! 852 ! check that the psi0 are localized and you have all the centers 853 IF (.NOT. linres_control%localized_psi0) & 854 CALL cp_warn(__LOCATION__, 'To get indirect spin-spin coupling parameters within '// & 855 'PBC you need to localize zero order orbitals') 856 ! 857 ! 858 ! read terms need to be calculated 859 ! FC 860 CALL section_vals_val_get(issc_section, "DO_FC", l_val=issc_env%do_fc) 861 ! SD 862 CALL section_vals_val_get(issc_section, "DO_SD", l_val=issc_env%do_sd) 863 ! PSO 864 CALL section_vals_val_get(issc_section, "DO_PSO", l_val=issc_env%do_pso) 865 ! DSO 866 CALL section_vals_val_get(issc_section, "DO_DSO", l_val=issc_env%do_dso) 867 ! 868 ! 869 ! read the list of atoms on which the issc need to be calculated 870 CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", n_rep_val=n_rep) 871 ! 872 ! 873 NULLIFY (issc_env%issc_on_atom_list) 874 n = 0 875 DO ir = 1, n_rep 876 NULLIFY (list) 877 CALL section_vals_val_get(issc_section, "ISSC_ON_ATOM_LIST", i_rep_val=ir, i_vals=list) 878 IF (ASSOCIATED(list)) THEN 879 CALL reallocate(issc_env%issc_on_atom_list, 1, n + SIZE(list)) 880 DO ini = 1, SIZE(list) 881 issc_env%issc_on_atom_list(ini + n) = list(ini) 882 ENDDO 883 n = n + SIZE(list) 884 ENDIF 885 ENDDO 886 ! 887 IF (.NOT. ASSOCIATED(issc_env%issc_on_atom_list)) THEN 888 ALLOCATE (issc_env%issc_on_atom_list(natom), STAT=istat) 889 CPASSERT(istat .EQ. 0) 890 DO iatom = 1, natom 891 issc_env%issc_on_atom_list(iatom) = iatom 892 ENDDO 893 ENDIF 894 issc_env%issc_natms = SIZE(issc_env%issc_on_atom_list) 895 ! 896 ! 897 ! Initialize the issc tensor 898 ALLOCATE (issc_env%issc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4), & 899 issc_env%issc_loc(3, 3, issc_env%issc_natms, issc_env%issc_natms, 4), & 900 STAT=istat) 901 CPASSERT(istat == 0) 902 issc_env%issc(:, :, :, :, :) = 0.0_dp 903 issc_env%issc_loc(:, :, :, :, :) = 0.0_dp 904 ! 905 ! allocation 906 ALLOCATE (issc_env%efg_psi0(nspins, 6), issc_env%pso_psi0(nspins, 3), issc_env%fc_psi0(nspins), & 907 issc_env%psi1_efg(nspins, 6), issc_env%psi1_pso(nspins, 3), issc_env%psi1_fc(nspins), & 908 issc_env%dso_psi0(nspins, 3), issc_env%psi1_dso(nspins, 3), & 909 STAT=istat) 910 CPASSERT(istat == 0) 911 DO ispin = 1, nspins 912 !mo_coeff => current_env%psi0_order(ispin)%matrix 913 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 914 CALL cp_fm_get_info(mo_coeff, ncol_global=m, nrow_global=nao) 915 916 NULLIFY (tmp_fm_struct) 917 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, & 918 ncol_global=m, & 919 context=mo_coeff%matrix_struct%context) 920 DO idir = 1, 6 921 NULLIFY (issc_env%psi1_efg(ispin, idir)%matrix, issc_env%efg_psi0(ispin, idir)%matrix) 922 CALL cp_fm_create(issc_env%psi1_efg(ispin, idir)%matrix, tmp_fm_struct) 923 CALL cp_fm_create(issc_env%efg_psi0(ispin, idir)%matrix, tmp_fm_struct) 924 ENDDO 925 DO idir = 1, 3 926 NULLIFY (issc_env%psi1_pso(ispin, idir)%matrix, issc_env%pso_psi0(ispin, idir)%matrix, & 927 issc_env%dso_psi0(ispin, idir)%matrix) 928 CALL cp_fm_create(issc_env%psi1_pso(ispin, idir)%matrix, tmp_fm_struct) 929 CALL cp_fm_create(issc_env%pso_psi0(ispin, idir)%matrix, tmp_fm_struct) 930 CALL cp_fm_create(issc_env%psi1_dso(ispin, idir)%matrix, tmp_fm_struct) 931 CALL cp_fm_create(issc_env%dso_psi0(ispin, idir)%matrix, tmp_fm_struct) 932 ENDDO 933 NULLIFY (issc_env%psi1_fc(ispin)%matrix, issc_env%fc_psi0(ispin)%matrix) 934 CALL cp_fm_create(issc_env%psi1_fc(ispin)%matrix, tmp_fm_struct) 935 CALL cp_fm_create(issc_env%fc_psi0(ispin)%matrix, tmp_fm_struct) 936 CALL cp_fm_struct_release(tmp_fm_struct) 937 ENDDO 938 ! 939 ! prepare for allocation 940 ALLOCATE (first_sgf(natom)) 941 ALLOCATE (last_sgf(natom)) 942 CALL get_particle_set(particle_set, qs_kind_set, & 943 first_sgf=first_sgf, & 944 last_sgf=last_sgf) 945 ALLOCATE (row_blk_sizes(natom)) 946 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf) 947 DEALLOCATE (first_sgf) 948 DEALLOCATE (last_sgf) 949 950 ! 951 ! efg, pso and fc operators 952 CALL dbcsr_allocate_matrix_set(issc_env%matrix_efg, 6) 953 ALLOCATE (issc_env%matrix_efg(1)%matrix) 954 CALL dbcsr_create(matrix=issc_env%matrix_efg(1)%matrix, & 955 name="efg (3xx-rr)/3", & 956 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 957 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 958 nze=0, mutable_work=.TRUE.) 959 CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_efg(1)%matrix, sab_orb) 960 961 ALLOCATE (issc_env%matrix_efg(2)%matrix, & 962 issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(4)%matrix, & 963 issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(6)%matrix) 964 CALL dbcsr_copy(issc_env%matrix_efg(2)%matrix, issc_env%matrix_efg(1)%matrix, & 965 'efg xy') 966 CALL dbcsr_copy(issc_env%matrix_efg(3)%matrix, issc_env%matrix_efg(1)%matrix, & 967 'efg xz') 968 CALL dbcsr_copy(issc_env%matrix_efg(4)%matrix, issc_env%matrix_efg(1)%matrix, & 969 'efg (3yy-rr)/3') 970 CALL dbcsr_copy(issc_env%matrix_efg(5)%matrix, issc_env%matrix_efg(1)%matrix, & 971 'efg yz') 972 CALL dbcsr_copy(issc_env%matrix_efg(6)%matrix, issc_env%matrix_efg(1)%matrix, & 973 'efg (3zz-rr)/3') 974 975 CALL dbcsr_set(issc_env%matrix_efg(1)%matrix, 0.0_dp) 976 CALL dbcsr_set(issc_env%matrix_efg(2)%matrix, 0.0_dp) 977 CALL dbcsr_set(issc_env%matrix_efg(3)%matrix, 0.0_dp) 978 CALL dbcsr_set(issc_env%matrix_efg(4)%matrix, 0.0_dp) 979 CALL dbcsr_set(issc_env%matrix_efg(5)%matrix, 0.0_dp) 980 CALL dbcsr_set(issc_env%matrix_efg(6)%matrix, 0.0_dp) 981 ! 982 ! PSO 983 CALL dbcsr_allocate_matrix_set(issc_env%matrix_pso, 3) 984 ALLOCATE (issc_env%matrix_pso(1)%matrix) 985 CALL dbcsr_create(matrix=issc_env%matrix_pso(1)%matrix, & 986 name="pso x", & 987 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, & 988 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 989 nze=0, mutable_work=.TRUE.) 990 CALL cp_dbcsr_alloc_block_from_nbl(issc_env%matrix_pso(1)%matrix, sab_orb) 991 992 ALLOCATE (issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(3)%matrix) 993 CALL dbcsr_copy(issc_env%matrix_pso(2)%matrix, issc_env%matrix_pso(1)%matrix, & 994 'pso y') 995 CALL dbcsr_copy(issc_env%matrix_pso(3)%matrix, issc_env%matrix_pso(1)%matrix, & 996 'pso z') 997 CALL dbcsr_set(issc_env%matrix_pso(1)%matrix, 0.0_dp) 998 CALL dbcsr_set(issc_env%matrix_pso(2)%matrix, 0.0_dp) 999 CALL dbcsr_set(issc_env%matrix_pso(3)%matrix, 0.0_dp) 1000 ! 1001 ! DSO 1002 CALL dbcsr_allocate_matrix_set(issc_env%matrix_dso, 3) 1003 ALLOCATE (issc_env%matrix_dso(1)%matrix, issc_env%matrix_dso(2)%matrix, issc_env%matrix_dso(3)%matrix) 1004 CALL dbcsr_copy(issc_env%matrix_dso(1)%matrix, issc_env%matrix_efg(1)%matrix, & 1005 'dso x') 1006 CALL dbcsr_copy(issc_env%matrix_dso(2)%matrix, issc_env%matrix_efg(1)%matrix, & 1007 'dso y') 1008 CALL dbcsr_copy(issc_env%matrix_dso(3)%matrix, issc_env%matrix_efg(1)%matrix, & 1009 'dso z') 1010 CALL dbcsr_set(issc_env%matrix_dso(1)%matrix, 0.0_dp) 1011 CALL dbcsr_set(issc_env%matrix_dso(2)%matrix, 0.0_dp) 1012 CALL dbcsr_set(issc_env%matrix_dso(3)%matrix, 0.0_dp) 1013 ! 1014 ! FC 1015 CALL dbcsr_allocate_matrix_set(issc_env%matrix_fc, 1) 1016 ALLOCATE (issc_env%matrix_fc(1)%matrix) 1017 CALL dbcsr_copy(issc_env%matrix_fc(1)%matrix, issc_env%matrix_efg(1)%matrix, & 1018 'fc') 1019 CALL dbcsr_set(issc_env%matrix_fc(1)%matrix, 0.0_dp) 1020 1021 DEALLOCATE (row_blk_sizes) 1022 ! 1023 ! Conversion factors 1024 IF (output_unit > 0) THEN 1025 WRITE (output_unit, "(T2,A,T60,I4,A)")& 1026 & "ISSC| spin-spin coupling computed for ", issc_env%issc_natms, ' atoms' 1027 ENDIF 1028 1029 CALL cp_print_key_finished_output(output_unit, logger, lr_section,& 1030 & "PRINT%PROGRAM_RUN_INFO") 1031 1032 CALL timestop(handle) 1033 1034 END SUBROUTINE issc_env_init 1035 1036! ************************************************************************************************** 1037!> \brief Deallocate the issc environment 1038!> \param issc_env ... 1039!> \par History 1040! ************************************************************************************************** 1041 SUBROUTINE issc_env_cleanup(issc_env) 1042 1043 TYPE(issc_env_type) :: issc_env 1044 1045 CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_env_cleanup', & 1046 routineP = moduleN//':'//routineN 1047 1048 INTEGER :: idir, ispin 1049 1050 issc_env%ref_count = issc_env%ref_count - 1 1051 IF (issc_env%ref_count == 0) THEN 1052 IF (ASSOCIATED(issc_env%issc_on_atom_list)) THEN 1053 DEALLOCATE (issc_env%issc_on_atom_list) 1054 ENDIF 1055 IF (ASSOCIATED(issc_env%issc)) THEN 1056 DEALLOCATE (issc_env%issc) 1057 ENDIF 1058 IF (ASSOCIATED(issc_env%issc_loc)) THEN 1059 DEALLOCATE (issc_env%issc_loc) 1060 ENDIF 1061 ! 1062 !efg_psi0 1063 IF (ASSOCIATED(issc_env%efg_psi0)) THEN 1064 DO idir = 1, SIZE(issc_env%efg_psi0, 2) 1065 DO ispin = 1, SIZE(issc_env%efg_psi0, 1) 1066 CALL cp_fm_release(issc_env%efg_psi0(ispin, idir)%matrix) 1067 ENDDO 1068 ENDDO 1069 DEALLOCATE (issc_env%efg_psi0) 1070 ENDIF 1071 ! 1072 !pso_psi0 1073 IF (ASSOCIATED(issc_env%pso_psi0)) THEN 1074 DO idir = 1, SIZE(issc_env%pso_psi0, 2) 1075 DO ispin = 1, SIZE(issc_env%pso_psi0, 1) 1076 CALL cp_fm_release(issc_env%pso_psi0(ispin, idir)%matrix) 1077 ENDDO 1078 ENDDO 1079 DEALLOCATE (issc_env%pso_psi0) 1080 ENDIF 1081 ! 1082 !dso_psi0 1083 IF (ASSOCIATED(issc_env%dso_psi0)) THEN 1084 DO idir = 1, SIZE(issc_env%dso_psi0, 2) 1085 DO ispin = 1, SIZE(issc_env%dso_psi0, 1) 1086 CALL cp_fm_release(issc_env%dso_psi0(ispin, idir)%matrix) 1087 ENDDO 1088 ENDDO 1089 DEALLOCATE (issc_env%dso_psi0) 1090 ENDIF 1091 ! 1092 !fc_psi0 1093 IF (ASSOCIATED(issc_env%fc_psi0)) THEN 1094 DO ispin = 1, SIZE(issc_env%fc_psi0, 1) 1095 CALL cp_fm_release(issc_env%fc_psi0(ispin)%matrix) 1096 ENDDO 1097 DEALLOCATE (issc_env%fc_psi0) 1098 ENDIF 1099 ! 1100 !psi1_efg 1101 IF (ASSOCIATED(issc_env%psi1_efg)) THEN 1102 DO idir = 1, SIZE(issc_env%psi1_efg, 2) 1103 DO ispin = 1, SIZE(issc_env%psi1_efg, 1) 1104 CALL cp_fm_release(issc_env%psi1_efg(ispin, idir)%matrix) 1105 ENDDO 1106 ENDDO 1107 DEALLOCATE (issc_env%psi1_efg) 1108 ENDIF 1109 ! 1110 !psi1_pso 1111 IF (ASSOCIATED(issc_env%psi1_pso)) THEN 1112 DO idir = 1, SIZE(issc_env%psi1_pso, 2) 1113 DO ispin = 1, SIZE(issc_env%psi1_pso, 1) 1114 CALL cp_fm_release(issc_env%psi1_pso(ispin, idir)%matrix) 1115 ENDDO 1116 ENDDO 1117 DEALLOCATE (issc_env%psi1_pso) 1118 ENDIF 1119 ! 1120 !psi1_dso 1121 IF (ASSOCIATED(issc_env%psi1_dso)) THEN 1122 DO idir = 1, SIZE(issc_env%psi1_dso, 2) 1123 DO ispin = 1, SIZE(issc_env%psi1_dso, 1) 1124 CALL cp_fm_release(issc_env%psi1_dso(ispin, idir)%matrix) 1125 ENDDO 1126 ENDDO 1127 DEALLOCATE (issc_env%psi1_dso) 1128 ENDIF 1129 ! 1130 !psi1_fc 1131 IF (ASSOCIATED(issc_env%psi1_fc)) THEN 1132 DO ispin = 1, SIZE(issc_env%psi1_fc, 1) 1133 CALL cp_fm_release(issc_env%psi1_fc(ispin)%matrix) 1134 ENDDO 1135 DEALLOCATE (issc_env%psi1_fc) 1136 ENDIF 1137 ! 1138 ! cubes 1139 !IF(ASSOCIATED(issc_env%list_cubes)) THEN 1140 ! DEALLOCATE(issc_env%list_cubes,STAT=istat) 1141 ! CPPostcondition(istat==0,cp_failure_level,routineP,failure) 1142 !ENDIF 1143 ! 1144 !matrix_efg 1145 IF (ASSOCIATED(issc_env%matrix_efg)) THEN 1146 CALL dbcsr_deallocate_matrix_set(issc_env%matrix_efg) 1147 ENDIF 1148 ! 1149 !matrix_pso 1150 IF (ASSOCIATED(issc_env%matrix_pso)) THEN 1151 CALL dbcsr_deallocate_matrix_set(issc_env%matrix_pso) 1152 ENDIF 1153 ! 1154 !matrix_dso 1155 IF (ASSOCIATED(issc_env%matrix_dso)) THEN 1156 CALL dbcsr_deallocate_matrix_set(issc_env%matrix_dso) 1157 ENDIF 1158 ! 1159 !matrix_fc 1160 IF (ASSOCIATED(issc_env%matrix_fc)) THEN 1161 CALL dbcsr_deallocate_matrix_set(issc_env%matrix_fc) 1162 ENDIF 1163 1164 ENDIF ! ref count 1165 1166 END SUBROUTINE issc_env_cleanup 1167 1168END MODULE qs_linres_issc_utils 1169