# ALPS 2 Tutorial:MC-09 Directed Worm Algorithm:Intensive Example 2

Jump to: navigation, search

# Description

This example is about error budgeting.

This section is under construction.

# Topic 1: Error budgeting in total particle number

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the total particle number.

## A first reference to start the business

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Parameters:

```tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)
```
```params=[]
```
```params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 't'               : 1.
, 'U'               : 8.11
, 'T'               : 1.
, 'mu_homogeneous'  : 4.05
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(tof_phase)
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)
```

Write:

```h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.N",params);
```

Or read:

```h5_infiles = pyalps.getInputH5Files(prefix = "parm.errorbudget.N");
```

A trial run of course:

```pyalps.runApplication('dwa', h5_infiles[0]);
```

## Pinning chemical potential

### 5% fluctuation in chemical potential?

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Write from scratch:

```tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)

params=[]
for ratio in [0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05]:
params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 't'               : 1.
, 'U'               : 8.11
, 'T'               : 1.
, 'mu_homogeneous'  : ratio * 4.05
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(tof_phase)
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)

h5_infiles = pyalps.writeInputH5Files('parm.errorbudget.N',params);
```

Read from archive:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.N');
```

Check for yourself, for instance:

```pyalps.getParameters(h5_infiles[0]);
```

### Automated processing

```def parallel_automated_processing(taskfile):
pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
batch_run_script = taskfile + '.thermalize.script' ,
batch_next_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
batch_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.startRunScript(
batch_run_script= taskfile + '.thermalize.script' ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

### Observing on your own

First the headers again:

```import matplotlib;
matplotlib.use('Agg');
```
```import pyalps;
import pyalps.plot;
import pyalps.dwa;
```

We collect useful data:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number')
```

then plot it:

```fig = matplotlib.pyplot.figure();
pyalps.plot.plot(N);
fig.savefig('parm.errorbudget.N.N.jpg');
fig.show();
```

or send it to yourself:

```pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients
, attachment='parm.errorbudget.N.N.jpg'
, subject='Picture: parm.errorbudget.N.N.jpg'
);
```

### Correcting chemical potential for fixed N

First the headers again:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Easily done like this:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', fixedY=280000*numpy.array([0.96, 0.97, 0.98, 0.99, 1., 1.01, 1.02, 1.03, 1.04]));
```

Prepare new input files:

```h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.N.1", params);
```

Run them:

```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

That's really convenient, don't you agree?

## Visibility

### Turn on respective measurements

Gather the files:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.N.1')
```

Switch on local density measurements:

```for h5_infile in h5_infiles:
pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Local Density]', 1);
pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Local Density]', 1);
```

Switch on green function measurements:

```for h5_infile in h5_infiles:
pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Green Function]', 1);
pyalps.dwa.switchParameter(pyalps.input2output(h5_infile), 'MEASURE[Green Function]', 1);
```

### Automated processing

```def parallel_automatic_processing(taskfile):
pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
batch_run_script = taskfile + '.measure1.script' ,
batch_next_run_script = taskfile + '.measure.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
n = 20 ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
batch_run_script = taskfile + '.measure.script' ,
batch_next_run_script = taskfile + '.sendmail.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has finished measurements for 20 iterations.')" ,
cmd_lang = 'python' ,
loc = locals() ,
batch_submit = True ,
batch_run_script = taskfile + '.sendmail.script' ,
batch_noRun = True
);

pyalps.dwa.startRunScript(
batch_run_script= taskfile + '.measure1.script' ,
batch_cmd_prefix = 'bsub -J jobN -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automatic_processing(h5_infile);
```

### Extracting final data(plot)

First, the headers again:

```import matplotlib;
matplotlib.use('Agg');
```
```import matplotlib.pyplot;
from mpl_toolkits.mplot3d.axes3d import Axes3D;
import numpy;
import numpy.fft;
import pyalps;
import pyalps.dwa;
```

Then, evaluate visibility:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N.1');

