diff --git a/src/input_cp2k_mixed.F b/src/input_cp2k_mixed.F index 6ab2e8d72f..0fad51423f 100644 --- a/src/input_cp2k_mixed.F +++ b/src/input_cp2k_mixed.F @@ -615,6 +615,17 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, "RECURSIVE_DIAGONALIZATION", & + description="Perform block diagonalization recursively until only two blocks remain. "// & + "For example, if the elements of a 8x8 matrix are first collected into 4 four blocks "// & + "(using keyword BLOCK), this keyword will block diagonalize the resulting 4x4 matrix "// & + "yielding a 2x2 matrix. In this example, the new blocks would be assembled from the "// & + " first and last 2 blocks of the previous matrix.", & + usage="RECURSIVE_DIAGONALIZATION TRUE", type_of_var=logical_t, & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_mixed_cdft_block_section END MODULE input_cp2k_mixed diff --git a/src/mixed_cdft_methods.F b/src/mixed_cdft_methods.F index a0fecb3b6e..4809f02ef5 100644 --- a/src/mixed_cdft_methods.F +++ b/src/mixed_cdft_methods.F @@ -85,14 +85,11 @@ MODULE mixed_cdft_methods mixed_cdft_type,& mixed_cdft_type_create,& mixed_cdft_work_type_release - USE mixed_cdft_utils, ONLY: hfun_zero,& - map_permutation_to_states,& - mixed_cdft_init_structures,& - mixed_cdft_parse_settings,& - mixed_cdft_print_couplings,& - mixed_cdft_redistribute_arrays,& - mixed_cdft_release_work,& - mixed_cdft_transfer_settings + USE mixed_cdft_utils, ONLY: & + hfun_zero, map_permutation_to_states, mixed_cdft_assemble_block_diag, & + mixed_cdft_diagonalize_blocks, mixed_cdft_get_blocks, mixed_cdft_init_structures, & + mixed_cdft_parse_settings, mixed_cdft_print_couplings, mixed_cdft_read_block_diag, & + mixed_cdft_redistribute_arrays, mixed_cdft_release_work, mixed_cdft_transfer_settings USE mixed_environment_types, ONLY: get_mixed_env,& mixed_environment_type,& set_mixed_env @@ -1785,6 +1782,8 @@ CONTAINS !> \param force_env the force_env that holds the CDFT states !> \par History !> 11.17 created [Nico Holmberg] +!> 01.18 added recursive diagonalization +!> split to subroutines [Nico Holmberg] ! ************************************************************************************************** SUBROUTINE mixed_cdft_block_diag(force_env) TYPE(force_env_type), POINTER :: force_env @@ -1792,24 +1791,15 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_block_diag', & routineP = moduleN//':'//routineN - CHARACTER(LEN=20) :: ilabel, jlabel - CHARACTER(LEN=3) :: tmp - INTEGER :: handle, i, icol, info, iounit, & - ipermutation, irow, j, k, l, nblk, & - nforce_eval, npermutations, & - work_array_size - INTEGER, DIMENSION(:), POINTER :: tmplist - LOGICAL :: explicit, has_duplicates, ignore_excited - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, H_mat_copy, H_offdiag, S_mat, & - S_mat_copy, S_offdiag + INTEGER :: handle, i, iounit, irecursion, j, n, & + nblk, nforce_eval, nrecursion + LOGICAL :: ignore_excited TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: H_block, S_block TYPE(cp_logger_type), POINTER :: logger TYPE(mixed_cdft_type), POINTER :: mixed_cdft - TYPE(section_vals_type), POINTER :: block_section, force_env_section, & - print_section + TYPE(section_vals_type), POINTER :: force_env_section, print_section EXTERNAL :: dsygv @@ -1823,269 +1813,75 @@ CONTAINS logger => cp_get_default_logger() CALL timeset(routineN, handle) - CALL force_env_get(force_env=force_env, & - force_env_section=force_env_section) - print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO") - iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog') - block_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%BLOCK_DIAGONALIZE") - CALL section_vals_get(block_section, explicit=explicit) - IF (.NOT. explicit) & - CALL cp_abort(__LOCATION__, & - "Block diagonalization of CDFT Hamiltonian was requested, but the "// & - "corresponding input section is missing!") CPASSERT(ALLOCATED(mixed_cdft%results%S)) CPASSERT(ALLOCATED(mixed_cdft%results%H)) nforce_eval = SIZE(mixed_cdft%results%S, 1) + + CALL force_env_get(force_env=force_env, & + force_env_section=force_env_section) + print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO") + iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog') ! Read block definitions from input - CALL section_vals_val_get(block_section, "BLOCK", n_rep_val=nblk) - ALLOCATE (blocks(nblk)) - DO i = 1, nblk - NULLIFY (blocks(i)%array) - CALL section_vals_val_get(block_section, "BLOCK", i_rep_val=i, i_vals=tmplist) - IF (SIZE(tmplist) < 1) & - CPABORT("Each BLOCK must contain at least 1 state.") - ALLOCATE (blocks(i)%array(SIZE(tmplist))) - blocks(i)%array(:) = tmplist(:) - END DO - CALL section_vals_val_get(block_section, "IGNORE_EXCITED", l_val=ignore_excited) - ! Check that the requested states exist - DO i = 1, nblk - DO j = 1, SIZE(blocks(i)%array) - IF (blocks(i)%array(j) < 1 .OR. blocks(i)%array(j) > nforce_eval) & - CPABORT("Requested state does not exist.") - END DO - END DO - ! Check for duplicates - has_duplicates = .FALSE. - DO i = 1, nblk - ! Within same block - DO j = 1, SIZE(blocks(i)%array) - DO k = j+1, SIZE(blocks(i)%array) - IF (blocks(i)%array(j) == blocks(i)%array(k)) has_duplicates = .TRUE. + CALL mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion) + nblk = SIZE(blocks) + ! Start block diagonalization + DO irecursion = 1, nrecursion + ! Print block definitions + IF (iounit > 0 .AND. irecursion == 1) THEN + WRITE (iounit, '(/,T3,A)') '-------------------------- CDFT BLOCK DIAGONALIZATION ------------------------' + WRITE (iounit, '(T3,A)') 'Block diagonalizing the mixed CDFT Hamiltonian' + WRITE (iounit, '(T3,A,I3)') 'Number of blocks:', nblk + WRITE (iounit, '(T3,A,L3)') 'Ignoring excited states within blocks:', ignore_excited + WRITE (iounit, '(/,T3,A)') 'List of CDFT states for each block' + DO i = 1, nblk + WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array)) END DO - END DO - ! Within different blocks - DO j = i+1, nblk - DO k = 1, SIZE(blocks(i)%array) - DO l = 1, SIZE(blocks(j)%array) - IF (blocks(i)%array(k) == blocks(j)%array(l)) has_duplicates = .TRUE. - END DO - END DO - END DO - END DO - IF (has_duplicates) CPABORT("Duplicate states are not allowed.") - ! Print block definitions - IF (iounit > 0) THEN - WRITE (iounit, '(/,T3,A)') '-------------------------- CDFT BLOCK DIAGONALIZATION ------------------------' - WRITE (iounit, '(T3,A)') 'Block diagonalizing the mixed CDFT Hamiltonian' - WRITE (iounit, '(T3,A,I3)') 'Number of blocks:', nblk - WRITE (iounit, '(T3,A,L3)') 'Ignoring excited states within blocks:', ignore_excited - WRITE (iounit, '(/,T3,A)') 'List of CDFT states for each block' - DO i = 1, nblk - WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array)) - END DO - END IF - ! Get the Hamiltonian and overlap matrices of each block - ALLOCATE (H_block(nblk), S_block(nblk)) - DO i = 1, nblk - NULLIFY (H_block(i)%array) - NULLIFY (S_block(i)%array) - ALLOCATE (H_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) - ALLOCATE (S_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) - icol = 0 - DO j = 1, SIZE(blocks(i)%array) - irow = 0 - icol = icol+1 - DO k = 1, SIZE(blocks(i)%array) - irow = irow+1 - H_block(i)%array(irow, icol) = mixed_cdft%results%H(blocks(i)%array(k), blocks(i)%array(j)) - S_block(i)%array(irow, icol) = mixed_cdft%results%S(blocks(i)%array(k), blocks(i)%array(j)) - END DO - END DO - ! Check that none of the interaction energies is repulsive - IF (ANY(H_block(i)%array .GE. 0.0_dp)) & - CALL cp_abort(__LOCATION__, & - "At least one of the interaction energies within block "//TRIM(ADJUSTL(cp_to_string(i)))// & - " is repulsive.") - END DO - ! Diagonalize blocks - ALLOCATE (eigenvalues(nblk)) - DO i = 1, nblk - NULLIFY (eigenvalues(i)%array) - ALLOCATE (eigenvalues(i)%array(SIZE(blocks(i)%array))) - eigenvalues(i)%array = 0.0_dp - ! Workspace query - ALLOCATE (work(1)) - info = 0 - ALLOCATE (H_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) - ALLOCATE (S_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) - H_mat_copy(:, :) = H_block(i)%array(:, :) ! Need explicit copies because dsygv destroys original values - S_mat_copy(:, :) = S_block(i)%array(:, :) - CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_mat_copy, SIZE(blocks(i)%array), & - S_mat_copy, SIZE(blocks(i)%array), eigenvalues(i)%array, work, -1, info) - work_array_size = NINT(work(1)) - DEALLOCATE (H_mat_copy, S_mat_copy) - ! Allocate work array - DEALLOCATE (work) - ALLOCATE (work(work_array_size)) - work = 0.0_dp - ! Solve Hc = eSc - info = 0 - CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_block(i)%array, SIZE(blocks(i)%array), & - S_block(i)%array, SIZE(blocks(i)%array), eigenvalues(i)%array, work, work_array_size, info) - IF (info /= 0) THEN - IF (info > nforce_eval) THEN - CPABORT("Matrix S is not positive definite") - ELSE - CPABORT("Diagonalization of H matrix failed.") - END IF END IF - DEALLOCATE (work) - END DO - ! Assemble the block diagonalized matrices - IF (ignore_excited) THEN - ALLOCATE (H_mat(nblk, nblk), S_mat(nblk, nblk)) - ELSE - ALLOCATE (H_mat(nforce_eval, nforce_eval), S_mat(nforce_eval, nforce_eval)) - END IF - ! The diagonal contains the eigenvalues of each block - IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Eigenvalues of the block diagonalized states" - H_mat(:, :) = 0.0_dp - S_mat(:, :) = 0.0_dp - k = 1 - DO i = 1, nblk - IF (iounit > 0) WRITE (iounit, '(T6,A,I3)') "Block", i - DO j = 1, SIZE(eigenvalues(i)%array) - H_mat(k, k) = eigenvalues(i)%array(j) - S_mat(k, k) = 1.0_dp - k = k+1 + ! Recursive diagonalization: update counters and references + IF (irecursion > 1) THEN + nblk = nblk/2 + ALLOCATE (blocks(nblk)) + j = 1 + DO i = 1, nblk + NULLIFY (blocks(i)%array) + ALLOCATE (blocks(i)%array(2)) + blocks(i)%array = (/j, j+1/) + j = j+2 + END DO + ! Print info IF (iounit > 0) THEN - IF (j == 1) THEN - WRITE (iounit, '(T9,A,T58,(3X,F20.14))') 'Ground state energy:', eigenvalues(i)%array(j) - ELSE - WRITE (iounit, '(T9,A,I2,A,T58,(3X,F20.14))') & - 'Excited state (', j-1, ' ) energy:', eigenvalues(i)%array(j) - END IF + WRITE (iounit, '(/, T3,A)') 'Recursive block diagonalization of the mixed CDFT Hamiltonian' + WRITE (iounit, '(T6,A)') 'Block diagonalization is continued until only two matrix blocks remain.' + WRITE (iounit, '(T6,A)') 'The new blocks are formed by collecting pairs of blocks from the previous' + WRITE (iounit, '(T6,A)') 'block diagonalized matrix in ascending order.' + WRITE (iounit, '(/,T3,A,I3,A,I3)') 'Recursion step:', irecursion-1, ' of ', nrecursion-1 + WRITE (iounit, '(/,T3,A)') 'List of old block indices for each new block' + DO i = 1, nblk + WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array)) + END DO END IF - IF (ignore_excited .AND. j == 1) EXIT - END DO - END DO - ! Transform the off-diagonal blocks using the eigenvectors of each block - npermutations = nblk*(nblk-1)/2 - IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Interactions between block diagonalized states" - DO ipermutation = 1, npermutations - CALL map_permutation_to_states(nblk, ipermutation, i, j) - ! Get the untransformed off-diagonal block - ALLOCATE (H_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array))) - ALLOCATE (S_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array))) - icol = 0 - DO k = 1, SIZE(blocks(j)%array) - irow = 0 - icol = icol+1 - DO l = 1, SIZE(blocks(i)%array) - irow = irow+1 - H_offdiag(irow, icol) = mixed_cdft%results%H(blocks(i)%array(l), blocks(j)%array(k)) - S_offdiag(irow, icol) = mixed_cdft%results%S(blocks(i)%array(l), blocks(j)%array(k)) - END DO - END DO - ! Check that none of the interaction energies is repulsive - IF (ANY(H_offdiag .GE. 0.0_dp)) & - CALL cp_abort(__LOCATION__, & - "At least one of the interaction energies between blocks "//TRIM(ADJUSTL(cp_to_string(i)))// & - " and "//TRIM(ADJUSTL(cp_to_string(j)))//" is repulsive.") - ! Now transform: C_i^T * H * C_j - H_offdiag(:, :) = MATMUL(H_offdiag, H_block(j)%array) - H_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), H_offdiag) - S_offdiag(:, :) = MATMUL(S_offdiag, H_block(j)%array) - S_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), S_offdiag) - ! Make sure the transformation preserves the sign of elements in the S and H matrices - ! The S/H matrices contain only positive/negative values so that any sign flipping occurs in the - ! same elements in both matrices - ! Check for sign flipping using the S matrix - IF (ANY(S_offdiag .LT. 0.0_dp)) THEN - DO l = 1, SIZE(S_offdiag, 2) - DO k = 1, SIZE(S_offdiag, 1) - IF (S_offdiag(k, l) .LT. 0.0_dp) THEN - S_offdiag(k, l) = -1.0_dp*S_offdiag(k, l) - H_offdiag(k, l) = -1.0_dp*H_offdiag(k, l) - END IF - END DO - END DO END IF + ! Get the Hamiltonian and overlap matrices of each block + CALL mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block) + ! Diagonalize blocks + CALL mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues) + ! Assemble the block diagonalized matrices IF (ignore_excited) THEN - H_mat(i, j) = H_offdiag(1, 1) - H_mat(j, i) = H_mat(i, j) - S_mat(i, j) = S_offdiag(1, 1) - S_mat(j, i) = S_mat(i, j) + n = nblk ELSE - irow = 1 - icol = 1 - DO k = 1, i-1 - irow = irow+SIZE(blocks(k)%array) - END DO - DO k = 1, j-1 - icol = icol+SIZE(blocks(k)%array) - END DO - H_mat(irow:irow+SIZE(H_offdiag, 1)-1, icol:icol+SIZE(H_offdiag, 2)-1) = H_offdiag(:, :) - H_mat(icol:icol+SIZE(H_offdiag, 2)-1, irow:irow+SIZE(H_offdiag, 1)-1) = TRANSPOSE(H_offdiag) - S_mat(irow:irow+SIZE(H_offdiag, 1)-1, icol:icol+SIZE(H_offdiag, 2)-1) = S_offdiag(:, :) - S_mat(icol:icol+SIZE(H_offdiag, 2)-1, irow:irow+SIZE(H_offdiag, 1)-1) = TRANSPOSE(S_offdiag) + n = nforce_eval END IF - IF (iounit > 0) THEN - WRITE (iounit, '(/,T3,A)') REPEAT('#', 39) - WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### Blocks I =', i, ' and J = ', j, ' ######' - WRITE (iounit, '(T3,A)') REPEAT('#', 39) - WRITE (iounit, '(T3,A)') 'Interaction energies' - DO irow = 1, SIZE(H_offdiag, 1) - ilabel = "(ground state)" - IF (irow > 1) THEN - IF (ignore_excited) EXIT - WRITE (tmp, '(I3)') irow-1 - ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" - END IF - DO icol = 1, SIZE(H_offdiag, 2) - jlabel = "(ground state)" - IF (icol > 1) THEN - IF (ignore_excited) EXIT - WRITE (tmp, '(I3)') icol-1 - jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" - END IF - WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', H_offdiag(irow, icol) - END DO - END DO - WRITE (iounit, '(T3,A)') 'Overlaps' - DO irow = 1, SIZE(H_offdiag, 1) - ilabel = "(ground state)" - IF (irow > 1) THEN - IF (ignore_excited) EXIT - ilabel = "(excited state)" - WRITE (tmp, '(I3)') irow-1 - ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" - END IF - DO icol = 1, SIZE(H_offdiag, 2) - jlabel = "(ground state)" - IF (icol > 1) THEN - IF (ignore_excited) EXIT - WRITE (tmp, '(I3)') icol-1 - jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" - END IF - WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', S_offdiag(irow, icol) - END DO - END DO - END IF - DEALLOCATE (H_offdiag, S_offdiag) - END DO - CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat, S=S_mat) - ! Deallocate work - DEALLOCATE (H_mat, S_mat) - DO i = 1, nblk - DEALLOCATE (H_block(i)%array) - DEALLOCATE (S_block(i)%array) - DEALLOCATE (eigenvalues(i)%array) - DEALLOCATE (blocks(i)%array) - END DO - DEALLOCATE (H_block, S_block, eigenvalues, blocks) + CALL mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, n, iounit) + ! Deallocate work + DO i = 1, nblk + DEALLOCATE (H_block(i)%array) + DEALLOCATE (S_block(i)%array) + DEALLOCATE (eigenvalues(i)%array) + DEALLOCATE (blocks(i)%array) + END DO + DEALLOCATE (H_block, S_block, eigenvalues, blocks) + END DO ! recursion IF (iounit > 0) & WRITE (iounit, '(T3,A)') & '------------------------------------------------------------------------------' diff --git a/src/mixed_cdft_utils.F b/src/mixed_cdft_utils.F index 27571e1e77..7f46b781fa 100644 --- a/src/mixed_cdft_utils.F +++ b/src/mixed_cdft_utils.F @@ -12,7 +12,9 @@ MODULE mixed_cdft_utils USE atomic_kind_types, ONLY: atomic_kind_type USE cell_types, ONLY: cell_type - USE cp_array_utils, ONLY: cp_1d_r_p_type + USE cp_array_utils, ONLY: cp_1d_i_p_type,& + cp_1d_r_p_type,& + cp_2d_r_p_type USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& cp_blacs_env_retain,& @@ -61,6 +63,7 @@ MODULE mixed_cdft_utils USE input_constants, ONLY: becke_cutoff_element,& shape_function_gaussian USE input_section_types, ONLY: section_vals_duplicate,& + section_vals_get,& section_vals_get_subs_vals,& section_vals_release,& section_vals_type,& @@ -73,6 +76,7 @@ MODULE mixed_cdft_utils mp_wait,& mp_waitall USE mixed_cdft_types, ONLY: mixed_cdft_result_type_release,& + mixed_cdft_result_type_set,& mixed_cdft_settings_type,& mixed_cdft_type,& mixed_cdft_work_type_init @@ -111,7 +115,9 @@ MODULE mixed_cdft_utils PUBLIC :: mixed_cdft_parse_settings, mixed_cdft_transfer_settings, & mixed_cdft_init_structures, mixed_cdft_redistribute_arrays, & mixed_cdft_print_couplings, map_permutation_to_states, hfun_zero, & - mixed_cdft_release_work + mixed_cdft_release_work, mixed_cdft_read_block_diag, & + mixed_cdft_get_blocks, mixed_cdft_diagonalize_blocks, & + mixed_cdft_assemble_block_diag CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_utils' @@ -1497,4 +1503,379 @@ CONTAINS END SUBROUTINE hfun_zero +! ************************************************************************************************** +!> \brief Read input section related to block diagonalization of the mixed CDFT Hamiltonian matrix. +!> \param force_env the force_env that holds the CDFT states +!> \param blocks list of CDFT states defining the matrix blocks +!> \param ignore_excited flag that determines if excited states resulting from the block +!> diagonalization process should be ignored +!> \param nrecursion integer that determines how many steps of recursive block diagonalization +!> is performed (1 if disabled) +!> \par History +!> 01.18 created [Nico Holmberg] +! ************************************************************************************************** + SUBROUTINE mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion) + TYPE(force_env_type), POINTER :: force_env + TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:), & + INTENT(OUT) :: blocks + LOGICAL, INTENT(OUT) :: ignore_excited + INTEGER, INTENT(OUT) :: nrecursion + + CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_read_block_diag', & + routineP = moduleN//':'//routineN + + INTEGER :: i, j, k, l, nblk, nforce_eval + INTEGER, DIMENSION(:), POINTER :: tmplist + LOGICAL :: do_recursive, explicit, has_duplicates + TYPE(section_vals_type), POINTER :: block_section, force_env_section + + EXTERNAL :: dsygv + + NULLIFY (force_env_section, block_section) + CPASSERT(ASSOCIATED(force_env)) + nforce_eval = SIZE(force_env%sub_force_env) + + CALL force_env_get(force_env=force_env, & + force_env_section=force_env_section) + block_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%BLOCK_DIAGONALIZE") + + CALL section_vals_get(block_section, explicit=explicit) + IF (.NOT. explicit) & + CALL cp_abort(__LOCATION__, & + "Block diagonalization of CDFT Hamiltonian was requested, but the "// & + "corresponding input section is missing!") + + CALL section_vals_val_get(block_section, "BLOCK", n_rep_val=nblk) + ALLOCATE (blocks(nblk)) + DO i = 1, nblk + NULLIFY (blocks(i)%array) + CALL section_vals_val_get(block_section, "BLOCK", i_rep_val=i, i_vals=tmplist) + IF (SIZE(tmplist) < 1) & + CPABORT("Each BLOCK must contain at least 1 state.") + ALLOCATE (blocks(i)%array(SIZE(tmplist))) + blocks(i)%array(:) = tmplist(:) + END DO + CALL section_vals_val_get(block_section, "IGNORE_EXCITED", l_val=ignore_excited) + CALL section_vals_val_get(block_section, "RECURSIVE_DIAGONALIZATION", l_val=do_recursive) + ! Check that the requested states exist + DO i = 1, nblk + DO j = 1, SIZE(blocks(i)%array) + IF (blocks(i)%array(j) < 1 .OR. blocks(i)%array(j) > nforce_eval) & + CPABORT("Requested state does not exist.") + END DO + END DO + ! Check for duplicates + has_duplicates = .FALSE. + DO i = 1, nblk + ! Within same block + DO j = 1, SIZE(blocks(i)%array) + DO k = j+1, SIZE(blocks(i)%array) + IF (blocks(i)%array(j) == blocks(i)%array(k)) has_duplicates = .TRUE. + END DO + END DO + ! Within different blocks + DO j = i+1, nblk + DO k = 1, SIZE(blocks(i)%array) + DO l = 1, SIZE(blocks(j)%array) + IF (blocks(i)%array(k) == blocks(j)%array(l)) has_duplicates = .TRUE. + END DO + END DO + END DO + END DO + IF (has_duplicates) CPABORT("Duplicate states are not allowed.") + nrecursion = 1 + IF (do_recursive) THEN + IF (MODULO(nblk, 2) /= 0) THEN + CALL cp_warn(__LOCATION__, & + "Number of blocks not divisible with 2. Recursive diagonalization not possible. "// & + "Calculation proceeds without.") + nrecursion = 1 + ELSE + nrecursion = nblk/2 + END IF + IF (nrecursion /= 1 .AND. .NOT. ignore_excited) & + CALL cp_abort(__LOCATION__, & + "Keyword IGNORE_EXCITED must be active for recursive diagonalization.") + END IF + + END SUBROUTINE mixed_cdft_read_block_diag + +! ************************************************************************************************** +!> \brief Assembles the matrix blocks from the mixed CDFT Hamiltonian. +!> \param mixed_cdft the env that holds the CDFT states +!> \param blocks list of CDFT states defining the matrix blocks +!> \param H_block list of Hamiltonian matrix blocks +!> \param S_block list of overlap matrix blocks +!> \par History +!> 01.18 created [Nico Holmberg] +! ************************************************************************************************** + SUBROUTINE mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block) + TYPE(mixed_cdft_type), POINTER :: mixed_cdft + TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks + TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:), & + INTENT(OUT) :: H_block, S_block + + CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_get_blocks', & + routineP = moduleN//':'//routineN + + INTEGER :: i, icol, irow, j, k, nblk + + EXTERNAL :: dsygv + + CPASSERT(ASSOCIATED(mixed_cdft)) + + nblk = SIZE(blocks) + ALLOCATE (H_block(nblk), S_block(nblk)) + DO i = 1, nblk + NULLIFY (H_block(i)%array) + NULLIFY (S_block(i)%array) + ALLOCATE (H_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) + ALLOCATE (S_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) + icol = 0 + DO j = 1, SIZE(blocks(i)%array) + irow = 0 + icol = icol+1 + DO k = 1, SIZE(blocks(i)%array) + irow = irow+1 + H_block(i)%array(irow, icol) = mixed_cdft%results%H(blocks(i)%array(k), blocks(i)%array(j)) + S_block(i)%array(irow, icol) = mixed_cdft%results%S(blocks(i)%array(k), blocks(i)%array(j)) + END DO + END DO + ! Check that none of the interaction energies is repulsive + IF (ANY(H_block(i)%array .GE. 0.0_dp)) & + CALL cp_abort(__LOCATION__, & + "At least one of the interaction energies within block "//TRIM(ADJUSTL(cp_to_string(i)))// & + " is repulsive.") + END DO + + END SUBROUTINE mixed_cdft_get_blocks + +! ************************************************************************************************** +!> \brief Diagonalizes each of the matrix blocks. +!> \param blocks list of CDFT states defining the matrix blocks +!> \param H_block list of Hamiltonian matrix blocks +!> \param S_block list of overlap matrix blocks +!> \param eigenvalues list of eigenvalues for each block +!> \par History +!> 01.18 created [Nico Holmberg] +! ************************************************************************************************** + SUBROUTINE mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues) + TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks + TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: H_block, S_block + TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:), & + INTENT(OUT) :: eigenvalues + + CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_diagonalize_blocks', & + routineP = moduleN//':'//routineN + + INTEGER :: i, info, nblk, work_array_size + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat_copy, S_mat_copy + + EXTERNAL :: dsygv + + nblk = SIZE(blocks) + ALLOCATE (eigenvalues(nblk)) + DO i = 1, nblk + NULLIFY (eigenvalues(i)%array) + ALLOCATE (eigenvalues(i)%array(SIZE(blocks(i)%array))) + eigenvalues(i)%array = 0.0_dp + ! Workspace query + ALLOCATE (work(1)) + info = 0 + ALLOCATE (H_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) + ALLOCATE (S_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array))) + H_mat_copy(:, :) = H_block(i)%array(:, :) ! Need explicit copies because dsygv destroys original values + S_mat_copy(:, :) = S_block(i)%array(:, :) + CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_mat_copy, SIZE(blocks(i)%array), & + S_mat_copy, SIZE(blocks(i)%array), eigenvalues(i)%array, work, -1, info) + work_array_size = NINT(work(1)) + DEALLOCATE (H_mat_copy, S_mat_copy) + ! Allocate work array + DEALLOCATE (work) + ALLOCATE (work(work_array_size)) + work = 0.0_dp + ! Solve Hc = eSc + info = 0 + CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_block(i)%array, SIZE(blocks(i)%array), & + S_block(i)%array, SIZE(blocks(i)%array), eigenvalues(i)%array, work, work_array_size, info) + IF (info /= 0) THEN + IF (info > SIZE(blocks(i)%array)) THEN + CPABORT("Matrix S is not positive definite") + ELSE + CPABORT("Diagonalization of H matrix failed.") + END IF + END IF + DEALLOCATE (work) + END DO + + END SUBROUTINE mixed_cdft_diagonalize_blocks + +! ************************************************************************************************** +!> \brief Assembles the new block diagonalized mixed CDFT Hamiltonian and overlap matrices. +!> \param mixed_cdft the env that holds the CDFT states +!> \param blocks list of CDFT states defining the matrix blocks +!> \param H_block list of Hamiltonian matrix blocks +!> \param eigenvalues list of eigenvalues for each block +!> \param n size of the new Hamiltonian and overlap matrices +!> \param iounit the output unit +!> \par History +!> 01.18 created [Nico Holmberg] +! ************************************************************************************************** + SUBROUTINE mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, & + n, iounit) + TYPE(mixed_cdft_type), POINTER :: mixed_cdft + TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks + TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: H_block + TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues + INTEGER :: n, iounit + + CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_assemble_block_diag', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=20) :: ilabel, jlabel + CHARACTER(LEN=3) :: tmp + INTEGER :: i, icol, ipermutation, irow, j, k, l, & + nblk, npermutations + LOGICAL :: ignore_excited + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, H_offdiag, S_mat, S_offdiag + + EXTERNAL :: dsygv + + ALLOCATE (H_mat(n, n), S_mat(n, n)) + nblk = SIZE(blocks) + ignore_excited = (nblk == n) + ! The diagonal contains the eigenvalues of each block + IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Eigenvalues of the block diagonalized states" + H_mat(:, :) = 0.0_dp + S_mat(:, :) = 0.0_dp + k = 1 + DO i = 1, nblk + IF (iounit > 0) WRITE (iounit, '(T6,A,I3)') "Block", i + DO j = 1, SIZE(eigenvalues(i)%array) + H_mat(k, k) = eigenvalues(i)%array(j) + S_mat(k, k) = 1.0_dp + k = k+1 + IF (iounit > 0) THEN + IF (j == 1) THEN + WRITE (iounit, '(T9,A,T58,(3X,F20.14))') 'Ground state energy:', eigenvalues(i)%array(j) + ELSE + WRITE (iounit, '(T9,A,I2,A,T58,(3X,F20.14))') & + 'Excited state (', j-1, ' ) energy:', eigenvalues(i)%array(j) + END IF + END IF + IF (ignore_excited .AND. j == 1) EXIT + END DO + END DO + ! Transform the off-diagonal blocks using the eigenvectors of each block + npermutations = nblk*(nblk-1)/2 + IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Interactions between block diagonalized states" + DO ipermutation = 1, npermutations + CALL map_permutation_to_states(nblk, ipermutation, i, j) + ! Get the untransformed off-diagonal block + ALLOCATE (H_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array))) + ALLOCATE (S_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array))) + icol = 0 + DO k = 1, SIZE(blocks(j)%array) + irow = 0 + icol = icol+1 + DO l = 1, SIZE(blocks(i)%array) + irow = irow+1 + H_offdiag(irow, icol) = mixed_cdft%results%H(blocks(i)%array(l), blocks(j)%array(k)) + S_offdiag(irow, icol) = mixed_cdft%results%S(blocks(i)%array(l), blocks(j)%array(k)) + END DO + END DO + ! Check that none of the interaction energies is repulsive + IF (ANY(H_offdiag .GE. 0.0_dp)) & + CALL cp_abort(__LOCATION__, & + "At least one of the interaction energies between blocks "//TRIM(ADJUSTL(cp_to_string(i)))// & + " and "//TRIM(ADJUSTL(cp_to_string(j)))//" is repulsive.") + ! Now transform: C_i^T * H * C_j + H_offdiag(:, :) = MATMUL(H_offdiag, H_block(j)%array) + H_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), H_offdiag) + S_offdiag(:, :) = MATMUL(S_offdiag, H_block(j)%array) + S_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), S_offdiag) + ! Make sure the transformation preserves the sign of elements in the S and H matrices + ! The S/H matrices contain only positive/negative values so that any sign flipping occurs in the + ! same elements in both matrices + ! Check for sign flipping using the S matrix + IF (ANY(S_offdiag .LT. 0.0_dp)) THEN + DO l = 1, SIZE(S_offdiag, 2) + DO k = 1, SIZE(S_offdiag, 1) + IF (S_offdiag(k, l) .LT. 0.0_dp) THEN + S_offdiag(k, l) = -1.0_dp*S_offdiag(k, l) + H_offdiag(k, l) = -1.0_dp*H_offdiag(k, l) + END IF + END DO + END DO + END IF + IF (ignore_excited) THEN + H_mat(i, j) = H_offdiag(1, 1) + H_mat(j, i) = H_mat(i, j) + S_mat(i, j) = S_offdiag(1, 1) + S_mat(j, i) = S_mat(i, j) + ELSE + irow = 1 + icol = 1 + DO k = 1, i-1 + irow = irow+SIZE(blocks(k)%array) + END DO + DO k = 1, j-1 + icol = icol+SIZE(blocks(k)%array) + END DO + H_mat(irow:irow+SIZE(H_offdiag, 1)-1, icol:icol+SIZE(H_offdiag, 2)-1) = H_offdiag(:, :) + H_mat(icol:icol+SIZE(H_offdiag, 2)-1, irow:irow+SIZE(H_offdiag, 1)-1) = TRANSPOSE(H_offdiag) + S_mat(irow:irow+SIZE(H_offdiag, 1)-1, icol:icol+SIZE(H_offdiag, 2)-1) = S_offdiag(:, :) + S_mat(icol:icol+SIZE(H_offdiag, 2)-1, irow:irow+SIZE(H_offdiag, 1)-1) = TRANSPOSE(S_offdiag) + END IF + IF (iounit > 0) THEN + WRITE (iounit, '(/,T3,A)') REPEAT('#', 39) + WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### Blocks I =', i, ' and J = ', j, ' ######' + WRITE (iounit, '(T3,A)') REPEAT('#', 39) + WRITE (iounit, '(T3,A)') 'Interaction energies' + DO irow = 1, SIZE(H_offdiag, 1) + ilabel = "(ground state)" + IF (irow > 1) THEN + IF (ignore_excited) EXIT + WRITE (tmp, '(I3)') irow-1 + ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" + END IF + DO icol = 1, SIZE(H_offdiag, 2) + jlabel = "(ground state)" + IF (icol > 1) THEN + IF (ignore_excited) EXIT + WRITE (tmp, '(I3)') icol-1 + jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" + END IF + WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', H_offdiag(irow, icol) + END DO + END DO + WRITE (iounit, '(T3,A)') 'Overlaps' + DO irow = 1, SIZE(H_offdiag, 1) + ilabel = "(ground state)" + IF (irow > 1) THEN + IF (ignore_excited) EXIT + ilabel = "(excited state)" + WRITE (tmp, '(I3)') irow-1 + ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" + END IF + DO icol = 1, SIZE(H_offdiag, 2) + jlabel = "(ground state)" + IF (icol > 1) THEN + IF (ignore_excited) EXIT + WRITE (tmp, '(I3)') icol-1 + jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")" + END IF + WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', S_offdiag(irow, icol) + END DO + END DO + END IF + DEALLOCATE (H_offdiag, S_offdiag) + END DO + CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat, S=S_mat) + ! Deallocate work + DEALLOCATE (H_mat, S_mat) + + END SUBROUTINE mixed_cdft_assemble_block_diag + END MODULE mixed_cdft_utils diff --git a/tests/QS/regtest-cdft-5/He2H-mixed-cdft-4.inp b/tests/QS/regtest-cdft-5/He2H-mixed-cdft-4.inp new file mode 100644 index 0000000000..46314dca0d --- /dev/null +++ b/tests/QS/regtest-cdft-5/He2H-mixed-cdft-4.inp @@ -0,0 +1,164 @@ +@SET RESTART_WFN TRUE +@SET WFN_FILE_1 He2H-cdft-state-1a-1_0.wfn +@SET WFN_FILE_1b He2H-cdft-state-1b-1_0.wfn +@SET WFN_FILE_2 He2H-cdft-state-2a-1_0.wfn +@SET WFN_FILE_2b He2H-cdft-state-2b-1_0.wfn +@SET PROJECT_NAME He2H-mixed-cdft-4 +@SET NAME ${PROJECT_NAME} +@SET WRITE_WFN 0 +@SET CHARGE 1 +@SET WRITE_CUBE FALSE +@SET XYZFILE He2H.xyz + +@SET BECKE_ACTIVE TRUE +@SET BECKE_FRAGMENT FALSE +@SET MAX_SCF 5 +! He+ [H He]0 +@SET BECKE_TARGET_1 1.2 +@SET BECKE_STR_1 0.561093953718 +@SET BECKE_TARGET_1b 0.8 +@SET BECKE_STR_1b 2.744257349834 +! He0 [H He]+ +@SET BECKE_TARGET_2 2.4 +@SET BECKE_STR_2 -3.951144435823 +@SET BECKE_TARGET_2b 2.3 +@SET BECKE_STR_2b -0.538906046282 + +! Use slightly perturbed values to avoid linear dependencies +! He+ [H He]0 +@SET BECKE_TARGET_1 1.2 +@SET BECKE_STR_1P 0.461093953718 +@SET BECKE_TARGET_1b 0.8 +@SET BECKE_STR_1bP 2.644257349834 +! He0 [H He]+ +@SET BECKE_TARGET_2 2.4 +@SET BECKE_STR_2P -3.851144435823 +@SET BECKE_TARGET_2b 2.3 +@SET BECKE_STR_2bP -0.438906046282 + +@SET BECKE_GLOBAL_CUTOFF FALSE +@SET BECKE_CUTOFF_ELEMENT TRUE + +@SET BECKE_ADJUST_SIZE TRUE + +@SET BECKE_ATOMIC_CHARGES FALSE + +@SET BECKE_CAVITY_CONFINE TRUE +@SET BECKE_CAVITY_SHAPE VDW +@SET BECKE_CAVITY_PRINT FALSE + +@SET BECKE_SHOULD_SKIP TRUE + +@SET BECKE_IN_MEMORY TRUE + +&GLOBAL + PROJECT ${PROJECT_NAME} + RUN_TYPE ENERGY + PRINT_LEVEL MEDIUM +&END GLOBAL + +&MULTIPLE_FORCE_EVALS + FORCE_EVAL_ORDER 2 3 4 5 6 7 8 9 + MULTIPLE_SUBSYS F +&END +&FORCE_EVAL + METHOD MIXED + &MIXED + MIXING_TYPE MIXED_CDFT + NGROUPS 1 + &MIXED_CDFT + LAMBDA 1.0 + COUPLING 1 + CI TRUE + NONORTHO_COUPLING TRUE + BLOCK_DIAGONALIZE TRUE + &BLOCK_DIAGONALIZE + BLOCK 1 2 + BLOCK 3 4 + BLOCK 5 6 + BLOCK 7 8 + ! Recursive diagonalization: + ! Blocks 1+2 and 3+4 get combined into new blocks + ! and the block diagonalization is repeated + RECURSIVE_DIAGONALIZATION TRUE + &END BLOCK_DIAGONALIZE + &END MIXED_CDFT + &PRINT + &PROGRAM_RUN_INFO + &END + &END PRINT + &END MIXED + @include subsys.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_1} + @SET BECKE_TARGET ${BECKE_TARGET_1} + @SET PROJECT_NAME ${NAME}-state-1a + @SET WFN_FILE ${WFN_FILE_1} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_1b} + @SET BECKE_TARGET ${BECKE_TARGET_1b} + @SET PROJECT_NAME ${NAME}-state-1b + @SET WFN_FILE ${WFN_FILE_1b} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_1P} + @SET BECKE_TARGET ${BECKE_TARGET_1} + @SET PROJECT_NAME ${NAME}-state-1aP + @SET WFN_FILE ${WFN_FILE_1} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_1bP} + @SET BECKE_TARGET ${BECKE_TARGET_1b} + @SET PROJECT_NAME ${NAME}-state-1bP + @SET WFN_FILE ${WFN_FILE_1b} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_2} + @SET BECKE_TARGET ${BECKE_TARGET_2} + @SET PROJECT_NAME ${NAME}-state-2a + @SET WFN_FILE ${WFN_FILE_2} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_2b} + @SET BECKE_TARGET ${BECKE_TARGET_2b} + @SET PROJECT_NAME ${NAME}-state-2b + @SET WFN_FILE ${WFN_FILE_2b} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_2P} + @SET BECKE_TARGET ${BECKE_TARGET_2} + @SET PROJECT_NAME ${NAME}-state-2aP + @SET WFN_FILE ${WFN_FILE_2} + @include dft-common-params.inc +&END FORCE_EVAL + +&FORCE_EVAL + METHOD QS + @SET BECKE_STR ${BECKE_STR_2bP} + @SET BECKE_TARGET ${BECKE_TARGET_2b} + @SET PROJECT_NAME ${NAME}-state-2bP + @SET WFN_FILE ${WFN_FILE_2b} + @include dft-common-params.inc +&END FORCE_EVAL diff --git a/tests/QS/regtest-cdft-5/TEST_FILES b/tests/QS/regtest-cdft-5/TEST_FILES index 70dddac939..6308bc1eb0 100644 --- a/tests/QS/regtest-cdft-5/TEST_FILES +++ b/tests/QS/regtest-cdft-5/TEST_FILES @@ -11,4 +11,5 @@ He2H-cdft-state-2b.inp 71 2e-9 He2H-mixed-cdft-1.inp 77 5e-8 -6.84813719910651 He2H-mixed-cdft-2.inp 77 2e-7 -7.90368522838573 He2H-mixed-cdft-3.inp 77 5e-8 -5.78009324749503 +He2H-mixed-cdft-4.inp 77 5e-8 -6.73674737139200 #EOF