Gabriel WlazłowskiAll input and output files for this example are available on Zenodo.
- Target
- Step 1: Create Working Directories
- Step 2: Static Calculation Setup (st-jj-example)
- Step 3: Compile & Run (Static)
- Step 4: Analyze Static Results
- Step 5: Time-Dependent Setup (td-jj-example)
- Step 6: Compile & Run (Time-Dependent)
- Step 7: Analyze Time Evolution
Target
Simulate quasi-1D Josephson dynamics in a unitary Fermi gas confined by
V^{(ext)}_{\sigma}(x) = \frac{1}{2}\omega_x^2 x^2 + V_0e^{-2x^2/w^2} + k x,- Harmonic trap: confinement.
- Gaussian barrier: splits the cloud into two reservoirs.
- Linear tilt: generates an initial population imbalance.
Step 1: Create Working Directories
Copy the project templates:
cp -r $WSLDA/st-project-template ./st-jj-example
cp -r $WSLDA/td-project-template ./td-jj-example
Step 2: Static Calculation Setup (st-jj-example)
Edit predefines.h
Define lattice and functional:
/**
* Define lattice size and lattice spacing.
* */
#define NX (256)
#define NY (16)
#define NZ (16)
// ...
#define FUNCTIONAL SLDAE
/**
* Sets the effective mass to be equal to 1.
* */
#define SLDA_FORCE_A1
// ...Adjust other tags if needed.
Edit problem-definition.h
Define the external potential:
#include "../extensions/wslda_utils.h"
// ...
double v_ext(int ix, int iy, int iz, int it, int spin, double *params, size_t extra_data_size, void *extra_data)
{
double x = DX*(ix-NX/2);
double V_ho = harmonic_oscillator_smooth_edges(x, params[11], params[2]*LX*0.5, params[3]*LX*0.5, 1.0);
double V_barrier = params[4] * exp(-2.0*pow(x/params[5],2));
double V_tilt = params[6]*x/(LX/2); // linear tilt potential
return V_ho + V_barrier + V_tilt;
}
// ...Human-Readable Parameters
Instead of specifying \omega_x directly, define the Thomas–Fermi radius:
\frac{1}{2}\omega_x^2 x_{TF}^2 = \mu
\quad \Rightarrow \quad
\omega_x = \frac{\sqrt{2\mu}}{x_{TF}}.Barrier parameters are given in units of \varepsilon_F and k_F.
Conversion is handled in process_params(...):
void process_params(double *params, double *kF, double *mu, size_t extra_data_size, void *extra_data)
{
// if akF is active, then overwrite the value of input->sclgth
// value of kF is generated by the function referencekF(...)
if(input->akF!=0.0) input->sclgth = input->akF/kF[0];
// PROCESS INPUT FILE PARAMETERS
// Definition of R_TF: 0.5*m*omega^2 * R_TF^2 = mu
// => omega = sqrt(2*mu/m)/R_TF
params[11] = sqrt(2.*mu[SPINA])/ (params[1]*LX*0.5); // omega_x
double eF = kF[0]*kF[0]/2.0;
params[4] = params[4]*eF; // convert barrier height to internal units
params[5] = params[5]/kF[0]; // convert barrier width to internal units
}Edit input.txt
Key user parameters:
# -------------- USER DEFINED PARAMETERS --------------
# Data flow: [Read params from input file] -> [execute process_params( )]
# -> [pass params to functions]
params[1] = 0.80 # Thomas-Fermi radius for x direction, in units of LX/2
params[2] = 0.90 # Smoothing parameter x1 for harmonic_oscillator_smooth_edges
params[3] = 0.97 # Smoothing parameter x2 for harmonic_oscillator_smooth_edges
params[4] = 0.24 # height of the barrier, in eF units
params[5] = 5.0 # width of barrier, in 1/kF units
params[6] = 0.015 # tilt potential (value at the edge of the trap)
# ...Other input parameters relevant for this run:
# ...
# ------------------- INITIALIZATION ------------------
inittype 0 # create uniform solution and start from it,
# ...
# ------------------- INPUT/OUTPUT ------------------
outprefix jj-r1 # all output files will start with this prefix
writewf 1 # write wf at the end of computation: yes=1, no=0
# we will evolve them by td-wslda-1d code
# variables to write
writevar rho delta j V_ext
sclgth -10000.0 # scattering length in units of lattice spacing
referencekF 1.0 # in this case, I can define the value of kF=1
# ...
# ----------------- SELF-CONSISTENT LOOP -----------------
energyconveps 1.0e-6 # energy convergence epsilon
npartconveps 1.0e+6 # large number means that this condition will always be satisfied
muchange 0.0 # do not change chemical potential,
# so computation will be for fixed chemical potential
maxiters 200 # maximum number of iterations
spinsymmetry 1 # to impose Na=Nb
# ...
# ----------- IN CASE OF: inittype=0 ----------------
# See: Wiki -> Initialization of the solver
init0na 0.01689 # required density for component a
init0nb 0.01689 # and component b, so corresponding kF=1
init0muchange 0.5 # change rate of chemical potentials, default muchange
init0eps 1.0e-9 # epsilon for convergence for uniform solution
init0maxiter 1000 # maximum number of iterations when searching for a uniform solution
# ...Step 3: Compile & Run (Static)
To generate the binary for a quasi-1D simulation, execute the command:
make 1dIf the compilation completes successfully, a binary named st-wslda-1d will be created.
This program should be executed within an MPI environment.
mpirun -np 4 ./st-wslda-1d input.txt
# START OF THE MAIN FUNCTION
# CODE: ST-WSLDA-1D
...
# ALGORITHM CONVERGED!
...
# EXTRA SAVING ITERATION DONE.See here for full stdout. The file produced by the logger is here.
Step 4: Analyze Static Results
Example: density profile and population imbalance z=(N_L-N_R)/(N_L+N_R), where N_L and N_R are particle numbers in left and right reservoirs
import numpy as np
import matplotlib.pyplot as plt
# to get wdata lib use: pip install wdata
from wdata.io import WData, Var
# Load data
data = WData.load("./st-jj-example/jj-r1.wtxt")
x=data.xyz[0]
rho=data.rho_a[-1,:]+data.rho_b[-1,:] # take last iteration
# compute initial population imbalance
NL = np.sum(rho[x<0])
NR = np.sum(rho[x>0])
z=(NL-NR)/(NL+NR) # initial population imbalance
print(z)
# Plot rho(x)
plt.plot(x,rho)
plt.show()The result is
For a more advanced script for plotting static results, see Zenodo.
Step 5: Time-Dependent Setup (td-jj-example)
Workflow is analogous.
Edit predefines.h
Set compilation flags to be compatible with your problem and with data produced by static calcs.
Edit problem-definition.h
Edit the function that defines the external potential. You can also add time dependence to your external potential function if you wish.
Edit input.txt
In the input file, the key parameters for this example are:
# -------------- USER DEFINED PARAMETERS --------------
# ...
params[6] = 0.0 # no tilt
# ...
# ------------------- INITIALIZATION ------------------
inittype 1 # start from st-wslda-1d solution,
# inprefix points to data from static calculations
# ...
# ------------------- INPUT/OUTPUT ------------------
outprefix jj-r2 # all output files will start with this prefix
inprefix ../st-jj-example/jj-r1 # prefix of static result
# variables to write
writevar rho delta j V_ext
sclgth -10000.0 # scattering length in units of lattice spacing
# ...
# ------------------- TIME EVOLUTION ----------------
dt 0.0035 # integration time step, in units of eF^-1
measurements 200 # number of requested measurements
timesteps 1430 # number of time steps between each measurement
# measurements are written with time resolution: timesteps*dt
# the total trajectory length is: measurements*timesteps*dt
selfstart 1 # use Taylor expansion for first 5 steps?
# ...Step 6: Compile & Run (Time-Dependent)
GPU resources required.
To generate the binary for a quasi-1D simulation, execute the command:
make 1dIf the compilation completes successfully, a binary named td-wslda-1d will be created.
This program should be executed within an MPI environment.
mpirun -np 1 ./td-wslda-1d input.txt
# START OF THE MAIN FUNCTION
# CODE: TD-WSLDA-1D
...
# REFERENCE VALUES: kF=1.000000, eF=0.500000, Effg=300.968458
time*eF Na Nb Na+Nb ETOT EKIN EPOT EPAIR ECURRENT EPOTEXT EPAIREXT EVELEXT rt
0.0000 501.61409628 501.61409628 1003.22819255 0.42101872 1.56019953 -0.38908978 -0.87549813 0.00000000 0.12540710 0.00000000 0.00000000
# SELFSTART: EXECUTING TAYLOR EXPANSION OF THE EVOLUTION OPERATOR.
# SELFSTART: DONE.
0.0000 501.61409627 501.61409627 1003.22819253 0.42104851 1.56019923 -0.38908978 -0.87546805 0.00000000 0.12540710 0.00000000 0.00000000
...
1001.0000 501.61409631 501.61409631 1003.22819262 0.42108490 1.55974541 -0.38904358 -0.87506952 0.00025957 0.12519302 0.00000000 0.00000000 118.64
# CHECKPOINT INFO: MODE=WRITE: DATA SIZE= 0.26 GB
# CHECKPOINT INFO: MPI_NP_PER_IO_GROUP=24.
# CHECKPOINT INFO: WRITE TIME= 1.74 sec
# CHECKPOINT INFO: WRITE SPEED= 0.148 GB/sec
# CREATING CHECK STAMP: `jj-r2_check.stamp`See here for full stdout. The file produced by the logger is here.
Step 7: Analyze Time Evolution
Plot population imbalance and phase difference as a function of time.
import numpy as np
import matplotlib.pyplot as plt
# to get wdata lib use: pip install wdata
from wdata.io import WData, Var
# Load data
data = WData.load("./td-jj-example/jj-r2.wtxt")
x=data.xyz[0]
rho=data.rho_a+data.rho_b # total density for all time frames..
delta=data.delta/data.eF # normalized order parameter for all time frames..
# compute initial population imbalance for each time frame...
NL = np.sum(rho[:,x<0], axis=1)
NR = np.sum(rho[:,x>0], axis=1)
z=(NL-NR)/(NL+NR) # population imbalance
# ... and phase difference
phase_diff = np.array([np.angle(delta[frame, 128 + 10] / delta[frame, 128 - 10]) / np.pi for frame in range(delta.shape[0])])
# plot population imbalance and phase difference as a function of time
# make two panels
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))
ax1.plot(z, color='blue')
ax2.plot(phase_diff, color='red')
plt.show()
For a more advanced script for plotting time-dependent results, see Zenodo.