From 0bf0ab03676e3f86b3fa8e5dfea73fcc5cc632b9 Mon Sep 17 00:00:00 2001 From: annahehn <36073704+annahehn@users.noreply.github.com> Date: Sun, 24 May 2026 15:27:23 +0200 Subject: [PATCH] =?UTF-8?q?Printing=20Cartesian=20overlap=20matrix=20in=20?= =?UTF-8?q?line=20with=20Cartesian=20MOs=20such=20tha=E2=80=A6=20(#5271)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Thomas D. Kuehne --- src/aobasis/basis_set_types.F | 84 +++++++++++++++++++++-- src/cp_dbcsr_output.F | 42 +++++++++--- src/input_cp2k_print_dft.F | 5 ++ src/particle_methods.F | 23 ++++++- src/qs_mo_io.F | 81 ++++++++++++++++++++--- src/qs_overlap.F | 92 ++++++++++++++++++++------ src/qs_scf_output.F | 13 ++-- tests/QS/regtest-gpw-1/TEST_FILES.toml | 2 + tests/QS/regtest-gpw-1/h2test.inp | 73 ++++++++++++++++++++ 9 files changed, 361 insertions(+), 54 deletions(-) create mode 100644 tests/QS/regtest-gpw-1/h2test.inp diff --git a/src/aobasis/basis_set_types.F b/src/aobasis/basis_set_types.F index 3820ca5c71..3b2b3555fa 100644 --- a/src/aobasis/basis_set_types.F +++ b/src/aobasis/basis_set_types.F @@ -85,7 +85,7 @@ MODULE basis_set_types m => NULL(), ncgf_set => NULL(), & npgf => NULL(), nsgf_set => NULL(), nshell => NULL() REAL(KIND=dp), DIMENSION(:, :), POINTER :: cphi => NULL(), pgf_radius => NULL(), sphi => NULL(), & - scon => NULL(), zet => NULL() + scon => NULL(), zet => NULL(), ccon => NULL() INTEGER, DIMENSION(:, :), POINTER :: first_cgf => NULL(), first_sgf => NULL(), l => NULL(), & last_cgf => NULL(), last_sgf => NULL(), n => NULL() REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc => NULL() @@ -193,6 +193,7 @@ CONTAINS IF (ASSOCIATED(gto_basis_set%pgf_radius)) DEALLOCATE (gto_basis_set%pgf_radius) IF (ASSOCIATED(gto_basis_set%sphi)) DEALLOCATE (gto_basis_set%sphi) IF (ASSOCIATED(gto_basis_set%scon)) DEALLOCATE (gto_basis_set%scon) + IF (ASSOCIATED(gto_basis_set%ccon)) DEALLOCATE (gto_basis_set%ccon) IF (ASSOCIATED(gto_basis_set%zet)) DEALLOCATE (gto_basis_set%zet) IF (ASSOCIATED(gto_basis_set%first_cgf)) DEALLOCATE (gto_basis_set%first_cgf) IF (ASSOCIATED(gto_basis_set%first_sgf)) DEALLOCATE (gto_basis_set%first_sgf) @@ -254,9 +255,16 @@ CONTAINS basis_set_out%nsgf_set = basis_set_in%nsgf_set maxco = SIZE(basis_set_in%cphi, 1) ALLOCATE (basis_set_out%cphi(maxco, ncgf), basis_set_out%sphi(maxco, nsgf), basis_set_out%scon(maxco, nsgf)) + ALLOCATE (basis_set_out%ccon(maxco, ncgf)) basis_set_out%cphi = basis_set_in%cphi basis_set_out%sphi = basis_set_in%sphi basis_set_out%scon = basis_set_in%scon + basis_set_out%ccon = 0.0_dp + IF (ASSOCIATED(basis_set_in%ccon)) THEN + IF ((SIZE(basis_set_in%ccon, 1) == maxco) .AND. (SIZE(basis_set_in%ccon, 2) == ncgf)) THEN + basis_set_out%ccon = basis_set_in%ccon + END IF + END IF maxpgf = MAXVAL(basis_set_in%npgf) ALLOCATE (basis_set_out%pgf_radius(maxpgf, nset), basis_set_out%zet(maxpgf, nset)) basis_set_out%pgf_radius = basis_set_in%pgf_radius @@ -443,6 +451,8 @@ CONTAINS pbasis%sphi = 0.0_dp ALLOCATE (pbasis%scon(maxco, ncgf)) pbasis%scon = 0.0_dp + ALLOCATE (pbasis%ccon(maxco, ncgf)) + pbasis%ccon = 0.0_dp ALLOCATE (pbasis%set_radius(nset)) ALLOCATE (pbasis%pgf_radius(mpgf, nset)) pbasis%pgf_radius = 0.0_dp @@ -570,6 +580,7 @@ CONTAINS CALL reallocate(basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(basis_set%pgf_radius, 1, maxpgf, 1, nset) END SUBROUTINE combine_basis_sets @@ -622,12 +633,13 @@ CONTAINS !> \param maxder ... !> \param short_kind_radius ... !> \param npgf_seg_sum number of primitives in "segmented contraction format" +!> \param ccon ... ! ************************************************************************************************** SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, & nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, & m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, & last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, & - npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum) + npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon) ! Get informations about a Gaussian-type orbital (GTO) basis set. @@ -655,6 +667,7 @@ CONTAINS INTEGER, INTENT(IN), OPTIONAL :: maxder REAL(KIND=dp), INTENT(OUT), OPTIONAL :: short_kind_radius INTEGER, INTENT(OUT), OPTIONAL :: npgf_seg_sum + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: ccon INTEGER :: iset, nder @@ -684,6 +697,7 @@ CONTAINS IF (PRESENT(pgf_radius)) pgf_radius => gto_basis_set%pgf_radius IF (PRESENT(sphi)) sphi => gto_basis_set%sphi IF (PRESENT(scon)) scon => gto_basis_set%scon + IF (PRESENT(ccon)) ccon => gto_basis_set%ccon IF (PRESENT(zet)) zet => gto_basis_set%zet IF (PRESENT(first_cgf)) first_cgf => gto_basis_set%first_cgf IF (PRESENT(first_sgf)) first_sgf => gto_basis_set%first_sgf @@ -789,7 +803,7 @@ CONTAINS END SELECT ! Initialise the transformation matrices "pgf" -> "cgf" - CALL init_cphi_and_sphi(gto_basis_set) + CALL init_cphi_and_sphi(gto_basis_set, .FALSE.) CALL timestop(handle) @@ -798,8 +812,9 @@ CONTAINS ! ************************************************************************************************** !> \brief ... !> \param gto_basis_set ... +!> \param lccon ... ! ************************************************************************************************** - SUBROUTINE init_cphi_and_sphi(gto_basis_set) + SUBROUTINE init_cphi_and_sphi(gto_basis_set, lccon) ! Initialise the matrices for the transformation of primitive Cartesian ! Gaussian-type functions to contracted Cartesian (cphi) and spherical @@ -808,6 +823,7 @@ CONTAINS ! - Creation (20.09.2000,MK) TYPE(gto_basis_set_type), INTENT(INOUT) :: gto_basis_set + LOGICAL, INTENT(IN), OPTIONAL :: lccon CHARACTER(len=*), PARAMETER :: routineN = 'init_cphi_and_sphi' @@ -815,7 +831,10 @@ CONTAINS ipgf, iset, ishell, l, last_sgf, lmax, & lmin, n, n1, n2, ncgf, nn, nn1, nn2, & npgf, nsgf + LOGICAL :: my_lccon + my_lccon = .FALSE. + IF (PRESENT(lccon)) my_lccon = lccon ! ------------------------------------------------------------------------- ! Build the Cartesian transformation matrix "cphi" @@ -896,6 +915,36 @@ CONTAINS END DO END IF + IF (my_lccon) THEN + IF (.NOT. ASSOCIATED(gto_basis_set%ccon)) THEN + CALL reallocate(gto_basis_set%ccon, 1, SIZE(gto_basis_set%cphi, 1), 1, gto_basis_set%ncgf) + ELSE IF ((SIZE(gto_basis_set%ccon, 1) /= SIZE(gto_basis_set%cphi, 1)) .OR. & + (SIZE(gto_basis_set%ccon, 2) /= gto_basis_set%ncgf)) THEN + CALL reallocate(gto_basis_set%ccon, 1, SIZE(gto_basis_set%cphi, 1), 1, gto_basis_set%ncgf) + END IF + n = SIZE(gto_basis_set%ccon, 1) + gto_basis_set%ccon = 0.0_dp + IF (n > 0) THEN + DO iset = 1, gto_basis_set%nset + lmin = gto_basis_set%lmin(iset) + lmax = gto_basis_set%lmax(iset) + npgf = gto_basis_set%npgf(iset) + nn = ncoset(lmax) - ncoset(lmin - 1) + DO ishell = 1, gto_basis_set%nshell(iset) + first_sgf = gto_basis_set%first_cgf(ishell, iset) + last_sgf = gto_basis_set%last_cgf(ishell, iset) + DO ipgf = 1, npgf + nn1 = (ipgf - 1)*ncoset(lmax) + ncoset(lmin - 1) + 1 + nn2 = ipgf*ncoset(lmax) + n1 = (ipgf - 1)*nn + 1 + n2 = ipgf*nn + gto_basis_set%ccon(n1:n2, first_sgf:last_sgf) = gto_basis_set%cphi(nn1:nn2, first_sgf:last_sgf) + END DO + END DO + END DO + END IF + END IF + CALL timestop(handle) END SUBROUTINE init_cphi_and_sphi @@ -1118,7 +1167,7 @@ CONTAINS ! Initialise the transformation matrices "pgf" -> "cgf" - CALL init_cphi_and_sphi(gto_basis_set) + CALL init_cphi_and_sphi(gto_basis_set, .TRUE.) CALL timestop(handle) @@ -1255,6 +1304,7 @@ CONTAINS CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(gto_basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%lx, 1, ncgf) CALL reallocate(gto_basis_set%ly, 1, ncgf) CALL reallocate(gto_basis_set%lz, 1, ncgf) @@ -1410,6 +1460,7 @@ CONTAINS CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(gto_basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%lx, 1, ncgf) CALL reallocate(gto_basis_set%ly, 1, ncgf) CALL reallocate(gto_basis_set%lz, 1, ncgf) @@ -1561,6 +1612,7 @@ CONTAINS CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(gto_basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%lx, 1, ncgf) CALL reallocate(gto_basis_set%ly, 1, ncgf) CALL reallocate(gto_basis_set%lz, 1, ncgf) @@ -1713,6 +1765,7 @@ CONTAINS CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(gto_basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%lx, 1, ncgf) CALL reallocate(gto_basis_set%ly, 1, ncgf) CALL reallocate(gto_basis_set%lz, 1, ncgf) @@ -1793,13 +1846,14 @@ CONTAINS !> \param n ... !> \param gcc ... !> \param short_kind_radius ... +!> \param ccon ... !> \author MK ! ************************************************************************************************** SUBROUTINE set_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, & nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, & lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, & cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, & - last_cgf, last_sgf, n, gcc, short_kind_radius) + last_cgf, last_sgf, n, gcc, short_kind_radius, ccon) TYPE(gto_basis_set_type), INTENT(INOUT) :: gto_basis_set CHARACTER(LEN=default_string_length), INTENT(IN), & @@ -1818,6 +1872,7 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & POINTER :: gcc REAL(KIND=dp), INTENT(IN), OPTIONAL :: short_kind_radius + REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: ccon IF (PRESENT(name)) gto_basis_set%name = name IF (PRESENT(aliases)) gto_basis_set%aliases = aliases @@ -1845,6 +1900,7 @@ CONTAINS IF (PRESENT(pgf_radius)) gto_basis_set%pgf_radius(:, :) = pgf_radius(:, :) IF (PRESENT(sphi)) gto_basis_set%sphi(:, :) = sphi(:, :) IF (PRESENT(scon)) gto_basis_set%scon(:, :) = scon(:, :) + IF (PRESENT(ccon)) gto_basis_set%ccon(:, :) = ccon(:, :) IF (PRESENT(zet)) gto_basis_set%zet(:, :) = zet(:, :) IF (PRESENT(first_cgf)) gto_basis_set%first_cgf(:, :) = first_cgf(:, :) IF (PRESENT(first_sgf)) gto_basis_set%first_sgf(:, :) = first_sgf(:, :) @@ -2050,6 +2106,8 @@ CONTAINS WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%sphi WRITE (UNIT=output_unit, FMT="(A1)") "SCON" WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%scon + WRITE (UNIT=output_unit, FMT="(A1)") "CCON" + WRITE (UNIT=output_unit, FMT="(12F10.5)") gto_basis_set%ccon END IF @@ -2536,6 +2594,7 @@ CONTAINS CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf) CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf) + CALL reallocate(gto_basis_set%ccon, 1, maxco, 1, ncgf) CALL reallocate(gto_basis_set%lx, 1, ncgf) CALL reallocate(gto_basis_set%ly, 1, ncgf) CALL reallocate(gto_basis_set%lz, 1, ncgf) @@ -2711,11 +2770,12 @@ CONTAINS INTEGER, ALLOCATABLE, DIMENSION(:) :: sort_index INTEGER, ALLOCATABLE, DIMENSION(:, :) :: icgf_set, isgf_set INTEGER, DIMENSION(:), POINTER :: lx, ly, lz, m, npgf + LOGICAL :: ccon_available REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp REAL(dp), DIMENSION(:), POINTER :: set_radius REAL(dp), DIMENSION(:, :), POINTER :: zet REAL(KIND=dp), DIMENSION(:), POINTER :: norm_cgf - REAL(KIND=dp), DIMENSION(:, :), POINTER :: cphi, scon, sphi + REAL(KIND=dp), DIMENSION(:, :), POINTER :: ccon, cphi, scon, sphi NULLIFY (set_radius, zet) @@ -2798,6 +2858,13 @@ CONTAINS sphi = 0.0_dp ALLOCATE (scon(SIZE(basis_set%scon, 1), SIZE(basis_set%scon, 2))) scon = 0.0_dp + ALLOCATE (ccon(SIZE(basis_set%cphi, 1), SIZE(basis_set%cphi, 2))) + ccon = 0.0_dp + ccon_available = ASSOCIATED(basis_set%ccon) + IF (ccon_available) THEN + ccon_available = (SIZE(basis_set%ccon, 1) == SIZE(ccon, 1)) .AND. & + (SIZE(basis_set%ccon, 2) == SIZE(ccon, 2)) + END IF ALLOCATE (sgf_symbol(SIZE(basis_set%sgf_symbol))) ALLOCATE (m(SIZE(basis_set%m))) @@ -2814,6 +2881,7 @@ CONTAINS ly(icgf_new) = basis_set%ly(icgf_old) lz(icgf_new) = basis_set%lz(icgf_old) cphi(:, icgf_new) = basis_set%cphi(:, icgf_old) + IF (ccon_available) ccon(:, icgf_new) = basis_set%ccon(:, icgf_old) cgf_symbol(icgf_new) = basis_set%cgf_symbol(icgf_old) END DO DO is = 1, is_max @@ -2843,6 +2911,8 @@ CONTAINS basis_set%sphi => sphi DEALLOCATE (basis_set%scon) basis_set%scon => scon + IF (ASSOCIATED(basis_set%ccon)) DEALLOCATE (basis_set%ccon) + basis_set%ccon => ccon DEALLOCATE (basis_set%m) basis_set%m => m diff --git a/src/cp_dbcsr_output.F b/src/cp_dbcsr_output.F index 45f6b5c4bd..2bb78a1201 100644 --- a/src/cp_dbcsr_output.F +++ b/src/cp_dbcsr_output.F @@ -35,7 +35,8 @@ MODULE cp_dbcsr_output USE machine, ONLY: m_flush USE mathlib, ONLY: symmetrize_matrix USE message_passing, ONLY: mp_para_env_type - USE orbital_pointers, ONLY: nso + USE orbital_pointers, ONLY: nco,& + nso USE particle_methods, ONLY: get_particle_set USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& @@ -156,10 +157,11 @@ CONTAINS !> \param scale ... !> \param output_unit ... !> \param omit_headers Write only the matrix data, not the row/column headers +!> \param cartesian_basis Use Cartesian instead of spherical basis labels ! ************************************************************************************************** SUBROUTINE cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, & first_row, last_row, first_col, last_col, scale, & - output_unit, omit_headers) + output_unit, omit_headers, cartesian_basis) TYPE(dbcsr_type) :: sparse_matrix INTEGER, INTENT(IN) :: before, after @@ -168,11 +170,12 @@ CONTAINS INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, first_col, last_col REAL(dp), INTENT(IN), OPTIONAL :: scale INTEGER, INTENT(IN) :: output_unit - LOGICAL, INTENT(IN), OPTIONAL :: omit_headers + LOGICAL, INTENT(IN), OPTIONAL :: omit_headers, cartesian_basis CHARACTER(LEN=default_string_length) :: matrix_name INTEGER :: col1, col2, dim_col, dim_row, row1, row2 - LOGICAL :: my_omit_headers, print_sym + LOGICAL :: my_cartesian_basis, my_omit_headers, & + print_sym REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set @@ -233,11 +236,14 @@ CONTAINS ELSE my_omit_headers = .FALSE. END IF + my_cartesian_basis = .FALSE. + IF (PRESENT(cartesian_basis)) my_cartesian_basis = cartesian_basis CALL dbcsr_get_info(sparse_matrix, name=matrix_name) IF (print_sym) THEN CALL write_matrix_sym(matrix, matrix_name, before, after, qs_env, para_env, & - row1, row2, col1, col2, output_unit, my_omit_headers) + row1, row2, col1, col2, output_unit, my_omit_headers, & + cartesian_basis=my_cartesian_basis) ELSE CALL write_matrix_gen(matrix, matrix_name, before, after, para_env, & row1, row2, col1, col2, output_unit, my_omit_headers) @@ -328,10 +334,11 @@ CONTAINS !> \param last_col ... !> \param output_unit ... !> \param omit_headers Write only the matrix data, not the row/column headers +!> \param cartesian_basis Use Cartesian instead of spherical basis labels !> \author Creation (01.07.2003,MK) ! ************************************************************************************************** SUBROUTINE write_matrix_sym(matrix, matrix_name, before, after, qs_env, para_env, & - first_row, last_row, first_col, last_col, output_unit, omit_headers) + first_row, last_row, first_col, last_col, output_unit, omit_headers, cartesian_basis) REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix CHARACTER(LEN=*), INTENT(IN) :: matrix_name @@ -341,7 +348,9 @@ CONTAINS INTEGER, INTENT(IN) :: first_row, last_row, first_col, & last_col, output_unit LOGICAL, INTENT(IN) :: omit_headers + LOGICAL, INTENT(IN), OPTIONAL :: cartesian_basis + CHARACTER(LEN=12), DIMENSION(:), POINTER :: cgf_symbol CHARACTER(LEN=2) :: element_symbol CHARACTER(LEN=25) :: fmtstr1 CHARACTER(LEN=35) :: fmtstr2 @@ -353,6 +362,7 @@ CONTAINS INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf INTEGER, DIMENSION(:), POINTER :: nshell INTEGER, DIMENSION(:, :), POINTER :: lshell + LOGICAL :: my_cartesian_basis TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(gto_basis_set_type), POINTER :: orb_basis_set TYPE(particle_type), DIMENSION(:), POINTER :: particle_set @@ -368,6 +378,9 @@ CONTAINS natom = SIZE(particle_set) + my_cartesian_basis = .FALSE. + IF (PRESENT(cartesian_basis)) my_cartesian_basis = cartesian_basis + CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf) ALLOCATE (first_sgf(natom)) @@ -425,20 +438,27 @@ CONTAINS CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & - nset=nset, nshell=nshell, l=lshell, sgf_symbol=sgf_symbol) + nset=nset, nshell=nshell, l=lshell, & + cgf_symbol=cgf_symbol, sgf_symbol=sgf_symbol) isgf = 1 DO iset = 1, nset DO ishell = 1, nshell(iset) l = lshell(ishell, iset) - DO iso = 1, nso(l) + DO iso = 1, MERGE(nco(l), nso(l), my_cartesian_basis) IF ((irow >= first_row) .AND. (irow <= last_row)) THEN IF (omit_headers) THEN WRITE (UNIT=output_unit, FMT=fmtstr2) & (matrix(irow, jcol), jcol=from, to) ELSE - WRITE (UNIT=output_unit, FMT=fmtstr2) & - irow, iatom, element_symbol, sgf_symbol(isgf), & - (matrix(irow, jcol), jcol=from, to) + IF (my_cartesian_basis) THEN + WRITE (UNIT=output_unit, FMT=fmtstr2) & + irow, iatom, element_symbol, cgf_symbol(isgf), & + (matrix(irow, jcol), jcol=from, to) + ELSE + WRITE (UNIT=output_unit, FMT=fmtstr2) & + irow, iatom, element_symbol, sgf_symbol(isgf), & + (matrix(irow, jcol), jcol=from, to) + END IF END IF END IF isgf = isgf + 1 diff --git a/src/input_cp2k_print_dft.F b/src/input_cp2k_print_dft.F index 53037bbfa2..c6292557b5 100644 --- a/src/input_cp2k_print_dft.F +++ b/src/input_cp2k_print_dft.F @@ -549,6 +549,11 @@ CONTAINS default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="CARTESIAN_OVERLAP", & + description="Print the Cartesian overlap matrix in line with Cartesian MO coefficients.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="ENERGIES", & variants=s2a("EIGENVALUES", "EIGVALS"), & description="Print the MO energies (eigenvalues).", & diff --git a/src/particle_methods.F b/src/particle_methods.F index 8cdaa9dbe9..60858eec45 100644 --- a/src/particle_methods.F +++ b/src/particle_methods.F @@ -108,6 +108,7 @@ CONTAINS !> \param nsgf ... !> \param nmao ... !> \param basis ... +!> \param ncgf ... !> \date 14.01.2002 !> \par History !> - particle type cleaned (13.10.2003,MK) @@ -116,12 +117,13 @@ CONTAINS !> \version 1.0 ! ************************************************************************************************** SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, & - nmao, basis) + nmao, basis, ncgf) TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: first_sgf, last_sgf, nsgf, nmao TYPE(gto_basis_set_p_type), DIMENSION(:), OPTIONAL :: basis + INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: ncgf INTEGER :: ikind, iparticle, isgf, nparticle, ns @@ -140,6 +142,9 @@ CONTAINS IF (PRESENT(nmao)) THEN CPASSERT(SIZE(nmao) >= nparticle) END IF + IF (PRESENT(ncgf)) THEN + CPASSERT(SIZE(ncgf) >= nparticle) + END IF IF (PRESENT(first_sgf) .OR. PRESENT(last_sgf) .OR. PRESENT(nsgf)) THEN isgf = 0 @@ -161,6 +166,22 @@ CONTAINS END DO END IF + IF (PRESENT(ncgf)) THEN + DO iparticle = 1, nparticle + CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind) + IF (PRESENT(basis)) THEN + IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, ncgf=ns) + ELSE + ns = 0 + END IF + ELSE + CALL get_qs_kind(qs_kind_set(ikind), ncgf=ns) + END IF + ncgf(iparticle) = ns + END DO + END IF + IF (PRESENT(first_sgf)) THEN IF (SIZE(first_sgf) > nparticle) first_sgf(nparticle + 1) = isgf + 1 END IF diff --git a/src/qs_mo_io.F b/src/qs_mo_io.F index e95b5dd2f2..31ab05dcea 100644 --- a/src/qs_mo_io.F +++ b/src/qs_mo_io.F @@ -22,6 +22,7 @@ MODULE qs_mo_io USE atomic_kind_types, ONLY: get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& + gto_basis_set_p_type,& gto_basis_set_type USE cp_dbcsr_api, ONLY: dbcsr_binary_write,& dbcsr_create,& @@ -30,7 +31,9 @@ MODULE qs_mo_io dbcsr_type USE cp_dbcsr_contrib, ONLY: dbcsr_checksum USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& - copy_fm_to_dbcsr + copy_fm_to_dbcsr,& + dbcsr_deallocate_matrix_set + USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix USE cp_files, ONLY: close_file,& open_file USE cp_fm_types, ONLY: cp_fm_get_info,& @@ -68,13 +71,20 @@ MODULE qs_mo_io USE qs_density_matrices, ONLY: calculate_density_matrix USE qs_dftb_types, ONLY: qs_dftb_atom_type USE qs_dftb_utils, ONLY: get_dftb_atom_param + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type USE qs_kind_types, ONLY: get_qs_kind,& get_qs_kind_set,& qs_kind_type + USE qs_ks_types, ONLY: qs_ks_env_type USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues USE qs_mo_occupation, ONLY: set_mo_occupation USE qs_mo_types, ONLY: get_mo_set,& mo_set_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,& + release_neighbor_list_sets + USE qs_neighbor_lists, ONLY: setup_neighbor_list + USE qs_overlap, ONLY: build_overlap_matrix_simple #include "./base/base_uses.f90" IMPLICIT NONE @@ -995,6 +1005,7 @@ CONTAINS !> \param cpart ... !> \param sim_step ... !> \param umo_set ... +!> \param qs_env ... !> \date 15.05.2001 !> \par History: !> - Optionally print Cartesian MOs (20.04.2005, MK) @@ -1007,7 +1018,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, & dft_section, before, kpoint, final_mos, spin, & - solver_method, rtp, cpart, sim_step, umo_set) + solver_method, rtp, cpart, sim_step, umo_set, qs_env) TYPE(mo_set_type), INTENT(IN) :: mo_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set @@ -1020,6 +1031,7 @@ CONTAINS LOGICAL, INTENT(IN), OPTIONAL :: rtp INTEGER, INTENT(IN), OPTIONAL :: cpart, sim_step TYPE(mo_set_type), INTENT(IN), OPTIONAL :: umo_set + TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env CHARACTER(LEN=12) :: symbol CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol @@ -1036,22 +1048,27 @@ CONTAINS CHARACTER(LEN=40) :: fmtstr3 CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol INTEGER :: after, first_mo, from, homo, iatom, icgf, ico, icol, ikind, imo, irow, iset, & - isgf, ishell, iso, iw, jcol, last_mo, left, lmax, lshell, nao, natom, ncgf, ncol, nmo, & - nset, nsgf, numo, right, scf_step, to, width + isgf, ishell, iso, iw, jcol, last_mo, left, lmax, lshell, nao, natom, ncgf, ncol, nkind, & + nmo, nset, nsgf, numo, right, scf_step, to, width INTEGER, DIMENSION(:), POINTER :: mo_index_range, nshell INTEGER, DIMENSION(:, :), POINTER :: l - LOGICAL :: ionode, my_final, my_rtp, & - print_cartesian, print_eigvals, & - print_eigvecs, print_occup, & - should_output + LOGICAL :: ionode, my_final, my_rtp, omit_headers, print_cartesian, print_cartesian_overlap, & + print_eigvals, print_eigvecs, print_occup, should_output REAL(KIND=dp) :: gap, maxocc REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mo_eigenvalues, mo_occupation_numbers REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers TYPE(cp_fm_type), POINTER :: mo_coeff, umo_coeff TYPE(cp_logger_type), POINTER :: logger - TYPE(gto_basis_set_type), POINTER :: orb_basis_set + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sro + TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: orb_basis_set_list + TYPE(gto_basis_set_type), POINTER :: orb_basis_set, orbbasis + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sro_list TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter + TYPE(qs_kind_type), POINTER :: qs_kind + TYPE(qs_ks_env_type), POINTER :: ks_env NULLIFY (bcgf_symbol) NULLIFY (bsgf_symbol) @@ -1068,6 +1085,7 @@ CONTAINS CALL section_vals_val_get(dft_section, "PRINT%MO%CARTESIAN", l_val=print_cartesian) CALL section_vals_val_get(dft_section, "PRINT%MO%MO_INDEX_RANGE", i_vals=mo_index_range) CALL section_vals_val_get(dft_section, "PRINT%MO%NDIGITS", i_val=after) + CALL section_vals_val_get(dft_section, "PRINT%MO%CARTESIAN_OVERLAP", l_val=print_cartesian_overlap) after = MIN(MAX(after, 1), 16) ! Do we print the final MO information after SCF convergence is reached (default: no) @@ -1161,6 +1179,49 @@ CONTAINS END IF END IF + IF (PRESENT(qs_env)) THEN + IF (ASSOCIATED(qs_env) .AND. my_final .AND. print_cartesian_overlap) THEN + NULLIFY (qs_kind_set) + CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set) + nkind = SIZE(qs_kind_set) + + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN + ALLOCATE (orb_basis_set_list(nkind)) + DO ikind = 1, nkind + qs_kind => qs_kind_set(ikind) + NULLIFY (orb_basis_set_list(ikind)%gto_basis_set) + NULLIFY (orbbasis) + CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB") + IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis + END DO + NULLIFY (sro_list) + CALL setup_neighbor_list(sro_list, orb_basis_set_list, qs_env=qs_env) + NULLIFY (sro) + NULLIFY (para_env) + CALL get_qs_env(qs_env, ks_env=ks_env, para_env=para_env) + CALL build_overlap_matrix_simple(ks_env, sro, & + orb_basis_set_list, orb_basis_set_list, sro_list, .TRUE.) + CALL release_neighbor_list_sets(sro_list) + + iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", & + extension=".Log") + CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) + CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) + after = MIN(MAX(after, 1), 16) + IF (ASSOCIATED(sro)) THEN + CALL cp_dbcsr_write_sparse_matrix(sro(1)%matrix, 4, after, qs_env, para_env, & + output_unit=iw, omit_headers=omit_headers, & + cartesian_basis=.TRUE.) + END IF + CALL cp_print_key_finished_output(iw, logger, qs_env%input, & + "DFT%PRINT%AO_MATRICES/OVERLAP") + IF (ASSOCIATED(sro)) CALL dbcsr_deallocate_matrix_set(sro) + DEALLOCATE (orb_basis_set_list) + END IF + END IF + END IF + iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO", & ignore_should_output=should_output, & extension=".MOLog") @@ -1441,10 +1502,10 @@ CONTAINS ! Release work storage - DEALLOCATE (smatrix) IF (print_cartesian) THEN DEALLOCATE (cmatrix) END IF + DEALLOCATE (smatrix) ELSE IF (print_occup .OR. print_eigvals) THEN diff --git a/src/qs_overlap.F b/src/qs_overlap.F index ba80693bf5..3bc52c5249 100644 --- a/src/qs_overlap.F +++ b/src/qs_overlap.F @@ -286,11 +286,11 @@ CONTAINS IF (dokp) THEN CALL dbcsr_allocate_matrix_set(matrixkp_s, maxder, nimg) CALL create_sab_matrix(ks_env, matrixkp_s, matrix_name, basis_set_list_a, basis_set_list_b, & - sab_nl, do_symmetric) + sab_nl, do_symmetric, lcart=.FALSE.) ELSE CALL dbcsr_allocate_matrix_set(matrix_s, maxder) CALL create_sab_matrix(ks_env, matrix_s, matrix_name, basis_set_list_a, basis_set_list_b, & - sab_nl, do_symmetric) + sab_nl, do_symmetric, lcart=.FALSE.) END IF maxs = maxder @@ -563,6 +563,7 @@ CONTAINS !> \param basis_set_list_a basis set list to be used for bra in !> \param basis_set_list_b basis set list to be used for ket in !> \param sab_nl pair list (must be consistent with basis sets!) +!> \param lcart ... !> \date 11.03.2016 !> \par History !> Simplified version of build_overlap_matrix @@ -570,13 +571,14 @@ CONTAINS !> \version 1.0 ! ************************************************************************************************** SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, & - basis_set_list_a, basis_set_list_b, sab_nl) + basis_set_list_a, basis_set_list_b, sab_nl, lcart) TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_nl + LOGICAL, OPTIONAL :: lcart CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_matrix_simple' @@ -587,13 +589,14 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & npgfb, nsgfa, nsgfb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb - LOGICAL :: do_symmetric, found, trans + LOGICAL :: do_symmetric, found, ldocart, trans REAL(KIND=dp) :: dab, rab2 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint REAL(KIND=dp), DIMENSION(3) :: rab REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b - REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb + REAL(KIND=dp), DIMENSION(:, :), POINTER :: ccona, cconb, rpgfa, rpgfb, scon_a, & + scon_b, zeta, zetb TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint TYPE(dft_control_type), POINTER :: dft_control @@ -606,6 +609,8 @@ CONTAINS !$ INTEGER(KIND=int_8) :: iatom8 !$ INTEGER, PARAMETER :: nlock = 501 + ldocart = .FALSE. + IF (PRESENT(lcart)) ldocart = lcart NULLIFY (dft_control) CALL timeset(routineN, handle) @@ -624,8 +629,13 @@ CONTAINS nkind = SIZE(qs_kind_set) CALL dbcsr_allocate_matrix_set(matrix_s, 1) - CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, & - sab_nl, do_symmetric) + IF (.NOT. ldocart) THEN + CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, & + sab_nl, do_symmetric, lcart=ldocart) + ELSE + CALL create_sab_matrix(ks_env, matrix_s, "Cartesian Overlap Matrix", basis_set_list_a, basis_set_list_b, & + sab_nl, do_symmetric, lcart=ldocart) + END IF ldsab = 0 DO ikind = 1, nkind @@ -639,12 +649,12 @@ CONTAINS !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED (ldsab,sab_nl,do_symmetric,ncoset,natom,& -!$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks) & +!$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks,ldocart) & !$OMP PRIVATE (oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,& !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, & !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, & !$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, & -!$OMP slot, lock_num, hash, hash1, hash2, iatom8 ) +!$OMP slot, lock_num, hash, hash1, hash2, iatom8, ccona, cconb ) !$OMP SINGLE !$ ALLOCATE (locks(nlock)) @@ -684,6 +694,11 @@ CONTAINS rpgfa => basis_set_a%pgf_radius set_radius_a => basis_set_a%set_radius scon_a => basis_set_a%scon + IF (ldocart) THEN + first_sgfa => basis_set_a%first_cgf + nsgfa => basis_set_a%ncgf_set + ccona => basis_set_a%ccon + END IF zeta => basis_set_a%zet ! basis jkind first_sgfb => basis_set_b%first_sgf @@ -695,6 +710,11 @@ CONTAINS rpgfb => basis_set_b%pgf_radius set_radius_b => basis_set_b%set_radius scon_b => basis_set_b%scon + IF (ldocart) THEN + first_sgfb => basis_set_b%first_cgf + nsgfb => basis_set_b%ncgf_set + cconb => basis_set_b%ccon + END IF zetb => basis_set_b%zet IF (do_symmetric) THEN @@ -741,8 +761,13 @@ CONTAINS lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & rab, sab=oint(:, :, 1)) ! Contraction - CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), & - cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans) + IF (.NOT. ldocart) THEN + CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), & + cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans) + ELSE + CALL contraction(oint(:, :, 1), owork, ca=ccona(:, sgfa:), na=n1, ma=nsgfa(iset), & + cb=cconb(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans) + END IF !$OMP CRITICAL(blockadd) CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, & sgfa, sgfb, trans=trans) @@ -1028,9 +1053,10 @@ CONTAINS !> \param basis_set_list_b Basis set used for |b> !> \param sab_nl Overlap neighbor list !> \param symmetric Is symmetry used in the neighbor list? +!> \param lcart ... ! ************************************************************************************************** SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, & - basis_set_list_a, basis_set_list_b, sab_nl, symmetric) + basis_set_list_a, basis_set_list_b, sab_nl, symmetric, lcart) TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s @@ -1039,6 +1065,7 @@ CONTAINS TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_nl LOGICAL, INTENT(IN) :: symmetric + LOGICAL, INTENT(IN), OPTIONAL :: lcart CHARACTER(LEN=12) :: cgfsym CHARACTER(LEN=32) :: symmetry_string @@ -1049,10 +1076,14 @@ CONTAINS TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + LOGICAL:: my_lcart + CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, & qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist) natom = SIZE(particle_set) + my_lcart = .FALSE. + IF (PRESENT(lcart)) my_lcart = lcart IF (PRESENT(matrix_name)) THEN mname = matrix_name @@ -1064,10 +1095,17 @@ CONTAINS ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) - CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & - basis=basis_set_list_a) - CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & - basis=basis_set_list_b) + IF (.NOT. my_lcart) THEN + CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & + basis=basis_set_list_a) + CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & + basis=basis_set_list_b) + ELSE + CALL get_particle_set(particle_set, qs_kind_set, ncgf=row_blk_sizes, & + basis=basis_set_list_a) + CALL get_particle_set(particle_set, qs_kind_set, ncgf=col_blk_sizes, & + basis=basis_set_list_b) + END IF ! prepare for allocation IF (symmetric) THEN @@ -1118,9 +1156,10 @@ CONTAINS !> \param basis_set_list_b Basis set used for |b> !> \param sab_nl Overlap neighbor list !> \param symmetric Is symmetry used in the neighbor list? +!> \param lcart ... ! ************************************************************************************************** SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, & - basis_set_list_a, basis_set_list_b, sab_nl, symmetric) + basis_set_list_a, basis_set_list_b, sab_nl, symmetric, lcart) TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s @@ -1129,12 +1168,14 @@ CONTAINS TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_nl LOGICAL, INTENT(IN) :: symmetric + LOGICAL, INTENT(in), OPTIONAL :: lcart CHARACTER(LEN=12) :: cgfsym CHARACTER(LEN=32) :: symmetry_string CHARACTER(LEN=default_string_length) :: mname, name INTEGER :: i1, i2, natom INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes + LOGICAL :: my_lcart TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set @@ -1143,6 +1184,8 @@ CONTAINS qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist) natom = SIZE(particle_set) + my_lcart = .FALSE. + IF (PRESENT(lcart)) my_lcart = lcart IF (PRESENT(matrix_name)) THEN mname = matrix_name @@ -1152,10 +1195,17 @@ CONTAINS ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) - CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & - basis=basis_set_list_a) - CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & - basis=basis_set_list_b) + IF (.NOT. my_lcart) THEN + CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & + basis=basis_set_list_a) + CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & + basis=basis_set_list_b) + ELSE + CALL get_particle_set(particle_set, qs_kind_set, ncgf=row_blk_sizes, & + basis=basis_set_list_a) + CALL get_particle_set(particle_set, qs_kind_set, ncgf=col_blk_sizes, & + basis=basis_set_list_b) + END IF ! prepare for allocation IF (symmetric) THEN diff --git a/src/qs_scf_output.F b/src/qs_scf_output.F index 8a4f61fede..292cd37b00 100644 --- a/src/qs_scf_output.F +++ b/src/qs_scf_output.F @@ -220,6 +220,7 @@ CONTAINS TYPE(mp_para_env_type), POINTER :: para_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(preconditioner_type), POINTER :: local_preconditioner + TYPE(qs_environment_type), POINTER :: cart_overlap_qs_env TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(scf_control_type), POINTER :: scf_control TYPE(section_vals_type), POINTER :: dft_section, input @@ -434,6 +435,8 @@ CONTAINS END IF ! OT is used ! Print MO information + NULLIFY (cart_overlap_qs_env) + IF ((ikp == 1) .AND. (ispin == 1)) cart_overlap_qs_env => qs_env IF (nspin > 1) THEN SELECT CASE (ispin) CASE (1) @@ -446,19 +449,21 @@ CONTAINS IF (ASSOCIATED(umo_set)) THEN CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, & final_mos=final_mos, spin=TRIM(spin), solver_method=solver_method, & - umo_set=umo_set) + umo_set=umo_set, qs_env=cart_overlap_qs_env) ELSE CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, & - final_mos=final_mos, spin=TRIM(spin), solver_method=solver_method) + final_mos=final_mos, spin=TRIM(spin), solver_method=solver_method, & + qs_env=cart_overlap_qs_env) END IF ELSE IF (ASSOCIATED(umo_set)) THEN CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, & final_mos=final_mos, solver_method=solver_method, & - umo_set=umo_set) + umo_set=umo_set, qs_env=cart_overlap_qs_env) ELSE CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, & - final_mos=final_mos, solver_method=solver_method) + final_mos=final_mos, solver_method=solver_method, & + qs_env=cart_overlap_qs_env) END IF END IF diff --git a/tests/QS/regtest-gpw-1/TEST_FILES.toml b/tests/QS/regtest-gpw-1/TEST_FILES.toml index be210bb4ce..8356dca72e 100644 --- a/tests/QS/regtest-gpw-1/TEST_FILES.toml +++ b/tests/QS/regtest-gpw-1/TEST_FILES.toml @@ -72,4 +72,6 @@ "h2q.inp" = [{matcher="E_total", tol=1.0E-10, ref=-0.76162210587786}] # Mol Dipole Voronoi "moldip_voronoi.inp" = [{matcher="E_total", tol=1.0E-10, ref=-41.88176215391840}] +# Test printing of cartesian overlap matrix with header +"h2test.inp" = [{matcher="E_total", tol=3e-13, ref=-1.06531398809803}] #EOF diff --git a/tests/QS/regtest-gpw-1/h2test.inp b/tests/QS/regtest-gpw-1/h2test.inp new file mode 100644 index 0000000000..d2802d479c --- /dev/null +++ b/tests/QS/regtest-gpw-1/h2test.inp @@ -0,0 +1,73 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT h2 + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME EMSL_BASIS_SETS + BASIS_SET_FILE_NAME BASIS_MOLOPT_UCL + &MGRID + CUTOFF 400 + REL_CUTOFF 60 + &END MGRID + &POISSON + PERIODIC NONE + PSOLVER WAVELET + &END POISSON + &PRINT + &AO_MATRICES + OVERLAP + &END AO_MATRICES + &MO + CARTESIAN .TRUE. + CARTESIAN_OVERLAP .TRUE. + COEFFICIENTS .True. + FILENAME ./mos + MO_INDEX_RANGE 0 200 + &END MO + &END PRINT + &QS + EPS_DEFAULT 1.0E-14 + LMAXN0 4 + METHOD gapw + &END QS + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-6 + MAX_SCF 300 + SCF_GUESS RESTART + &DIAGONALIZATION + &END DIAGONALIZATION + &END SCF + &XC + &XC_FUNCTIONAL pbe + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6. 6. 6. + PERIODIC NONE + &END CELL + &COORD + h 0.00000000000000 0.00000000000000 0.58312867712598 + h 0.00000000000000 0.00000000000000 0.08687132338553 + &END COORD + &KIND H + BASIS_SET Ahlrichs-VDZ + ELEMENT H + HARD_EXP_RADIUS 1.00 + LEBEDEV_GRID 50 + POTENTIAL ALL + RADIAL_GRID 50 + &END KIND + &TOPOLOGY + &CENTER_COORDINATES T + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL