joselado / dmrgpy

DMRGPy is a Python library to compute quasi-one-dimensional spin chains and fermionic systems using matrix product states with DMRG as implemented in ITensor. Most of the computations can be performed both with DMRG and exact diagonalization for small systems, which allows one to benchmark the results.
GNU General Public License v3.0
87 stars 20 forks source link

New folder in dmrgpy #19

Closed HamidArianZad closed 1 year ago

HamidArianZad commented 1 year ago

Dear Jose,

I installed your great package dmrgpy in ubuntu 2022. It works when I run the predefined python scripts in folder examples in dmrgpy. However, I need to simulate the magnetization of Heisenberg chains and ladders for different values of exchange interactions so I need to create separated folder for each one to have their source codes. Hence, I created a separated folder with a different name in dmrgpy folder that includes a main.py (for test, I just copied your script from folder 1dIsing) . But whenever I run the code I get below error:

Traceback (most recent call last): File "/home/hamid/dmrgpy/examples/DSL Project_dmrgpy/DSL_dmrgpy_Delta1J1_S40/main.py", line 9, in <module> sc = spinchain.Spin_Chain(spins) # create the spin chain File "/home/hamid/anaconda3/lib/python3.9/site-packages/dmrgpy/spinchain.py", line 52, in __init__ Many_Body_Chain.__init__(self,sites) File "/home/hamid/anaconda3/lib/python3.9/site-packages/dmrgpy/manybodychain.py", line 88, in __init__ self.initialize() File "/home/hamid/anaconda3/lib/python3.9/site-packages/dmrgpy/manybodychain.py", line 93, in initialize self.execute(lambda: write_sites(self)) # write the different sites File "/home/hamid/anaconda3/lib/python3.9/site-packages/dmrgpy/manybodychain.py", line 415, in execute self.to_folder() # go to folder File "/home/hamid/anaconda3/lib/python3.9/site-packages/dmrgpy/manybodychain.py", line 113, in to_folder os.chdir(self.path) # go to calculation folder FileNotFoundError: [Errno 2] No such file or directory: '/home/hamid/dmrgpy/examples/DSL Project_dmrgpy/DSL_dmrgpy_Delta1J1_S40/.mpsfolder/' While, the same main.py script in folder 1dIsing and also other folders runs normally. What should I do to have my own folder in dmrgpy folder?

joselado commented 1 year ago

Dear HamidArianZad ,

thank you for your message. To make the script work in any folder you need to replace the line

import os ; import sys ; sys.path.append(os.getcwd()+'/../../src')

by

import os ; import sys ; sys.path.append(YOUR_PATH_TO_DMRGPY_SRC)

where YOUR_PATH_TO_DMRGPY_SRC is the direction to the dmrgpy library. In the original examples, that direction is written as a relative path, so that if you put the example elsewhere you need to put the right path for the library.

I hope that this solves the issue you had.

Best regards, Jose

HamidArianZad commented 1 year ago

Thank you for your quick respond, However, with changing the path, I still get the same error.

In fact I changed the path: import os ; import sys ; sys.path.append(os.getcwd()+'/../../src') to sys.path.append('/home/hamid/dmrgpy/examples/DSLProject_dmrgpy')

in which /home/hamid/dmrgpy/examples/DSLProject_dmrgpy is the path of my own main.py file. In fact, I just created an extra folder DSLProject_dmrgpy (including a main.py) file in the 'examples' folder.

joselado commented 1 year ago

Dear HamidArianZad,

Try replacing the line

sys.path.append('/home/hamid/dmrgpy/examples/DSLProject_dmrgpy')

by

sys.path.append('/home/hamid/dmrgpy/src')

The path should the one of the dmrgpy libray, not the main script. I hope this helps.

Best regards, Jose

HamidArianZad commented 1 year ago

Thank you very much, It works!!

HamidArianZad commented 1 year ago

Dear Jose,

Thank you for your great package dmrgpy, I could go through most of the files related to various spin models, fermionic and bosonic systems. It is surprising that this package also supports mixed-spin models. Even more interesting, I could solve a particular Heisenberg spin model with cyclic four-spin exchange interaction by implementing dmrg, that is great.

