StatesOverlap.py

import numpy as np
import scipy as sp
import sys
import argparse
import time
import math
import HofstadterThreeBody as HHModule
import GenericModule as gm
    
def SingleOverlap(firstState, secondState):
    """
    Optimize version of CalculateMatrixElements() using CreateLGMatrixOptimized()
    """
   
    overlap = ( np.vdot(firstState,secondState) )
        
    #print(deltaEnergies)

    return overlap
    
def MultipleOverlap(firstState, eigenstates, nbrSecondStates, threshold=0.2):
    """
    """
    overlaps = []
    largeOverlaps = []
    indexLargeOverlaps = []
    for n in np.arange(0, nbrSecondStates):
        currOverlap = SingleOverlap(firstState, eigenstates[:,n])
        overlaps.append(currOverlap)
        if np.abs(currOverlap) > threshold:
            largeOverlaps.append(np.abs(currOverlap))
            indexLargeOverlaps.append(n)
        
    return np.abs(overlaps), largeOverlaps, indexLargeOverlaps
    

# System parameters
J = 1.
U = 0.
U3 = 0.
FluxDensity = 0.2
trapConf = 0.
threshold = 0.2

# Handle the parameters parsed by the command line
parser = argparse.ArgumentParser()
parser.add_argument('-N', type=int, help='number of particles')
parser.add_argument('-L', type=int, help='side of the square lattice of size LxL')
parser.add_argument('-J', type=float, help='tunneling energy')
parser.add_argument('-U', type=float, help='two-body onsite interaction (only in softcore mode)')
parser.add_argument('-Utwo', type=float, help='two-body onsite interaction (only in softcore mode) for the second state')
parser.add_argument('-U3', type=float, help='three-body onsite interaction (only in softcore mode)')
parser.add_argument('-U3two', type=float, help='three-body onsite interaction (only in softcore mode) for the second state')
parser.add_argument('--conf', type=float, help='harmonic trap confinement strength (v0) as v0 * r^2')
parser.add_argument('--gamma', type=float, default=2, help='trap steepness (g) as v0 * (r)^g (default=2)')
parser.add_argument('--alpha', type=float, help='magnetic flux density as alpha=p/q')
parser.add_argument('--alphatwo', type=float, help='magnetic flux density as alpha=p/q for the second state')
parser.add_argument('--hardcore', type=int, nargs='?', const=1, default=0, help='hardcore bosons mode')
parser.add_argument('-n1', type=int, default=0, help='eigenstate number of the first state (e.g. n=0 is the ground state)')
parser.add_argument('-n2', type=int, default=0, help='eigenstate number of the second state (e.g. n=0 is the ground state)')
# Multiple overlap
parser.add_argument('--multipleoverlap', type=int, nargs='?', const=1, default=0, help='to calculate matrix elements with quenched states from H = H_0 + eps*O_l (laser parameters to be specified)')
parser.add_argument('--nbrsecondstates', type=int, help='the number of excited states to do the overlap with <n1|...>')
# Multiple overlaps varying simultaneously U,U3 parameters
parser.add_argument('--uoverlap', type=int, nargs='?', const=1, default=0, help='enable the fidelity check spanning U,U3 parameters')
parser.add_argument('--Ustep', type=float, help='increasing/decreasing step for U/U3 (initial values are --Utwo and --U3two)')
# Fidelity check changing the flux density alpha
parser.add_argument('--overlapflux', type=int, nargs='?', const=1, default=0, help='enable the fidelity check spanning U,U3 parameters')
parser.add_argument('--alphainit', type=float, help='increasing/decreasing step for U/U3 (initial values are --Utwo and --U3two)')
parser.add_argument('--alphafinal', type=float, help='increasing/decreasing step for U/U3 (initial values are --Utwo and --U3two)')
parser.add_argument('--alphastep', type=float, help='increasing/decreasing step for U/U3 (initial values are --Utwo and --U3two)')
# Fidelity defined as F_alpha = |<psi_n1(alpha)|psi_n1(alpha - dAlpha)>|ˆ2
parser.add_argument('--fidelityflux', type=int, nargs='?', const=1, default=0, help='calculate F_alpha = |<psi_n1(alpha)|psi_n1(alpha - dAlpha)>|^2 spanning different fluxes (need to specify -n1 and --alphainit, --alphafinal, --alphastep)')
# Fidelity defined as F_alpha = |<psi_n1(U,U3)|psi_n1(U-dU,U3-dU3)>|ˆ2
parser.add_argument('--ufidelity', type=int, nargs='?', const=1, default=0, help='calculate F_alpha = |<psi_n1(U,U3)|psi_n1(U-dU,U3-dU3)>|ˆ2 spanning different U,U3 (need to specify -n1 and --Ustep)')
parser.add_argument('--ufidelityadaptive', type=int, nargs='?', const=1, default=0, help='calculate F_alpha = |<psi_n(U,U3)|psi_np(U-dU,U3-dU3)>|ˆ2 where the indices n,np adapts to the maximum overlaps generated by StatesOverlap.py')
args = parser.parse_args()