for resultFile in [resultFiles[0]]:
L = pyalps.getParameters(pyalps.input2output(h5_infile))[0]['L'];
green_tof = pyalps.getMeasurements(pyalps.input2output(h5_infile), observable='Green Function:TOF')[0]['mean']['value'].reshape([L,L,L]);
momentum_density = numpy.fft.fftn(green_tof).real;

V0   = numpy.array([8.805, 8. ,  8. ]);
wlen = numpy.array([765., 843., 843.]);
a    = 101;
m    = 86.99;

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
```

This section is under construction.

# Topic 2: Error budgeting in scattering length

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the scattering length.

## Varying scattering length at fixed chemical potential

Let's first fix the chemical potential while varying scattering length, and see what happens out of it.

### Parameters

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Write from scratch:

```tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)

ps = [];
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
V0   = numpy.array([8.805, 8. ,  8. ]);
wlen = numpy.array([765., 843., 843.]);
a    = ratio*101;
m    = 86.99;
L    = 160;

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
p = {};
p.update({'ratio_a' : ratio , 'U' : band.Ut()[2]})
ps.append(p);

params=[]
for p in ps:
params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 'ratio_a'         : p['ratio_a']
, 't'               : 1.
, 'U'               : p['U']
, 'T'               : 1.
, 'mu_homogeneous'  : 4.0265
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(tof_phase)
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)

h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a",params);
```

Read from archive:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a');
```

Check for yourself, for instance:

```pyalps.getParameters(h5_infiles[4]);
```

Have a first taste of simulation:

```pyalps.runApplication('dwa', h5_infiles[4], T=2700);
```

### Automated processing

```def parallel_automated_processing(taskfile):
pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
batch_run_script = taskfile + '.thermalize.script' ,
batch_next_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
batch_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.startRunScript(
batch_run_script= taskfile + '.thermalize.script' ,
batch_cmd_prefix = 'bsub -J joba -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

### Observing on your own

First the headers again:

```import matplotlib;
matplotlib.use('Agg');
```
```import pyalps;
import pyalps.plot;
import pyalps.dwa;
```

We collect useful data:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
N = pyalps.collectXY(data, x='ratio_a', y='Total Particle Number')
```

then plot it:

```fig = matplotlib.pyplot.figure();
pyalps.plot.plot(N);
fig.savefig('parm.errorbudget.a.N.jpg');
fig.show();
```

or send it to yourself:

```pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients
, attachment='parm.errorbudget.a.N.jpg'
, subject='Picture: parm.errorbudget.a.N.jpg'
);
```

## Varying scattering length at fixed N

### Parameters

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Write from scratch:

```tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)

ps = [];
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
V0   = numpy.array([8.805, 8. ,  8. ]);
wlen = numpy.array([765., 843., 843.]);
a    = ratio*101;
m    = 86.99;
L    = 160;

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
p = {};
p.update({'ratio_a' : ratio , 'U' : band.Ut()[2]})
ps.append(p);

params=[]
for p in ps:
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 'ratio_a'         : p['ratio_a']
, 't'               : 1.
, 'U'               : p['U']
, 'T'               : 1.
, 'mu_homogeneous'  : ratio*4.0265
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(tof_phase)
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)

h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a",params);
```

Read from archive:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a');
```

### Observing on your own

First the headers again:

```import matplotlib;
matplotlib.use('Agg');
```
```import pyalps;
import pyalps.plot;
import pyalps.dwa;
```

We collect useful data:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_a'])
```

then plot it:

```fig = matplotlib.pyplot.figure();
pyalps.plot.plot(N);
fig.savefig('parm.errorbudget.a.N.jpg');
fig.show();
```

or send it to yourself:

```pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients
, attachment='parm.errorbudget.a.N.jpg'
, subject='Picture: parm.errorbudget.a.N.jpg'
);
```

### Fine tuning chemical potential

First the headers again:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Easily done like this:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_a'], fixedY=280000);
```

Prepare new input files:

```h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.a.1", params);
```

Run them:

```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

That's really convenient, don't you agree?

This section is under construction.

# Topic 3: Error budgeting in lattice potential

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the lattice potential along the x-direction.

### Parameters

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Write from scratch:

