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

From ALPS
Jump to: navigation, search

Waist correction

Starting the business

Bandstructure

Headers:

import numpy;
import pyalps;
import pyalps.dwa;

Try:

V0   = numpy.array([8.35, 8.35, 8.35]);
wlen = numpy.array([843., 843., 843.]);
a = 101
m = 86.99
L = 160

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

See:

band

Consider a parabolic trapping of 10.5 Hz:

omega = numpy.array([10.5, 10.5, 10.5]);
VT = pyalps.dwa.trappingFrequency(omega=omega, mass=m, wavelength=wlen) / band.t();

Parameters (no waist)

Parameters:

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

params=[]

for seed in range(10,210,10):
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : 2.          # Attempt to fix at N=280000 , lower bound: 2 , upper bound: 3        
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    #, 'V0'              : (8.35 * (154.89/4.38277)) 
    #, 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100
    , 'SEED'            : seed                 
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.nowaist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.nowaist");

Parameters (waist)

Parameters:

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

params=[]

for seed in range(10,210,10):
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : 2.          # Attempt to fix at N=280000 , lower bound: 2 , upper bound: 3        
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    , 'V0'              : (8.35 * (154.89/4.38277))
    , 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100
    , 'SEED'            : seed                 
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.waist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.waist");

A preliminary run

A trial run like this:

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

Automated processing (no waist)

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', tolerance=0.001, simplified=True)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -J jobno -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 jobno -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 jobno -n1 -W60' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);

Automated processing (waist)

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', tolerance=0.001, simplified=True)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -J jobyes -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 jobyes -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 jobyes -n1 -W60' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);

Tuning chemical potential

Parameters (no waist)

Headers:

import numpy;
import pyalps;
import pyalps.dwa;

Parameters:

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

params=[]

for mu in [2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9]:
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : mu          # Attempt to fix at N=280000 , lower bound: 2 , upper bound: 3        
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    #, 'V0'              : (8.35 * (154.89/4.38277)) 
    #, 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100                
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.nowaist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.nowaist");

Parameters (waist)

Headers:

import numpy;
import pyalps;
import pyalps.dwa;

Parameters:

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

params=[]

for mu in [2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9]:
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : mu          # Attempt to fix at N=280000 , lower bound: 2 , upper bound: 3        
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    , 'V0'              : (8.35 * (154.89/4.38277))
    , 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100               
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.waist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.waist");

Observing on your own

First the headers again:

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

We gather result files:

resultFiles = pyalps.getResultFiles(prefix='parm.nowaist');

or

resultFiles = pyalps.getResultFiles(prefix='parm.waist');

and collect data:

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.N.jpg');
fig.show();

or send it to yourself:

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

Unexpected norun

Sometimes, due to certain node failure, jobs may not even start. The following is the way to fix this:

import pyalps;
import pyalps.dwa;
all_h5_infiles = pyalps.getInputH5Files(prefix='parm.nowaist');
h5_infiles = [];
for h5_infile in all_h5_infiles:
   ar = pyalps.hdf5.h5ar(h5_infile);
   if not ar.is_group('/simulation/worldlines'):
      h5_infiles.append(h5_infile);
for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);

Checking thermalization

Gather result files:

resultFiles = pyalps.getResultFiles(prefix='parm.nowaist');

or:

resultFiles = pyalps.getResultFiles(prefix='parm.waist');

Check now:

for resultFile in resultFiles:
  mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
  N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
  thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', tolerance=0.001, simplified=True);
  print(resultFile + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

Checking convergence

Gather result files:

resultFiles = pyalps.getResultFiles(prefix='parm.nowaist');

or:

resultFiles = pyalps.getResultFiles(prefix='parm.waist');

Check now:

for resultFile in resultFiles:
  mu=pyalps.getParameters(resultFile)[0]['mu_homogeneous'];
  N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
  converged=pyalps.dwa.converged(resultFile, 'Total Particle Number', simplified=True);
  print(resultFile + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , converged : ' + str(converged));

Rethermalization

Gather the files:

h5_infiles = pyalps.getInputH5Files(prefix='parm.nowaist');

or:

h5_infiles = pyalps.getInputH5Files(prefix='parm.waist');

Copy worldlines over:

for h5_infile in h5_infiles:
  pyalps.dwa.extract_worldlines(pyalps.input2output(h5_infile), h5_infile);
  print('Done : ' + h5_infile);

Rethermalize:

for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);

Fine tuning the chemical potential

First the headers again:

import numpy;
import pyalps;
import pyalps.dwa;

Easily done like this:

resultFiles = pyalps.getResultFiles(prefix='parm.nowaist');

or

resultFiles = pyalps.getResultFiles(prefix='parm.waist');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', fixedY=280000);

Prepare new input files:

h5_infiles = pyalps.writeInputH5Files("parm.nowaist.1", params);

or

h5_infiles = pyalps.writeInputH5Files("parm.waist.1", params);

Run them:

for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);


Into the business

Parameters (no waist)

Parameters:

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

params=[]

for seed in range(10,210,10):
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : 2.54          # Attempt to fix at N=280000        
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    #, 'V0'              : (8.35 * (154.89/4.38277)) 
    #, 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100
    , 'SEED'            : seed                 
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.nowaist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.nowaist");

Parameters (waist)

Parameters:

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

params=[]

