Bulk helium

pathvis++ Helium.0.h5

We now see a screen like this:

The simulation box can be rotated by left-clicking and dragging the mouse. You may zoom with the scroll wheel or by middle-clicking and dragging left (zoom out) or right (zoom in). At the bottom of the screen is a slider that allows us to change through frames of each successive Monte Carlo step that was dumped to the output file. You should notice here that the particle paths are are much smaller than the inter-particle spacing. This compactness essentially indicates that at all points through imaginary time, the particle is in nearly the same place. In this regime, the helium atoms would be reasonably well approximated by treating them as classical particles.

Let's now lower the temperature and see what happens when quantum effects become important. We will lower the temperature to approximately 4.2 K. Change tau (Section: System) to 0.0125 and find an appropriate value to change NumTimeSlices (Section: System) to in order to get the desired temperature. (Caveat: Typically, we wouldn't change tau significantly during a study because the time step error would change but we are doing it for the tutorial to keep the number of time slices small for different temperatures). Now, run your simulation and again view the simulation using pathvis++. You should find paths much like this:

Notice that these paths are significantly larger than the paths at 64 K. Their non-negligible extent is a manifestation of zero point motion. In this regime, the system would be poorly represented by a classical potential (which "point" in space along the path would you even use to represent the classical particle?). There are many simulations in which zero point motion is important (for example ambient water), even though exchange statistics aren't critical. You can tell that statistics aren't relevant in this simulation because the paths aren't permuting onto each other. The inter-particle separation is sufficiently large that a permutation would be very unfavorable in terms of the energetics of the springs. Notice that in the simulations at both 6 K and 64 K, the winding number is 0.

Now, let's run a simulation where the bosonic nature of Helium is important. Set the temperature to 2.0K (let tau be 0.025). Run this simulation and view the results using pathvis++. As you push the slider and advance through Monte Carlo steps we see the paths in our system expand and begin to permute onto each other. The existence of permuted paths indicates the importance of boson exchange statistics in this simulation. At some point, you will see a series of permuting paths wind across the entire box (to see this, it is sometimes helpful to toggle to Wrap mode with the button at top. This forces the paths to stay within the box whereas the other mode tries to enforce continuity of the paths). This wrapping around the box is an indication that the system is superfluid. See mathematically how superfluidity manifests itself within Path Integral Monte Carlo . Winding paths are indicated in yellow. Of course, the jumping around of these paths sometimes obfuscates their global structure which makes it hard to see what's going on. If we click "Smooth", these paths will be smoothed by truncating their Fourier spectrum. The "Detail" slider controls how many Fourier components are retained, i.e. you can adjust this slider to tune the desired amount of detail in the path. Finally, once you have a configuration and perspective you like, you can click on tubes to add depth cues to the paths. At this point, you should have a visualization similar to

In summary, we see that as the temperature decreases, quantum effects become increasingly important, specifically zero point motion and then quantum statistics.

Summary

64 K

4.2 K

1 K

Before we move on, we should make sure we are using simulation parameters (time step and system size) that lead to a well-converged result (i.e. introduce negligible systematic errors). We will also want to tune the algorithm (i.e., you may notice that you are spending too much time computing observables and not enough time actually performing the Monte Carlo moves of the simulation).

1. Checking Convergence
2. Algorithm Tuning

Of course, as we change other parameters, we would want to re-optimize these settings. However, an initial adjustment is adequate for our purposes.

[edit] Superfluid Transition

Let us now proceed to calculate the superfluid transition temperature of liquid helium. We can start to get a rough estimate of this value by measuring the superfluid fraction at a variety of different temperatures. Let us do a series of calculations from . The temperature is set as where is set by the product of the input variables tau (Section: System) and NumTimeSlices (Section: System). Let's make sure our time step is 0.025.

