Gabriel WlazłowskiAll the files (inputs and outputs) for this example are accessible via Zenodo.
Target
Simulate quasi-1D Josephson dynamics. Initialize the simulation of a unitary Fermi gas confined in an external potential
V^{(ext)}_{\sigma}(x) = \frac{1}{2}\omega_x^2 x^2 + V_0e^{-2x^2/w^2} + k x,where the first term corresponds to harmonic confinement, the second describes the potential barrier separating the cloud into two reservoirs, and the last introduces a small linear tilt to generate a population imbalance between the left and right reservoirs.
Step 1: Creating working folders
We start by copying template folders for static (st) and time-dependent (td) simulations:
cp -r $WSLDA/st-project-template ./st-jj-example
cp -r $WSLDA/td-project-template ./td-jj-exampleStep 2: Editing files for the static calculations
All edits should be done in the working folder st-jj-example.
Edit predefines.h
Select the lattice size and functional, and inspect all other tags that can be relevant for your problem:
/**
* 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
// ...Edit problem-definition.h
In this file, you provide a formula for the external potential. We will parametrize the function so it can be controlled via an input file.
#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;
}
// ...It is more convenient to provide parameter values that define the potential in human-readable units. For example, rather than specifying the trapping frequency \omega_x directly, user provides the desired Thomas-Fermi radius x_{TF}, which is related to the chemical potential through
\frac{1}{2}\omega_x^2 x_{TF}^2 = \mu
\quad \Rightarrow \quad
\omega_x = \frac{\sqrt{2\mu}}{x_{TF}}.Another example concerns the barrier height V_0 and width w, which are more conveniently expressed in terms of the Fermi energy \varepsilon_F and Fermi wave vector k_F. To do a conversion from human-readable units to code units, we can use the process_params(...) function.
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 logger.h
In this file user can customize the output. In this example, we use the default logger form.
Edit input.txt
The input file provides run-time parameters. The most important in this case are the parameters that define the potential
# -------------- 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 are
# ...
# ------------------- 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: Compiling and running the static code
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: Analyzing the results of the static run
As a simple example, we will plot the density profile and compute the population imbalance, defined as z=(N_L-N_R)/(N_L+N_R), where N_L and N_R are particle numbers in left and right reservoirs, using Python
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: Editing files for the time-dependent calculations
All edits should be done in the working folder td-jj-example.
The workflow is analogous to that for static calculations. Many of the functions can be copied & pasted from static files.
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 logger.h
In this file user can customize the output. In this example, we use the default logger form.
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: Compiling and running the time-dependent code
We do this part in an analogous way to the case of static calculations. Note that the time-dependent run requires GPUs to run.
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: Analyzing the results of the time-dependent run
As an example of analysis, we will plot population imbalance and phase differnce accors the junction as a function of time, using simple python script
# todo