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

From ALPS
Jump to: navigation, search

Description

This example is about experimental validation.

Example parameters

The usual headers:

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                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 160                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 8.11               
  , 'T'               : 2.5                 
  , 'mu_homogeneous'  : 3.92               
  , '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'          : 1000000             
  , 'SKIP'            : 1000                
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 160                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 8.11               
  , 'T'               : 4.                
  , 'mu_homogeneous'  : 3.00               
  , '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'          : 10000000             
  , 'SKIP'            : 10000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 160                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 8.11               
  , 'T'               : 5.6                 
  , 'mu_homogeneous'  : 1.675               
  , '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'          : 1000000000             
  , 'SKIP'            : 1000000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               :  200                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 8.11               
  , 'T'               : 6.7                 
  , 'mu_homogeneous'  : -0.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'          : 1000000000             
  , 'SKIP'            : 1000000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 200                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 8.11               
  , 'T'               : 10.                 
  , 'mu_homogeneous'  : -6.60                
  , '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'          : 1000000000             
  , 'SKIP'            : 1000000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 120                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 1.                 
  , 'mu_homogeneous'  : 11.72                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 1000000             
  , 'SKIP'            : 1000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 120                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 2.25                 
  , 'mu_homogeneous'  : 11.57                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 10000000             
  , 'SKIP'            : 10000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 120                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 2.63                 
  , 'mu_homogeneous'  : 11.476                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 10000000             
  , 'SKIP'            : 10000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 120                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 3.33                 
  , 'mu_homogeneous'  : 11.15                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 100000000             
  , 'SKIP'            : 100000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 160                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 5.                 
  , 'mu_homogeneous'  : 9.95                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 100000000             
  , 'SKIP'            : 100000                 
  }
)
params.append(
  { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
  , 'MODEL'           : 'boson Hubbard'                        
  , 'L'               : 200                             
  , 'Nmax'            : 20                 
  , 't'               : 1.                 
  , 'U'               : 27.5               
  , 'T'               : 7.69                 
  , 'mu_homogeneous'  : 7.2                
  , 'mu'              : 'mu_homogeneous - (0.02395*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.01248*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.01326*(z-(L-1)/2.)*(z-(L-1)/2.))'
  , 'tof_phase'       : str(tof_phase) 
  , 'SWEEPS'          : 200000000             
  , 'SKIP'            : 100000                 
  }
)

Write parameters into file:

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

Retrieve them if previously exist:

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

Have a preliminary taste:

pyalps.runApplication('dwa', h5_infiles);

Band structure calculations

First, the headers:

import numpy;
import pyalps.dwa;

There are 2 sets of experiments being mimicked here:

a. U/t = 8.11

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

b. U/t = 27.5

V0   = numpy.array([12.64 , 11.75 , 11.75]);
wlen = numpy.array([765., 843., 843.]);
a    = 101;
m    = 86.99;
L    = 200;

Finally:

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

Simulation

Thermalization

Laptop mode

Thermalizing task 0 on your laptop for instance:

pyalps.dwa.recursiveRun( 
  "pyalps.runApplication('dwa', h5_infiles[0], T=2700)" ,  cmd_lang = 'python' ,
  follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(h5_infiles[0]), outfile=h5_infiles[0])" ,
  break_if = "pyalps.dwa.thermalized(pyalps.input2output(h5_infiles[0]), 'Total Particle Number', simplified=True)" ,
  write_status = "pyalps.dwa.write_status(pyalps.input2output(h5_infiles[0]), 'Thermalizing')" ,
  loc = locals()
);

Cluster mode

Complete thermalization in the cluster:

