Difference between revisions of "ALPS 2 Tutorial:MC-09 Directed Worm Algorithm:Intensive Example 2"

From ALPS
Jump to: navigation, search
(Checking thermalization)
(Checking thermalization)
Line 1,002: Line 1,002:
  
 
  resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.wlen');
 
  resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.wlen');
 +
h5_infiles = [];
 
  for resultFile in resultFiles:
 
  for resultFile in resultFiles:
 
   ratio_wlenx = pyalps.getParameters(resultFile)[0]['ratio_wlenx'];
 
   ratio_wlenx = pyalps.getParameters(resultFile)[0]['ratio_wlenx'];
Line 1,007: Line 1,008:
 
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
 
   N=pyalps.getMeasurements(resultFile,observable='Total Particle Number')[0]['mean']['value'];
 
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
 
   thermalized=pyalps.dwa.thermalized(resultFile, 'Total Particle Number', simplified=True);
   print(resultFile + ' , ratio_wlenx : ' + str(ratio_wlenx) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));
+
   if not thermalized:
 +
    h5_infiles.append(resultFile);
 +
    print(resultFile + ' , ratio_wlenx : ' + str(ratio_wlenx) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));
  
 
=== Observing on your own ===
 
=== Observing on your own ===

Revision as of 02:00, 19 June 2013

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

To check whether they are thermalized:

resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.N.1');
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', simplified=True);
  print('mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

Condensate fraction

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');
results = []; 
for resultFile in resultFiles:
  L = int(pyalps.getParameters(resultFile)[0]['L']);
  green_tof = pyalps.getMeasurements(resultFile, 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;
  
  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);
    
  N = pyalps.getMeasurements(resultFile, observable='Total Particle Number')[0]['mean']['value'];
  
  results.append([N, condensate_fraction]);

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

To check whether they are thermalized:

resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.a.1');
for resultFile in resultFiles:
  ratio_a = pyalps.getParameters(resultFile)[0]['ratio_a'];
  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', simplified=True);
  print('ratio_a : ' + str(ratio_a) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

Condensate fraction

Turn on respective measurements

Gather the files:

h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.a.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 joba -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 joba -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 joba -n1 -W60' ,
    loc = locals()
  );
for h5_infile in h5_infiles:
  parallel_automatic_processing(h5_infile);

Extracting final data(plot)

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.8, 0.82, 0.84, 0.88, 0.9, 0.92, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.08, 1.10, 1.12, 1.14, 1.16, 1.18, 1.2]:
    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);

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.errorbudget.V0');
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);

Unexpected worldlines configuration corruption

You may check worldlines configuration like the following:

import pyalps;
import pyalps.dwa;

Checking infiles:

all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
h5_infiles = [];
for h5_infile in all_h5_infiles:
   ar = pyalps.hdf5.h5ar(h5_infile);
   if ar.is_group('/simulation/worldlines'):
     wl = pyalps.dwa.extract_worldlines(h5_infile);
     valid = wl.is_valid(20);
     if not valid:
       h5_infiles.append(h5_infile);
       print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking outfiles:

all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.V0');
h5_infiles = [];
for h5_infile in all_h5_infiles:
   ar = pyalps.hdf5.h5ar(h5_infile);
   if ar.is_group('/simulation/worldlines'):
     wl = pyalps.dwa.extract_worldlines(pyalps.input2output(h5_infile));
     valid = wl.is_valid(20);
     if not valid:
       h5_infiles.append(h5_infile);
       print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking thermalization

resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.V0');
for resultFile in resultFiles:
  ratio_V0x = pyalps.getParameters(resultFile)[0]['ratio_V0x'];
  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', simplified=True);
  print(resultFile + ' , ratio_V0x : ' + str(ratio_V0x) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

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

First the headers again:

import numpy;
import pyalps;
import pyalps.dwa;

Easily done like this:

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

Prepare new input files:

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

Run them:

for h5_infile in h5_infiles:
  parallel_automated_processing(h5_infile);

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.8, 0.82, 0.84, 0.88, 0.9, 0.92, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.08, 1.10, 1.12, 1.14, 1.16, 1.18, 1.2]:
    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);

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.errorbudget.wlen');
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);

Unexpected worldlines configuration corruption

You may check worldlines configuration like the following:

import pyalps;
import pyalps.dwa;

Checking infiles:

all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
h5_infiles = [];
for h5_infile in all_h5_infiles:
   ar = pyalps.hdf5.h5ar(h5_infile);
   if ar.is_group('/simulation/worldlines'):
     wl = pyalps.dwa.extract_worldlines(h5_infile);
     valid = wl.is_valid(20);
     if not valid:
       h5_infiles.append(h5_infile);
       print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking outfiles:

all_h5_infiles = pyalps.getInputH5Files(prefix='parm.errorbudget.wlen');
h5_infiles = [];
for h5_infile in all_h5_infiles:
   ar = pyalps.hdf5.h5ar(h5_infile);
   if ar.is_group('/simulation/worldlines'):
     wl = pyalps.dwa.extract_worldlines(pyalps.input2output(h5_infile));
     valid = wl.is_valid(20);
     if not valid:
       h5_infiles.append(h5_infile);
       print('h5_infile : ' + h5_infile + ' , valid : ' + str(valid));

Checking thermalization

resultFiles = pyalps.getResultFiles(prefix='parm.errorbudget.wlen');
h5_infiles = [];
for resultFile in resultFiles:
  ratio_wlenx = pyalps.getParameters(resultFile)[0]['ratio_wlenx'];
  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', simplified=True);
  if not thermalized:
    h5_infiles.append(resultFile);
    print(resultFile + ' , ratio_wlenx : ' + str(ratio_wlenx) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));

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.wlen');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
N = pyalps.collectXY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_wlenx'])

then plot it:

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

or send it to yourself:

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

Fine tuning the chemical potential

This section is under construction