-
Notifications
You must be signed in to change notification settings - Fork 103
Fortran Simulation Quantum
The program qmc_walk_sho solves the diffusion equation in imaginary time
corresponding to the Schrodinger equation,
for a single simple harmonic oscillator.
Atomic units are chosen so that the effective diffusion coefficient is D=1/2.
A few hundred independent systems, or walkers, are simulated using a simple random walk
and a crude creation/destruction scheme based on the difference between the potential energy
and the trial energy.
The scheme is described in the text.
The value of the trial energy et is updated regularly,
and the hope is that, after convergence,
it will be equal to the correct ground-state energy for the system which, in this case, is 1/2.
The updating scheme, and several of the default parameters,
are taken from the following paper
Reasonable results for the energy and the ground-state wavefunction,
which is accumulated as a histogram of walker positions,
should be obtained using the default input values,
with an empty namelist &nml /;
these defaults include setting et initially to the exact ground state energy.
Other values such as &nml et=0.6 / may be supplied through the namelist in the usual way.
This type of simulation is sensitive to the initial value,
and quite noisy:
possible improvements are discussed in general terms in the text.
The program qmc_pi_sho carries out a path integral Monte Carlo simulation
for a single simple harmonic oscillator,
at a specified temperature and ring-polymer size P.
Larger values of P give energies closer to the exact quantum mechanical canonical ensemble average.
For this simple model,
exact results can also be calculated for the finite values of P used in the simulation
- KS Schweizer, RM Stratt, D Chandler, PG Wolynes, J Chem Phys, 75, 1347 (1981),
- M Takahashi, M Imada, J Phys Soc Japan, 53, 3765 (1984),
and a routine to evaluate these is included in the example. No special techniques are used to accelerate the simulation; standard Metropolis moves are employed. Default parameters correspond to P=8, T=0.2. The table below is for test runs at various values of P, keeping the same temperature, which generates a range of average energies between the classical limit E=0.2 and the quantum limit E=0.506784; in each case we compare with the exactly-known value for the same P.
| P | E (MC) | E (exact) |
|---|---|---|
| 2 | 0.3218(3) | 0.321951 |
| 3 | 0.3933(4) | 0.392308 |
| 4 | 0.4312(3) | 0.431618 |
| 5 | 0.4543(4) | 0.454545 |
| 6 | 0.4694(6) | 0.468708 |
| 7 | 0.4778(6) | 0.477941 |
| 8 | 0.4846(9) | 0.484244 |
The program qmc_pi_lj applies the path-integral method to the Lennard-Jones fluid.
The simplest, primitive, algorithm is used,
together with the crudest estimators for energy and pressure.
The program uses
single-particle Monte Carlo moves for the individual beads in the ring polymer,
along with translations of the centre-of-mass of each polymer.
As mentioned in the text, there are many improvements of all these aspects
of the algorithm, which are recommended for production work.
The program takes in configuration files cnf##.inp where the ## reflects
the polymer bead number, in the range 1 to P.
These files have the same format as the classical Lennard-Jones configurations.
They may be prepared in the same way as usual,
from the initialize program, from a run of mc_nvt_lj
at the appropriate density and temperature,
or from a run of qmc_pi_lj for a different value of P.
It does no harm if these starting configurations are simply duplicates of each other,
provided that a preliminary run is conducted to allow the polymers to equilibrate,
after which all the output files cnf##.out may be renamed to cnf##.inp.
For testing, we compare with a set of simulations of neon,
which are mainly based on an empirical pair potential, but include selected results for Lennard-Jones for the case P=32. The LJ parameters for neon are ε=36.8K, σ=0.2789nm, atomic mass m=20.18u, and hence a reduced de Boer parameter λ=0.092σ, which is the default value in the program. We choose their lowest-temperature state point, (T,ρ)=(25.8K,36.28nm-3)=(0.701087,0.787069) in reduced LJ units. We use N=108 atoms, compared with Neumann and Zoppi's N=256, and our runs are five times shorter (10 blocks of 10000 sweeps); these differences should not be critical. The maximum centre-of-mass displacement is 0.1σ. Because the intra-polymer springs increase in strength with P, the maximum displacement parameter for individual bead moves is reduced from 0.06σ for P=2 down to 0.02σ for P=32. These choices give reasonable acceptance ratios for both kinds of move. In the table below we report: the rms radius R of the ring polymers, the average spring potential E(spring) per atom, the kinetic energy KE per atom, total energy E per atom, and pressure p (the last two including the standard long-range correction for the applied cutoff Rc=2.5σ).
| P | R | E(spring) | KE | E | p |
|---|---|---|---|---|---|
| 1 † | 0.0 | 0.0 | 1.052 | -4.692(1) | -0.756(5) |
| 2 | 0.04043(1) | 0.8963(7) | 1.2070(7) | -4.454(1) | -0.408(6) |
| 4 | 0.04942(1) | 2.906(1) | 1.301(1) | -4.320(3) | -0.231(9) |
| 8 | 0.05181(2) | 7.073(3) | 1.340(3) | -4.261(3) | -0.150(8) |
| 16 | 0.05247(2) | 15.479(4) | 1.347(4) | -4.252(4) | -0.148(5) |
| 32 | 0.05261(4) | 32.287(8) | 1.365(8) | -4.233(8) | -0.140(6) |
| 32 ‡ | 0.053 | 32.3008 | 1.352 | -4.227 | -0.039 |
† For completeness the P=1 runs were performed using mc_nvt_lj.
‡ These results are taken from Table I of Neumann and Zoppi (2002), converted into LJ reduced units.
A drawback of this state point is that the pressure is negative, suggesting instability with respect to phase separation. Nonetheless we have seen no evidence of crystalline structure in the simulated configurations, and assume that the liquid phase is metastable.