1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Distribution of the spin orbit integral matrix. 8!> \par History 9!> \author VW (27.02.2009) 10! ************************************************************************************************** 11MODULE qs_spin_orbit 12 USE ai_spin_orbit, ONLY: pso 13 USE basis_set_types, ONLY: gto_basis_set_p_type,& 14 gto_basis_set_type 15 USE block_p_types, ONLY: block_p_type 16 USE cell_types, ONLY: cell_type,& 17 pbc 18 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 19 USE cp_log_handling, ONLY: cp_get_default_logger,& 20 cp_logger_type 21 USE cp_output_handling, ONLY: cp_p_file,& 22 cp_print_key_finished_output,& 23 cp_print_key_should_output,& 24 cp_print_key_unit_nr 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 27 dbcsr_p_type 28 USE input_section_types, ONLY: section_vals_val_get 29 USE kinds, ONLY: dp 30 USE orbital_pointers, ONLY: init_orbital_pointers,& 31 ncoset 32 USE particle_types, ONLY: particle_type 33 USE qs_environment_types, ONLY: get_qs_env,& 34 qs_environment_type 35 USE qs_kind_types, ONLY: get_qs_kind,& 36 get_qs_kind_set,& 37 qs_kind_type 38 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 39 neighbor_list_iterate,& 40 neighbor_list_iterator_create,& 41 neighbor_list_iterator_p_type,& 42 neighbor_list_iterator_release,& 43 neighbor_list_set_p_type 44#include "./base/base_uses.f90" 45 46 IMPLICIT NONE 47 48 PRIVATE 49 50! *** Global parameters *** 51 52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_spin_orbit' 53 54! *** Public subroutines *** 55 56 PUBLIC :: build_pso_matrix 57 58CONTAINS 59 60! ************************************************************************************************** 61!> \brief Calculation of the paramagnetic spin orbit matrix over 62!> Cartesian Gaussian functions. 63!> \param qs_env ... 64!> \param matrix_so ... 65!> \param rc ... 66!> \date 27.02.2009 67!> \author VW 68!> \version 1.0 69! ************************************************************************************************** 70 71 SUBROUTINE build_pso_matrix(qs_env, matrix_so, rc) 72 73 TYPE(qs_environment_type), POINTER :: qs_env 74 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_so 75 REAL(dp), DIMENSION(3), INTENT(IN) :: rc 76 77 CHARACTER(len=*), PARAMETER :: routineN = 'build_pso_matrix', & 78 routineP = moduleN//':'//routineN 79 80 INTEGER :: after, handle, i, iatom, icol, ikind, inode, irow, iset, iw, jatom, jkind, jset, & 81 last_jatom, ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, & 82 nseta, nsetb, sgfa, sgfb 83 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 84 npgfb, nsgfa, nsgfb 85 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 86 LOGICAL :: found, new_atom_b, omit_headers 87 REAL(KIND=dp) :: dab, rab2 88 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 89 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr_work, soab 90 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc 91 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 92 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 93 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: soint 94 TYPE(cell_type), POINTER :: cell 95 TYPE(cp_logger_type), POINTER :: logger 96 TYPE(cp_para_env_type), POINTER :: para_env 97 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 98 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 99 TYPE(neighbor_list_iterator_p_type), & 100 DIMENSION(:), POINTER :: nl_iterator 101 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 102 POINTER :: sab_orb 103 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 104 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 105 TYPE(qs_kind_type), POINTER :: qs_kind 106 107 CALL timeset(routineN, handle) 108 109 NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env) 110 NULLIFY (logger) 111 112 logger => cp_get_default_logger() 113 114 CALL get_qs_env(qs_env=qs_env, & 115 qs_kind_set=qs_kind_set, & 116 particle_set=particle_set, & 117 neighbor_list_id=neighbor_list_id, & 118 para_env=para_env, & 119 sab_orb=sab_orb, & 120 cell=cell) 121 122 nkind = SIZE(qs_kind_set) 123 natom = SIZE(particle_set) 124 125! *** Allocate work storage *** 126 127 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 128 maxco=maxco, & 129 maxlgto=maxlgto, & 130 maxsgf=maxsgf) 131 132 ldai = ncoset(maxlgto + 1) 133 CALL init_orbital_pointers(ldai) 134 135 ALLOCATE (rr_work(0:2*maxlgto + 2, ldai, ldai)) 136 ALLOCATE (soab(maxco, maxco, 3)) 137 ALLOCATE (work(maxco, maxsgf)) 138 ALLOCATE (soint(3)) 139 140 rr_work(:, :, :) = 0.0_dp 141 soab(:, :, :) = 0.0_dp 142 work(:, :) = 0.0_dp 143 144 ALLOCATE (basis_set_list(nkind)) 145 DO ikind = 1, nkind 146 qs_kind => qs_kind_set(ikind) 147 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 148 IF (ASSOCIATED(basis_set_a)) THEN 149 basis_set_list(ikind)%gto_basis_set => basis_set_a 150 ELSE 151 NULLIFY (basis_set_list(ikind)%gto_basis_set) 152 END IF 153 END DO 154 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 155 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 156 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 157 iatom=iatom, jatom=jatom, r=rab) 158 basis_set_a => basis_set_list(ikind)%gto_basis_set 159 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 160 basis_set_b => basis_set_list(jkind)%gto_basis_set 161 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 162 ra = pbc(particle_set(iatom)%r, cell) 163 ! basis ikind 164 first_sgfa => basis_set_a%first_sgf 165 la_max => basis_set_a%lmax 166 la_min => basis_set_a%lmin 167 npgfa => basis_set_a%npgf 168 nseta = basis_set_a%nset 169 nsgfa => basis_set_a%nsgf_set 170 rpgfa => basis_set_a%pgf_radius 171 set_radius_a => basis_set_a%set_radius 172 sphi_a => basis_set_a%sphi 173 zeta => basis_set_a%zet 174 ! basis jkind 175 first_sgfb => basis_set_b%first_sgf 176 lb_max => basis_set_b%lmax 177 lb_min => basis_set_b%lmin 178 npgfb => basis_set_b%npgf 179 nsetb = basis_set_b%nset 180 nsgfb => basis_set_b%nsgf_set 181 rpgfb => basis_set_b%pgf_radius 182 set_radius_b => basis_set_b%set_radius 183 sphi_b => basis_set_b%sphi 184 zetb => basis_set_b%zet 185 186 IF (inode == 1) last_jatom = 0 187 188 rb = rab + ra 189 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 190 dab = SQRT(rab2) 191 rac = pbc(ra, rc, cell) 192 rbc = rac - rab 193 194 IF (jatom /= last_jatom) THEN 195 new_atom_b = .TRUE. 196 last_jatom = jatom 197 ELSE 198 new_atom_b = .FALSE. 199 END IF 200 201 IF (new_atom_b) THEN 202 IF (iatom <= jatom) THEN 203 irow = iatom 204 icol = jatom 205 ELSE 206 irow = jatom 207 icol = iatom 208 END IF 209 210 DO i = 1, 3 211 NULLIFY (soint(i)%block) 212 CALL dbcsr_get_block_p(matrix=matrix_so(i)%matrix, & 213 row=irow, col=icol, BLOCK=soint(i)%block, found=found) 214 ENDDO 215 END IF 216 217 DO iset = 1, nseta 218 219 ncoa = npgfa(iset)*ncoset(la_max(iset)) 220 sgfa = first_sgfa(1, iset) 221 222 DO jset = 1, nsetb 223 224 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 225 226 ncob = npgfb(jset)*ncoset(lb_max(jset)) 227 sgfb = first_sgfb(1, jset) 228 229 ! *** Calculate the primitive fermi contact integrals *** 230 231 CALL pso(la_max(iset), la_min(iset), npgfa(iset), & 232 rpgfa(:, iset), zeta(:, iset), & 233 lb_max(jset), lb_min(jset), npgfb(jset), & 234 rpgfb(:, jset), zetb(:, jset), & 235 rac, rbc, rab, soab, SIZE(rr_work, 1), SIZE(rr_work, 2), rr_work) 236 237 ! *** Contraction step *** 238 239 DO i = 1, 3 240 241 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 242 1.0_dp, soab(1, 1, i), SIZE(soab, 1), & 243 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 244 0.0_dp, work(1, 1), SIZE(work, 1)) 245 246 IF (iatom <= jatom) THEN 247 248 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 249 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 250 work(1, 1), SIZE(work, 1), & 251 1.0_dp, soint(i)%block(sgfa, sgfb), & 252 SIZE(soint(i)%block, 1)) 253 254 ELSE 255 256 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 257 -1.0_dp, work(1, 1), SIZE(work, 1), & 258 sphi_a(1, sgfa), SIZE(sphi_a, 1), & 259 1.0_dp, soint(i)%block(sgfb, sgfa), & 260 SIZE(soint(i)%block, 1)) 261 ENDIF 262 263 ENDDO 264 265 ENDDO 266 267 ENDDO 268 269 END DO 270 CALL neighbor_list_iterator_release(nl_iterator) 271 272 ! *** Release work storage *** 273 274 DEALLOCATE (basis_set_list) 275 276 DEALLOCATE (soab) 277 278 DEALLOCATE (work) 279 280 DEALLOCATE (soint) 281 282! *** Print the spin orbit matrix, if requested *** 283 284 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 285 qs_env%input, "DFT%PRINT%AO_MATRICES/PSO"), cp_p_file)) THEN 286 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/PSO", & 287 extension=".Log") 288 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 289 after = MIN(MAX(after, 1), 16) 290 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 291 CALL cp_dbcsr_write_sparse_matrix(matrix_so(1)%matrix, 4, after, qs_env, & 292 para_env, output_unit=iw, omit_headers=omit_headers) 293 CALL cp_dbcsr_write_sparse_matrix(matrix_so(2)%matrix, 4, after, qs_env, & 294 para_env, output_unit=iw, omit_headers=omit_headers) 295 CALL cp_dbcsr_write_sparse_matrix(matrix_so(3)%matrix, 4, after, qs_env, & 296 para_env, output_unit=iw, omit_headers=omit_headers) 297 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 298 "DFT%PRINT%AO_MATRICES/PSO") 299 END IF 300 301 CALL timestop(handle) 302 303 END SUBROUTINE build_pso_matrix 304 305! ************************************************************************************************** 306 307END MODULE qs_spin_orbit 308 309