def parallel_automatic_processing_thermalize(taskfile):
  pyalps.dwa.switchParameter(taskfile, 'MEASURE', 0);
  
  pyalps.dwa.recursiveRun(
    "pyalps.runApplication('dwa', taskfile, T=2700)" , cmd_lang = 'python' ,
    follow_up_script = "pyalps.dwa.extract_worldlines(pyalps.input2output(taskfile), taskfile)" ,
    end_script = "pyalps.dwa.switchParameter(taskfile, 'MEASURE', 1)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Speed thermalizing')" ,
    loc = locals() ,
    batch_submit = True , 
    batch_cmd_prefix = 'bsub -n1 -W60' ,
    batch_run_script = taskfile + '.speed.script' ,
    batch_next_run_script = taskfile + '.thermalize1.script' ,
    batch_noRun = True
  )
  
  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)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing-1')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    batch_run_script = taskfile + '.thermalize1.script' ,
    batch_next_run_script = taskfile + '.thermalize2.script' ,
    batch_noRun = True
  );
  
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', taskfile, T=27200)" ,  cmd_lang = 'python' ,
    follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
    end_script = "pyalps.dwa.write_benchmark(pyalps.input2output(taskfile), 'Total Particle Number')" ,
    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-2')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    batch_run_script = taskfile + '.thermalize2.script' ,
    batch_next_run_script = taskfile + '.sendmail.script' ,
    batch_noRun = True
  );
  
  pyalps.dwa.recursiveRun(
     "pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has thermalized.')" ,
     cmd_lang = 'python' ,
     loc = locals() ,
     batch_submit = True ,
     batch_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
  ); 
  
  pyalps.dwa.startRunScript(
    batch_run_script= taskfile + '.speed.script' ,
    batch_cmd_prefix = 'bsub -n1 -W60' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing_thermalize(h5_infile);

Checking for errors

Always do checks every now and then, for instance:

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

If you need rethermalization

If you are forced to kill the jobs before completion, you can start the script from somewhere, say thermalize2.script like the following:

def parallel_automatic_processing_rethermalize(taskfile):
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', taskfile, T=27200)" ,  cmd_lang = 'python' ,
    follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
    end_script = "pyalps.dwa.write_benchmark(pyalps.input2output(taskfile), 'Total Particle Number')" ,
    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-2')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    batch_run_script = taskfile + '.thermalize2.script' ,
    batch_next_run_script = taskfile + '.sendmail.script' ,
    batch_noRun = True
  );
  
  pyalps.dwa.recursiveRun(
     "pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has rethermalized.')" ,
     cmd_lang = 'python' ,
     loc = locals() ,
     batch_submit = True ,
     batch_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
  ); 
  
  pyalps.dwa.startRunScript(
    batch_run_script= taskfile + '.thermalize2.script' ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing_rethermalize(h5_infile);

Measurements

Add in any new measurements

Gather the files:

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

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);

Switch off local density measurements:

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

Switch off green function measurements:

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

Laptop mode

On your laptop:

pyalps.runApplication('dwa', h5_infiles[0] , T=2700);
pyalps.dwa.recursiveRun( 
  "pyalps.runApplication('dwa', pyalps.input2output(h5_infiles[0]), T=2700)" ,  cmd_lang = 'python' ,
  n = 10 ,
  write_status = "pyalps.dwa.write_status(pyalps.input2output(h5_infiles[0]), 'Converging')" , 
  loc = locals()
);

Cluster mode (1st submission)

def parallel_automatic_processing_24hr(taskfile):
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', taskfile, T=27200)" ,  cmd_lang = 'python' ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    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=27200)" ,  cmd_lang = 'python' ,
    n = 2 ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    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 24 hr.')" ,
     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 -n1 -W480' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing_24hr(h5_infile);

Cluster mode (resubmission)

def parallel_automatic_processing_resubmit(taskfile):
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=27200)" ,  cmd_lang = 'python' ,
    n = 5 ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Measuring')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    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 another 5 cycles.')" ,
     cmd_lang = 'python' ,
     loc = locals() ,
     batch_submit = True ,
     batch_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
  ); 
  
  pyalps.dwa.startRunScript(
    batch_run_script= taskfile + '.measure.script' ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing_resubmit(h5_infile);

Converge data further

List of physical observables to be converged:

obs = [ 'Green Function:0'
      , 'Momentum Distribution:0'
      , 'Momentum Distribution:TOF:0'
      ]

Check the autocorrelations:

for h5_infile in h5_infiles:
  print(pyalps.input2output(h5_infile) + ' : ' + str( pyalps.dwa.tau( pyalps.input2output(h5_infile) , obs ) ));

You may either increase skip:

pyalps.dwa.increase_skip('parm.validation.task1.in.h5');

or decrease skip:

pyalps.dwa.decrease_skip('parm.validation.task1.in.h5');

Now, converge data:

def parallel_automatic_processing_converge(taskfile):
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', pyalps.input2output(taskfile), T=27200)" ,  cmd_lang = 'python' ,
    break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile) , [ 'Green Function:0' , 'Momentum Distribution:0' , 'Momentum Distribution:TOF:0' ] , simplified=True)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
    loc = locals() ,
    batch_submit = True ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    batch_run_script = taskfile + '.converge.script' ,
    batch_next_run_script = taskfile + '.sendmail.script' ,
    batch_noRun = True
  ); 
  
  pyalps.dwa.recursiveRun(
     "pyalps.sendmail('pingnang@phys.ethz.ch', message= taskfile + ' has converged.')" ,
     cmd_lang = 'python' ,
     loc = locals() ,
     batch_submit = True ,
     batch_run_script = taskfile + '.sendmail.script' ,
     batch_noRun = True
  ); 
  
  pyalps.dwa.startRunScript(
    batch_run_script= taskfile + '.converge.script' ,
    batch_cmd_prefix = 'bsub -n1 -W480' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing_converge(h5_infile);

Density profile

The headers:

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

To extract the useful density profile in task1, say:

h5_infile = 'parm.validation.task1.in.h5'
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'
               );