```tof_phase = pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[765,843,843], mass=86.99)

ps = [];
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
V0   = numpy.array([ratio*8.805, 8. ,  8. ]);
wlen = numpy.array([765., 843., 843.]);
a    = 101;
m    = 86.99;
L    = 160;

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
p = {};
p.update({'ratio_V0x' : ratio , 'U' : band.Ut()[2] , 'tx_t' : band.t()[0]/band.t()[2]})
ps.append(p);

params=[]
for p in ps:
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 'ratio_V0x'       : p['ratio_V0x']
, 't'               : 1.
, 'tx_t'            : p['tx_t']
, 'U'               : p['U']
, 'T'               : 1.
, 'mu_homogeneous'  : ratio*4.0265
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(tof_phase)
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)

h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.V0",params);
```

Read from archive:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
```

### Automated processing

```def parallel_automated_processing(taskfile):
pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
batch_run_script = taskfile + '.thermalize.script' ,
batch_next_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
batch_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.startRunScript(
batch_run_script= taskfile + '.thermalize.script' ,
batch_cmd_prefix = 'bsub -J jobV0 -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

### Observing on your own

First the headers again:

```import matplotlib;
matplotlib.use('Agg');
```
```import pyalps;
import pyalps.plot;
import pyalps.dwa;
```

We collect useful data:

```resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.V0');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_V0x'])
```

then plot it:

```fig = matplotlib.pyplot.figure();
pyalps.plot.plot(N);
fig.savefig('parm.errorbudget.V0.N.jpg');
fig.show();
```

or send it to yourself:

```pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients
, attachment='parm.errorbudget.V0.N.jpg'
, subject='Picture: parm.errorbudget.V0.N.jpg'
);
```

### Fine tuning the chemical potential

This section is under construction

# Topic 4: Error budgeting in laser wavelength

We want to investigate how the condensate fraction and visibility varies with up to 5% fluctuation in the laser wavelength along the x-direction.

### Parameters

The usual headers:

```import numpy;
import pyalps;
import pyalps.dwa;
```

Write from scratch:

```ps = [];
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
V0   = numpy.array([ratio*ratio*8.805, 8. ,  8. ]);
wlen = numpy.array([ratio*765., 843., 843.]);
a    = 101;
m    = 86.99;
L    = 160;

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
p = {};
p.update({'ratio_wlenx' : ratio , 'U' : band.Ut()[2] , 'tx_t' : band.t()[0]/band.t()[2] , 'tof_phase' : pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[ratio*765,843,843], mass=86.99) })
ps.append(p);

params=[]
for p in ps:
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
params.append(
{ 'LATTICE'         : 'inhomogeneous simple cubic lattice'
, 'MODEL'           : 'boson Hubbard'
, 'L'               : 160
, 'Nmax'            : 20
, 'ratio_wlenx'     : p['ratio_wlenx']
, 't'               : 1.
, 'tx_t'            : p['tx_t']
, 'U'               : p['U']
, 'T'               : 1.
, 'mu_homogeneous'  : ratio*4.0265
, 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
, 'tof_phase'       : str(p['tof_phase'])
, 'SWEEPS'          : 100000
, 'SKIP'            : 100
}
)

h5_infiles = pyalps.writeInputH5Files("parm.errorbudget.wlen",params);
```

Read from archive:

```h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
```

### Automated processing

```def parallel_automated_processing(taskfile):
pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', taskfile, T=2700)" ,  cmd_lang = 'python' ,
follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
batch_run_script = taskfile + '.thermalize.script' ,
batch_next_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.recursiveRun(
"pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=2700)" ,  cmd_lang = 'python' ,
break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
loc = locals() ,
batch_submit = True ,
batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
batch_run_script = taskfile + '.converge.script' ,
batch_noRun = True
);

pyalps.dwa.startRunScript(
batch_run_script= taskfile + '.thermalize.script' ,
batch_cmd_prefix = 'bsub -J jobwlen -n1 -W60' ,
loc = locals()
);
```
```for h5_infile in h5_infiles:
parallel_automated_processing(h5_infile);
```

### Fine tuning the chemical potential

This section is under construction