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

Tamama's project

Headers:

import numpy;
import numpy.fft;
import pyalps;
import pyalps.dwa;

Fluctuation in N:

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

Fluctuation in a:

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

Fluctuation in V0:

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

Fluctuation in wlen:

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

Code:

obs = []; 
for resultFile in resultFiles:
  quantity = {};
  quantity['params'] = pyalps.getParameters(resultFile)[0];
  quantity['N'] = pyalps.getMeasurements(resultFile, observable='Total Particle Number')[0]['mean']['value']
  quantity['E'] = pyalps.getMeasurements(resultFile, observable='Energy')[0]['mean']['value']
  
  L = int(pyalps.getParameters(resultFile)[0]['L']);
  
  density = pyalps.getMeasurements(resultFile, observable='Local Density')[0]['mean']['value'].reshape([L,L,L]);
  density_sum2 = numpy.sum(density, axis=2);
  quantity['density:0'] = density[L/2,L/2,L/2];
  quantity['density_sum2:0'] = density_sum2[L/2,L/2]; 
  
  quantity['green:tof:0'] = pyalps.getMeasurements(resultFile, observable='Green Function:0')[0]['mean']['value'];
  
  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])
  
  if nk_min < 0:
    nk_min = 0;
  
  condensate_fraction = nk_max;
  visibility = (nk_max - nk_min) / (nk_max + nk_min);
  
  quantity['condensate fraction'] = condensate_fraction;
  quantity['visibility'] = visibility;
  
  obs.append(quantity);

Quantities:

N = numpy.array([item['N'] for item in obs]);

E = numpy.array([item['E'] for item in obs]);

n0 = numpy.array([item['density:0'] for item in obs]);

cn0 = numpy.array([item['density_sum2:0'] for item in obs]); 

fc = numpy.array([item['condensate fraction'] for item in obs]);

v = numpy.array([item['visibility'] for item in obs]);

ratio = numpy.array([item['params']['ratio_a'] for item in obs])

ratio = numpy.array([item['params']['ratio_V0x'] for item in obs])

ratio = numpy.array([item['params']['ratio_wlenx'] for item in obs])

Percentage changes (N):

N0 = 280000.;

pc = {};
pc['Energy'] = abs(0.5*(numpy.max(E) - numpy.min(E)) / E[abs(N-N0)/N0 < 1e-3]) ;
pc['density:0'] = abs(0.5*(numpy.max(n0) - numpy.min(n0)) / n0[abs(N-N0)/N0 < 1e-3]) ; 
pc['density_sum2:0'] = abs(0.5*(numpy.max(cn0) - numpy.min(cn0)) / cn0[abs(N-N0)/N0 < 1e-3]) ;  
pc['condensate fraction'] = abs(0.5*(numpy.max(fc) - numpy.min(fc)) / fc[abs(N-N0)/N0 < 1e-3]) ;
pc['visibility'] = abs(0.5*(numpy.max(v) - numpy.min(v)) / v[abs(N-N0)/N0 < 1e-3]) ;

Percentage changes (a):

pc = {};
pc['Energy'] = abs(0.5*(numpy.max(E) - numpy.min(E)) / E[abs(ratio - 1.) < 1e-3]) ;
pc['density:0'] = abs(0.5*(numpy.max(n0) - numpy.min(n0)) / n0[abs(ratio - 1.) < 1e-3]) ; 
pc['density_sum2:0'] = abs(0.5*(numpy.max(cn0) - numpy.min(cn0)) / cn0[abs(ratio - 1.) < 1e-3]) ;  
pc['condensate fraction'] = abs(0.5*(numpy.max(fc) - numpy.min(fc)) / fc[abs(ratio - 1.) < 1e-3]) ;
pc['visibility'] = abs(0.5*(numpy.max(v) - numpy.min(v)) / v[abs(ratio -1.) < 1e-3]) ;

Plotting:

matplotlib.pyplot.plot(N,E,'o');
matplotlib.pyplot.show();

matplotlib.pyplot.plot(N,n0,'o');
matplotlib.pyplot.show();

matplotlib.pyplot.plot(N,cn0,'o');
matplotlib.pyplot.show();

matplotlib.pyplot.plot(N,fc,'o');
matplotlib.pyplot.show();

matplotlib.pyplot.plot(N,v,'o');
matplotlib.pyplot.show();