XAS_TDP| switch to proper DBCSR symmetry types + remove some tmp arrays

src/xas_tdp_atom.F

@@ -37,7 +37,8 @@ MODULE xas_tdp_atom
         dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
         dbcsr_get_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
-        dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_replicate_all, dbcsr_type
+        dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_replicate_all, dbcsr_type, &
+        dbcsr_type_antisymmetric, dbcsr_type_symmetric
    USE dbcsr_tensor_api,                ONLY: dbcsr_t_destroy,&
                                               dbcsr_t_get_block,&
                                               dbcsr_t_iterator_blocks_left,&
@@ -491,9 +492,9 @@ SUBROUTINE compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
                      factor = orbtramat(l)%s2c(iso, ico) &
                               *SQRT(4.0_dp*pi/dfac(2*l + 1)*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
 
-                     CALL daxpy(nsgf_set(iset), factor, &
-                                sphi(start_c + ncoset(l - 1) + ico, sgfi:sgfi + nsgf_set(iset) - 1), 1, &
-                                sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1), 1)
+                     sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) = &
+                        sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) &
+                        + factor*sphi(start_c + ncoset(l - 1) + ico, sgfi:sgfi + nsgf_set(iset) - 1)
 
                   END DO !ico
                END DO !iso
@@ -735,8 +736,8 @@ SUBROUTINE get_exat_ri_sinv(ri_sinv, whole_s, neighbors, idx_to_nb, basis_set_ri
       CALL dbcsr_distribution_new(sinv_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
                                   col_dist=col_dist)
 
-      CALL dbcsr_create(matrix=ri_sinv, name="RI_SINV", matrix_type="S", dist=sinv_dist, &
-                        row_blk_size=nsgf, col_blk_size=nsgf)
+      CALL dbcsr_create(matrix=ri_sinv, name="RI_SINV", matrix_type=dbcsr_type_symmetric, &
+                        dist=sinv_dist, row_blk_size=nsgf, col_blk_size=nsgf)
       !reserving the blocks in the correct pattern
       DO inb = 1, nnb
          ri = pbc(particle_set(neighbors(inb))%r, cell)
@@ -1086,11 +1087,11 @@ SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
 
                   IF (sinv_found) THEN
                      DO i = 1, blk_size_ri(katom)
-                        CALL daxpy(SIZE(work2), factor*work1(i, 1), sinv_block(i, :), 1, work2(:), 1)
+                        work2(:) = work2(:) + factor*work1(i, 1)*sinv_block(i, :)
                      END DO
                   ELSE
                      DO i = 1, blk_size_ri(katom)
-                        CALL daxpy(SIZE(work2), factor*work1(i, 1), sinv_blockt(:, i), 1, work2(:), 1)
+                        work2(:) = work2(:) + factor*work1(i, 1)*sinv_blockt(:, i)
                      END DO
                   END IF
 
@@ -1172,7 +1173,7 @@ SUBROUTINE put_density_on_atomic_grid(rho_set, ri_dcoeff, atom_kind, do_gga, bat
                                                             nsoi, nspins, sgfi, starti
       INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nsgf_set
       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf
-      REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: so, work1, work2
+      REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: so
       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :)          :: dso
       REAL(dp), DIMENSION(:, :), POINTER                 :: dgr1, dgr2, ga, gr, ri_sphi_so, zet
       REAL(dp), DIMENSION(:, :, :), POINTER              :: dga1, dga2, rhoa, rhob
@@ -1237,16 +1238,12 @@ SUBROUTINE put_density_on_atomic_grid(rho_set, ri_dcoeff, atom_kind, do_gga, bat
             sgfi = first_sgf(1, iset)
             nsoi = npgf(iset)*nsoset(lmax(iset))
             nsgfi = nsgf_set(iset)
-            ALLOCATE (work1(nsoi, 1), work2(nsgfi, 1))
 
-            work2(:, 1) = ri_dcoeff(ispin)%array(sgfi:sgfi + nsgfi - 1)
-            CALL dgemm('N', 'N', nsoi, 1, nsgfi, 1.0_dp, ri_sphi_so(1:nsoi, sgfi:sgfi + nsgfi - 1), &
-                       nsoi, work2, nsgfi, 0.0_dp, work1, nsoi)
-            ri_dcoeff_so(ispin)%array((iset - 1)*maxso + 1:(iset - 1)*maxso + nsoi) = work1(:, 1)
+            CALL dgemv('N', nsoi, nsgfi, 1.0_dp, ri_sphi_so(1:nsoi, sgfi:sgfi + nsgfi - 1), nsoi, &
+                       ri_dcoeff(ispin)%array(sgfi:sgfi + nsgfi - 1), 1, 0.0_dp, &
+                       ri_dcoeff_so(ispin)%array((iset - 1)*maxso + 1:(iset - 1)*maxso + nsoi), 1)
 
-            DEALLOCATE (work1, work2)
          END DO
-
       END DO
 
       !allocate space to store the spherical orbitals on the grid
@@ -3424,7 +3421,7 @@ SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env)
 !  Build the matrix_soc based on the matrix_s (but anti-symmetric)
       DO i = 1, 3
          CALL dbcsr_create(matrix_soc(i)%matrix, name="SOC MATRIX", template=matrix_s(1)%matrix, &
-                           matrix_type="A")
+                           matrix_type=dbcsr_type_antisymmetric)
       END DO
 
 !  Iterate over its diagonal blocks and fill=it

