All 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-example

Step 2: Editing files for the static calculations

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 form of logger.

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 1d

If 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.