Time-of-flight images

The headers:

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;

To extract the useful 2D column-integrated time-of-flight image from task1, say:

h5_infile = 'parm.validation.task1.in.h5'
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);

q_z   = numpy.array(band.q(2));
wk2_z = numpy.array(band.wk2(2));
wk2_z = numpy.array([q_z, wk2_z]).transpose();
wk2_z0= wk2_z[wk2_z[:,0] == 0.][0][1];
wk2_z = numpy.transpose(wk2_z[wk2_z[:,1]/wk2_z0 > 1e-4]);
q_z   = wk2_z[0]
wk2_z = wk2_z[1]

momentum_density = numpy.tile(momentum_density, reps=(1,1,2*((q_z[q_z[:] >= 0].size / L)+1)));

dummy = numpy.zeros(momentum_density.shape[2]);
dummy[dummy.size/2:dummy.size/2 + q_z[q_z[:] >= 0].size] = wk2_z[q_z[:] >= 0]
dummy[dummy.size/2-q_z[q_z[:] < 0].size:dummy.size/2]  = wk2_z[q_z[:] < 0]    

momentum_density = numpy.tensordot(momentum_density, dummy, axes=([2],[0]))
momentum_density = numpy.tile(momentum_density, reps=(4,4));

q_x   = numpy.array(band.q(0));
wk2_x = numpy.array(band.wk2(0));
wk2_x = wk2_x[(q_x[:] >= -2.)*(q_x[:] < 2.)];
q_x   = q_x[(q_x[:] >= -2.)*(q_x[:] < 2.)]

q_y   = numpy.array(band.q(1));
wk2_y = numpy.array(band.wk2(1));
wk2_y = wk2_y[(q_y[:] >= -2.)*(q_y[:] < 2.)];
q_y   = q_y[(q_y[:] >= -2.)*(q_y[:] < 2.)]
 
wk2 = numpy.outer(wk2_x, wk2_y);

q_x = numpy.array([q_x] * wk2.shape[1], float);
q_y = numpy.array([q_y] * wk2.shape[0], float).transpose();
tof_image = wk2 * momentum_density;
#tof_image = 0.5 *(tof_image[:,:] + pyalps.dwa.cyclic_shift(tof_image[::-1,::-1],by=(1,1)))

To say extract condensate fraction and visibility:

selection_max = numpy.outer([(q_x[0,:] >= -0.05) * (q_x[0,:] <= 0.05)] , [(q_y[:,0] >= -0.05) * (q_y[:,0] <= 0.05)])
selection_min = numpy.outer([(q_x[0,:] >= 0.45) * (q_x[0,:] <= 0.55)] , [(q_y[:,0] >= 0.45) * (q_y[:,0] <= 0.55)])                   

nk_max = sum(tof_image[selection_max])
nk_min = sum(tof_image[selection_min])

condensate_fraction = nk_max;
visibility = (nk_max - nk_min) / (nk_max + nk_min);

To obtain the plot:

#Coarse the grid...
mag = 16;
q_x = q_x[0:q_x.shape[0]:mag, 0:q_x.shape[1]:mag];
q_y = q_y[0:q_y.shape[0]:mag, 0:q_y.shape[1]:mag];
tof_image = tof_image[1:tof_image.shape[0]:mag, 0:tof_image.shape[1]:mag] * mag * mag
tof_image[tof_image < 0] = 0      
fig = matplotlib.pyplot.figure();
ax = fig.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( q_y, q_x, tof_image, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=10);
fig.savefig(h5_infile + '.tof_image.jpg');
fig.show();

To see 1D cross section:

matplotlib.pyplot.figure();
matplotlib.pyplot.plot(q_x[q_x.shape[0]/2,:], tof_image[tof_image.shape[0]/2,:]);
matplotlib.pyplot.show();