` n = 8 spins = ["S=1/2" for i in range(n)] # spin 1/2 heisenberg chain sc = spinchain.Spin_Chain(spins) # create the spin chain h = 0 mz_t = np.zeros(n) # global magnetization E_GSt = np.zeros(n) # ground-state energy h_GS = np.zeros(n) # magnetic field

mz_t[0] = 0
h_GS[0] = 0
E_GSt[0] = 0

for i in range(1, n): mz_t[i] = i/n # global magnetization per site h = h + sc.Sx[i-1] sc.Sx[i] h = h + sc.Sy[i-1] sc.Sy[i] h = h + sc.Sz[i-1] * sc.Sz[i] sc.set_hamiltonian(h) E_GSt[i] = sc.gs_energy() # compute the ground-state energy of the model in each site print("Finish job ", i, "/", n-1)

print("Energy", E_GSt)

for i in range(1, n): h_GS[i] = E_GSt[i] - E_GSt[i-1] # computing magnetic field per site i from subtracting local GS_energies

print("h_GS", h_GS)

inds = range(n) # indexes of the sites plt.plot(h_GS, mz_t, marker="o", label="TotMag vs field")

plt.figlegend()

plt.show() `

But I get different result from DMRG results of ALPS package. I checked all folders of dmrgpy to address this problem, but in all files you defined magnetic field b as a separated parameter. How can I compute the global magnetization, ground-state energy per site E_GSt[i] and local magnetic field b[i] using dmrgpy?

joselado commented 1 year ago

Dear HamidArianZad,

thank you for the nice feedback.

You can compute the global magnetization and local magnetic field with

mzi = [sc.vev(szi)) for szi in sc.Sz] # magnetization in all the sites mz = sum(mzi) # total magnetization

regarding the "ground-state energy per site E_GSt[i]", this would depend on how you would like to define this quantity in terms of expectation values. If you would define the "local energy" with a certain operator Hi, you could just compute

ei = sc.vev(Hi)

Let me know if this answers your question.

Best regards, Jose

HamidArianZad commented 1 year ago

Dear Jose,

I used your notation to obtain the corresponding global magnetization vs local magnetic field. I still get different result from DMRG of ALPS and ED method of QuSpin. Actually, I used below code for obtaining the global magnetization and the ground-state energy of a spin-1/2 model using ED method implemented by QuSpin:


L = 20                      # number of spin-1/2 particles
sz_min = 0          
sz_max = 10

E0_Exact = np.zeros(sz_max+1)     # ground-state energy
Mag = np.zeros(sz_max+1)            # global magnetization 

for j in range(sz_min, sz_max+1):
    Mag[j] = j          # global magnetization
    # compute spin-1/2 basis
        basis = spin_basis_1d(L, Nup=sz_max+j)    # Nup is the number of up spins for which  sz_max<= Nup <= L.
        H_XXZ = hamiltonian(static, dynamic, basis=basis, dtype=np.float64)
        E_exact, V_exact = H_XXZ.eigsh(k=1, which='SA')    # ground-state energy of the system when Nup = sz_max+j
        E0_Exact[j] = E_exact

b_exact = np.zeros(sz_max+1)    # local magnetic field
b_exact[0] = 0

for i in range(sz_min, sz_max+1):
       b_exact[i] = E0_Exact[i] - E0_Exact[i-1]

In above code, sz_max=L//2 and sz_min = 0 that are eigenvalues of the operator SzT. In this code in each step j I considered the Hamiltonian of the system when L//2+j number of spins are up. Hence, the global magnetization and the corresponding ground-state energy can be derived for each step j.

I used the same procedure in dmrgpy as following code:


n = 20                      # number of spin-1/2 particles
sz_min = 0          
sz_max = 10

spins = ["S=1/2" for i in range(n)]  # spin 1/2 heisenberg chain
sc = spinchain.Spin_Chain(spins)  # create the spin chain
h = 0

E_GSt = np.zeros(sz_max+1)    # Ground-state energy
Mz_t = np.zeros(sz_max+1)

E_GSt[0] = 0    
Mz_t[0] = 0         # total magnetization

for i in range(1, n):
    h = h + sc.Sx[i-1] * sc.Sx[i]
    h = h + sc.Sy[i-1] * sc.Sy[i]
    h = h + sc.Sz[i-1] * sc.Sz[i]

sc.set_hamiltonian(h)

for i in range(1, sz_max+1):
    sc.Nup = sz_max + i        # selecting sz_max+i number of up spins
    mz_t[i] = i / n  # total magnetization per site
    E_GSt[i] = sc.vev(h)   # ground-state energy of the system

b_GS = np.zeros(sz_max+1)    # magnetic field
b_GS[0] = 0       

for i in range(1, sz_max+1):
    b_GS[i] = E_GSt[i] - E_GSt[i-1]  # computing magnetic field per site

