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

From ALPS
Jump to: navigation, search
(Fine tuning chemical potential)
(Fine tuning chemical potential)
Line 583: Line 583:
 
   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('ratio_a : ' + str(ratio_a) + 'mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));
+
   print('ratio_a : ' + str(ratio_a) + ' , mu : ' + str(mu) + ' , N : ' + str(N) + ' , thermalized : ' + str(thermalized));
  
 
== Condensate fraction ==
 
== Condensate fraction ==

Revision as of 23:05, 17 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.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