analyze_full_wavefunction_sym Subroutine

public subroutine analyze_full_wavefunction_sym(sym_labels, ilut_list_opt)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(out), allocatable :: sym_labels(:)
integer(kind=n_int), intent(in), optional :: ilut_list_opt(:,:)

Contents


Source Code

    subroutine analyze_full_wavefunction_sym(sym_labels, ilut_list_opt)
        ! give the symmetry eigenvalues of a wavefunction in ilut_format
        ! of certain symmetry operations, depending on the lattice point group
        ! e.g. for 2D: rot90, rot180, rot270, m_x, m_y, m_d, m_o and the S^2 EV
        ! on can either provide a wavfunction in ilut_list, or otherwise
        ! it is obtained from the FCIQMC wavefunction depending on the
        ! input
        real(dp), intent(out), allocatable :: sym_labels(:)
        integer(n_int), intent(in), optional :: ilut_list_opt(:, :)
        character(*), parameter :: this_routine = "analyze_full_wavefunction_sym"
        integer(n_int), allocatable :: ilut_list(:, :), trans_pg_wf(:, :, :), &
                                       ilut_spin(:, :)
        integer :: n_syms, i, ms, nI(nel)

        ASSERT(associated(lat))

        if (lat%get_ndim() == 2) then
            ! for 2D it is the 8 fold point group symmetry and S^2 EV..
            n_syms = 10
            ! activate all the symmetries
            t_symmetry_mirror = .true.
            t_symmetry_rotation = .true.
            t_symmetry_inversion = .true.
        else
            call stop_all(this_routine, "not yet implemented!")
        end if

        allocate(sym_labels(n_syms))
        sym_labels = 0.0_dp

        if (present(ilut_list_opt)) then
            allocate(ilut_list(0:niftot, size(ilut_list_opt, 2)), &
                      source=ilut_list_opt)

        else
            call init_symmetry_states()
            allocate(ilut_list(0:NIfTot, size(symmetry_states_ilut, 2)), &
                      source=symmetry_states_ilut)

        end if

        if (lat%get_ndim() == 2) then
            trans_pg_wf = apply_2D_point_group(ilut_list)
        else
            call stop_all(this_routine, "not yet implemented!")
        end if

        ! now we have to calculate <y|y'> to get the symmetry EV
        ! but dont do the S^2 yet..
        do i = 1, n_syms - 1
            sym_labels(i) = calc_overlap(ilut_list, trans_pg_wf(:, :, i))
        end do

        do i = 1, size(ilut_list, 2)
            sym_labels(n_syms) = sym_labels(n_syms) + &
                                 ssquared_contrib(ilut_list(:, i), .false., size(ilut_list, 2), &
                                                  ilut_list)
        end do

        ! and we need the S_z^2 contribution:
        call decode_bit_det(nI, ilut_list(:, 1))
        ms = sum(get_spin_pn(nI))
        sym_labels(n_syms) = sym_labels(n_syms) + real(ms * (ms + 2) / 4.0, dp)

    end subroutine analyze_full_wavefunction_sym