if args.N is not None: N = args.N
if args.L is not None: L = args.L
if args.J is not None: J = args.J

if args.U is not None: U = args.U

if args.Utwo is not None:
    Utwo = args.Utwo
else:
    Utwo = 0

if args.U3 is not None: U3 = args.U3

if args.U3two is not None:
    U3two = args.U3two
else:
    U3two = 0

if args.alpha is not None: FluxDensity = args.alpha

if args.alphatwo is not None:
    FluxDensityTwo = args.alphatwo
else:
    FluxDensityTwo = FluxDensity

if args.conf is not None: trapConf = args.conf
if args.nbrsecondstates is not None: nbrSecondStates = args.nbrsecondstates

if args.n1 is not None: n1 = args.n1
if args.n2 is not None: n2 = args.n2

gamma = args.gamma

if args.hardcore == 0:
    hardcore = False
elif args.hardcore == 1:
    hardcore = True

Ns = L*L
c = HHModule.FindCenter(L)

Dim = int(gm.HilbertDim(N,Ns,hardcore))

fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N)
print(f'Loading the first state --> {fileName}')
firstState = gm.LoadVector(fileName)

if args.multipleoverlap == 1:
    # Load the eigenvectors for the multiple overlaps
    
    print('Loading all the eigenstates for the overlap...')
    eigVec = np.zeros((Dim,nbrSecondStates), dtype=complex)
    for n in np.arange(0,nbrSecondStates):
        fileName = gm.GenFilename(hardcore, L, J, Utwo, trapConf, gamma, n, U3=U3two, alpha=FluxDensityTwo, N=N)
        print(f'--> {fileName}')
        eigVec[:,n] = gm.LoadVector(fileName)

    print('Calculating the overlaps...')
    overlaps, largeOverlaps, idxLargeOverlaps = MultipleOverlap(firstState, eigVec, nbrSecondStates, threshold=threshold)
    print(f'Large overlaps (>{threshold}) and associated eigenstate index: ')
    for n in np.arange(0,len(largeOverlaps)):
        print(f'n={idxLargeOverlaps[n]}, overlap={largeOverlaps[n]}')
        
if args.uoverlap == 1:
    # file saving
    U3vals = []
    maxOverlapIndex = []
    maxOverlaps = []

    dU = args.Ustep
    nbrDecSteps = int(Utwo / dU)
    eigVec = np.zeros((Dim,nbrSecondStates), dtype=complex)
    for m in np.arange(0,nbrDecSteps):
        print(f'Loading all the eigenstates for U={round(Utwo,2)}, U3={round(U3two,2)}...')
        for n in np.arange(0,nbrSecondStates):
            fileName = gm.GenFilename(hardcore, L, J, round(Utwo,2), trapConf, gamma, n, U3=round(U3two,2), alpha=FluxDensityTwo, N=N)
            print(f'--> {fileName}')
            eigVec[:,n] = gm.LoadVector(fileName)
        overlaps, largeOverlaps, idxLargeOverlaps = MultipleOverlap(firstState, eigVec, nbrSecondStates, threshold=threshold)
        
        print(f'Large overlaps (>{threshold}) and associated eigenstate index: ')
        print(f'n={idxLargeOverlaps}, overlap={largeOverlaps}')
        U3vals.append(U3two)
        idxMax = np.argmax(largeOverlaps)
        maxOverlapIndex.append(idxLargeOverlaps[idxMax])
        maxOverlaps.append(np.max(largeOverlaps))
        
        # Here we save all the overlaps for a given interaction strenght
        fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N) + '_alloverlaps'
        gm.AppendArraysTwoColFile(fileName, np.full((len(overlaps)),U3two), overlaps)
            
        Utwo = Utwo - dU
        U3two = U3two + dU
            
    fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N, maxOverlap=True)
    gm.SaveArraysThreeColFile(fileName, U3vals, maxOverlapIndex, maxOverlaps)
    print(f'Maximum overlaps saved in --> {fileName}')
        
if args.overlapflux == 1:
    alphaList = []
    maxOverlapIndex = []
    maxOverlaps = []

    dFlux = args.alphastep
    alphaInit = args.alphainit
    alphaFinal = args.alphafinal
    
    eigVec = np.zeros((Dim,nbrSecondStates), dtype=complex)
    for a in np.arange(alphaInit,alphaFinal,dFlux):
        print(f'Loading all the eigenstates for alpha={round(a,4)}...')
        for n in np.arange(0,nbrSecondStates):
            fileName = gm.GenFilename(hardcore, L, J, round(Utwo,2), trapConf, gamma, n, U3=round(U3two,2), alpha=round(a,4), N=N)
            print(f'--> {fileName}')
            eigVec[:,n] = gm.LoadVector(fileName)
        overlaps, largeOverlaps, idxLargeOverlaps = MultipleOverlap(firstState, eigVec, nbrSecondStates, threshold=threshold)
        
        print(f'Large overlaps (>{threshold}) and associated eigenstate index: ')
        print(f'n={idxLargeOverlaps}, overlap={largeOverlaps}')
        print(f'{overlaps}')
        alphaList.append(a)
        idxMax = np.argmax(largeOverlaps)
        maxOverlapIndex.append(idxLargeOverlaps[idxMax])
        maxOverlaps.append(np.max(largeOverlaps))
        
    fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N, maxOverlap=True)
    gm.SaveArraysThreeColFile(fileName, alphaList, maxOverlapIndex, maxOverlaps)
    print(f'Maximum overlaps saved in --> {fileName}')
    