src/xas_tdp_kernel.F

@@ -20,7 +20,7 @@ MODULE xas_tdp_kernel
         dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
         dbcsr_multiply, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_reserve_block2d, &
-        dbcsr_set, dbcsr_transposed, dbcsr_type
+        dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
    USE dbcsr_tensor_api,                ONLY: dbcsr_t_get_block,&
                                               dbcsr_t_iterator_blocks_left,&
                                               dbcsr_t_iterator_next_block,&
@@ -231,7 +231,7 @@ SUBROUTINE coulomb(coul_ker, contr1_int, dist, blk_size, xas_tdp_env, xas_tdp_co
       PQ => xas_tdp_env%ri_inv_coul
 
       ! Create a normal type work matrix
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       ! Compute the product (aI|P) * (P|Q)^-1 * (Q|Jb) = (aI|Jb)
@@ -252,7 +252,7 @@ SUBROUTINE coulomb(coul_ker, contr1_int, dist, blk_size, xas_tdp_env, xas_tdp_co
       CALL dbcsr_finalize(work_mat)
 
       !Create the symmetric kernel matrix and redistribute work_mat into it
-      CALL dbcsr_create(coul_ker, name="COULOMB KERNEL", matrix_type="S", dist=dist, &
+      CALL dbcsr_create(coul_ker, name="COULOMB KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
       CALL dbcsr_complete_redistribute(work_mat, coul_ker)
 
@@ -300,7 +300,7 @@ SUBROUTINE sc_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_
       ndo_mo = donor_state%ndo_mo
       ri_atom = donor_state%at_index
       !normal type work matrix such that distribution of all spin quadrants match
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       rhs_int => contr1_int_PQ ! contains [ (aI|P)*(P|Q)^-1 ]^T
@@ -377,7 +377,7 @@ SUBROUTINE sc_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_
       CALL dbcsr_finalize(work_mat)
 
       ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
-      CALL dbcsr_create(xc_ker, name="SC OS XC KERNEL", matrix_type="S", dist=dist, &
+      CALL dbcsr_create(xc_ker, name="SC OS XC KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
       CALL dbcsr_complete_redistribute(work_mat, xc_ker)
 
@@ -428,7 +428,7 @@ SUBROUTINE ondiag_sf_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, x
       ndo_mo = donor_state%ndo_mo
       ri_atom = donor_state%at_index
       !normal type work matrix such that distribution of all spin quadrants match
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       !Create a lhs_int, in which the whole (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) will be put
@@ -468,8 +468,8 @@ SUBROUTINE ondiag_sf_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, x
       CALL dbcsr_finalize(work_mat)
 
       ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
-      CALL dbcsr_create(xc_ker, name="ON-DIAG SF OS XC KERNEL", matrix_type="S", dist=dist, &
-                        row_blk_size=blk_size, col_blk_size=blk_size)
+      CALL dbcsr_create(xc_ker, name="ON-DIAG SF OS XC KERNEL", matrix_type=dbcsr_type_symmetric, &
+                        dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
       CALL dbcsr_complete_redistribute(work_mat, xc_ker)
 
       !clean-up
@@ -523,7 +523,7 @@ SUBROUTINE rcs_xc(sg_xc_ker, tp_xc_ker, contr1_int_PQ, dist, blk_size, donor_sta
 
       ! Work structures
       ALLOCATE (fxc(nsgfp, nsgfp))
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       ! Case study: singlet and/or triplet ?
@@ -544,8 +544,8 @@ SUBROUTINE rcs_xc(sg_xc_ker, tp_xc_ker, contr1_int_PQ, dist, blk_size, donor_sta
          CALL dbcsr_finalize(work_mat)
 
          !Create the symmetric kernel matrix and redistribute work_mat into it
-         CALL dbcsr_create(sg_xc_ker, name="XC SINGLET KERNEL", matrix_type="S", dist=dist, &
-                           row_blk_size=blk_size, col_blk_size=blk_size)
+         CALL dbcsr_create(sg_xc_ker, name="XC SINGLET KERNEL", matrix_type=dbcsr_type_symmetric, &
+                           dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
          CALL dbcsr_complete_redistribute(work_mat, sg_xc_ker)
 
       END IF
@@ -567,8 +567,8 @@ SUBROUTINE rcs_xc(sg_xc_ker, tp_xc_ker, contr1_int_PQ, dist, blk_size, donor_sta
          CALL dbcsr_finalize(work_mat)
 
          !Create the symmetric kernel matrix and redistribute work_mat into it
-         CALL dbcsr_create(tp_xc_ker, name="XC TRIPLET KERNEL", matrix_type="S", dist=dist, &
-                           row_blk_size=blk_size, col_blk_size=blk_size)
+         CALL dbcsr_create(tp_xc_ker, name="XC TRIPLET KERNEL", matrix_type=dbcsr_type_symmetric, &
+                           dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
          CALL dbcsr_complete_redistribute(work_mat, tp_xc_ker)
 
       END IF
@@ -724,7 +724,7 @@ SUBROUTINE ondiag_ex(ondiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas
 
       CALL dbcsr_release(mats_desymm)
 
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       ! Loop over donor spin-orbitals. End matrix is symmetric => span only upper half
@@ -789,8 +789,8 @@ SUBROUTINE ondiag_ex(ondiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas
       END DO !iso
 
       CALL dbcsr_finalize(work_mat)
-      CALL dbcsr_create(ondiag_ex_ker, name="ONDIAG EX KERNEL", matrix_type="S", dist=dist, &
-                        row_blk_size=blk_size, col_blk_size=blk_size)
+      CALL dbcsr_create(ondiag_ex_ker, name="ONDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
+                        dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
       CALL dbcsr_complete_redistribute(work_mat, ondiag_ex_ker)
 
       !Clean-up
@@ -852,7 +852,7 @@ SUBROUTINE offdiag_ex_sc(offdiag_ex_ker, contr1_int, dist, blk_size, donor_state
       !the transpose of the later, and put the result in the correct spin quadrants
 
       !Create a normal type work matrix
-      CALL dbcsr_create(work_mat, name="WORK", matrix_type="N", dist=dist, &
+      CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
                         row_blk_size=blk_size, col_blk_size=blk_size)
 
       !Case study on closed-shell, ROKS or UKS
@@ -882,8 +882,8 @@ SUBROUTINE offdiag_ex_sc(offdiag_ex_ker, contr1_int, dist, blk_size, donor_state
       CALL dbcsr_finalize(work_mat)
 
       !Create the symmetric kernel matrix and redistribute work_mat into it
-      CALL dbcsr_create(offdiag_ex_ker, name="OFFDIAG EX KERNEL", matrix_type="S", dist=dist, &
-                        row_blk_size=blk_size, col_blk_size=blk_size)
+      CALL dbcsr_create(offdiag_ex_ker, name="OFFDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
+                        dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
       CALL dbcsr_complete_redistribute(work_mat, offdiag_ex_ker)
 
       !clean-up
@@ -1035,10 +1035,10 @@ SUBROUTINE contract2_AO_to_doMO(contr_int, op_type, donor_state, xas_tdp_env, xa
       CALL cp_dbcsr_dist2d_to_dist(opt_dist2d, opt_dbcsr_dist)
 
       ALLOCATE (aI_P, P_Ib, work, matrices(2))
-      CALL dbcsr_create(aI_P, dist=opt_dbcsr_dist, matrix_type="N", name="(aI|P)", &
+      CALL dbcsr_create(aI_P, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(aI|P)", &
                         row_blk_size=std_blk_size, col_blk_size=ri_blk_size)
 
-      CALL dbcsr_create(P_Ib, dist=opt_dbcsr_dist, matrix_type="N", name="(P|Ib)", &
+      CALL dbcsr_create(P_Ib, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(P|Ib)", &
                         row_blk_size=ri_blk_size, col_blk_size=std_blk_size)
 
       !reserve the blocks (needed for 3c contraction routines)
@@ -1050,8 +1050,9 @@ SUBROUTINE contract2_AO_to_doMO(contr_int, op_type, donor_state, xas_tdp_env, xa
       ALLOCATE (contr_int(ndo_so))
       DO i = 1, ndo_so
          ALLOCATE (contr_int(i)%matrix)
-         CALL dbcsr_create(matrix=contr_int(i)%matrix, template=matrix_s(1)%matrix, matrix_type="N", &
-                           row_blk_size=std_blk_size, col_blk_size=ri_blk_size)
+         CALL dbcsr_create(matrix=contr_int(i)%matrix, template=matrix_s(1)%matrix, &
+                           matrix_type=dbcsr_type_no_symmetry, row_blk_size=std_blk_size, &
+                           col_blk_size=ri_blk_size)
       END DO
 
       ! Only take the coeffs for atom on which MOs I,J are localized
@@ -1183,7 +1184,7 @@ SUBROUTINE contract2_AO_to_doMO_low(ab_Q, vec, mat_aIb, mat_bIa, atom_k)
       CHARACTER(LEN=*), PARAMETER :: routineN = 'contract2_AO_to_doMO_low', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: blk, handle, i, iatom, ind(3), jatom, &
+      INTEGER                                            :: blk, handle, i, iatom, ind(3), j, jatom, &
                                                             katom, n, s1, s2
       INTEGER, DIMENSION(:), POINTER                     :: atom_blk_size
       LOGICAL                                            :: found, t_found
@@ -1219,8 +1220,10 @@ SUBROUTINE contract2_AO_to_doMO_low(ab_Q, vec, mat_aIb, mat_bIa, atom_k)
             CALL dbcsr_get_block_p(matrix=mat_aIb, row=iatom, col=jatom, BLOCK=pblock, found=found)
 
             IF (found) THEN
-               DO i = 1, n
-                  CALL daxpy(s1*s2, vec(i), iabc(:, i, :), 1, pblock(:, :), 1)
+               DO j = 1, s2
+                  DO i = 1, n
+                     pblock(:, j) = pblock(:, j) + vec(i)*iabc(:, i, j)
+                  END DO
                END DO
             END IF
          END IF ! jatom == atom_k
@@ -1363,7 +1366,7 @@ SUBROUTINE ri_int_product(kernel, lhs_int, rhs_int, quadrants, qs_env, eps_filte
       CPASSERT(ANY(quadrants))
       ndo_so = SIZE(lhs_int)
       CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=nblk)
-      CALL dbcsr_create(prod, template=matrix_s(1)%matrix, matrix_type="N")
+      CALL dbcsr_create(prod, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
       my_mt = .FALSE.
       IF (PRESENT(mo_transpose)) my_mt = mo_transpose
 

src/xas_tdp_methods.F

@@ -47,14 +47,9 @@ MODULE xas_tdp_methods
                                               cp_print_key_unit_nr,&
                                               debug_print_level
    USE cp_para_types,                   ONLY: cp_para_env_type
-   USE dbcsr_api,                       ONLY: dbcsr_add_on_diag,&
-                                              dbcsr_copy,&
-                                              dbcsr_create,&
-                                              dbcsr_finalize,&
-                                              dbcsr_get_info,&
-                                              dbcsr_multiply,&
-                                              dbcsr_p_type,&
-                                              dbcsr_set
+   USE dbcsr_api,                       ONLY: &
+        dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_finalize, dbcsr_get_info, &
+        dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_type_no_symmetry
    USE input_constants,                 ONLY: &
         do_loc_none, do_potential_coulomb, do_potential_id, do_potential_short, &
         do_potential_truncated, op_loc_berry, state_loc_list, tddfpt_singlet, tddfpt_spin_cons, &
@@ -759,11 +754,11 @@ SUBROUTINE xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
       ALLOCATE (xas_tdp_env%q_projector(nspins))
       ALLOCATE (xas_tdp_env%q_projector(1)%matrix)
       CALL dbcsr_create(xas_tdp_env%q_projector(1)%matrix, name="Q PROJECTOR ALPHA", &
-                        template=matrix_s(1)%matrix, matrix_type="N")
+                        template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
       IF (do_os) THEN
          ALLOCATE (xas_tdp_env%q_projector(2)%matrix)
          CALL dbcsr_create(xas_tdp_env%q_projector(2)%matrix, name="Q PROJECTOR BETA", &
-                           template=matrix_s(1)%matrix, matrix_type="N")
+                           template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
       END IF
 
 !     In the case of spin-restricted calculations, rho_ao includes double occupency => 0.5 prefactor