for seed in range(10,210,10):
  params.append(  
    { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
    , 'MODEL'           : 'boson Hubbard'                        
    , 'L'               : 160                             
    , 'Nmax'            : 20                 
    , 't'               : 1.                 
    , 'U'               : 8.1               
    , 'T'               : 1.                 
    , 'mu_homogeneous'  : 2.45          # Attempt to fix at N=280000         
    , 'mu'              : 'mu_homogeneous - (0.00369*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.00369*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.00369*(z-(L-1)/2.)*(z-(L-1)/2.))'
    , 'tof_phase'       : str(tof_phase)
    , 'V0'              : (8.35 * (154.89/4.38277))
    , 'waist'           : 355.87
    , 'SWEEPS'          : 100000             
    , 'SKIP'            : 100
    , 'SEED'            : seed                 
    }
  )

Write:

h5_infiles = pyalps.writeInputH5Files("parm.waist",params);

Or read:

h5_infiles = pyalps.getInputH5Files(prefix = "parm.waist");


Turn on respective measurements

Gather the files:

h5_infiles = pyalps.getInputH5Files(prefix='parm.nowaist')

or:

h5_infiles = pyalps.getInputH5Files(prefix='parm.waist')

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 (no waist)

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 jobno -n1 -W60 -r' ,
    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 jobno -n1 -W60 -r' ,
    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 jobno -n1 -W60 -r' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing(h5_infile);

Automated processing (waist)

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 jobyes -n1 -W60 -r' ,
    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 jobyes -n1 -W60 -r' ,
    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 jobyes -n1 -W60 -r' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing(h5_infile);


Local density

The headers:

import matplotlib; 
matplotlib.use('Agg');
import matplotlib.pyplot;
from mpl_toolkits.mplot3d.axes3d import Axes3D;
import numpy;
import pyalps;
import pyalps.dwa;
h5_infiles = pyalps.getInputH5Files(prefix='parm.nowaist');
h5_infile = h5_infiles[0]
L = int(pyalps.getParameters(pyalps.input2output(h5_infile))[0]['L']);
density = pyalps.getMeasurements(pyalps.input2output(h5_infile), observable='Local Density')[0]['mean']['value'].reshape([L,L,L]);
cross_section_density = density[:,:,L/2];
column_integrated_density = numpy.sum(density, axis=2);
x = numpy.array([range(L)] * L, float) - (L-1)/2.;
y = numpy.transpose(x);

To plot the cross section density along a line:

matplotlib.pyplot.figure();
matplotlib.pyplot.plot(x[L/2,:], cross_section_density[L/2,:]);
matplotlib.pyplot.show();

or column integrated density along a line:

matplotlib.pyplot.figure();
matplotlib.pyplot.plot(y[:,L/2], column_integrated_density[:,L/2]);
matplotlib.pyplot.show();

To see a cross section density on a plane:

fig = matplotlib.pyplot.figure();
ax = fig.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( x, y, cross_section_density, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=10);
fig.savefig(h5_infile + '.cross-section-density.jpg');
fig.show();

To see a column integrated density on a plane:

fig = matplotlib.pyplot.figure();
ax = fig.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( x, y, column_integrated_density, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=10);
fig.savefig(h5_infile + '.column-integrated-density.jpg');
fig.show();

Send the plot to yourself:

pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment=h5_infile +'.cross-section-density.jpg'
               , subject='Picture: ' + h5_infile + '.cross-section-density.jpg'
               );
pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment=h5_infile + '.column-integrated-density.jpg'
               , subject='Picture: ' + h5_infile + '.column-integrated-density.jpg'
               );


For this project

The headers:

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

This:

h5_infiles = pyalps.getInputH5Files(prefix='parm.nowaist');

or this:

h5_infiles = pyalps.getInputH5Files(prefix='parm.waist');
obs = []; 
for h5_infile in h5_infiles:
  L = int(pyalps.getParameters(pyalps.input2output(h5_infile))[0]['L']);
  density = pyalps.getMeasurements(pyalps.input2output(h5_infile), observable='Local Density')[0]['mean']['value'].reshape([L,L,L]);
  cross_section_1d_density = (density[:,L/2,L/2] + density[L/2,:,L/2] + density[L/2,L/2,:])/3.;
  cross_section_1d_density = (cross_section_1d_density + cross_section_1d_density[::-1])/2.;
  x = numpy.array(range(L), float) - (L-1)/2.;
  obs.append({ 'x' : x , 'obs' : cross_section_1d_density});
  print('Done : ' + h5_infile);
x = obs[0]['x'];
obs = numpy.array([item['obs'] for item in obs]);
value = obs.mean(axis=0);
error = obs.std(axis=0)/math.sqrt(obs.shape[0]);

Data

Saving:

numpy.savetxt('parm.nowaist.cross-section-1d-density.dat', numpy.array([x,value,error]).transpose())

or:

numpy.savetxt('parm.waist.cross-section-1d-density.dat', numpy.array([x,value,error]).transpose())

Email:

Send the plot to yourself:

pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment='parm.nowaist.cross-section-1d-density.dat'
               , subject='Picture: parm.nowaist.cross-section-1d-density.dat'
               );

or

pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment='parm.waist.cross-section-1d-density.dat'
               , subject='Picture: parm.nowaist.cross-section-1d-density.dat'
               );

Plotting figure

To plot the cross section density along a line:

fig = matplotlib.pyplot.figure();
matplotlib.pyplot.errorbar(x, value, error);
fig.savefig('parm.nowaist.cross-section-1d-density.jpg');

or

fig.savefig('parm.waist.cross-section-1d-density.jpg');

Send the plot to yourself:

pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment='parm.nowaist.cross-section-1d-density.jpg'
               , subject='Picture: parm.nowaist.cross-section-1d-density.jpg'
               );

or

pyalps.sendmail( 'pingnang@phys.ethz.ch'    # email address of recipients 
               , attachment='parm.waist.cross-section-1d-density.jpg'
               , subject='Picture: parm.nowaist.cross-section-1d-density.jpg'
               );