"""
Fortran module VariationalOps: July 2017 (dj, updated)
Containing the variational algorithms for tensor networks.
**Authors**
* D. Jaschke
* M. L. Wall
**Details**
The following subroutines / functions are defined for the
symmetric tensors.
+-----------------------------------------+-------------+---------+---------+
| procedure | include.f90 | mpi.f90 | omp.f90 |
+=========================================+=============+=========+=========+
| find_groundstate_two | X | | |
+-----------------------------------------+-------------+---------+---------+
| grow_system | X | | |
+-----------------------------------------+-------------+---------+---------+
| Fit_HPsi | X | | |
+-----------------------------------------+-------------+---------+---------+
| FitMPS | X | | |
+-----------------------------------------+-------------+---------+---------+
| orth | X | | |
+-----------------------------------------+-------------+---------+---------+
"""
[docs]def find_excited_two_mps_mpo():
"""
fortran-subroutine - ?? (mlw)
Find the excited state Psi (represented as an MPS) of the Hamiltonian H
(represented as an MPO) variationally. Psi should contain the initial
guess on input.
**Arguments**
Psi : TYPE(mps), inout
On entry the initial guess, on exit the excited state.
Psibs : TYPE(mps), inout
Energetically lower lying states which are required to
orthogonalize against those states.
ne : INTEGER, in
Excited state to be found.
Hts : TYPE(sr_matrix_tensor)(\*), inout
List of all nearest-neighbor two-site MPO matrices.
Ham : TYPE(mpo), inout
Hamiltonian for the system, represented as MPO.
LR : TYPE(tensorlist)(\*), inout
Left right overlaps with Hamiltonian / MPO.
LRAB : TYPE(tensorlist)(\*), inout
Left right overlaps with energetically lower lying states.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the excited state
search.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_excited_two_mps_mpo(Psi, Psibs, ne, Hts, &
Ham, Cp, energy, converged, conv_val, pbc, errst)
type(mps), intent(inout) :: Psi
type(mps), pointer, intent(inout) :: Psibs(:)
integer, intent(in) :: ne
type(sr_matrix_tensor), dimension(:), intent(inout) :: Hts
type(mpo), intent(inout) :: Ham
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! MPO squared
type(mpo) :: Hsq
! Overlap state-Hamiltonian and overlap state-state
type(tensorlist), dimension(:), allocatable :: LR, LRAB
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_excited_two_mps_'//&
! 'mpo: variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:134', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site excited state.'
call setuplr(LR, Psi, Ham, Psi, Psi%oc)
allocate(LRAB(ne))
do ii = 1, ne
call setuplr(LRAB(ii), Psi, Psibs(ii), Psi%oc)
end do
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print - (energy: there was no iMPS growth before)
energy = 0.0_rKind
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'sweep', jj, &
'bd', chi, 'energy', energy
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call find_excited_two_step_mps_mpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mps_'//&
! 'mpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:194', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call find_excited_two_step_mps_mpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mps_'//&
! 'mpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:205', &
! errst=errst)) return
end do
Psi%oc = 2
sense = 1
do kk = 1, (k1 - 1)
call find_excited_two_step_mps_mpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mps_'//&
! 'mpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:217', &
! errst=errst)) return
end do
Psi%oc = k1
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'bond dimension, variance'//&
', last variance, difference in variance, discarded '//&
'weight', chi, conv_val, varianceSave, &
abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
if((verbose > 0) .and. (.not. converged)) write(slog, *) &
'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
if((verbose > 0) .and. converged) write(slog, *) &
'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
do ii = 1, size(LRAB, 1)
do jj = 1, size(LRAB(ii)%Li)
call destroy(LRAB(ii)%Li(jj))
end do
deallocate(LRAB(ii)%Li)
end do
deallocate(LR, LRAB)
if(Cp%conv_method == 'S') deallocate(firstlambda)
end subroutine find_excited_two_mps_mpo
"""
return
[docs]def find_excited_two_mpsc_mpoc():
"""
fortran-subroutine - ?? (mlw)
Find the excited state Psi (represented as an MPS) of the Hamiltonian H
(represented as an MPO) variationally. Psi should contain the initial
guess on input.
**Arguments**
Psi : TYPE(mpsc), inout
On entry the initial guess, on exit the excited state.
Psibs : TYPE(mpsc), inout
Energetically lower lying states which are required to
orthogonalize against those states.
ne : INTEGER, in
Excited state to be found.
Hts : TYPE(sr_matrix_tensorc)(\*), inout
List of all nearest-neighbor two-site MPO matrices.
Ham : TYPE(mpoc), inout
Hamiltonian for the system, represented as MPO.
LR : TYPE(tensorlistc)(\*), inout
Left right overlaps with Hamiltonian / MPO.
LRAB : TYPE(tensorlistc)(\*), inout
Left right overlaps with energetically lower lying states.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the excited state
search.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_excited_two_mpsc_mpoc(Psi, Psibs, ne, Hts, &
Ham, Cp, energy, converged, conv_val, pbc, errst)
type(mpsc), intent(inout) :: Psi
type(mpsc), pointer, intent(inout) :: Psibs(:)
integer, intent(in) :: ne
type(sr_matrix_tensorc), dimension(:), intent(inout) :: Hts
type(mpoc), intent(inout) :: Ham
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! MPO squared
type(mpoc) :: Hsq
! Overlap state-Hamiltonian and overlap state-state
type(tensorlistc), dimension(:), allocatable :: LR, LRAB
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_excited_two_mpsc_'//&
! 'mpoc: variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:134', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site excited state.'
call setuplr(LR, Psi, Ham, Psi, Psi%oc)
allocate(LRAB(ne))
do ii = 1, ne
call setuplr(LRAB(ii), Psi, Psibs(ii), Psi%oc)
end do
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print - (energy: there was no iMPS growth before)
energy = 0.0_rKind
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'sweep', jj, &
'bd', chi, 'energy', energy
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call find_excited_two_step_mpsc_mpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mpsc_'//&
! 'mpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:194', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call find_excited_two_step_mpsc_mpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mpsc_'//&
! 'mpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:205', &
! errst=errst)) return
end do
Psi%oc = 2
sense = 1
do kk = 1, (k1 - 1)
call find_excited_two_step_mpsc_mpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_mpsc_'//&
! 'mpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:217', &
! errst=errst)) return
end do
Psi%oc = k1
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'bond dimension, variance'//&
', last variance, difference in variance, discarded '//&
'weight', chi, conv_val, varianceSave, &
abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
if((verbose > 0) .and. (.not. converged)) write(slog, *) &
'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
if((verbose > 0) .and. converged) write(slog, *) &
'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
do ii = 1, size(LRAB, 1)
do jj = 1, size(LRAB(ii)%Li)
call destroy(LRAB(ii)%Li(jj))
end do
deallocate(LRAB(ii)%Li)
end do
deallocate(LR, LRAB)
if(Cp%conv_method == 'S') deallocate(firstlambda)
end subroutine find_excited_two_mpsc_mpoc
"""
return
[docs]def find_excited_two_qmps_qmpo():
"""
fortran-subroutine - ?? (mlw)
Find the excited state Psi (represented as an MPS) of the Hamiltonian H
(represented as an MPO) variationally. Psi should contain the initial
guess on input.
**Arguments**
Psi : TYPE(qmps), inout
On entry the initial guess, on exit the excited state.
Psibs : TYPE(qmps), inout
Energetically lower lying states which are required to
orthogonalize against those states.
ne : INTEGER, in
Excited state to be found.
Hts : TYPE(sr_matrix_qtensor)(\*), inout
List of all nearest-neighbor two-site MPO matrices.
Ham : TYPE(qmpo), inout
Hamiltonian for the system, represented as MPO.
LR : TYPE(qtensorlist)(\*), inout
Left right overlaps with Hamiltonian / MPO.
LRAB : TYPE(qtensorlist)(\*), inout
Left right overlaps with energetically lower lying states.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the excited state
search.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_excited_two_qmps_qmpo(Psi, Psibs, ne, Hts, &
Ham, Cp, energy, converged, conv_val, pbc, errst)
type(qmps), intent(inout) :: Psi
type(qmps), pointer, intent(inout) :: Psibs(:)
integer, intent(in) :: ne
type(sr_matrix_qtensor), dimension(:), intent(inout) :: Hts
type(qmpo), intent(inout) :: Ham
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! MPO squared
type(qmpo) :: Hsq
! Overlap state-Hamiltonian and overlap state-state
type(qtensorlist), dimension(:), allocatable :: LR, LRAB
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_excited_two_qmps_'//&
! 'qmpo: variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:134', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site excited state.'
call setuplr(LR, Psi, Ham, Psi, Psi%oc)
allocate(LRAB(ne))
do ii = 1, ne
call setuplr(LRAB(ii), Psi, Psibs(ii), Psi%oc)
end do
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print - (energy: there was no iMPS growth before)
energy = 0.0_rKind
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'sweep', jj, &
'bd', chi, 'energy', energy
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call find_excited_two_step_qmps_qmpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmps_'//&
! 'qmpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:194', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call find_excited_two_step_qmps_qmpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmps_'//&
! 'qmpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:205', &
! errst=errst)) return
end do
Psi%oc = 2
sense = 1
do kk = 1, (k1 - 1)
call find_excited_two_step_qmps_qmpo(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmps_'//&
! 'qmpo: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:217', &
! errst=errst)) return
end do
Psi%oc = k1
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'bond dimension, variance'//&
', last variance, difference in variance, discarded '//&
'weight', chi, conv_val, varianceSave, &
abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
if((verbose > 0) .and. (.not. converged)) write(slog, *) &
'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
if((verbose > 0) .and. converged) write(slog, *) &
'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
do ii = 1, size(LRAB, 1)
do jj = 1, size(LRAB(ii)%Li)
call destroy(LRAB(ii)%Li(jj))
end do
deallocate(LRAB(ii)%Li)
end do
deallocate(LR, LRAB)
if(Cp%conv_method == 'S') deallocate(firstlambda)
end subroutine find_excited_two_qmps_qmpo
"""
return
[docs]def find_excited_two_qmpsc_qmpoc():
"""
fortran-subroutine - ?? (mlw)
Find the excited state Psi (represented as an MPS) of the Hamiltonian H
(represented as an MPO) variationally. Psi should contain the initial
guess on input.
**Arguments**
Psi : TYPE(qmpsc), inout
On entry the initial guess, on exit the excited state.
Psibs : TYPE(qmpsc), inout
Energetically lower lying states which are required to
orthogonalize against those states.
ne : INTEGER, in
Excited state to be found.
Hts : TYPE(sr_matrix_qtensorc)(\*), inout
List of all nearest-neighbor two-site MPO matrices.
Ham : TYPE(qmpoc), inout
Hamiltonian for the system, represented as MPO.
LR : TYPE(qtensorclist)(\*), inout
Left right overlaps with Hamiltonian / MPO.
LRAB : TYPE(qtensorclist)(\*), inout
Left right overlaps with energetically lower lying states.
Cp : TYPE(ConvParam), in
Contains the convergence parameters for the excited state
search.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_excited_two_qmpsc_qmpoc(Psi, Psibs, ne, Hts, &
Ham, Cp, energy, converged, conv_val, pbc, errst)
type(qmpsc), intent(inout) :: Psi
type(qmpsc), pointer, intent(inout) :: Psibs(:)
integer, intent(in) :: ne
type(sr_matrix_qtensorc), dimension(:), intent(inout) :: Hts
type(qmpoc), intent(inout) :: Ham
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! MPO squared
type(qmpoc) :: Hsq
! Overlap state-Hamiltonian and overlap state-state
type(qtensorclist), dimension(:), allocatable :: LR, LRAB
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_excited_two_qmpsc_'//&
! 'qmpoc: variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:134', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site excited state.'
call setuplr(LR, Psi, Ham, Psi, Psi%oc)
allocate(LRAB(ne))
do ii = 1, ne
call setuplr(LRAB(ii), Psi, Psibs(ii), Psi%oc)
end do
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print - (energy: there was no iMPS growth before)
energy = 0.0_rKind
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'sweep', jj, &
'bd', chi, 'energy', energy
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call find_excited_two_step_qmpsc_qmpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmpsc_'//&
! 'qmpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:194', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call find_excited_two_step_qmpsc_qmpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmpsc_'//&
! 'qmpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:205', &
! errst=errst)) return
end do
Psi%oc = 2
sense = 1
do kk = 1, (k1 - 1)
call find_excited_two_step_qmpsc_qmpoc(Psi, Psibs, &
ne, Hts, Ham, LR, LRAB, kk, sense, Cp, local_tol, &
energy, innererr, errst=errst)
!if(prop_error('find_excited_two_qmpsc_'//&
! 'qmpoc: find_excited_two_step '//&
! 'failed.', 'VariationalOps_include.f90:217', &
! errst=errst)) return
end do
Psi%oc = k1
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'bond dimension, variance'//&
', last variance, difference in variance, discarded '//&
'weight', chi, conv_val, varianceSave, &
abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
if((verbose > 0) .and. (.not. converged)) write(slog, *) &
'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
if((verbose > 0) .and. converged) write(slog, *) &
'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance (or sing val diff)=', conv_val
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
do ii = 1, size(LRAB, 1)
do jj = 1, size(LRAB(ii)%Li)
call destroy(LRAB(ii)%Li(jj))
end do
deallocate(LRAB(ii)%Li)
end do
deallocate(LR, LRAB)
if(Cp%conv_method == 'S') deallocate(firstlambda)
end subroutine find_excited_two_qmpsc_qmpoc
"""
return
[docs]def find_groundstate_two_mpo():
"""
fortran-subroutine - July 2017 (dj, updated)
Find the ground state psi (represented as an MPS) of an Hamiltonian H.
**Arguments**
Psi : TYPE(mps), inout
Find the ground state based on the initial guess Psi. On exit
the optimized wave function.
Hts : TYPE(sr_matrix_tensor)(\*), inout
All nearest neighbor two-site MPO matrices.
Ham : TYPE(mpo), inout
The Hamiltonian represented as MPO.
LR : TYPE(tensorlist), inout
The left-right overlap of the wave-function Psi with the
Hamiltonian H.
max_chi : INTEGER, in
Maximal bond dimension.
max_lanczos_iter : INTEGER, in
Maximal number of Lanczos vectors generated.
lanczos_tol : REAL, in
Tolerance for the Lanczos algorithm as soft cut-off before
reaching max_lanczos_iter.
local_tol : REAL, in
Defining the soft cut-off for the singular values.
conv_tol : REAL, in
The tolerance criterion, either for the variance or the
maximum relative difference for the first singular values.
conv_method : CHARACTER, in
'V' for variance and 'S' for relative difference between
the first singular values.
energy : REAL, out
Energy of Psi.
converged : LOGICAL, out
Flag for convergence of the ground state search.
conv_val : REAL, out
The actual value of the convergence criterion achieved,
either the variance or the relative difference for the
first singular values.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_groundstate_two_mpo(Psi, Hts, Ham, LR, &
Cp, energy, converged, conv_val, pbc, errst)
type(mps), intent(inout) :: Psi
type(sr_matrix_tensor), dimension(:), intent(inout) :: Hts
type(mpo), intent(inout) :: Ham
type(tensorlist), dimension(:), allocatable :: LR
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! Square of the Hamiltonian for variance
type(mpo) :: Hsq
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_groundstate_two_mpo:'//&
! ' variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:568', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site algorithm'
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
firstlambda = 0.0_rKind
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'inner', jj, &
'bond dimension chi', chi
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call step_groundstate_two_mps_mpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpo :'//&
! 'step_groundstate_two (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call step_groundstate_two_mps_mpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpo :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call step_groundstate_two_mps_mpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpo :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'energy bond dimension'//&
', variance, last variance, difference in variance'//&
', discarded weight', energy, chi, conv_val, &
varianceSave, abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
! Deallocate stuff
if(Cp%conv_method == 'S') then
deallocate(firstlambda)
end if
! Print results
if((verbose > 0) .and. (converged)) then
write(slog, *) 'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
elseif(verbose > 0) then
write(slog, *) 'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
end if
end subroutine find_groundstate_two_mpo
"""
return
[docs]def find_groundstate_two_mpoc():
"""
fortran-subroutine - July 2017 (dj, updated)
Find the ground state psi (represented as an MPS) of an Hamiltonian H.
**Arguments**
Psi : TYPE(mpsc), inout
Find the ground state based on the initial guess Psi. On exit
the optimized wave function.
Hts : TYPE(sr_matrix_tensorc)(\*), inout
All nearest neighbor two-site MPO matrices.
Ham : TYPE(mpoc), inout
The Hamiltonian represented as MPO.
LR : TYPE(tensorlistc), inout
The left-right overlap of the wave-function Psi with the
Hamiltonian H.
max_chi : INTEGER, in
Maximal bond dimension.
max_lanczos_iter : INTEGER, in
Maximal number of Lanczos vectors generated.
lanczos_tol : REAL, in
Tolerance for the Lanczos algorithm as soft cut-off before
reaching max_lanczos_iter.
local_tol : REAL, in
Defining the soft cut-off for the singular values.
conv_tol : REAL, in
The tolerance criterion, either for the variance or the
maximum relative difference for the first singular values.
conv_method : CHARACTER, in
'V' for variance and 'S' for relative difference between
the first singular values.
energy : REAL, out
Energy of Psi.
converged : LOGICAL, out
Flag for convergence of the ground state search.
conv_val : REAL, out
The actual value of the convergence criterion achieved,
either the variance or the relative difference for the
first singular values.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_groundstate_two_mpoc(Psi, Hts, Ham, LR, &
Cp, energy, converged, conv_val, pbc, errst)
type(mpsc), intent(inout) :: Psi
type(sr_matrix_tensorc), dimension(:), intent(inout) :: Hts
type(mpoc), intent(inout) :: Ham
type(tensorlistc), dimension(:), allocatable :: LR
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! Square of the Hamiltonian for variance
type(mpoc) :: Hsq
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_groundstate_two_mpoc:'//&
! ' variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:568', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site algorithm'
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
firstlambda = 0.0_rKind
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'inner', jj, &
'bond dimension chi', chi
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call step_groundstate_two_mpsc_mpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpoc :'//&
! 'step_groundstate_two (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call step_groundstate_two_mpsc_mpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpoc :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call step_groundstate_two_mpsc_mpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_mpoc :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'energy bond dimension'//&
', variance, last variance, difference in variance'//&
', discarded weight', energy, chi, conv_val, &
varianceSave, abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
! Deallocate stuff
if(Cp%conv_method == 'S') then
deallocate(firstlambda)
end if
! Print results
if((verbose > 0) .and. (converged)) then
write(slog, *) 'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
elseif(verbose > 0) then
write(slog, *) 'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
end if
end subroutine find_groundstate_two_mpoc
"""
return
[docs]def find_groundstate_two_qmpo():
"""
fortran-subroutine - July 2017 (dj, updated)
Find the ground state psi (represented as an MPS) of an Hamiltonian H.
**Arguments**
Psi : TYPE(qmps), inout
Find the ground state based on the initial guess Psi. On exit
the optimized wave function.
Hts : TYPE(sr_matrix_qtensor)(\*), inout
All nearest neighbor two-site MPO matrices.
Ham : TYPE(qmpo), inout
The Hamiltonian represented as MPO.
LR : TYPE(qtensorlist), inout
The left-right overlap of the wave-function Psi with the
Hamiltonian H.
max_chi : INTEGER, in
Maximal bond dimension.
max_lanczos_iter : INTEGER, in
Maximal number of Lanczos vectors generated.
lanczos_tol : REAL, in
Tolerance for the Lanczos algorithm as soft cut-off before
reaching max_lanczos_iter.
local_tol : REAL, in
Defining the soft cut-off for the singular values.
conv_tol : REAL, in
The tolerance criterion, either for the variance or the
maximum relative difference for the first singular values.
conv_method : CHARACTER, in
'V' for variance and 'S' for relative difference between
the first singular values.
energy : REAL, out
Energy of Psi.
converged : LOGICAL, out
Flag for convergence of the ground state search.
conv_val : REAL, out
The actual value of the convergence criterion achieved,
either the variance or the relative difference for the
first singular values.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_groundstate_two_qmpo(Psi, Hts, Ham, LR, &
Cp, energy, converged, conv_val, pbc, errst)
type(qmps), intent(inout) :: Psi
type(sr_matrix_qtensor), dimension(:), intent(inout) :: Hts
type(qmpo), intent(inout) :: Ham
type(qtensorlist), dimension(:), allocatable :: LR
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! Square of the Hamiltonian for variance
type(qmpo) :: Hsq
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_groundstate_two_qmpo:'//&
! ' variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:568', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site algorithm'
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
firstlambda = 0.0_rKind
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'inner', jj, &
'bond dimension chi', chi
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call step_groundstate_two_qmps_qmpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpo :'//&
! 'step_groundstate_two (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call step_groundstate_two_qmps_qmpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpo :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call step_groundstate_two_qmps_qmpo(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpo :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'energy bond dimension'//&
', variance, last variance, difference in variance'//&
', discarded weight', energy, chi, conv_val, &
varianceSave, abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
! Deallocate stuff
if(Cp%conv_method == 'S') then
deallocate(firstlambda)
end if
! Print results
if((verbose > 0) .and. (converged)) then
write(slog, *) 'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
elseif(verbose > 0) then
write(slog, *) 'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
end if
end subroutine find_groundstate_two_qmpo
"""
return
[docs]def find_groundstate_two_qmpoc():
"""
fortran-subroutine - July 2017 (dj, updated)
Find the ground state psi (represented as an MPS) of an Hamiltonian H.
**Arguments**
Psi : TYPE(qmpsc), inout
Find the ground state based on the initial guess Psi. On exit
the optimized wave function.
Hts : TYPE(sr_matrix_qtensorc)(\*), inout
All nearest neighbor two-site MPO matrices.
Ham : TYPE(qmpoc), inout
The Hamiltonian represented as MPO.
LR : TYPE(qtensorclist), inout
The left-right overlap of the wave-function Psi with the
Hamiltonian H.
max_chi : INTEGER, in
Maximal bond dimension.
max_lanczos_iter : INTEGER, in
Maximal number of Lanczos vectors generated.
lanczos_tol : REAL, in
Tolerance for the Lanczos algorithm as soft cut-off before
reaching max_lanczos_iter.
local_tol : REAL, in
Defining the soft cut-off for the singular values.
conv_tol : REAL, in
The tolerance criterion, either for the variance or the
maximum relative difference for the first singular values.
conv_method : CHARACTER, in
'V' for variance and 'S' for relative difference between
the first singular values.
energy : REAL, out
Energy of Psi.
converged : LOGICAL, out
Flag for convergence of the ground state search.
conv_val : REAL, out
The actual value of the convergence criterion achieved,
either the variance or the relative difference for the
first singular values.
pbc : LOGICAL, in
If PBC are used in any rule set. There is a check before,
but the debugging mode should ensure for now that there
are no calls with such a check.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine find_groundstate_two_qmpoc(Psi, Hts, Ham, LR, &
Cp, energy, converged, conv_val, pbc, errst)
type(qmpsc), intent(inout) :: Psi
type(sr_matrix_qtensorc), dimension(:), intent(inout) :: Hts
type(qmpoc), intent(inout) :: Ham
type(qtensorclist), dimension(:), allocatable :: LR
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: energy
logical, intent(out) :: converged
real(KIND=rKind), intent(out) :: conv_val
logical, intent(in) :: pbc
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! saving the cumulated error of splits (SVD) in an inner loop
real(KIND=rKind) :: innererr
! for convergence with first singular value
real(KIND=rKind) :: lambdakk
real(KIND=rKind), dimension(:), allocatable :: firstlambda
! for convergence with variance
real(KIND=rKind) :: variancesave, desratio
logical :: chifail, stuck
! Square of the Hamiltonian for variance
type(qmpoc) :: Hsq
!if(present(errst)) errst = 0
!if(pbc) then
! errst = raise_error('find_groundstate_two_qmpoc:'//&
! ' variational does not have PBC.', &
! 99, 'VariationalOps_include.f90:568', errst=errst)
! return
!end if
if(verbose > 0) write(slog, *) 'Beginning two-site algorithm'
if(Cp%local_tol < 0.0_rKind) then
local_tol = Cp%conv_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%local_tol
end if
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
if(Cp%conv_method == 'V') then
! Use variance for convergence check
variancesave = 10000.0_rKind
elseif(Cp%conv_method == 'S') then
! Use difference singular values for convergence check
allocate(firstlambda(Psi%ll + 1))
firstlambda = 0.0_rKind
do kk = 1, (Psi%ll + 1)
if(Psi%haslambda(kk)) then
firstlambda(kk) = maxvalue(Psi%Lambda(kk))
end if
end do
else
stop 'Bad convergence method.'
end if
! Set for first print
chi = maxchi(Psi)
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0
if(verbose > 0) write(slog, *) 'outer', ii, 'inner', jj, &
'bond dimension chi', chi
sense = 1
do kk = k1, max(Psi%ll - 2, 1)
call step_groundstate_two_qmpsc_qmpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpoc :'//&
! 'step_groundstate_two (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call step_groundstate_two_qmpsc_qmpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpoc :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call step_groundstate_two_qmpsc_qmpoc(Psi, Hts, Ham, &
LR, kk, sense, Cp, local_tol, energy, innererr, &
errst=errst)
!if(prop_error('find_groundstate_two_qmpoc :'//&
! 'step_groundstate_two (2) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
chi = maxchi(Psi)
if(Cp%conv_method == 'V') then
! Use variance
call shiftandsquare(Hsq, Ham, energy)
call meas_mpo(conv_val, Hsq, Psi)
call destroy(Hsq)
conv_val = abs(conv_val)
converged = (conv_val < Cp%conv_tol)
if(verbose > 0) write(slog, *) 'energy bond dimension'//&
', variance, last variance, difference in variance'//&
', discarded weight', energy, chi, conv_val, &
varianceSave, abs(varianceSave - conv_val), innererr
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(variancesave - conv_val) < innererr)
variancesave = conv_val
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = conv_val
end if
elseif(Cp%conv_method == 'S') then
! Use difference singular values
conv_val = 0.0_rKind
! Leave away boundaries while looping - Lambdas should
! be all set
do kk = 2, Psi%ll
lambdakk = maxvalue(Psi%Lambda(kk))
conv_val = max(conv_val, &
abs(firstlambda(kk) - lambdakk) &
/ abs(lambdakk))
firstlambda(kk) = lambdakk
end do
converged = (conv_val <= Cp%conv_tol)
end if
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%conv_tol / conv_val
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
! Deallocate stuff
if(Cp%conv_method == 'S') then
deallocate(firstlambda)
end if
! Print results
if((verbose > 0) .and. (converged)) then
write(slog, *) 'Converged in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
elseif(verbose > 0) then
write(slog, *) 'Failed to converge in', sweeps_used, &
'sweeps. Energy:', energy, &
' with variance/Lambda_ratio=', conv_val
end if
end subroutine find_groundstate_two_qmpoc
"""
return
[docs]def grow_system_mpo():
"""
fortran-subroutine - July 2017 (dj, updated)
Use the infinite size algorithm to grow the system up to a size L.
**Arguments**
Psi : TYPE(mps), inout
On exit, the initial guess for the ground state built with
the iMPS algorithm.
ll : INTEGER, in
Number of sites in the MPS.
Ham : TYPE(mpo), in
This is the MPO representing the Hamiltonian where the ground state
should be found.
LR : TYPE(tensorlist), inout
One exit, the left-right overlap of the state and the MPO.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
The quantum numbers of both symmetries. Only accessed for
for simulations with quantum numbers.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file. Only accessed for simulations
with quantum numbers.
Cp : TYPE(convparam), in
Contains the convergence parameters for the algorithms.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine grow_system_mpo(Psi, ll, Ham, LR, nqs, qs, IOObj, Cp, &
hasexcited, errst)
type(mps), intent(inout) :: Psi
integer, intent(in) :: ll
type(mpo), intent(inout) :: Ham
type(tensorlist), dimension(:), allocatable, intent(inout) :: LR
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
type(convparam), intent(in) :: Cp
logical, intent(in), optional :: hasexcited
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Iterating over the splittings
integer :: nn
! flag for single-site variational step in the case of odd L
logical :: singlesite
! number of two-site variational steps
integer :: niter
! Local dimension of the Hilbert space
integer :: locdim
! variable for error
real(KIND=rKind) :: err
! Contains the eigevalue for each added two sites
type(tensor) :: Eigvals
! Two site tensor
type(tensor) :: Theta
! Two site MPO matrix
type(sr_matrix_tensor) :: Hts
!if(present(errst)) errst = 0
! Set the number of two-site variational steps and set single-site
! step for odd
singlesite = .not. (mod(ll, 2) == 0)
niter = floor(0.5_rKind * ll)
if(singlesite) then
call create(Eigvals, [niter + 1])
else
call create(Eigvals, [niter])
end if
! Allocate the MPS
Psi%ll = ll
allocate(Psi%Aa(ll), Psi%haslambda(ll + 1), Psi%lambda(ll + 1), Psi%can(ll))
Psi%haslambda = .false.
Psi%can = 'o'
! Allocate the transfer matrices
allocate(LR(ll))
! Initialization of the algorithm (nn = 0)
! ===============================
nn = 0
! Create first random guess
call state_first_mps(Theta, locdim, Ham%Ws(1), nqs, qs, IOObj)
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(1), Hts, LR(ll), Eigvals%elem(nn + 1), &
Theta, .true., .true., &
Cp%max_num_lanczos_iter, Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpo: Lanczos (1).', &
! 'VariationalOps_include.f90:967', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err)
if(ll - nn /= nn + 2) then
call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2), errst=errst)
!if(prop_error('grow_system_mpo: copy failed.', &
! 'VariationalOps_include.f90:977', errst=errst)) return
end if
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at n-th'//&
'iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! Iterations growing the system
! =============================
do nn = 1, (niter - 1)
call state_iter_mps(Theta, Psi%Aa(nn), Psi%Aa(ll - nn + 1), &
Psi%Lambda(nn + 1), Psi%Lambda(nn), &
locdim, IOObj%qstateunit, nn, errst=errst)
!if(prop_error('grow_system_mpo: state_iter_mps '// &
! 'failed.', 'VariationalOps_include.f90:999', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn), errst=errst)
!if(prop_error('grow_system_mpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1003', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn + 1), errst=errst)
!if(prop_error('grow_system_mpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1007', errst=errst)) return
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_mpo: ptm_right_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1012', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), &
errst=errst)
!if(prop_error('grow_system_mpo: ptm_left_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1019', errst=errst)) return
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(nn), Hts, LR(ll - nn + 1), Eigvals%elem(nn + 1), &
Theta, .false., .false., Cp%max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpo: Lanczos (2).', &
! 'VariationalOps_include.f90:1026', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err, errst=errst)
!if(prop_error('grow_system_mpo: split failed.', &
! 'VariationalOps_include.f90:1034', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn + 1), errst=errst)
!if(prop_error('grow_system_mpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1038', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn), errst=errst)
!if(prop_error('grow_system_mpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1042', errst=errst)) return
if(nn /= niter - 1) call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2))
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
end do
! Finalize algorithm depending on the case odd/even
! -------------------------------------------------
! That is the value carried over from the loop, but to be explicit
nn = niter
if(singlesite) then
! Copy of singular values was not included in the if case above
call copy(Psi%Lambda(ll - nn + 1), Psi%Lambda(nn + 1))
call state_single_mps(Psi%Aa(nn + 1), Psi%Aa(nn), &
Psi%Aa(ll - nn + 1), locdim, &
IOObj%qstateunit)
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_mpo: ptm_right_mpo.', &
! 'VariationalOps_include.f90:1075', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_mpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1081', errst=errst)) return
call Lanczos(LR(nn), Ham%Ws(nn + 1), LR(nn + 2), &
Eigvals%elem(nn + 1), Psi%Aa(nn + 1), &
.false., .false., Cp%Max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpo: Lanczos (3).', &
! 'VariationalOps_include.f90:1088', errst=errst)) return
Psi%can(nn + 1) = 'c'
Psi%oc = nn + 1
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! We need a dummy transfer matrix for this site
call copy(LR(nn + 1), LR(nn))
elseif(ll == 2) then
! Some things are not initialized as expected for a larger
! system
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.true., errst=errst)
!if(prop_error('grow_system_mpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1106', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
else
! We begin on the left of the two MPS, and so we only add
! the right overlap
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.false., Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_mpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1124', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
end if
call state_finalize_mps(IOObj%qstateunit, hasexcited)
call destroy(Eigvals)
end subroutine grow_system_mpo
"""
return
[docs]def grow_system_mpoc():
"""
fortran-subroutine - July 2017 (dj, updated)
Use the infinite size algorithm to grow the system up to a size L.
**Arguments**
Psi : TYPE(mpsc), inout
On exit, the initial guess for the ground state built with
the iMPS algorithm.
ll : INTEGER, in
Number of sites in the MPS.
Ham : TYPE(mpoc), in
This is the MPO representing the Hamiltonian where the ground state
should be found.
LR : TYPE(tensorlistc), inout
One exit, the left-right overlap of the state and the MPO.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
The quantum numbers of both symmetries. Only accessed for
for simulations with quantum numbers.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file. Only accessed for simulations
with quantum numbers.
Cp : TYPE(convparam), in
Contains the convergence parameters for the algorithms.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine grow_system_mpoc(Psi, ll, Ham, LR, nqs, qs, IOObj, Cp, &
hasexcited, errst)
type(mpsc), intent(inout) :: Psi
integer, intent(in) :: ll
type(mpoc), intent(inout) :: Ham
type(tensorlistc), dimension(:), allocatable, intent(inout) :: LR
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
type(convparam), intent(in) :: Cp
logical, intent(in), optional :: hasexcited
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Iterating over the splittings
integer :: nn
! flag for single-site variational step in the case of odd L
logical :: singlesite
! number of two-site variational steps
integer :: niter
! Local dimension of the Hilbert space
integer :: locdim
! variable for error
real(KIND=rKind) :: err
! Contains the eigevalue for each added two sites
type(tensor) :: Eigvals
! Two site tensor
type(tensorc) :: Theta
! Two site MPO matrix
type(sr_matrix_tensorc) :: Hts
!if(present(errst)) errst = 0
! Set the number of two-site variational steps and set single-site
! step for odd
singlesite = .not. (mod(ll, 2) == 0)
niter = floor(0.5_rKind * ll)
if(singlesite) then
call create(Eigvals, [niter + 1])
else
call create(Eigvals, [niter])
end if
! Allocate the MPS
Psi%ll = ll
allocate(Psi%Aa(ll), Psi%haslambda(ll + 1), Psi%lambda(ll + 1), Psi%can(ll))
Psi%haslambda = .false.
Psi%can = 'o'
! Allocate the transfer matrices
allocate(LR(ll))
! Initialization of the algorithm (nn = 0)
! ===============================
nn = 0
! Create first random guess
call state_first_mpsc(Theta, locdim, Ham%Ws(1), nqs, qs, IOObj)
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(1), Hts, LR(ll), Eigvals%elem(nn + 1), &
Theta, .true., .true., &
Cp%max_num_lanczos_iter, Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpoc: Lanczos (1).', &
! 'VariationalOps_include.f90:967', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err)
if(ll - nn /= nn + 2) then
call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2), errst=errst)
!if(prop_error('grow_system_mpoc: copy failed.', &
! 'VariationalOps_include.f90:977', errst=errst)) return
end if
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at n-th'//&
'iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! Iterations growing the system
! =============================
do nn = 1, (niter - 1)
call state_iter_mpsc(Theta, Psi%Aa(nn), Psi%Aa(ll - nn + 1), &
Psi%Lambda(nn + 1), Psi%Lambda(nn), &
locdim, IOObj%qstateunit, nn, errst=errst)
!if(prop_error('grow_system_mpoc: state_iter_mpsc '// &
! 'failed.', 'VariationalOps_include.f90:999', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn), errst=errst)
!if(prop_error('grow_system_mpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1003', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn + 1), errst=errst)
!if(prop_error('grow_system_mpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1007', errst=errst)) return
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_mpoc: ptm_right_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1012', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), &
errst=errst)
!if(prop_error('grow_system_mpoc: ptm_left_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1019', errst=errst)) return
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(nn), Hts, LR(ll - nn + 1), Eigvals%elem(nn + 1), &
Theta, .false., .false., Cp%max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpoc: Lanczos (2).', &
! 'VariationalOps_include.f90:1026', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err, errst=errst)
!if(prop_error('grow_system_mpoc: split failed.', &
! 'VariationalOps_include.f90:1034', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn + 1), errst=errst)
!if(prop_error('grow_system_mpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1038', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn), errst=errst)
!if(prop_error('grow_system_mpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1042', errst=errst)) return
if(nn /= niter - 1) call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2))
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
end do
! Finalize algorithm depending on the case odd/even
! -------------------------------------------------
! That is the value carried over from the loop, but to be explicit
nn = niter
if(singlesite) then
! Copy of singular values was not included in the if case above
call copy(Psi%Lambda(ll - nn + 1), Psi%Lambda(nn + 1))
call state_single_mpsc(Psi%Aa(nn + 1), Psi%Aa(nn), &
Psi%Aa(ll - nn + 1), locdim, &
IOObj%qstateunit)
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_mpoc: ptm_right_mpo.', &
! 'VariationalOps_include.f90:1075', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_mpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1081', errst=errst)) return
call Lanczos(LR(nn), Ham%Ws(nn + 1), LR(nn + 2), &
Eigvals%elem(nn + 1), Psi%Aa(nn + 1), &
.false., .false., Cp%Max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_mpoc: Lanczos (3).', &
! 'VariationalOps_include.f90:1088', errst=errst)) return
Psi%can(nn + 1) = 'c'
Psi%oc = nn + 1
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! We need a dummy transfer matrix for this site
call copy(LR(nn + 1), LR(nn))
elseif(ll == 2) then
! Some things are not initialized as expected for a larger
! system
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.true., errst=errst)
!if(prop_error('grow_system_mpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1106', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
else
! We begin on the left of the two MPS, and so we only add
! the right overlap
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.false., Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_mpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1124', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
end if
call state_finalize_mpsc(IOObj%qstateunit, hasexcited)
call destroy(Eigvals)
end subroutine grow_system_mpoc
"""
return
[docs]def grow_system_qmpo():
"""
fortran-subroutine - July 2017 (dj, updated)
Use the infinite size algorithm to grow the system up to a size L.
**Arguments**
Psi : TYPE(qmps), inout
On exit, the initial guess for the ground state built with
the iMPS algorithm.
ll : INTEGER, in
Number of sites in the MPS.
Ham : TYPE(qmpo), in
This is the MPO representing the Hamiltonian where the ground state
should be found.
LR : TYPE(qtensorlist), inout
One exit, the left-right overlap of the state and the MPO.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
The quantum numbers of both symmetries. Only accessed for
for simulations with quantum numbers.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file. Only accessed for simulations
with quantum numbers.
Cp : TYPE(convparam), in
Contains the convergence parameters for the algorithms.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine grow_system_qmpo(Psi, ll, Ham, LR, nqs, qs, IOObj, Cp, &
hasexcited, errst)
type(qmps), intent(inout) :: Psi
integer, intent(in) :: ll
type(qmpo), intent(inout) :: Ham
type(qtensorlist), dimension(:), allocatable, intent(inout) :: LR
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
type(convparam), intent(in) :: Cp
logical, intent(in), optional :: hasexcited
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Iterating over the splittings
integer :: nn
! flag for single-site variational step in the case of odd L
logical :: singlesite
! number of two-site variational steps
integer :: niter
! Local dimension of the Hilbert space
integer :: locdim
! variable for error
real(KIND=rKind) :: err
! Contains the eigevalue for each added two sites
type(tensor) :: Eigvals
! Two site tensor
type(qtensor) :: Theta
! Two site MPO matrix
type(sr_matrix_qtensor) :: Hts
!if(present(errst)) errst = 0
! Set the number of two-site variational steps and set single-site
! step for odd
singlesite = .not. (mod(ll, 2) == 0)
niter = floor(0.5_rKind * ll)
if(singlesite) then
call create(Eigvals, [niter + 1])
else
call create(Eigvals, [niter])
end if
! Allocate the MPS
Psi%ll = ll
allocate(Psi%Aa(ll), Psi%haslambda(ll + 1), Psi%lambda(ll + 1), Psi%can(ll))
Psi%haslambda = .false.
Psi%can = 'o'
! Allocate the transfer matrices
allocate(LR(ll))
! Initialization of the algorithm (nn = 0)
! ===============================
nn = 0
! Create first random guess
call state_first_qmps(Theta, locdim, Ham%Ws(1), nqs, qs, IOObj)
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(1), Hts, LR(ll), Eigvals%elem(nn + 1), &
Theta, .true., .true., &
Cp%max_num_lanczos_iter, Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpo: Lanczos (1).', &
! 'VariationalOps_include.f90:967', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err)
if(ll - nn /= nn + 2) then
call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2), errst=errst)
!if(prop_error('grow_system_qmpo: copy failed.', &
! 'VariationalOps_include.f90:977', errst=errst)) return
end if
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at n-th'//&
'iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! Iterations growing the system
! =============================
do nn = 1, (niter - 1)
call state_iter_qmps(Theta, Psi%Aa(nn), Psi%Aa(ll - nn + 1), &
Psi%Lambda(nn + 1), Psi%Lambda(nn), &
locdim, IOObj%qstateunit, nn, errst=errst)
!if(prop_error('grow_system_qmpo: state_iter_qmps '// &
! 'failed.', 'VariationalOps_include.f90:999', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn), errst=errst)
!if(prop_error('grow_system_qmpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1003', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn + 1), errst=errst)
!if(prop_error('grow_system_qmpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1007', errst=errst)) return
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_qmpo: ptm_right_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1012', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), &
errst=errst)
!if(prop_error('grow_system_qmpo: ptm_left_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1019', errst=errst)) return
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(nn), Hts, LR(ll - nn + 1), Eigvals%elem(nn + 1), &
Theta, .false., .false., Cp%max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpo: Lanczos (2).', &
! 'VariationalOps_include.f90:1026', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err, errst=errst)
!if(prop_error('grow_system_qmpo: split failed.', &
! 'VariationalOps_include.f90:1034', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn + 1), errst=errst)
!if(prop_error('grow_system_qmpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1038', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn), errst=errst)
!if(prop_error('grow_system_qmpo: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1042', errst=errst)) return
if(nn /= niter - 1) call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2))
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
end do
! Finalize algorithm depending on the case odd/even
! -------------------------------------------------
! That is the value carried over from the loop, but to be explicit
nn = niter
if(singlesite) then
! Copy of singular values was not included in the if case above
call copy(Psi%Lambda(ll - nn + 1), Psi%Lambda(nn + 1))
call state_single_qmps(Psi%Aa(nn + 1), Psi%Aa(nn), &
Psi%Aa(ll - nn + 1), locdim, &
IOObj%qstateunit)
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_qmpo: ptm_right_mpo.', &
! 'VariationalOps_include.f90:1075', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_qmpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1081', errst=errst)) return
call Lanczos(LR(nn), Ham%Ws(nn + 1), LR(nn + 2), &
Eigvals%elem(nn + 1), Psi%Aa(nn + 1), &
.false., .false., Cp%Max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpo: Lanczos (3).', &
! 'VariationalOps_include.f90:1088', errst=errst)) return
Psi%can(nn + 1) = 'c'
Psi%oc = nn + 1
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! We need a dummy transfer matrix for this site
call copy(LR(nn + 1), LR(nn))
elseif(ll == 2) then
! Some things are not initialized as expected for a larger
! system
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.true., errst=errst)
!if(prop_error('grow_system_qmpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1106', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
else
! We begin on the left of the two MPS, and so we only add
! the right overlap
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.false., Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_qmpo: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1124', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
end if
call state_finalize_qmps(IOObj%qstateunit, hasexcited)
call destroy(Eigvals)
end subroutine grow_system_qmpo
"""
return
[docs]def grow_system_qmpoc():
"""
fortran-subroutine - July 2017 (dj, updated)
Use the infinite size algorithm to grow the system up to a size L.
**Arguments**
Psi : TYPE(qmpsc), inout
On exit, the initial guess for the ground state built with
the iMPS algorithm.
ll : INTEGER, in
Number of sites in the MPS.
Ham : TYPE(qmpoc), in
This is the MPO representing the Hamiltonian where the ground state
should be found.
LR : TYPE(qtensorclist), inout
One exit, the left-right overlap of the state and the MPO.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
The quantum numbers of both symmetries. Only accessed for
for simulations with quantum numbers.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file. Only accessed for simulations
with quantum numbers.
Cp : TYPE(convparam), in
Contains the convergence parameters for the algorithms.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine grow_system_qmpoc(Psi, ll, Ham, LR, nqs, qs, IOObj, Cp, &
hasexcited, errst)
type(qmpsc), intent(inout) :: Psi
integer, intent(in) :: ll
type(qmpoc), intent(inout) :: Ham
type(qtensorclist), dimension(:), allocatable, intent(inout) :: LR
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
type(convparam), intent(in) :: Cp
logical, intent(in), optional :: hasexcited
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Iterating over the splittings
integer :: nn
! flag for single-site variational step in the case of odd L
logical :: singlesite
! number of two-site variational steps
integer :: niter
! Local dimension of the Hilbert space
integer :: locdim
! variable for error
real(KIND=rKind) :: err
! Contains the eigevalue for each added two sites
type(tensor) :: Eigvals
! Two site tensor
type(qtensorc) :: Theta
! Two site MPO matrix
type(sr_matrix_qtensorc) :: Hts
!if(present(errst)) errst = 0
! Set the number of two-site variational steps and set single-site
! step for odd
singlesite = .not. (mod(ll, 2) == 0)
niter = floor(0.5_rKind * ll)
if(singlesite) then
call create(Eigvals, [niter + 1])
else
call create(Eigvals, [niter])
end if
! Allocate the MPS
Psi%ll = ll
allocate(Psi%Aa(ll), Psi%haslambda(ll + 1), Psi%lambda(ll + 1), Psi%can(ll))
Psi%haslambda = .false.
Psi%can = 'o'
! Allocate the transfer matrices
allocate(LR(ll))
! Initialization of the algorithm (nn = 0)
! ===============================
nn = 0
! Create first random guess
call state_first_qmpsc(Theta, locdim, Ham%Ws(1), nqs, qs, IOObj)
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(1), Hts, LR(ll), Eigvals%elem(nn + 1), &
Theta, .true., .true., &
Cp%max_num_lanczos_iter, Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpoc: Lanczos (1).', &
! 'VariationalOps_include.f90:967', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err)
if(ll - nn /= nn + 2) then
call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2), errst=errst)
!if(prop_error('grow_system_qmpoc: copy failed.', &
! 'VariationalOps_include.f90:977', errst=errst)) return
end if
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at n-th'//&
'iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! Iterations growing the system
! =============================
do nn = 1, (niter - 1)
call state_iter_qmpsc(Theta, Psi%Aa(nn), Psi%Aa(ll - nn + 1), &
Psi%Lambda(nn + 1), Psi%Lambda(nn), &
locdim, IOObj%qstateunit, nn, errst=errst)
!if(prop_error('grow_system_qmpoc: state_iter_qmpsc '// &
! 'failed.', 'VariationalOps_include.f90:999', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn), errst=errst)
!if(prop_error('grow_system_qmpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1003', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn + 1), errst=errst)
!if(prop_error('grow_system_qmpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1007', errst=errst)) return
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_right_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1012', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), &
errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_left_mpo (1) '// &
! 'failed.', 'VariationalOps_include.f90:1019', errst=errst)) return
call sdot(Hts, Ham%Ws(nn + 1), Ham%Ws(ll - nn))
call Lanczos(LR(nn), Hts, LR(ll - nn + 1), Eigvals%elem(nn + 1), &
Theta, .false., .false., Cp%max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpoc: Lanczos (2).', &
! 'VariationalOps_include.f90:1026', errst=errst)) return
call destroy(Hts)
call split(Psi%Aa(nn + 1), Psi%Lambda(nn + 2), Psi%Aa(ll - nn), &
Theta, [1, 2], [3, 4], trunc=Cp%warmup_tol, &
ncut=Cp%warmup_chi, renorm='M', err=err, errst=errst)
!if(prop_error('grow_system_qmpoc: split failed.', &
! 'VariationalOps_include.f90:1034', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(nn + 1), errst=errst)
!if(prop_error('grow_system_qmpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1038', errst=errst)) return
!call check_qnum_tensdim(Psi%Aa(ll - nn), errst=errst)
!if(prop_error('grow_system_qmpoc: check_qnum_tensdim'// &
! ' failed.', 'VariationalOps_include.f90:1042', errst=errst)) return
if(nn /= niter - 1) call copy(Psi%Lambda(ll - nn), Psi%Lambda(nn + 2))
Psi%haslambda(nn + 2) = .true.
Psi%haslambda(ll - nn) = .true.
Psi%can(nn + 1) = 'l'
Psi%can(ll - nn) = 'r'
call destroy(Theta)
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
end do
! Finalize algorithm depending on the case odd/even
! -------------------------------------------------
! That is the value carried over from the loop, but to be explicit
nn = niter
if(singlesite) then
! Copy of singular values was not included in the if case above
call copy(Psi%Lambda(ll - nn + 1), Psi%Lambda(nn + 1))
call state_single_qmpsc(Psi%Aa(nn + 1), Psi%Aa(nn), &
Psi%Aa(ll - nn + 1), locdim, &
IOObj%qstateunit)
call ptm_right_mpo(LR(nn), Psi%Aa(nn), Ham%Ws(nn), Psi%Aa(nn), &
(nn == 1), Matin=LR(max(1, nn - 1)), errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_right_mpo.', &
! 'VariationalOps_include.f90:1075', errst=errst)) return
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
(nn == 1), Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1081', errst=errst)) return
call Lanczos(LR(nn), Ham%Ws(nn + 1), LR(nn + 2), &
Eigvals%elem(nn + 1), Psi%Aa(nn + 1), &
.false., .false., Cp%Max_num_lanczos_iter, &
Cp%lanczos_tol, errst=errst)
!if(prop_error('grow_system_qmpoc: Lanczos (3).', &
! 'VariationalOps_include.f90:1088', errst=errst)) return
Psi%can(nn + 1) = 'c'
Psi%oc = nn + 1
if(verbose > 0) write(slog, *) 'Energy and energy density at'//&
' n-th iteration of IMPS:', nn, Eigvals%elem(nn + 1), &
Eigvals%elem(nn + 1) / (2 * (nn + 1))
! We need a dummy transfer matrix for this site
call copy(LR(nn + 1), LR(nn))
elseif(ll == 2) then
! Some things are not initialized as expected for a larger
! system
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.true., errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1106', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
else
! We begin on the left of the two MPS, and so we only add
! the right overlap
call ptm_left_mpo(LR(ll - nn + 1), Psi%Aa(ll - nn + 1), &
Ham%Ws(ll - nn + 1), Psi%Aa(ll - nn + 1), &
.false., Matin=LR(min(ll, ll - nn + 2)), errst=errst)
!if(prop_error('grow_system_qmpoc: ptm_left_mpo.', &
! 'VariationalOps_include.f90:1124', errst=errst)) return
! Contract singular values to install orthogonality center
call contr_diag(Psi%Aa(ll / 2), Psi%Lambda(ll / 2 + 1), 3)
Psi%can(ll / 2) = 'c'
Psi%oc = ll / 2
! We have to get any transfer matrix for the last site (will
! be deleted, but must be there)
call copy(LR(nn), LR(ll - nn + 1))
end if
call state_finalize_qmpsc(IOObj%qstateunit, hasexcited)
call destroy(Eigvals)
end subroutine grow_system_qmpoc
"""
return
[docs]def Fit_Hpsi_mps_mpo():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(mps), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(mpo), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(mps), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_mps_mpo(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(mps), intent(inout) :: Psi
type(mpo), intent(inout) :: Ham
type(mps), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(tensorlist), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
real(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(mpo) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_mps_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_mps_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_mps_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mps_mpo : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_mps_mpo
"""
return
[docs]def Fit_Hpsi_mpsc_mpo():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(mpo), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(mpsc), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_mpsc_mpo(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(mpsc), intent(inout) :: Psi
type(mpo), intent(inout) :: Ham
type(mpsc), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(tensorlistc), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
complex(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(mpo) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_mpsc_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_mpsc_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_mpsc_mpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpo : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_mpsc_mpo
"""
return
[docs]def Fit_Hpsi_mpsc_mpoc():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(mpsc), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(mpoc), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(mpsc), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_mpsc_mpoc(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(mpsc), intent(inout) :: Psi
type(mpoc), intent(inout) :: Ham
type(mpsc), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(tensorlistc), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
complex(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(mpoc) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_mpsc_mpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_mpsc_mpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_mpsc_mpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_mpsc_mpoc : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_mpsc_mpoc
"""
return
[docs]def Fit_Hpsi_qmps_qmpo():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(qmps), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(qmpo), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(qmps), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_qmps_qmpo(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(qmps), intent(inout) :: Psi
type(qmpo), intent(inout) :: Ham
type(qmps), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(qtensorlist), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
real(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(qmpo) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_qmps_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_qmps_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_qmps_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmps_qmpo : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_qmps_qmpo
"""
return
[docs]def Fit_Hpsi_qmpsc_qmpo():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(qmpo), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(qmpsc), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_qmpsc_qmpo(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(qmpsc), intent(inout) :: Psi
type(qmpo), intent(inout) :: Ham
type(qmpsc), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(qtensorclist), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
complex(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(qmpo) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_qmpsc_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_qmpsc_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_qmpsc_qmpo(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpo : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_qmpsc_qmpo
"""
return
[docs]def Fit_Hpsi_qmpsc_qmpoc():
"""
fortran-subroutine - April 2016 (updated, dj)
Fit p to :math:`| \psi \\rangle \\approx H \cdot | \phi \\rangle` using
a two-site variational algorithm.
**Arguments**
Psi : TYPE(qmpsc), inout
The new fitted MPS representing :math:`H | r \\rangle`.
Ham : TYPE(qmpoc), inout
MPO to be applied, e.g. propagator.
Phi : TYPE(qmpsc), inout
Apply the MPO to this MPS and try to target this state H p in fit.
Cp : TYPE(ConvParam), in
Convergence parameter for the fit.
err : REAL, out
Summing up the local errors done during the procedure.
converged : LOGICAL, out
Flag if fit converged according the convergence parameters.
initmethod : CHARACTER, in
'R' for randomization of whole MPS, 'F' for adding fluctuations, else
initial guess is a copy of p. (Currently methods 'R' and 'F' are equal).
convmethod : CHARACTER, in
'V' check on convergence by checking on the variance measuring the
MPO squared. 'E' uses the error for the estimate of convergence.
renorm : CHARACTER, in
Flag if the singular values should be renormalized when splitting
two sites. 'N' : no normalization; 'M' normalize for MPS.
**Details**
(template defined in VariationalOps_transfer.f90)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine Fit_Hpsi_qmpsc_qmpoc(Psi, Ham, Phi, Cp, err, converged, &
initmethod, convmethod, renorm, errst)
type(qmpsc), intent(inout) :: Psi
type(qmpoc), intent(inout) :: Ham
type(qmpsc), intent(inout) :: Phi
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
character, intent(in) :: initmethod, convmethod
character, intent(in) :: renorm
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! save orthogonalization for randomization
integer :: ocsave
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Error of one inner sweep
real(KIND=rKind) :: innererr
! The overalp between Psi, Ham, and Phi
type(qtensorclist), dimension(:), allocatable :: LR(:)
! for the convergence
real(KIND=rKind) :: variance, variancesave, desratio, truenorm, mynorm
logical :: chifail, stuck
! For overlap <Psi | H | Phi>
complex(KIND=rKind) :: overlap
real(KIND=rKind) :: roverlap
! square of Ham
type(qmpoc) :: Hsq
if(verbose > 0) write(slog, *) 'Beginning fitting Psi to H Phi'
if(Cp%hlocal_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%hpsi_tol / (4.0_rKind * Phi%ll)
else
local_tol = Cp%hlocal_tol
end if
! Initialize guess to be p or randomized MPS
call copy(Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : copy (1) failed.', &
! errst=errst)) return
if((initmethod == 'R') .or. (initmethod == 'F')) then
ocsave = Psi%oc
call randomize(Psi)
call orthonormalize(Psi, ocsave, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : orhtonormalize '//&
! 'failed.', errst=errst)) return
end if
! Setup the left-right overlaps for
call setuplr(LR, Psi, Ham, Phi, Phi%oc, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : setuplr '//&
! 'failed.', errst=errst)) return
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
varianceSave = 10000.0_rKind
err = 0.0_rKind
if(convmethod == 'V') then
! Get the norm of H|\psi>
call shiftandsquare(Hsq, Ham, 0.0_rKind)
call meas_mpo(truenorm, Hsq, Phi)
call destroy(Hsq)
end if
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
innererr = 0.0_rKind
sense = 1
do kk = k1, (Psi%ll - 2)
call fit_hpsi_step_qmpsc_qmpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : '//&
! 'step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fit_hpsi_step_qmpsc_qmpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : '//&
! 'step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call fit_hpsi_step_qmpsc_qmpoc(Psi, Ham, LR, Phi, &
kk, sense, Cp, local_tol, renorm, innererr, &
errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : '//&
! 'step (3) failed.', &
! errst=errst)) return
end do
err = err + innererr
! Convergence checks
! ..................
chi = maxchi(Psi)
mynorm = norm(Psi)
if(convmethod == 'V') then
call meas_psi_mpo_phi(overlap, Ham, Psi, Phi, errst=errst)
!if(prop_error('Fit_Hpsi_qmpsc_qmpoc : '//&
! 'meas_psi_mpo_phi failed.', &
! errst=errst)) return
roverlap = 2.0_rKind * real(overlap, KIND=rKind)
!tolerance= <\phi|\phi>+<\psi|H^2|\psi>-2Re(<\phi|H\psi>)
variance = (truenorm - mynorm) / truenorm
converged = (abs(variance) < Cp%hpsi_tol * Psi%ll)
elseif(convmethod == 'E') then
converged = (abs(innererr) < Cp%hpsi_tol)
end if
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
if(convmethod == 'V') then
if(abs(varianceSave - abs(variance)) < innererr) then
stuck = .true.
end if
elseif(convmethod == 'E') then
stuck = .true.
end if
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
variancesave = abs(variance)
end if
! Set stuck true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
if(convmethod == 'V') desratio = Cp%hpsi_tol / variance
if(convmethod == 'E') desratio = Cp%hpsi_tol / err
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do ii = 1, size(LR, 1)
do jj = 1, size(LR(ii)%Li, 1)
call destroy(LR(ii)%Li(jj))
end do
deallocate(LR(ii)%Li)
end do
deallocate(LR)
end subroutine Fit_Hpsi_qmpsc_qmpoc
"""
return
[docs]def FitMPS_mps_real():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(mps), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_mps_real(Psi, weight, Phis, Cp, err, converged, errst)
type(mps) :: Psi
real(KIND=rKind), intent(in) :: weight(:)
type(mps), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
real(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(tensorlist), pointer :: LR(:)
! Two site tensors
type(tensor) :: Theta, Theta2
! temporary tensor
type(tensor) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_mps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_mps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_mps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_mps_real
"""
return
[docs]def FitMPS_mpsc_real():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(mpsc), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_mpsc_real(Psi, weight, Phis, Cp, err, converged, errst)
type(mpsc) :: Psi
real(KIND=rKind), intent(in) :: weight(:)
type(mpsc), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(tensorlistc), pointer :: LR(:)
! Two site tensors
type(tensorc) :: Theta, Theta2
! temporary tensor
type(tensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_mpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_mpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_mpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_mpsc_real
"""
return
[docs]def FitMPS_mpsc_complex():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(mpsc), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_mpsc_complex(Psi, weight, Phis, Cp, err, converged, errst)
type(mpsc) :: Psi
complex(KIND=rKind), intent(in) :: weight(:)
type(mpsc), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(tensorlistc), pointer :: LR(:)
! Two site tensors
type(tensorc) :: Theta, Theta2
! temporary tensor
type(tensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_mpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_mpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_mpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_mpsc_complex
"""
return
[docs]def FitMPS_qmps_real():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(qmps), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_qmps_real(Psi, weight, Phis, Cp, err, converged, errst)
type(qmps) :: Psi
real(KIND=rKind), intent(in) :: weight(:)
type(qmps), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
real(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(qtensorlist), pointer :: LR(:)
! Two site tensors
type(qtensor) :: Theta, Theta2
! temporary tensor
type(qtensor) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_qmps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_qmps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_qmps_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_qmps_real
"""
return
[docs]def FitMPS_qmpsc_real():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(qmpsc), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_qmpsc_real(Psi, weight, Phis, Cp, err, converged, errst)
type(qmpsc) :: Psi
real(KIND=rKind), intent(in) :: weight(:)
type(qmpsc), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(qtensorclist), pointer :: LR(:)
! Two site tensors
type(qtensorc) :: Theta, Theta2
! temporary tensor
type(qtensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_qmpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_qmpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_qmpsc_real(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_qmpsc_real
"""
return
[docs]def FitMPS_qmpsc_complex():
"""
fortran-subroutine - July 2017 (dj, updated)
Fit Psi to \sum_k weight(k)*Phis(k)
**Arguments**
r : TYPE(qmpsc), inout
??
weight : REAL_OR_COMPLEX(*), in
on input some scaling factor
Phis : , inout
??
maxbondD : INTEGER, in
maximal bond dimension for splitting up T-sites.
err : REAL, inout
error done during subroutine (exact definition ??)
converged : LOGICAL, inout
Convergence flag (exact definition?)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine FitMPS_qmpsc_complex(Psi, weight, Phis, Cp, err, converged, errst)
type(qmpsc) :: Psi
complex(KIND=rKind), intent(in) :: weight(:)
type(qmpsc), pointer, intent(inout) :: Phis(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Number of states in Phis
integer :: nstates
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! The overlap between each state in Phis and the state Psi
type(qtensorclist), pointer :: LR(:)
! Two site tensors
type(qtensorc) :: Theta, Theta2
! temporary tensor
type(qtensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
!if(present(errst)) errst = 0
if(verbose > 0) write(slog, *) 'Beginning fitting of Psi to Phis'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Setup the left-right overlaps for all Phis
nstates = size(weight, 1)
allocate(LR(nstates), dotvec(nstates))
do pp = 1, nstates
call setuplr(LR(pp), Psi, Phis(pp), Psi%oc)
end do
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call fitMPS_step_qmpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call fitMPS_step_qmpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
sense = 1
do kk = 1, (k1 - 1)
call fitMPS_step_qmpsc_complex(Psi, weight, &
Phis, LR, kk, sense, Cp, local_tol, err)
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = 1, nstates
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1])
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
end if
call destroy(Tmpe)
dotvec(pp) = weight(pp) * dot(Theta, Theta2)
call destroy(Theta2)
end do
call destroy(Theta)
!alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
alpc = real(sum(dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(1.0_rKind - abs(alpc)) < Cp%psi_tol)
if(ii >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(alpsave - abs(alpc)) < Cp%psi_tol)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, nstates
do ii = 1, Psi%ll
call destroy(LR(pp)%Li(ii))
end do
deallocate(LR(pp)%Li)
end do
deallocate(LR, dotvec)
end subroutine FitMPS_qmpsc_complex
"""
return
[docs]def orth_mps():
"""
fortran-subroutine - Orthogonalize Psi against other MPS Phis.
**Arguments**
Psi : TYP(mps), inout
Orthogonalize this MPS against the MPSs Phis.
Phis : , inout
The MPS Psi will be orthogonal to these MPSs.
minn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
maxn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
LR : TYPE(tensorlist), inout
The partial overlaps at the orthogonality center.
err : REAL, out
Sum of all errors from splitting.
converged : LOGICAL, out
Flag for convergence of the variational orthogonalization.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine orth_mps(Psi, Phis, minn, maxn, LR, Cp, err, converged, &
errst)
type(mps), intent(inout) :: Psi
type(mps), pointer, intent(inout) :: Phis(:)
integer, intent(in) :: minn, maxn
type(tensorlist), pointer, intent(inout) :: LR(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
real(KIND=rKind), dimension(:), allocatable :: dotvec
! Copy of the original incoming Psi
type(mps) :: Psp
! Variational overlap between original and new state
type(tensorlist) :: Lrv
! Two site tensors
type(tensor) :: Theta, Theta2
! temporary tensor
type(tensor) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
if(verbose > 0) write(slog, *) 'Beginning two-site orthogonalization'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Psp is Psi and Psi becomes the variational state
call copy(Psp, Psi)
! variational overlap
call setuplr(Lrv, Psi, Psp, Psi%oc)
! Get workspace for dot product between Psi and Phis
allocate(dotvec(maxn))
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call orth_step_mps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mps : step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call orth_step_mps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mps : step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call orth_step_mps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mps : step (3) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = minn, maxn
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1], &
errst=errst)
!if(prop_error('orth_mps : contr (1) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
!if(prop_error('orth_mps : contr (2) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Tmpe)
dotvec(pp) = dot(Theta, Theta2, errst=errst)
!if(prop_error('orth_mps : dot failed', &
! errst=errst)) return
call destroy(Theta2)
end do
call destroy(Theta)
alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(alpc) < Cp%psi_tol)
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(abs(alpsave) - alpc) < err)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, size(LRV%Li)
call destroy(LRV%Li(pp))
end do
deallocate(dotvec, LRV%Li)
call destroy(Psp)
end subroutine orth_mps
"""
return
[docs]def orth_mpsc():
"""
fortran-subroutine - Orthogonalize Psi against other MPS Phis.
**Arguments**
Psi : TYP(mpsc), inout
Orthogonalize this MPS against the MPSs Phis.
Phis : , inout
The MPS Psi will be orthogonal to these MPSs.
minn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
maxn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
LR : TYPE(tensorlistc), inout
The partial overlaps at the orthogonality center.
err : REAL, out
Sum of all errors from splitting.
converged : LOGICAL, out
Flag for convergence of the variational orthogonalization.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine orth_mpsc(Psi, Phis, minn, maxn, LR, Cp, err, converged, &
errst)
type(mpsc), intent(inout) :: Psi
type(mpsc), pointer, intent(inout) :: Phis(:)
integer, intent(in) :: minn, maxn
type(tensorlistc), pointer, intent(inout) :: LR(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! Copy of the original incoming Psi
type(mpsc) :: Psp
! Variational overlap between original and new state
type(tensorlistc) :: Lrv
! Two site tensors
type(tensorc) :: Theta, Theta2
! temporary tensor
type(tensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
if(verbose > 0) write(slog, *) 'Beginning two-site orthogonalization'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Psp is Psi and Psi becomes the variational state
call copy(Psp, Psi)
! variational overlap
call setuplr(Lrv, Psi, Psp, Psi%oc)
! Get workspace for dot product between Psi and Phis
allocate(dotvec(maxn))
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call orth_step_mpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mpsc : step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call orth_step_mpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mpsc : step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call orth_step_mpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_mpsc : step (3) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = minn, maxn
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1], &
errst=errst)
!if(prop_error('orth_mpsc : contr (1) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
!if(prop_error('orth_mpsc : contr (2) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Tmpe)
dotvec(pp) = dot(Theta, Theta2, errst=errst)
!if(prop_error('orth_mpsc : dot failed', &
! errst=errst)) return
call destroy(Theta2)
end do
call destroy(Theta)
alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(alpc) < Cp%psi_tol)
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(abs(alpsave) - alpc) < err)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, size(LRV%Li)
call destroy(LRV%Li(pp))
end do
deallocate(dotvec, LRV%Li)
call destroy(Psp)
end subroutine orth_mpsc
"""
return
[docs]def orth_qmps():
"""
fortran-subroutine - Orthogonalize Psi against other MPS Phis.
**Arguments**
Psi : TYP(qmps), inout
Orthogonalize this MPS against the MPSs Phis.
Phis : , inout
The MPS Psi will be orthogonal to these MPSs.
minn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
maxn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
LR : TYPE(qtensorlist), inout
The partial overlaps at the orthogonality center.
err : REAL, out
Sum of all errors from splitting.
converged : LOGICAL, out
Flag for convergence of the variational orthogonalization.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine orth_qmps(Psi, Phis, minn, maxn, LR, Cp, err, converged, &
errst)
type(qmps), intent(inout) :: Psi
type(qmps), pointer, intent(inout) :: Phis(:)
integer, intent(in) :: minn, maxn
type(qtensorlist), pointer, intent(inout) :: LR(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
real(KIND=rKind), dimension(:), allocatable :: dotvec
! Copy of the original incoming Psi
type(qmps) :: Psp
! Variational overlap between original and new state
type(qtensorlist) :: Lrv
! Two site tensors
type(qtensor) :: Theta, Theta2
! temporary tensor
type(qtensor) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
if(verbose > 0) write(slog, *) 'Beginning two-site orthogonalization'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Psp is Psi and Psi becomes the variational state
call copy(Psp, Psi)
! variational overlap
call setuplr(Lrv, Psi, Psp, Psi%oc)
! Get workspace for dot product between Psi and Phis
allocate(dotvec(maxn))
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call orth_step_qmps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmps : step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call orth_step_qmps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmps : step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call orth_step_qmps(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmps : step (3) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = minn, maxn
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1], &
errst=errst)
!if(prop_error('orth_qmps : contr (1) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
!if(prop_error('orth_qmps : contr (2) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Tmpe)
dotvec(pp) = dot(Theta, Theta2, errst=errst)
!if(prop_error('orth_qmps : dot failed', &
! errst=errst)) return
call destroy(Theta2)
end do
call destroy(Theta)
alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(alpc) < Cp%psi_tol)
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(abs(alpsave) - alpc) < err)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, size(LRV%Li)
call destroy(LRV%Li(pp))
end do
deallocate(dotvec, LRV%Li)
call destroy(Psp)
end subroutine orth_qmps
"""
return
[docs]def orth_qmpsc():
"""
fortran-subroutine - Orthogonalize Psi against other MPS Phis.
**Arguments**
Psi : TYP(qmpsc), inout
Orthogonalize this MPS against the MPSs Phis.
Phis : , inout
The MPS Psi will be orthogonal to these MPSs.
minn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
maxn : INTEGER, in
First MPS in Phis to be orthogonal to Psi.
LR : TYPE(qtensorclist), inout
The partial overlaps at the orthogonality center.
err : REAL, out
Sum of all errors from splitting.
converged : LOGICAL, out
Flag for convergence of the variational orthogonalization.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine orth_qmpsc(Psi, Phis, minn, maxn, LR, Cp, err, converged, &
errst)
type(qmpsc), intent(inout) :: Psi
type(qmpsc), pointer, intent(inout) :: Phis(:)
integer, intent(in) :: minn, maxn
type(qtensorclist), pointer, intent(inout) :: LR(:)
type(ConvParam), intent(in) :: Cp
real(KIND=rKind), intent(out) :: err
logical, intent(out) :: converged
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, jj, kk, pp
! first and final site
integer :: k1
! direction of looping
integer :: sense
! number of sweeps used
integer :: sweeps_used
! actual bond dimension of Psi
integer :: chi
! Adapted local tolerance
real(KIND=rKind) :: local_tol
! Calculate the overlap to each Phi
complex(KIND=rKind), dimension(:), allocatable :: dotvec
! Copy of the original incoming Psi
type(qmpsc) :: Psp
! Variational overlap between original and new state
type(qtensorclist) :: Lrv
! Two site tensors
type(qtensorc) :: Theta, Theta2
! temporary tensor
type(qtensorc) :: Tmpe
! for the convergence
real(KIND=rKind) :: alpsave, alpc, desratio
logical :: chifail, stuck
if(verbose > 0) write(slog, *) 'Beginning two-site orthogonalization'
if(Cp%psi_local_tol < 0.0_rKind) then
! One sweep = 2 * L local truncations
local_tol = Cp%psi_tol / (4.0_rKind * Psi%ll)
else
local_tol = Cp%psi_local_tol
end if
! Psp is Psi and Psi becomes the variational state
call copy(Psp, Psi)
! variational overlap
call setuplr(Lrv, Psi, Psp, Psi%oc)
! Get workspace for dot product between Psi and Phis
allocate(dotvec(maxn))
! Save orthogonality center as starting site and initialize variables
k1 = Psi%oc
converged = .false.
stuck = .false.
chifail = .false.
sweeps_used = 0
alpsave = 10.0_rKind
err = 0.0_rKind
outer : do ii = 1, Cp%max_outer_sweeps
inner : do jj = 1, Cp%max_inner_sweeps
sweeps_used = sweeps_used + 1
sense = 1
do kk = k1, (Psi%ll - 2)
call orth_step_qmpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmpsc : step (1) failed.', &
! errst=errst)) return
end do
sense = -1
do kk = (Psi%ll - 1), 2, (-1)
call orth_step_qmpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmpsc : step (2) failed.', &
! errst=errst)) return
end do
sense = 1
do kk = 1, (k1 - 1)
call orth_step_qmpsc(Psp, Psi, Phis, Lrv, LR, kk, sense, &
minn, maxn, dotvec, Cp, local_tol, err, errst=errst)
!if(prop_error('orth_qmpsc : step (3) failed.', &
! errst=errst)) return
end do
! Convergence checks
! ..................
! Calculate overlaps
dotvec = 0.0_rKind
call contr(Theta, Psi%Aa(k1), Psi%Aa(k1 + 1), [3], [1])
do pp = minn, maxn
call contr(Theta2, Phis(pp)%Aa(k1), Phis(pp)%Aa(k1 + 1), &
[3], [1])
if(k1 == Psi%ll - 1) then
call copy(Tmpe, Theta2)
else
call contr(Tmpe, Theta2, LR(pp)%Li(k1 + 2), [4], [1], &
errst=errst)
!if(prop_error('orth_qmpsc : contr (1) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Theta2)
if(k1 == 1) then
call copy(Theta2, Tmpe)
else
call contr(Theta2, LR(pp)%Li(k1 - 1), Tmpe, [2], [1])
!if(prop_error('orth_qmpsc : contr (2) '//&
! 'failed.', errst=errst)) return
end if
call destroy(Tmpe)
dotvec(pp) = dot(Theta, Theta2, errst=errst)
!if(prop_error('orth_qmpsc : dot failed', &
! errst=errst)) return
call destroy(Theta2)
end do
call destroy(Theta)
alpc = real(dot_product(dotvec, dotvec), KIND=rKind)
chi = maxchi(Psi)
converged = (abs(alpc) < Cp%psi_tol)
if(jj >= Cp%min_inner_sweeps) then
! Condition for being "stuck" (not globally converged
! but fully optimized for the given parameters)
stuck = (abs(abs(alpsave) - alpc) < err)
alpsave = abs(alpc)
! The bond dimension is not sufficient
chifail = (chi >= Cp%max_bond_dimension)
! Exit if converged or stuck
if(stuck .or. converged) exit
else
alpsave = abs(alpc)
end if
! Set stuck to true in case we exit inner loop
stuck = .true.
end do inner
! Check if we need more outer loops
if(converged) exit
! Not yet optimized, but not stuck
if(.not. stuck) cycle
! Ran out of bond dimension
if(chifail) exit
! If stuck, assume that variance is a linear function of
! truncation error to extrapolate to require truncation error
desratio = Cp%psi_tol / abs(1.0_rKind - abs(alpc))
local_tol = max(desratio * local_tol, numzero)
stuck = .false.
end do outer
do pp = 1, size(LRV%Li)
call destroy(LRV%Li(pp))
end do
deallocate(dotvec, LRV%Li)
call destroy(Psp)
end subroutine orth_qmpsc
"""
return
[docs]def state_first_mps():
"""
fortran-subroutine - July 2017 (dj, updated)
Generate an initial guess for the first two sites when growing an MPS.
**Arguments**
Theta : TYPE(tensor), inout
On exit the guess for the two sites.
locdim : INTEGER, in
Local dimension of a site.
Hw : TYPE(sr_matrix_tensor), in
MPO matrix to extract local dimension.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
Quantum number of the total system for both symmetries.
IOObj : TYPE(IOObject), in
Dummy for the interface.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_first_mps(Theta, locdim, Hw, nqs, qs, IOObj)
type(tensor), intent(inout) :: Theta
integer, intent(out) :: locdim
type(sr_matrix_tensor), intent(in) :: Hw
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
! Local variables
! Set local dimension
locdim = Hw%Row(1)%Op(1)%dl(1)
call create(Theta, [1, locdim, locdim, 1], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_first_mps
"""
return
[docs]def state_first_mpsc():
"""
fortran-subroutine - July 2017 (dj, updated)
Generate an initial guess for the first two sites when growing an MPS.
**Arguments**
Theta : TYPE(tensorc), inout
On exit the guess for the two sites.
locdim : INTEGER, in
Local dimension of a site.
Hw : TYPE(sr_matrix_tensorc), in
MPO matrix to extract local dimension.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
Quantum number of the total system for both symmetries.
IOObj : TYPE(IOObject), in
Dummy for the interface.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_first_mpsc(Theta, locdim, Hw, nqs, qs, IOObj)
type(tensorc), intent(inout) :: Theta
integer, intent(out) :: locdim
type(sr_matrix_tensorc), intent(in) :: Hw
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
! Local variables
! Set local dimension
locdim = Hw%Row(1)%Op(1)%dl(1)
call create(Theta, [1, locdim, locdim, 1], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_first_mpsc
"""
return
[docs]def state_iter_mps():
"""
fortran-subroutine - July 2017 (dj, updated)
Create an initial guess for two sites when growing the MPS during
the iterations.
**Arguments**
Theta : TYPE(tensor), inout
On exit, the initial guess for the two-site tensor.
Tl : TYPE(tensor), in
The tensor to the left of the new two-site tensor Theta.
Tr : TYPE(tensor), in
The tensor to the right of the new two-site tensor Theta.
Lam : TYPE(tensor), in
Contains the singular values needed to follow an iMPS type
of guess. These are the values of the last splitting.
Lamprev : TYPE(tensor), in
Contains the singular values needed to follow an iMPS type
of guess. These are the values of the second last splitting.
locdim : INTEGER, in
The local dimension of the Hilbert space of a single site.
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
nn : INTEGER, in
nn-th iteration in growing the system.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_iter_mps(Theta, Tl, Tr, Lam, Lamprev, locdim, unit, &
nn, errst)
type(tensor), intent(inout) :: Theta
type(tensor), intent(in) :: Tl, Tr
type(tensor), intent(in) :: Lam, Lamprev
integer, intent(in) :: locdim
integer, intent(in) :: unit, nn
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Temporary new left/right tensor
type(tensor) :: Nl, Nr
! inverse lambda
type(tensor) :: Laminv
! Have to catch the first iteration
if(nn == 1) then
call create(Theta, [Tl%dl(3), Tl%dl(2), Tr%dl(2), Tr%dl(1)], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
return
end if
! Create an educated guess for the two inner sites using McCulloch's
! prediction algorithm
! A-> Lambda[n-1] B[n-1]
! B-> Lambda[n-2]^{-1} A[n-1] Lambda[N-1]
! Left bond dimension is set by previous truncation, right bond dimension
! is set by current truncation
! Create left tensor
call copy(Nl, Tr)
call contr_diag(Nl, Lam, 1, errst=errst)
!if(prop_error('state_iter_mps: contr_diag (1) failed.', &
! errst=errst)) return
! Create right tensor
call copy(Nr, Tl)
call contr_diag(Nr, Lam, 3, errst=errst)
!if(prop_error('state_iter_mps: contr_diag (2) failed.', &
! errst=errst)) return
if(is_set(Lamprev)) then
call create(Laminv, Lamprev%dl)
Laminv%elem(:Laminv%dim) = 1.0_rKind / Lamprev%elem(:Lamprev%dim)
call contr_diag(Nr, Laminv, 1, errst=errst)
!if(prop_error('state_iter_mps: contr_diag (3) failed.', &
! errst=errst)) return
call destroy(Laminv)
!else
! Singular values is one - contraction does not change anything
end if
! Form theta two-site tensor
call contr(Theta, Nl, Nr, [3], [1], errst=errst)
!if(prop_error('state_iter_mps: contr failed.', &
! errst=errst)) return
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
call destroy(Nl)
call destroy(Nr)
end subroutine state_iter_mps
"""
return
[docs]def state_iter_mpsc():
"""
fortran-subroutine - July 2017 (dj, updated)
Create an initial guess for two sites when growing the MPS during
the iterations.
**Arguments**
Theta : TYPE(tensorc), inout
On exit, the initial guess for the two-site tensor.
Tl : TYPE(tensorc), in
The tensor to the left of the new two-site tensor Theta.
Tr : TYPE(tensorc), in
The tensor to the right of the new two-site tensor Theta.
Lam : TYPE(tensor), in
Contains the singular values needed to follow an iMPS type
of guess. These are the values of the last splitting.
Lamprev : TYPE(tensor), in
Contains the singular values needed to follow an iMPS type
of guess. These are the values of the second last splitting.
locdim : INTEGER, in
The local dimension of the Hilbert space of a single site.
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
nn : INTEGER, in
nn-th iteration in growing the system.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_iter_mpsc(Theta, Tl, Tr, Lam, Lamprev, locdim, unit, &
nn, errst)
type(tensorc), intent(inout) :: Theta
type(tensorc), intent(in) :: Tl, Tr
type(tensor), intent(in) :: Lam, Lamprev
integer, intent(in) :: locdim
integer, intent(in) :: unit, nn
integer, intent(out), optional :: errst
! Local variables
! ---------------
! Temporary new left/right tensor
type(tensorc) :: Nl, Nr
! inverse lambda
type(tensor) :: Laminv
! Have to catch the first iteration
if(nn == 1) then
call create(Theta, [Tl%dl(3), Tl%dl(2), Tr%dl(2), Tr%dl(1)], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
return
end if
! Create an educated guess for the two inner sites using McCulloch's
! prediction algorithm
! A-> Lambda[n-1] B[n-1]
! B-> Lambda[n-2]^{-1} A[n-1] Lambda[N-1]
! Left bond dimension is set by previous truncation, right bond dimension
! is set by current truncation
! Create left tensor
call copy(Nl, Tr)
call contr_diag(Nl, Lam, 1, errst=errst)
!if(prop_error('state_iter_mpsc: contr_diag (1) failed.', &
! errst=errst)) return
! Create right tensor
call copy(Nr, Tl)
call contr_diag(Nr, Lam, 3, errst=errst)
!if(prop_error('state_iter_mpsc: contr_diag (2) failed.', &
! errst=errst)) return
if(is_set(Lamprev)) then
call create(Laminv, Lamprev%dl)
Laminv%elem(:Laminv%dim) = 1.0_rKind / Lamprev%elem(:Lamprev%dim)
call contr_diag(Nr, Laminv, 1, errst=errst)
!if(prop_error('state_iter_mpsc: contr_diag (3) failed.', &
! errst=errst)) return
call destroy(Laminv)
!else
! Singular values is one - contraction does not change anything
end if
! Form theta two-site tensor
call contr(Theta, Nl, Nr, [3], [1], errst=errst)
!if(prop_error('state_iter_mpsc: contr failed.', &
! errst=errst)) return
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
call destroy(Nl)
call destroy(Nr)
end subroutine state_iter_mpsc
"""
return
[docs]def state_single_mps():
"""
fortran-subroutine - July 2017 (dj, updated)
Random guess for the last site of an odd system size.
**Arguments**
Theta : TYPE(tensor), inout
On exit, the initial guess for the single-site tensor.
Tl : TYPE(tensor), in
The tensor to the left of the new single-site tensor Theta.
Tr : TYPE(tensor), in
The tensor to the right of the new single-site tensor Theta.
locdim : INTEGER, in
The local dimension of the Hilbert space of a single site.
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_single_mps(Theta, Tl, Tr, locdim, unit)
type(tensor), intent(inout) :: Theta
type(tensor), intent(in) :: Tl, Tr
integer, intent(in) :: locdim
integer, intent(in) :: unit
! No local variables
! ------------------
call create(Theta, [Tl%dl(3), locdim, Tr%dl(1)], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_single_mps
"""
return
[docs]def state_single_mpsc():
"""
fortran-subroutine - July 2017 (dj, updated)
Random guess for the last site of an odd system size.
**Arguments**
Theta : TYPE(tensorc), inout
On exit, the initial guess for the single-site tensor.
Tl : TYPE(tensorc), in
The tensor to the left of the new single-site tensor Theta.
Tr : TYPE(tensorc), in
The tensor to the right of the new single-site tensor Theta.
locdim : INTEGER, in
The local dimension of the Hilbert space of a single site.
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_single_mpsc(Theta, Tl, Tr, locdim, unit)
type(tensorc), intent(inout) :: Theta
type(tensorc), intent(in) :: Tl, Tr
integer, intent(in) :: locdim
integer, intent(in) :: unit
! No local variables
! ------------------
call create(Theta, [Tl%dl(3), locdim, Tr%dl(1)], init='R')
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_single_mpsc
"""
return
[docs]def state_finalize_mps():
"""
fortran-subroutine - July 2017 (dj)
Finalize the state, meaning closing the open unit if symmetries are
used.
**Arguments**
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_finalize_mps(unit, hasexcited)
integer, intent(in) :: unit
logical, intent(in), optional :: hasexcited
! No local variables
! ------------------
end subroutine state_finalize_mps
"""
return
[docs]def state_finalize_mpsc():
"""
fortran-subroutine - July 2017 (dj)
Finalize the state, meaning closing the open unit if symmetries are
used.
**Arguments**
unit : INTEGER, in
Dummy variable for interface. Not accessed without symmetries.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_finalize_mpsc(unit, hasexcited)
integer, intent(in) :: unit
logical, intent(in), optional :: hasexcited
! No local variables
! ------------------
end subroutine state_finalize_mpsc
"""
return
[docs]def state_first_qmps():
"""
fortran-subroutine - July 2017 (dj, updated)
Random guess for systems with symmetries for the first two sites.
**Arguments**
Theta : TYPE(qtensor), inout
On exit the guess for the two sites.
locdim : INTEGER, in
Dummy variable for interface.
Hw : TYPE(sr_matrix_qtensor), in
Dummy variable for interface.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
Quantum number of the total system for both symmetries.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file.
**Details**
The subroutine can be used to create the first two-site tensor
in the qIMPS. Therefore the total incoming charge is set to zero (left
tensor will be first site at the end) and the total outgoing charge are
the total charges. (Right tensor will be last site at the end of the
intialization.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_first_qmps(Theta, locdim, Hw, nqs, qs, IOObj)
type(qtensor), intent(inout) :: Theta
integer, intent(in) :: locdim
type(sr_matrix_qtensor), intent(in) :: Hw
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
! Local variables
! ---------------
! for looping
integer :: ii
! number of blocks qith different quantum numbers
integer :: nblocks
! degeneracy dimension of block
integer :: dl, dr
! Number of conserved quantities (from conserved numbers or
! each symmetry)
integer :: snqs
! for string formatter
character(16) :: iwstring, specstring
open(unit=IOObj%qstateunit, file=IOObj%qstatename, &
action='read', status='old')
snqs = sum(nqs)
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') 2 * snqs
specstring = "("//trim(adjustl(iwstring))//"I16)"
! Number of blocks, i.e. irrep tuples (qi qj)
read(IOObj%qstateunit, '(1I16)') nblocks
call create(Theta, nqs, nblocks)
! Get q#s and degeneracy dimensions of each block
do ii = 1, nblocks
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(4 * snqs))
Theta%Data(Theta%nb)%qq = 0
read(IOObj%qstateunit, specstring) &
Theta%Data(Theta%nb)%qq(snqs + 1:2 * snqs), &
Theta%data(Theta%nb)%qq(2 * snqs + 1:3 * snqs)
Theta%Data(Theta%nb)%qq(3 * snqs + 1:4 * snqs) = qs
read(IOObj%qstateunit, '(2I16)') dl, dr
call create(Theta%Data(Theta%nb)%Tens, [1, dl, dr, 1], init='R')
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_first_qmps
"""
return
[docs]def state_first_qmpsc():
"""
fortran-subroutine - July 2017 (dj, updated)
Random guess for systems with symmetries for the first two sites.
**Arguments**
Theta : TYPE(qtensorc), inout
On exit the guess for the two sites.
locdim : INTEGER, in
Dummy variable for interface.
Hw : TYPE(sr_matrix_qtensorc), in
Dummy variable for interface.
nqs : INTEGER, in
Number of abelian and discrete symmetries.
qs : INTEGER(\*), in
Quantum number of the total system for both symmetries.
IOObj : TYPE(IOObject), in
Contains name of the file with the quantum numbers and
the unit to open the file.
**Details**
The subroutine can be used to create the first two-site tensor
in the qIMPS. Therefore the total incoming charge is set to zero (left
tensor will be first site at the end) and the total outgoing charge are
the total charges. (Right tensor will be last site at the end of the
intialization.)
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_first_qmpsc(Theta, locdim, Hw, nqs, qs, IOObj)
type(qtensorc), intent(inout) :: Theta
integer, intent(in) :: locdim
type(sr_matrix_qtensorc), intent(in) :: Hw
integer, dimension(2), intent(in) :: nqs
integer, dimension(:), intent(in) :: qs
type(IOObject), intent(in) :: IOObj
! Local variables
! ---------------
! for looping
integer :: ii
! number of blocks qith different quantum numbers
integer :: nblocks
! degeneracy dimension of block
integer :: dl, dr
! Number of conserved quantities (from conserved numbers or
! each symmetry)
integer :: snqs
! for string formatter
character(16) :: iwstring, specstring
open(unit=IOObj%qstateunit, file=IOObj%qstatename, &
action='read', status='old')
snqs = sum(nqs)
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') 2 * snqs
specstring = "("//trim(adjustl(iwstring))//"I16)"
! Number of blocks, i.e. irrep tuples (qi qj)
read(IOObj%qstateunit, '(1I16)') nblocks
call create(Theta, nqs, nblocks)
! Get q#s and degeneracy dimensions of each block
do ii = 1, nblocks
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(4 * snqs))
Theta%Data(Theta%nb)%qq = 0
read(IOObj%qstateunit, specstring) &
Theta%Data(Theta%nb)%qq(snqs + 1:2 * snqs), &
Theta%data(Theta%nb)%qq(2 * snqs + 1:3 * snqs)
Theta%Data(Theta%nb)%qq(3 * snqs + 1:4 * snqs) = qs
read(IOObj%qstateunit, '(2I16)') dl, dr
call create(Theta%Data(Theta%nb)%Tens, [1, dl, dr, 1], init='R')
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
end subroutine state_first_qmpsc
"""
return
[docs]def state_iter_qmps():
"""
fortran-subroutine - July 2017 (dj, updated)
Generate a random state for two site when growing the MPS.
**Arguments**
Theta : TYPE(qtensor), inout
On exit, the initial guess for the two-site tensor.
Tl : TYPE(qtensor), in
The tensor to the left of the new two-site tensor Theta.
Tr : TYPE(qtensor), in
The tensor to the right of the new two-site tensor Theta.
Lam : TYPE(qtensor), in
Dummy variable for interface. Not accessed with symmetries.
Lamprev : TYPE(qtensor), in
Dummy variable for interface. Not accessed with symmetries.
locdim : INTEGER, in
Dummy variable for interface. Not accessed with symmetries.
unit : INTEGER, in
On this unit the file with the quantum numbers for the
symmetries.
nn : INTEGER, in
nn-th iteration in growing the system. Not referenced in symmetric
version.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_iter_qmps(Theta, Tl, Tr, Lam, Lamprev, locdim, unit, &
nn, errst)
type(qtensor), intent(inout) :: Theta
type(qtensor), intent(inout) :: Tl, Tr
type(qtensor), intent(in) :: Lam, Lamprev
integer, intent(in) :: locdim
integer, intent(in) :: unit, nn
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, pp
! index of matching hash
integer :: jj
! representative for quantum numbers
integer :: ak, bk
! degeneracy dimension of block
integer :: dl, dr
! Number of conserved quantities (retrieved from array with quantum numbers)
integer :: snqs
! number of blocks
integer :: nblocks
! quantum number for block
integer, dimension(:), allocatable :: qi, qj
! for sorting hashes
integer :: nunia, nunib
integer, dimension(:), allocatable :: inda, indb, dega, degb
real(KIND=rKind), dimension(:), allocatable :: littlea, littleb
! string specifier
character(16) :: iwstring, specstring
!if(present(errst)) errst = 0
! Could use Tl%nqs here
snqs = sum(Tl%nqs)
! Find matching blocks of A and B (guaranteed to exist by construction)
! by finding the (possibly degenerate) q_{\beta} irreps of A
! and searching for matches in B
call set_hash(Tl, [3])
call set_hash(Tr, [1])
! Now sort them, they can be degenerate
allocate(inda(Tl%nb), dega(Tl%nb + 1), littlea(Tl%nb), &
indb(Tr%nb), degb(Tr%nb + 1), littleb(Tr%nb))
call ascending_hsort(Tl%hash(:Tl%nb), littlea, inda, nunia, dega)
call ascending_hsort(Tr%hash(:Tr%nb), littleb, indb, nunib, degb)
! Shift q_{\alpha} of B to physical value q_{\beta}-q_j
! Does not affect hash!
do ii = 1, Tr%nb
! Improve generality?
Tr%Data(ii)%qq(:snqs) = sumq(Tr%Data(ii)%qq(2 * snqs + 1:3 * snqs), &
- Tr%Data(ii)%qq(snqs + 1:2 * snqs), &
Tr%nqs)
end do
read(unit, '(1I16)') nblocks
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') 2 * snqs
specString = "("//trim(adjustl(iwstring))//"I16)"
allocate(qi(snqs), qj(snqs))
call create(Theta, Tl%nqs, nblocks * nunia * nunib)
do ii = 1, nblocks
read(unit, specstring) qi, qj
read(unit, '(2I16)') dl, dr
do pp = 1, nunib
jj = findtagindex(littleb(pp), littlea(1:nunia))
if(jj > 0) then
! Pick a representative of the A block and B block
ak = dega(jj) + 1
bk = degb(pp) + 1
! Add new tensor
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(4 * snqs))
Theta%Data(Theta%nb)%qq(1:snqs) = &
Tl%data(inda(ak))%qq(2 * snqs + 1:3*snqs)
Theta%data(Theta%nb)%qq(snqs+1:2*snqs) = qi
Theta%data(Theta%nb)%qq(2*snqs+1:3*snqs) = qj
Theta%data(Theta%nb)%qq(3*snqs+1:4*snqs) = &
Tr%Data(indb(bk))%qq(1:snqs)
call create(Theta%Data(Theta%nb)%Tens, &
[Tl%Data(inda(ak))%Tens%dl(3), dl, dr, &
Tr%Data(indb(bk))%Tens%dl(1)], init='R')
end if
end do
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
deallocate(qi, qj)
deallocate(inda, dega, littlea, indb, degb, littleb)
end subroutine state_iter_qmps
"""
return
[docs]def state_iter_qmpsc():
"""
fortran-subroutine - July 2017 (dj, updated)
Generate a random state for two site when growing the MPS.
**Arguments**
Theta : TYPE(qtensorc), inout
On exit, the initial guess for the two-site tensor.
Tl : TYPE(qtensorc), in
The tensor to the left of the new two-site tensor Theta.
Tr : TYPE(qtensorc), in
The tensor to the right of the new two-site tensor Theta.
Lam : TYPE(qtensor), in
Dummy variable for interface. Not accessed with symmetries.
Lamprev : TYPE(qtensor), in
Dummy variable for interface. Not accessed with symmetries.
locdim : INTEGER, in
Dummy variable for interface. Not accessed with symmetries.
unit : INTEGER, in
On this unit the file with the quantum numbers for the
symmetries.
nn : INTEGER, in
nn-th iteration in growing the system. Not referenced in symmetric
version.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_iter_qmpsc(Theta, Tl, Tr, Lam, Lamprev, locdim, unit, &
nn, errst)
type(qtensorc), intent(inout) :: Theta
type(qtensorc), intent(inout) :: Tl, Tr
type(qtensor), intent(in) :: Lam, Lamprev
integer, intent(in) :: locdim
integer, intent(in) :: unit, nn
integer, intent(out), optional :: errst
! Local variables
! ---------------
! for looping
integer :: ii, pp
! index of matching hash
integer :: jj
! representative for quantum numbers
integer :: ak, bk
! degeneracy dimension of block
integer :: dl, dr
! Number of conserved quantities (retrieved from array with quantum numbers)
integer :: snqs
! number of blocks
integer :: nblocks
! quantum number for block
integer, dimension(:), allocatable :: qi, qj
! for sorting hashes
integer :: nunia, nunib
integer, dimension(:), allocatable :: inda, indb, dega, degb
real(KIND=rKind), dimension(:), allocatable :: littlea, littleb
! string specifier
character(16) :: iwstring, specstring
!if(present(errst)) errst = 0
! Could use Tl%nqs here
snqs = sum(Tl%nqs)
! Find matching blocks of A and B (guaranteed to exist by construction)
! by finding the (possibly degenerate) q_{\beta} irreps of A
! and searching for matches in B
call set_hash(Tl, [3])
call set_hash(Tr, [1])
! Now sort them, they can be degenerate
allocate(inda(Tl%nb), dega(Tl%nb + 1), littlea(Tl%nb), &
indb(Tr%nb), degb(Tr%nb + 1), littleb(Tr%nb))
call ascending_hsort(Tl%hash(:Tl%nb), littlea, inda, nunia, dega)
call ascending_hsort(Tr%hash(:Tr%nb), littleb, indb, nunib, degb)
! Shift q_{\alpha} of B to physical value q_{\beta}-q_j
! Does not affect hash!
do ii = 1, Tr%nb
! Improve generality?
Tr%Data(ii)%qq(:snqs) = sumq(Tr%Data(ii)%qq(2 * snqs + 1:3 * snqs), &
- Tr%Data(ii)%qq(snqs + 1:2 * snqs), &
Tr%nqs)
end do
read(unit, '(1I16)') nblocks
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') 2 * snqs
specString = "("//trim(adjustl(iwstring))//"I16)"
allocate(qi(snqs), qj(snqs))
call create(Theta, Tl%nqs, nblocks * nunia * nunib)
do ii = 1, nblocks
read(unit, specstring) qi, qj
read(unit, '(2I16)') dl, dr
do pp = 1, nunib
jj = findtagindex(littleb(pp), littlea(1:nunia))
if(jj > 0) then
! Pick a representative of the A block and B block
ak = dega(jj) + 1
bk = degb(pp) + 1
! Add new tensor
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(4 * snqs))
Theta%Data(Theta%nb)%qq(1:snqs) = &
Tl%data(inda(ak))%qq(2 * snqs + 1:3*snqs)
Theta%data(Theta%nb)%qq(snqs+1:2*snqs) = qi
Theta%data(Theta%nb)%qq(2*snqs+1:3*snqs) = qj
Theta%data(Theta%nb)%qq(3*snqs+1:4*snqs) = &
Tr%Data(indb(bk))%qq(1:snqs)
call create(Theta%Data(Theta%nb)%Tens, &
[Tl%Data(inda(ak))%Tens%dl(3), dl, dr, &
Tr%Data(indb(bk))%Tens%dl(1)], init='R')
end if
end do
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
deallocate(qi, qj)
deallocate(inda, dega, littlea, indb, degb, littleb)
end subroutine state_iter_qmpsc
"""
return
[docs]def state_single_qmps():
"""
fortran-subroutine - July 2018 (dj, updated)
Generate the random guess for the last site for system with odd number
of sites and symmetry.
**Arguments**
Theta : TYPE(qtensor), inout
On exit, the initial guess for the single-site tensor.
Tl : TYPE(qtensor), inout
The tensor to the left of the new single-site tensor Theta.
Tr : TYPE(qtensor), inout
The tensor to the right of the new single-site tensor Theta.
locdim : INTEGER, in
Dummy variable for interface. Not accessed with symmetries.
unit : INTEGER, in
On this unit the file with the quantum numbers for the
symmetries.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_single_qmps(Theta, Tl, Tr, locdim, unit)
type(qtensor), intent(inout) :: Theta
type(qtensor), intent(inout) :: Tl, Tr
integer, intent(in) :: locdim
integer, intent(in) :: unit
! Local variables
! ---------------
! for looping
integer :: ii, pp
! index of matching hash
integer :: jj
! representative for quantum numbers
integer :: ak, bk
! degeneracy dimension of block
integer :: dl
! Number of conserved quantities (retrieved from array with quantum numbers)
integer :: snqs
! number of blocks
integer :: nblocks
! quantum number for block
integer, dimension(:), allocatable :: qi
! for sorting hashes
integer :: nunia, nunib
integer, dimension(:), allocatable :: inda, indb, dega, degb
real(KIND=rKind), dimension(:), allocatable :: littlea, littleb
! string specifier
character(16) :: iwstring, specstring
snqs = sum(Tl%nqs)
! Find matching blocks of A and B (guaranteed to exist by construction)
! by finding the (possibly degenerate) q_{\beta} irreps of A
! and searching for matches in B
call set_hash(Tl, [3])
call set_hash(Tr, [1])
! Now sort them, they can be degenerate
allocate(inda(Tl%nb), dega(Tl%nb + 1), littlea(Tl%nb), &
indb(Tr%nb), degb(Tr%nb + 1), littleb(Tr%nb))
call ascending_hsort(Tl%hash(:Tl%nb), littlea, inda, nunia, dega)
call ascending_hsort(Tr%hash(:Tr%nb), littleb, indb, nunib, degb)
! Shift q_{\alpha} of B to physical value q_{\beta}-q_j
! Does not affect hash!
do ii = 1, Tr%nb
! Improve generality?
Tr%Data(ii)%qq(:snqs) = sumq(Tr%Data(ii)%qq(2 * snqs + 1:3 * snqs), &
- Tr%Data(ii)%qq(snqs + 1:2 * snqs), &
Tr%nqs)
end do
read(unit, '(1I16)') nblocks
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') snqs
specString = "("//trim(adjustl(iwstring))//"I16)"
allocate(qi(snqs))
call create(Theta, Tl%nqs, nblocks * Tl%nb * Tr%nb)
do ii = 1, nblocks
read(unit, specString) qi
read(unit, '(1I16)') dl
! B may also be degenerate, but we scan the entire range
do pp = 1, nunib
jj = findtagindex(littleb(pp), littlea(:nunia))
if(jj > 0) then
! Pick a representative of the A block and B block
ak = dega(jj) + 1
bk = degb(pp) + 1
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(3 * snqs))
Theta%Data(Theta%nb)%qq(1:snqs) = &
Tl%Data(inda(ak))%qq(2 * snqs + 1:3 * snqs)
Theta%Data(Theta%nb)%qq(snqs + 1:2 * snqs) = qi
Theta%Data(Theta%nb)%qq(2 * snqs + 1:3 * snqs) = &
Tr%Data(indb(bk))%qq(1:snqs)
call create(Theta%Data(Theta%nb)%Tens, &
[Tl%Data(inda(ak))%Tens%dl(3), dl, &
Tr%Data(indb(bk))%Tens%dl(1)], init='R')
end if
end do
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
deallocate(qi)
deallocate(inda, dega, littlea, indb, degb, littleb)
end subroutine state_single_qmps
"""
return
[docs]def state_single_qmpsc():
"""
fortran-subroutine - July 2018 (dj, updated)
Generate the random guess for the last site for system with odd number
of sites and symmetry.
**Arguments**
Theta : TYPE(qtensorc), inout
On exit, the initial guess for the single-site tensor.
Tl : TYPE(qtensorc), inout
The tensor to the left of the new single-site tensor Theta.
Tr : TYPE(qtensorc), inout
The tensor to the right of the new single-site tensor Theta.
locdim : INTEGER, in
Dummy variable for interface. Not accessed with symmetries.
unit : INTEGER, in
On this unit the file with the quantum numbers for the
symmetries.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_single_qmpsc(Theta, Tl, Tr, locdim, unit)
type(qtensorc), intent(inout) :: Theta
type(qtensorc), intent(inout) :: Tl, Tr
integer, intent(in) :: locdim
integer, intent(in) :: unit
! Local variables
! ---------------
! for looping
integer :: ii, pp
! index of matching hash
integer :: jj
! representative for quantum numbers
integer :: ak, bk
! degeneracy dimension of block
integer :: dl
! Number of conserved quantities (retrieved from array with quantum numbers)
integer :: snqs
! number of blocks
integer :: nblocks
! quantum number for block
integer, dimension(:), allocatable :: qi
! for sorting hashes
integer :: nunia, nunib
integer, dimension(:), allocatable :: inda, indb, dega, degb
real(KIND=rKind), dimension(:), allocatable :: littlea, littleb
! string specifier
character(16) :: iwstring, specstring
snqs = sum(Tl%nqs)
! Find matching blocks of A and B (guaranteed to exist by construction)
! by finding the (possibly degenerate) q_{\beta} irreps of A
! and searching for matches in B
call set_hash(Tl, [3])
call set_hash(Tr, [1])
! Now sort them, they can be degenerate
allocate(inda(Tl%nb), dega(Tl%nb + 1), littlea(Tl%nb), &
indb(Tr%nb), degb(Tr%nb + 1), littleb(Tr%nb))
call ascending_hsort(Tl%hash(:Tl%nb), littlea, inda, nunia, dega)
call ascending_hsort(Tr%hash(:Tr%nb), littleb, indb, nunib, degb)
! Shift q_{\alpha} of B to physical value q_{\beta}-q_j
! Does not affect hash!
do ii = 1, Tr%nb
! Improve generality?
Tr%Data(ii)%qq(:snqs) = sumq(Tr%Data(ii)%qq(2 * snqs + 1:3 * snqs), &
- Tr%Data(ii)%qq(snqs + 1:2 * snqs), &
Tr%nqs)
end do
read(unit, '(1I16)') nblocks
! Generate the format type for an integer string with partDigs digits
write(iwstring, '(I4)') snqs
specString = "("//trim(adjustl(iwstring))//"I16)"
allocate(qi(snqs))
call create(Theta, Tl%nqs, nblocks * Tl%nb * Tr%nb)
do ii = 1, nblocks
read(unit, specString) qi
read(unit, '(1I16)') dl
! B may also be degenerate, but we scan the entire range
do pp = 1, nunib
jj = findtagindex(littleb(pp), littlea(:nunia))
if(jj > 0) then
! Pick a representative of the A block and B block
ak = dega(jj) + 1
bk = degb(pp) + 1
Theta%nb = Theta%nb + 1
allocate(Theta%Data(Theta%nb)%qq(3 * snqs))
Theta%Data(Theta%nb)%qq(1:snqs) = &
Tl%Data(inda(ak))%qq(2 * snqs + 1:3 * snqs)
Theta%Data(Theta%nb)%qq(snqs + 1:2 * snqs) = qi
Theta%Data(Theta%nb)%qq(2 * snqs + 1:3 * snqs) = &
Tr%Data(indb(bk))%qq(1:snqs)
call create(Theta%Data(Theta%nb)%Tens, &
[Tl%Data(inda(ak))%Tens%dl(3), dl, &
Tr%Data(indb(bk))%Tens%dl(1)], init='R')
end if
end do
end do
call scale(1.0_rKind / sqrt(norm(Theta)), Theta)
deallocate(qi)
deallocate(inda, dega, littlea, indb, degb, littleb)
end subroutine state_single_qmpsc
"""
return
[docs]def state_finalize_qmps():
"""
fortran-subroutine - June 2018 (dj, updated)
Finalize the state, meaning closing the open unit if symmetries are
used.
**Arguments**
unit : INTEGER, in
This is unit was open to read to quantum numbers. Close and
delete file.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_finalize_qmps(unit, hasexcited)
integer, intent(in) :: unit
logical, intent(in), optional :: hasexcited
! No local variables
! ------------------
if(present(hasexcited)) then
if(.not. hasexcited) then
close(unit, status='delete')
else
close(unit)
end if
else
close(unit, status='delete')
end if
end subroutine state_finalize_qmps
"""
return
[docs]def state_finalize_qmpsc():
"""
fortran-subroutine - June 2018 (dj, updated)
Finalize the state, meaning closing the open unit if symmetries are
used.
**Arguments**
unit : INTEGER, in
This is unit was open to read to quantum numbers. Close and
delete file.
hasexcited : LOGICAL, in
Intended use to prevent deleting file with initial state
in ground state search, where initial state is needed as
well for the excited state.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine state_finalize_qmpsc(unit, hasexcited)
integer, intent(in) :: unit
logical, intent(in), optional :: hasexcited
! No local variables
! ------------------
if(present(hasexcited)) then
if(.not. hasexcited) then
close(unit, status='delete')
else
close(unit)
end if
else
close(unit, status='delete')
end if
end subroutine state_finalize_qmpsc
"""
return