if args.fidelityflux == 1:
    dFlux = args.alphastep
    alphaInit = args.alphainit
    alphaFinal = args.alphafinal
    
    fidelity = []

    for a in np.arange(alphaInit,alphaFinal,dFlux):
        eigVec = np.zeros((Dim,2), dtype=complex)
        fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=round(a,6), N=N)
        eigVec[:,0] = gm.LoadVector(fileName)
        fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=round(a+dFlux,6), N=N)
        eigVec[:,1] = gm.LoadVector(fileName)

        print(f'Doing the overlap <psi_{n1}({round(a,6)})|psi_{n1}({round(a+dFlux,6)})>...')
        nextOverlap = SingleOverlap(eigVec[:,0], eigVec[:,1])
        fidelity.append(np.abs(nextOverlap)**2)
        
    print('Fidelities:')
    print(fidelity)
    
if args.ufidelity == 1:
    fidelity = []
    U3vals = []

    dU = args.Ustep
    nbrDecSteps = int(U / dU)
    for m in np.arange(0,nbrDecSteps):
        eigVec = np.zeros((Dim,2), dtype=complex)
        fileName = gm.GenFilename(hardcore, L, J, round(U,6), trapConf, gamma, n1, U3=round(U3,6), alpha=FluxDensity, N=N)
        eigVec[:,0] = gm.LoadVector(fileName)
        fileName = gm.GenFilename(hardcore, L, J, round(U-dU,6), trapConf, gamma, n1, U3=round(U3+dU,6), alpha=FluxDensity, N=N)
        eigVec[:,1] = gm.LoadVector(fileName)
        
        print(f'Doing the overlap <psi_{n1}({round(U,6)},{round(U3,6)})|psi_{n1}({round(U-dU,6)},{round(U3+dU,6)})>...')
        nextOverlap = SingleOverlap(eigVec[:,0], eigVec[:,1])
        fidelity.append(np.abs(nextOverlap)**2)
        U3vals.append(U3)
        
        U = U - dU
        U3 = U3 + dU

    fileName = gm.GenFilename(hardcore, L, J, 0, trapConf, gamma, n1, 0, alpha=FluxDensity, N=N) + '_ufidelity'
    gm.SaveArraysTwoColFile(fileName, U3vals, fidelity)
    print(f'Fidelity saved in ---> {fileName}')
    
    print('Fidelities:')
    print(fidelity)
    
    
if args.ufidelityadaptive == 1:
    fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N, maxOverlap=True)
    U3vals, maxOverlapIndex, maxOverlaps = gm.LoadFileThree(fileName)
    
    fidelity = []
    for m in np.arange(0,len(U3vals)-1):
        eigVec = np.zeros((Dim,2), dtype=complex)
        fileName = gm.GenFilename(hardcore, L, J, round(8.1-U3vals[m],6), trapConf, gamma, int(maxOverlapIndex[m]), U3=round(U3vals[m],6), alpha=FluxDensity, N=N)
        eigVec[:,0] = gm.LoadVector(fileName)
        fileName = gm.GenFilename(hardcore, L, J, round(8.1-U3vals[m+1],6), trapConf, gamma, int(maxOverlapIndex[m+1]), U3=round(U3vals[m+1],6), alpha=FluxDensity, N=N)
        eigVec[:,1] = gm.LoadVector(fileName)
        
        print(f'Doing the overlap <psi_{int(maxOverlapIndex[m])}({round(8.1-U3vals[m],6)},{round(U3vals[m],6)})|psi_{int(maxOverlapIndex[m+1])}({round(8.1-U3vals[m+1],6)},{round(U3vals[m+1],6)})>...')
        nextOverlap = SingleOverlap(eigVec[:,0], eigVec[:,1])
        fidelity.append(np.abs(nextOverlap)**2)
        
    print('Fidelities:')
    print(fidelity)
    
    fileName = gm.GenFilename(hardcore, L, J, U, trapConf, gamma, n1, U3=U3, alpha=FluxDensity, N=N) + '_ufidelityadaptive'
    gm.SaveArraysTwoColFile(fileName, [U3vals[i] for i in np.arange(0,len(fidelity))], fidelity)
    print(f'Fidelities saved in ---> {fileName}')