ALPS 2 Tutorial:MC-09 Directed Worm Algorithm:Intensive Example 3
From ALPS
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();