Mixed CDFT: Add recursive block diagonalization

svn-origin-rev: 18232
This commit is contained in:
Holmberg Nico 2018-01-25 13:31:44 +00:00
parent 4031ec9b8d
commit 368684fa46
5 changed files with 626 additions and 273 deletions

View file

@ -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

View file

@ -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)') &
'------------------------------------------------------------------------------'

View file

@ -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

View file

@ -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

View file

@ -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