Issue with Quantum Phase Transition for Ising Model in NetKet

[mahmudams]

I'm working on checking the quantum phase transition for the Ising model using NetKet. My goal is to get the plot below:

As an initial step, I computed the magnetization (defined as the structure factor in the code) for the case when h=0.

I utilized the exact sampler to compute the magnetization over 400,000 iterations. However, I'm consistently getting zero magnetization, which is clearly incorrect since the expected value should be 1.

Below is the code I'm using. Could anyone help me identify what might be causing this issue and how to resolve it?

Thank you in advance for your assistance!

# Define the 1D chain graph
g = nk.graph.Hypercube(length=L, n_dim=1, pbc=False)

# Define the Hilbert space of spin-1/2 particles on the graph
hi = nk.hilbert.Spin(s=0.5, N=g.n_nodes)

H = nk.operator.Ising(h=0.0, hilbert=hi, J=1, graph=g)

class JasShort(nn.Module):
    @nn.compact
    def __call__(self, x):
        # Define the two variational parameters J1 and J2
        j1 = self.param(
            "j1", nn.initializers.normal(), (1,), float
        )
        j2 = self.param(
            "j2", nn.initializers.normal(), (1,), float
        )

        # Compute the nearest-neighbor correlations
        corr1 = x * jnp.roll(x, -1, axis=-1)
        corr2 = x * jnp.roll(x, -2, axis=-1)

        # Sum the output
        return jnp.sum(j1 * corr1 + j2 * corr2, axis=-1)

model = JasShort()

# Create the local sampler on the Hilbert space
#sampler = nk.sampler.MetropolisLocal(hi)
sampler = nk.sampler.ExactSampler(hi)

vstate = nk.vqs.MCState(sampler, model, n_samples=1008)

optimizer = nk.optimizer.Sgd(learning_rate=0.01)

gs = nk.driver.VMC(H, optimizer, variational_state=vstate, preconditioner=nk.optimizer.SR(diag_shift=0.01))

# Construct the logger
test = nk.logging.RuntimeLog()

jas_energy = vstate.expect(H)
print(f"Optimized energy : {jas_energy}")
print(vstate.parameters)

# Define local operators for structure factor (magnetization)
structure_factor = nk.operator.LocalOperator(hi, dtype=complex)
for i in range(0, L):
  structure_factor += nk.operator.spin.sigmaz(hi, i)

# Run the optimization protocol
gs.run(out='test', n_iter=400000, obs={'Structure Factor': structure_factor})

# Load the data from the .log file
import json

data = json.load(open("test.log"))

# Check the structure of the data
print(data.keys())  # Print the keys to verify the structure

# Plot the magnetization vs iterations
iters = data['Energy']['iters']
energy = data['Energy']['Mean']
sf = data['Structure Factor']['Mean']

# Plot the energy and structure factor vs iterations
fig, ax1 = plt.subplots()
ax1.plot(iters, energy, color='blue', label='Energy')
ax1.set_ylabel('Energy')
ax1.set_xlabel('Iteration')

ax2 = ax1.twinx()
ax2.plot(iters, sf['real'], color='green', label='Structure Factor')
ax2.set_ylabel('Structure Factor')

ax1.legend(loc=2)
ax2.legend(loc=1)

plt.show()

print(f"Final optimized parameters: j1={vstate.parameters['j1']}, j2={vstate.parameters['j2']}")

[gcarleo]

On any finite system the magnetisation must be zero, it is non zero only in the thermodynamic limit. the way to check this is to look at correlators of the spin ( or the structure factor) but what you are computing is not a correlator, it's the sum of the spin sigmaz, which is then correctly giving you zero !

[gcarleo]

can you share the ED script used in that case ?

[mahmudams]

import netket as nk
import numpy as np
import matplotlib.pyplot as plt

def generate_basis_states(L):
    # Generate all possible basis states in binary representation
    return [format(i, f'0{L}b') for i in range(2**L)]

def calculate_magnetization(L, J, h_values):
    # Define the 1D chain graph
    g = nk.graph.Hypercube(length=L, n_dim=1, pbc=False)

    # Define the Hilbert space of spin-1/2 particles on the graph
    hi = nk.hilbert.Spin(s=0.5, N=g.n_nodes)

    magnetizations = []  # To store magnetization values for each h

    # Generate basis states as binary strings
    basis_states = generate_basis_states(L)

    # Loop over different values of the external magnetic field
    for h in h_values:
        # Define the Hamiltonian with transverse field and nearest-neighbor interaction
        H = sum([-h * nk.operator.spin.sigmax(hi, i) for i in range(L)])  # Transverse field
        H += sum([J * nk.operator.spin.sigmaz(hi, i) * nk.operator.spin.sigmaz(hi, (i + 1) % L) for i in range(L)])  # Nearest-neighbor interaction

        # Convert the Hamiltonian to dense matrix form
        H_dense = H.to_dense()

        # Find all eigenvalues and eigenvectors
        E, psi = np.linalg.eigh(H_dense)

        # Ground state energy and wavefunction
        ground_state_energy = E[0]
        ground_state_psi = psi[:, 0]

        # Magnetization calculation
        M = 0.0
        for i, bstate in enumerate(basis_states):
            # Convert basis state from binary string to integer
            bstate_int = int(bstate, 2)
            bstate_M = 0.0
            for bit in range(L):
                # Extract the spin value using bitwise operation
                spin = (bstate_int >> bit) & 1  # Get the bit at position bit
                spin_value = 1 if spin else -1  # Spin up (1) or down (-1)
                # Sum up contributions from each spin in the basis state
                bstate_M += (ground_state_psi[i] ** 2) * spin_value / L

            # Accumulate absolute magnetization for each basis state
            M += abs(bstate_M)

        magnetizations.append(M)

    return h_values, magnetizations

# Parameters for the simulation
L_values = [4, 6, 8, 10]  # List of lengths for the spin chain
J = -1.0                  # Interaction strength between nearest neighbors
h_values = np.logspace(-2, 2, 50)  # Magnetic field values from 0.01 to 100 in log scale

# Initialize plot
plt.figure(figsize=(12, 8))

# Loop over different values of L and calculate magnetizations
for L in L_values:
    h_vals, magnetizations = calculate_magnetization(L, J, h_values)
    plt.plot(h_vals, magnetizations, marker='o', linestyle='-', label=f'L={L}')

# Plot settings
plt.xscale('log')  # Logarithmic scale for the external magnetic field
plt.yscale('linear')  # Linear scale for the magnetization
plt.xlabel('External Magnetic Field (h)')
plt.ylabel('Magnetization |M|')
plt.title('Magnetization vs External Magnetic Field in the Transverse Field Ising Model')
plt.grid(True, which="both", ls="--")

# Check if h=J is within the plot range
if min(h_values) <= J <= max(h_values):
    plt.axvline(J, color='red', linestyle='--', label='Expected Transition Point h=J')

plt.legend()
plt.show()

[gcarleo]

the main difference is that in one case you are doing the absolute value of the magnetization, which indeed should reveal a transition (the ED case) whereas in the vmc case you do take the total magnetization, without the absolute value which indeed is supposed fo be zero. you can check that removing the abs in the ED code

[gcarleo]

the main limitation here is that you cannot take a function of a local operator (as far as I remember) so If you try to to take the abs of the magnetization in the Vmc code it will complain. this is in general because computing functions of operators is highly nontrivial, since they are generally off diagonal in the computational basis. A quick and dirty solution is to use the samples generated by the Variational State sampler and compute the abs of the magnetization manually on those samples + average them