There is a only a discrete set of temperatures (in our relevant temperature range) we can choose because must be an integral multiple of tau, whose values are constrained by the density matrix files that were generated. Set up 5 different simulations in this temperature range (by changing the NumTimeSlices (Section: System) variable and see if you can locate the Superfluid transition (remember, you will also need to change the name for the OutFileBase (Section: Observables) if you want it to write different output files). Once these runs have completed, use the python analysis script Report.py to examine them

Report.py theFileName

and then plot the Superfluid Fraction as a function of temperature. (5 minute runs should be sufficient to get a rough outline of the graph especially at higher temperature). A reasonably converged graph should look like this (the red points are from the PIMC++ simulation. The green lines are experimental results from Donnelly, 1967):

An example of what you may see after 5 minutes is here

Although we know that this transition must be sharp in the thermodynamic limit ( constant), it appears very smooth and not at any definitive temperature in our graph. At best we can only qualitatively locate the transition temperature. This is because 10 particles is very far from the thermodynamic limit. (Remember that phase transitions are only formally defined in the thermodynamic limit) We could obviously mitigate these finite size effects by going to a system of more than 10 particles (and we will do so in a moment). (An even better method is to do a variety of system sizes and then use a finite size scaling analysis to pinpoint the exact transition temperature, but an explanation of how this works is beyond the scope of the tutorial.) Before we embark on these larger systems, though, it behooves us understand how to use parallelization with PIMC++.

Because we want to calculate the superfluid fraction as quickly as possible, we will want to take advantage of different modes of parallelization within PIMC++. The code supports two modes of parallelism: cloning and time slice parallization.

[edit] Cloning Parallelization

We will skip this section of the tutorial, since the computer lab is not set up with a parallel run environment. You can skip to the next Section. You can read on if you are interested in parallelization in PIMC++. However, the lab is not set up to run PIMC++ in parallel.

The first mode of parallelism is to run many independent copies of the simulation. Because the error bars of our simulation depend on the number of independent samples we collect, Monte Carlo is "trivially parallel" in this sense. PIMC++ automatically takes advantage of this form of parallelism as its default behavior. If you run the code (using mpirun or the equivalent for your MPI implementation) on more than one processor, it will run independently on each processor, creating output files labeled basename.PROCNUMBER.h5. Let's see this now by running PIMC++ on multiple processors. For a preconfigured system (these machines aren't), you can run

runPIMC Helium.in -np 4

 Report.py Helium 

[edit] Time slice Parallelization

We will also skip this section, since the computer lab is not set up with a parallel environment.

In discussing time slice parallelization, we will refer to two important quantities. refers to the number of processors allocated to PIMC++ at runtime. For example, if we had used the generic MPI command %mpirun -np 4 pimc++ myInput.in or the Tungsten script %runPIMC nyInput.in -np 4, NP would be 4. The other important quantity is ProcsPerClone, specified in the PIMC++ input file. For the Cloning mode of parallelization discussed in the previous section, ProcsPerClone=1.

Unfortunately, for very large systems, the trivial (or natural!) parallelization of Monte Carlo may not be enough. This is because the amount of wall clock time it takes for the system to come into equilibrium may become very large. In this case, we can also parallelize over the time slices. In this case, the path of M slices is divided between n processors (i.e. processor 0 can be in charge of slices 0-20, processor 1 can be in charge of slices 20-40, processor 2 in charge of 40-60, etc.). The two modes of parallelism can be combined.

Effective use of time slice parallelism can make one's wallclock running time largely independent of temperature. We can utilize this mode of parallelism by editing the Section:Parallel of our input file. The relevant variable, ProcsPerClone, sets how many processors are used for each independent copy of the simulation. The time slices will then be distributed evenly over these processors. To choose the correct number of processors, we should keep the following considerations in mind. To begin with, it is a rule that each processor must have at least slices where l is the maximum level that is used in any Bisection Move throughout the simulation. On the other hand, you don't want to have too many slices per processor or we won't gain as much a speedup as desired. Also, the total number of processors must be a multiple of ProcsPerClone. If , PIMC++ will run clones, each of which has ProcsPerClone processors assigned to it, thereby using both methods of parallelism in the same simulation. When , a run will produce only one output file that looks as if the simulation had been run in series.

Let's now experiment with these modes of parallelism. We will do a system of 10 time slices split over 2 processors so each processor has 5 slices (practically you would more likely do 500 slices spread over 25 processors so each processor had 20 slices). Set ProcsPerClone (Section:Parallel) equal to 2 instead of 1 and the NumTimeSlices (Section: System) variable to 10. Now, run your job on 2 processors. On Tungsten,

runPIMC Helium.in -np 2

You'll notice that it produces a single output file. Run the analysis on this output file and examine it. If you had run this same temperature earlier, you should notice your error bars are approximately a factor of smaller than they were previously.

[edit] Superfluid Transition Temperature Calculation (larger system)

We will now try to pin down the superfluid transition temperature more accurately by doing a larger system. In the winter school, we will do this in a group setting (if you are following this tutorial yourself, acquiring all the data will involve performing several longer runs). We will do a system of 32 particles and again try to find the transition temperature. With 32 particles, the size of the box is now for the same saturated vapor pressure density. Everyone should make sure they make the changes:

• NumParticles (Section:System/Particles/Species) to 32 and
• Box (Section:System) to [11.361397, 11.361397, 11.361397]

Now, choose a temperature in the range from 1-3K (preferably something different from your neighbor) and adjust NumTimeSlices (Section: System) to reflect that.

Now, everyone should run their simulation using

pimc++ Helium.in

As it is progressing, run your analysis script (Report.py) to make sure it looks like it is producing reasonable data. We will collect these numbers on the board and try to better locate the transition temperature.

Continue now to H2 molecule