b_GS[sz_max] = E_GSt[sz_max] - E_GSt[sz_max-1]

But the result of the total magnetization per site mz_t versus magnetic field b_GS is incorrect and is different from the ED method used in QuSpin. The problem is that, here the dmrgpy can not recognize the selected number of up spins Nup.

joselado commented 1 year ago

Dear HamidArianZad,

dmrgpy works in the full Hilbert space, namely that one does not determine the number of up or down spins. The ground state obtained is the lowest energy state among all the manifolds.

Perhaps you can try comparing with a calculation where in those codes you do not constrain the number of up or down spins.

For example, you can compare how the magnetization VS magnetic field for Heisenberg model behaves for the different codes, taking expectation values always in the lowest energy state (dmrgpy does this automatically).

Best regards, Jose

HamidArianZad commented 1 year ago

Dear Jose, Thank you for your advise. I will calculate the magnetization of the spin model with respect to the magnetic field by adding: h = h + sc.Sz[i] * b.

n = 50
spins = ["S=1" for i in range(n+1)]  # spin 1 heisenberg chain
sc = spinchain.Spin_Chain(spins)     # create the spin chain

h = 0

J = 0.000276
beta = 0.00218

Si = [sc.Sx, sc.Sy, sc.Sz]  # store the three components

for i in range(1, n+1):
    h = h + J*sc.Sx[i-1] * sc.Sx[i]
    h = h + J*sc.Sy[i-1] * sc.Sy[i]
    h = h + J*sc.Sz[i-1] * sc.Sz[i]
    for S in Si:
        h = h + J*beta*S[i-1] * S[i] * S[i-1] * S[i]  # biquadratic

algebra.maxsize = 70000
sc.nsweeps = 10
sc.maxm = 100
sc.kpmmaxm = sc.maxm
sc.set_hamiltonian(h)
mzi = [sc.vev(szi) for szi in sc.Sz]  # magnetization in all the sites
Mz_t = sum(mzi)  # total magnetization

inds = range(n+1)  # indexes of the sites

plt.plot(inds, mzi,  color="b", linewidth=2, marker="o")

Below figures show the magnetization in sites for two sequential running the same code:

Sz_spin1_KS100_a Sz_spin1_KS100_b

. Whenever I run the code again and again, I get new results. Does it mean that the number of sweeps and maxm should be increased to converge the ground-state energy such that the same results are obtained during the successive running the code?

However I computed the magnetization in sites of the spin-1 modle with n=10, for sc.maxm = 300 and sc.nsweeps = 40, and again I get different results fpr successive running the code (see below figures):

Mag_BLBQ_S10_a Mag_BLBQ_S10_b


Could you also please tell me what is the cause of below error?

terminate called after throwing an instance of 'itensor::ITError'
  what():  Couldn't open file "psi_GS.mps" for writing
Aborted (core dumped)
Traceback (most recent call last):
  File "/home/hamid/dmrgpy/examples/DSLProject_dmrgpy/BLBQspin1_Mag_B.py", line 58, in <module>
    mz_dmrg = sc.vev(Mz, mode="DMRG").real

As you recommended, I computed the magnetization of a spin-1 chain with n=50 sites by introducing magnetic field b (i.e., b in np.arange(0, 1, 0.01)). I got the magnetization results for a few steps of field b (i.e., b=0,0.01, 0.02), but for higher values of the b the code is terminated by above error.

Thank you for your patience and your replies.

joselado commented 1 year ago

Dear HamidArianZad,

I think that the first point to keep in mind is that is good to normalize the energy scales, namely that the exchange coupling is 1. This is done for numerical stability as all the thresholds are defined with respect to the biggest energy scale.

Regarding the change in the results. The issue is that your ground state is two-fold degenerate in both cases (in the first because of topological degeneracy, in the second because you have an odd number of sites). Therefore, the ground state you obtain each time is a linear combination of the two ground states. If you compute expectation values that are invariant irrespective of which GS you take, then you would get the same result.

Best regards, Jose

HamidArianZad commented 1 year ago

Dear Jose,

Thank you very much. I could eventually overcome this problem by adding a very small magnetic field for i in range(0, n): h = h - 0.00005*sc.Sz[i] to break degeneracy and I achieved the lowest eigenenergy with pm 0.00005 error. This very small amount of magnetic field can be neglected or can be subtracted from the obtained ground-state energy. Now I can gain the same result for the average magnetization in sites, , during sequential running the code.