- Email: [email protected]
Validity of the low coupling approximation in the effective-potential Monte Carlo theory
18 March 19% PHYSICS ELSMER
LETTERS
A
Physics Letters A 212(1996) 161-166
Validity of the low coupling approximation in the effective-potential Monte Carlo theory Dominic Acocella a*‘,George K. Horton aq2,E. Roger Cowley b b Department
a Serin Physics Laboratory, Rutgers - the State University, Piscataway, NJ 08855-0849, USA of Physics, Camden College of Arts and Sciences, Rutgers - the State University, Camden, NJ 08102-1205,
Received 29 November
1995; accepted for publication Communicated by L.J. Sham
3 January
USA
1996
Abstract Recent fundamental advances in formulating a reliable all-temperature lattice dynamics contain an uncontrolled approximation - the low-coupling approximation - routinely used in implementing the effective-potential Monte Carlo theory. We find that this approximation is most reliable when used in conjunction with the improved effective-potential
Monte Carlo theory. PACS: 63.70. + h; 05.30. - d; 02.7O.Lq
There has recently been striking progress in lattice dynamics due to fundamental advances based on Feynman’s path-integral formulation of the partition function [ 11 and the subsequent effective-potential formulation of statistical mechanics [2,3]. Consequently, we now have the numerically efficient effective-potential Monte Carlo (EPMC) theory [4] and its successor, the improved effective-potential (IEP) theory [5], as well as the path-integral Monte Carlo (PIMC) theory [6]. The PIMC method is a numerical evaluation of the path integral and is in principle exact but rather time-consuming if a low-temperature specific heat or
’ Present address: Visual System Engineering. ics, Montreal, Canada H4L 4X4. * E-mail: [email protected]. Elsevier Science B.V. All rights reserved PII SO37S-9601(96)00037-O
CAE Electron-
elastic constant for a quantum crystal is required. No such work has yet been reported. The EPMC and IEP theories provide a fast and practical alternative. Where exact results are known, such as for a single particle in a potential well, EPMC theory is extremely accurate [7]. To make the EPMC and IEP theories computationally manageable, the low-coupling approximation (LCA) has routinely been made [4,8,9]. Without the LCA, EPMC and IEP simulations require computer times comparable to PIMC simulations. It is clearly unacceptable that the EPMC and IEP theories are subject to an uncontrolled approximation. An investigation of the validity of the LCA is essential to put the foundations of the new lattice dynamics on a sound footing, especially as the EPMC method is gaining in popularity [IO]. In this paper, we present the results of an EPMC simulation with-
162
D. Acocella et d/Physics
out making the LCA; we refer to this as an exactEPMC simulation, in contrast to an LCA-EPMC simulation in which the LCA is made. For our nearest-neighbor model of neon, we find that the LCA is not satisfactory. Our results suggest why this is so and show how to correct the error caused by the LCA. We recall that EPMC is a variational theory based on the path-integral representation of the partition function with a quadratic trial potential. The variations result in reducing the quantum partition function to one of classical form with an effective potential incorporating quantum effects. The full EPMC equations [4] are z=
(-!$.-)3n/2j &N
R
,-L-‘&R),
X+(R,-R,+x), t
D’%P = (v,,., - ~_&&JIB.~ - U,,.&
K IL? J/3
a2K( R) =
aRIm dR,, ’
(1)
where R is a 3N-dimensional coordinate. Since the partition function has the classical form, all classical effects are integrated exactly by the classical Monte Carlo simulation. So EPMC theory is exact in the classical limit and, furthermore, it reproduces the first term in the Wigner expansion [4]. Since the trial potential is quadratic, EPMC theory is exact for a harmonic oscillator. The purely quantum fluctuations of the atoms are taken into account by smearing the potential. Thus, the fully quantum harmonic term and all classical anharmonic terms are included exactly in the free energy, and
Letters A 212 (1996) 161-166
only the purely quantum anharmonic contributions are approximated. As a practical matter, a classical Monte Carlo simulation with these full EPMC equations takes much too much time. The effective potential V,,, must be solved iteratively, and each iteration involves diagonalization. The LCA speeds up the simulation dramatically by assuming that changes in a! and U, due to the deviation of the atoms from their equilibrium positions, are negligible. The result is that the smearing matrix DIJ and the frequencies o remain constant so that the only changes to V,,, come from the position dependence of 4(x). Thus, the iteration procedure need only be done once, at equi~ib~um, and the classical Monte Carlo simulation becomes an efficient way of evaluating the partition function. Although the LCA has never been investigated quantitatively, we can make several qualitative statements about it. The LCA is exact for a harmonic crystal. This is because the second derivative of the harmonic potential is constant and, thus, is unaffected by the smearing; this implies that U and w remain constant, as do CYand D”, which is precisely what is assumed in the LCA. Also, since the smearing represents the purely quantum fluctuations, the LCA can only affect the purely qu~tum contributions to the free energy and must be exact classically. In fact, the LCA even reproduces the first term of the Wigner expansion since this term has no frequency dependence [4]. Thus, the LCA preserves the high-temperature virtues of the full EPMC theory. Furthermore, the LCA reproduces exact-EPMC results at zero temperature, since only the equilibrium configuration contributes to the free energy in a classical Monte Carlo simulation at this temperature. Since the LCA is exact in both the high- and lowtemperature limits and is exact for a harmonic crystal, it has been assumed to be satisfacto~ for anharmanic crystals at all tem~ra~res. However, me error due to the LCA must be a maximum at some finite temperature, and we find that it is significant. To start our investigation of the LCA, we took a single particle in a quadratic-plus-qua&c potential well, for which exact results are easily obtained. Our purpose is only to acquire qualitative information on what kind of systematic errors are present in the EPMC formalism. This single-particle model has
D. Acocella et al./Physics
I
I
i
I
. 0.000
0
1
2
3
4
Temperature
(orb.
units)
5
Fig. 1. Difference between the exact free energy of a quartic oscillator and the exact-EPMC energy and the LCA-EPMC energy.
only one free parameter which was fixed by choosing the mass of the particle to correspond to that of neon. The temperature scale as well as the energy scale is thus fixed and we should not expect a comparison between the results of this single-particle model and those of a three-dimensional crystal to be quantitatively significant. It is the qualitative aspects of this single-particle model which are important. We compare, in Fig. 1, the exact free energy to that of an exact-EPMC calculation and to that of an LCA-EPMC calculation 3 by plotting the error (i.e., the difference with the exact free energy). Since the model potential has no natural energy scale, the axes are labeled with arbitrary units. Exact-EPMC theory at zero ~rn~ra~re is in error because it omits second-order (and higher-order) terms in the cumulam expansion. However, as the temperature increases, this error quickly goes to zero. The exact-
Leiters A 212 (1996) 161-166
EPMC theory is remarkably accurate at higher temperatures. On the other hand, the error in an LCAEPMC calculation initially grows, even though the size of the purely quantum fluctuations decreases. The error then slowly decreases, persisting well into the region where exact-EPMC calculations are virtually correct. It is also useful to know what systematic error the LCA introduces into LCA-EPMC calculations. To study this, we plot in Fig. 2 the difference between the LCA-EPMC free energy and the exact-EPMC free energy. The LCA-EPMC calculation gives the exact-EPMC result at zero degrees, as discussed earlier, and its error increases quickly to a maximum, and then decreases slowly. It is important to note that the size of this maximum error is of the same order as the exact-EPMC error at zero degrees (from Fig. 1). Although a single particle model cannot give any information on the volume dependence of the free energy, such as the pressure, we should expect similar behavior there. Before making a similar comparison for a threedimensional crystal, we discuss some computational aspects of an exact-EPMC simulation. It is. first and
0.006
I
1
I
I
‘;;; Ifi 0.005 5 .z e
0.004
z 5
0.003
z 0 $
0.002
z W
L= : 0.001
0.000
3 Technically, these arc not EPMC calculations since the classical integration is not done by Monte Carlo simulation. However, it is still the same effective-potential tbeory, and we call it EPMC theory to avoid confusion.
163
0
,
I
I
I
1
2
3
4
Temperature
(orb.
5
units)
Fig. 2. Difference between the LCA-EPMC and the exact-EPMC free energies.
164
D. Acocella et al./Physics
foremost, a classical Monte Carlo simulation with an effective potential. We use standard Metropolis sampling: one atom is moved; the change in the effective potential AV,, is calculated; if the change is negative then the move is accepted unconditionally and if it is positive then it is accepted with a probability e-a*“e~f. However, we note that V,,, is not easily calculated; Eqs. (1) are coupled and must be solved self-consistently. We start with an initial guess for U and w, solve for K,,, diagonalize it to get a new U and w,..., and iterate until convergence. Then we use these final lJ and w to calculate V,,,. This iterative calculation of V,,, must be done at every R, i.e., at every move. For each iteration step, the diagonalization procedure is of order N3 and, thus, is very time consuming. Optimization efforts, such as taking the largest convergence parameter possible and making the best initial guess possible, are necessary. In the end, each move requires over five minutes on a Sparcstation 2 for a lOO-atom crystal with periodic boundary conditions 4. To obtain acceptable accuracy for the energy and pressure, we required approximately 300000 moves, corresponding to 25000 hours of computation time for a single temperature and volume! We therefore used the 5 12-node CM-5 massively parallel supercomputer, which has a peak speed of 64 Gflops and which can theoretically reduce the total computation time of 300000 hours, for twelve different temperatures, to 21 hours. Of course, peak rates are never achieved and our total computation time on the CM-5 was equivalent to about 84 hours. However, we did not use the CM-5 as a single, very fast machine. Monte Carlo simulations are well suited to parallel architectures because each node can be used to simultaneously accumulate independent statistics. In fact, the CM-5 has four independent vector units per node, and each vector unit can accumulate independent statistics. Our program ran one temperature and volume on sixteen nodes, dividing the statistics into 64 independent blocks, one block per vector unit. Thus, twelve separate data points were run in parallel on 64 nodes, requiring approximately 600 hours.
4We use a Sparcstation 2 as a unit for comparison CM-5 computer is composed of Spare 2 processors.
because the
Letters A 212 (1996) 161-166 -4.20
-3.9
I
I
1
I
I ,a-
-4.0 -
-4.25
-4.1 -4.2 -
-4.45
’ 0.0
I
1
I
0.1
0.2
0.3
Temperature
(s/k)
Fig. 3. Internal energy at zero pressure versus temperature. Solid line: exact-EPMC; long-dashed line: LCA-EPMC; short-dashed line: IEP; solid circles: PIMC.
One problem with this procedure is that each of the 64 independent blocks must be thermalized separately, thus discarding 64 times more moves than would be necessary on a scalar machine. To avoid this problem, we started each block with a random thennalized configuration, which was obtained from an LCA-EPMC simulation. For a three-dimensional, nearest-neighbor Lenmud-Jones model of neon, we compare the results of a PIMC calculation 1111of the zero-pressure internal energy with those of an exact-EPMC calculation, an LCA-EPMC calculation [4], and with those of our recently published improved effective-potential Monte Carlo theory (IEP) [5]. The neon potential used was
where r = I x I, E = 72.09 X lo-l6 erg, and u = 2.7012 X 10m8 cm. We use this potential purely for the sake of comparison with earlier work. We will report on comparison with experiment using a realistic potential at a later date. We compare, in Fig. 3, the exact-EPMC internal energy with the LCA-EPMC internal energy. At zero
D. Acocella et al./Physics
temperature, as expected, they coincide since the LCA is exact at that temperature; and at T = 0.06, they are still rather close. However, at T = 0.14, differences arise, with the LCA-EPMC energy increasing, but the exact-EPMC energy decreasing and closely approaching the PIMC results. Even at the highest temperatures, near melting, the systematic error caused by the use of the LCA is significant. We note that the exact-EPMC energy varies unphysically with temperature, and the associated specific heat will be negative in that temperature range, In particular, it is impossible to avoid the negative specific heat behavior because EPMC theory is a variational theory, so the ground state energy must be higher than the exact energy. As the temperature increases, the exact-EPMC error decreases and the energy must, therefore, decrease. Thus, the negative specific heat is built into the formalism and is unavoidable. This illustrates the serious problems with the LCA-EPMC and exact-EPMC methods for our model of neon. Next, we compare our exact-EPMC results with the PIMC results. The size of the PIMC points corresponds to one standard deviation. It is gratifying that, at T = 0.25 and T = 0.41, the exact-EPMC results are indistinguishable from the PIMC results, as expected from the single particle case at higher temperatures. At lower temperatures, however, the exact-EPMC results are significantly in error, although, as the temperature increases, this error quickly goes to zero. Recall that second-order terms in the cumulant expansion are omitted in zero-temperature EPMC theory. In particular, the cubic term from the potential is of leading order in perturbation theory, which explains why the error in EPMC theory at zero temperature is rather large. According to the observations made when discussing the quartic oscillator, this means that the error in~oduced by the LCA will be of this order, and that it will persist to higher temperatures. This is exactly what happens. Including the cubic term in the EPMC theory would dramatically reduce the zero-temperature error and, thus, similarly reduce the error caused by the LCA at all temperatures. This is an important point and, in fact, is not specific to our particular illustration using Ne22. In EPMC theory in general, the odd terms in the potential do not contribute to the zero-point energy since the smearing is even. But, as the temperature in-
Letters A 212 (1996) 161-166
165
creases from zero, the atoms are moved away from their equilib~um sites in the classical Monte Carlo simulation, so that the odd terms do con~ibute at finite temperature [S]. In an exact-EPMC simulation, both the phonon terms (the frequency dependent terms in the effective potential) and the smeared potential 5 are configuration-dependent and, therefore, receive con~ibutions from all missing terms in the expansion of the potential. This is the reason why the exact-EPMC method becomes so accurate as the temperature increases. When using the LCA, however, the missing terms contribute to the smeared potential as the temperature increases, but neuer to the ~~un~n terms since these are constant. Thus, the terms which are missing at low temperature continue to be inadequately accounted for at higher temperatures. This explains why the error caused by the LCA is of the order of the exact-EPMC error at zero temperature and why it persists to higher temperatures. And since this is a general feature of the LCA, estimating the zero-temperature exact-EPMC error gives us an estimate, at all temperatures, of the systematic error introduced by the LCA. We now compare the exact-EPMC results with our IEP results. IEP theory adds the purely quantum cubic term as the leading correction to EPMC theory. This correction was modeled after the cubic correction introduced by the improved self-consistent phonon theory USC) to correct the deficiencies of the first-order self-consistent phonon theory [ 121. The form of the purely quantum cubic correction was chosen by using simple consistency arguments, but we could not be certain that these were legitimate. Keeping in mind our previous remarks, if the cubic contribution is included properly then the error introduced by the LCA must be very much reduced
at all tem~ra~res. The fact that this is what happens, i.e., that the IEP correction nearly completely cancels the errors due to the LCA, gives much
’ At zero temperature,the smeared potential (K in Eqs. (1)) is the zero-point potential energy and the phonon terms ((V,, - K) in Eqs. (1)) represent the zerepoint kinetic energy. But, at finite temperature, this is not the case and we refer to these quantities as the smeared potential and the phonon terms, rather than the potential energy and kinetic energy. The analogy is useful, however.
D. Acocella et ~l./Ph~si~s Letters A 212 (i996) 161-166
2
tem~ratures. This allows an estimation of the systematic error introduced by the LCA, for example by calculating the contribution of the cubic term. Furthermore, by using IEP theory, this error is not only estimated but also corrected for at all temperatures. An IEP simulation is as fast as an LCA-EPMC simulation and, for neon, reproduces PIMC results at all temperatures, as well as ISC results at the lowest temperatures [5] where no PIMC results are available. We conclude that, for solids in which exchange effects can be neglected, IEP theory is the best available numerically implementable, all-temperature lattice dynamical theory.
0.03 -
: E i n
0.02 -
zs ki
I=
0.01 0.00
. 0.0
0.2
0.3
Temperature
(e/k)
0.1
0.4
Fig. 4. Difference between the LCA-EPMC and the exact-EPMC internal energies at zero pressure versus temperature.
credibility to the IEP theory. The improvement in the ground-state energy shown by the IEP theory over the EPMC theory is approximately 0.7% and persists beyond half the melting ~rn~ra~re. This is a very si~i~c~t improvement. IEP theory was designed to be very reliable at zero temperature, and we now see that it is also very accurate at all temperatures, even for such a severe test as our model for Nez2. These results for the internal energy are typical and representative of other properties. Finally, we compare in Fig. 4 the LCA-EPMC internal energy at zero pressure with the exact-EPMC internal energy by plotting, in analogy with the single-particle case, their differences. Qualitatively, we see the same picture as in Fig. 2, although we realize that internal energies and free energies are different quantities, which justifies our use of a single-particle model for qualitative information. Again, the systematic error caused by the LCA rises from zero to a maximum and then falls. However, because neon melts at a relatively low temperature, we cannot follow the fall-off as we did in the single-particle model. Thus, we find that, although the LCA is necessary for practical reasons, LCA-EPMC theory is not satisfactory as a theory of solid neon. In general, the error the LCA introduces is of the size of the exactEPMC error at zero temperature and persists to high
This work was partially supported by the U.S. National Science Foundation under Grant No. DMR92-02907 and by the National Center for Supercomputing Applications under Grant No. DMR93-OOI5N and Grant No. DMR95-0005N.
References
[f] R.P. Feynman, Statistical mechanics (Addison-WesIey, Reading.. MA, 1988). [Z] R. Giachetti and V. Tognetti, Phys. Rev. Lett. 55 (1985) 912; Phys. Rev. B 33 (1986) 7647. [3] R.P. Feynman and H. Kleinert, Phys. Rev. A 34 (1986) 5080. [41 D. Acocella, G.K. Horton and E.R. Cowley, Phys. Rev. B 51 (1995) 11406. [51 D. Acocella, G.K. Horton and E.R. Cowley, Phys. Rev. Lett. 74 ( 1995) 4887. [6] K.E. Schmidt and D.M. Ceperly, in: The Monte Carlo method in condensed matter physics, ed. K. Binder (Springer, Berlin, 1992); A.R. McGum, in: Dynamical properties of solids, Vol. 7, eds. G.K. Horton and A.A. Maradudin (North-Holland, Amsterdam, 1995). [7] S. Srivastava and Vishwamittar, Phys. Rev. A 44 (1991) 8006. 181 S. Liu, G.K. Horton and E.R. Cowley, Phys. Rev. B 44 (1991) 11714. 191 A. Cuccoli, A. Macchi, V. Tognetti and R. Vaia, Phys. Rev. B 47 (1993) 14923. [lOI E.R. Cowley and G.K. Horton, in: Dynamical Properties of Solids, Vol. 7, eds. G.K. Horton and A.A. Maradudin (Noah-Hotl~d, Amsterdam, 1995). 1111 S. Liu, G.K. Horton, E.R. Cowley, A.R. McGum, A.A. M~adudin and R.F. Wallis, Phys. Rev. B 45 (1992) 9716. [I21 V.V. Gofdman, G.K. Horton and M.L. Klein, Phys. Rev. Lett. 21 (1%8) 1527.
LETTERS
A
Physics Letters A 212(1996) 161-166
Validity of the low coupling approximation in the effective-potential Monte Carlo theory Dominic Acocella a*‘,George K. Horton aq2,E. Roger Cowley b b Department
a Serin Physics Laboratory, Rutgers - the State University, Piscataway, NJ 08855-0849, USA of Physics, Camden College of Arts and Sciences, Rutgers - the State University, Camden, NJ 08102-1205,
Received 29 November
1995; accepted for publication Communicated by L.J. Sham
3 January
USA
1996
Abstract Recent fundamental advances in formulating a reliable all-temperature lattice dynamics contain an uncontrolled approximation - the low-coupling approximation - routinely used in implementing the effective-potential Monte Carlo theory. We find that this approximation is most reliable when used in conjunction with the improved effective-potential
Monte Carlo theory. PACS: 63.70. + h; 05.30. - d; 02.7O.Lq
There has recently been striking progress in lattice dynamics due to fundamental advances based on Feynman’s path-integral formulation of the partition function [ 11 and the subsequent effective-potential formulation of statistical mechanics [2,3]. Consequently, we now have the numerically efficient effective-potential Monte Carlo (EPMC) theory [4] and its successor, the improved effective-potential (IEP) theory [5], as well as the path-integral Monte Carlo (PIMC) theory [6]. The PIMC method is a numerical evaluation of the path integral and is in principle exact but rather time-consuming if a low-temperature specific heat or
’ Present address: Visual System Engineering. ics, Montreal, Canada H4L 4X4. * E-mail: [email protected]. Elsevier Science B.V. All rights reserved PII SO37S-9601(96)00037-O
CAE Electron-
elastic constant for a quantum crystal is required. No such work has yet been reported. The EPMC and IEP theories provide a fast and practical alternative. Where exact results are known, such as for a single particle in a potential well, EPMC theory is extremely accurate [7]. To make the EPMC and IEP theories computationally manageable, the low-coupling approximation (LCA) has routinely been made [4,8,9]. Without the LCA, EPMC and IEP simulations require computer times comparable to PIMC simulations. It is clearly unacceptable that the EPMC and IEP theories are subject to an uncontrolled approximation. An investigation of the validity of the LCA is essential to put the foundations of the new lattice dynamics on a sound footing, especially as the EPMC method is gaining in popularity [IO]. In this paper, we present the results of an EPMC simulation with-
162
D. Acocella et d/Physics
out making the LCA; we refer to this as an exactEPMC simulation, in contrast to an LCA-EPMC simulation in which the LCA is made. For our nearest-neighbor model of neon, we find that the LCA is not satisfactory. Our results suggest why this is so and show how to correct the error caused by the LCA. We recall that EPMC is a variational theory based on the path-integral representation of the partition function with a quadratic trial potential. The variations result in reducing the quantum partition function to one of classical form with an effective potential incorporating quantum effects. The full EPMC equations [4] are z=
(-!$.-)3n/2j &N
R
,-L-‘&R),
X+(R,-R,+x), t
D’%P = (v,,., - ~_&&JIB.~ - U,,.&
K IL? J/3
a2K( R) =
aRIm dR,, ’
(1)
where R is a 3N-dimensional coordinate. Since the partition function has the classical form, all classical effects are integrated exactly by the classical Monte Carlo simulation. So EPMC theory is exact in the classical limit and, furthermore, it reproduces the first term in the Wigner expansion [4]. Since the trial potential is quadratic, EPMC theory is exact for a harmonic oscillator. The purely quantum fluctuations of the atoms are taken into account by smearing the potential. Thus, the fully quantum harmonic term and all classical anharmonic terms are included exactly in the free energy, and
Letters A 212 (1996) 161-166
only the purely quantum anharmonic contributions are approximated. As a practical matter, a classical Monte Carlo simulation with these full EPMC equations takes much too much time. The effective potential V,,, must be solved iteratively, and each iteration involves diagonalization. The LCA speeds up the simulation dramatically by assuming that changes in a! and U, due to the deviation of the atoms from their equilibrium positions, are negligible. The result is that the smearing matrix DIJ and the frequencies o remain constant so that the only changes to V,,, come from the position dependence of 4(x). Thus, the iteration procedure need only be done once, at equi~ib~um, and the classical Monte Carlo simulation becomes an efficient way of evaluating the partition function. Although the LCA has never been investigated quantitatively, we can make several qualitative statements about it. The LCA is exact for a harmonic crystal. This is because the second derivative of the harmonic potential is constant and, thus, is unaffected by the smearing; this implies that U and w remain constant, as do CYand D”, which is precisely what is assumed in the LCA. Also, since the smearing represents the purely quantum fluctuations, the LCA can only affect the purely qu~tum contributions to the free energy and must be exact classically. In fact, the LCA even reproduces the first term of the Wigner expansion since this term has no frequency dependence [4]. Thus, the LCA preserves the high-temperature virtues of the full EPMC theory. Furthermore, the LCA reproduces exact-EPMC results at zero temperature, since only the equilibrium configuration contributes to the free energy in a classical Monte Carlo simulation at this temperature. Since the LCA is exact in both the high- and lowtemperature limits and is exact for a harmonic crystal, it has been assumed to be satisfacto~ for anharmanic crystals at all tem~ra~res. However, me error due to the LCA must be a maximum at some finite temperature, and we find that it is significant. To start our investigation of the LCA, we took a single particle in a quadratic-plus-qua&c potential well, for which exact results are easily obtained. Our purpose is only to acquire qualitative information on what kind of systematic errors are present in the EPMC formalism. This single-particle model has
D. Acocella et al./Physics
I
I
i
I
. 0.000
0
1
2
3
4
Temperature
(orb.
units)
5
Fig. 1. Difference between the exact free energy of a quartic oscillator and the exact-EPMC energy and the LCA-EPMC energy.
only one free parameter which was fixed by choosing the mass of the particle to correspond to that of neon. The temperature scale as well as the energy scale is thus fixed and we should not expect a comparison between the results of this single-particle model and those of a three-dimensional crystal to be quantitatively significant. It is the qualitative aspects of this single-particle model which are important. We compare, in Fig. 1, the exact free energy to that of an exact-EPMC calculation and to that of an LCA-EPMC calculation 3 by plotting the error (i.e., the difference with the exact free energy). Since the model potential has no natural energy scale, the axes are labeled with arbitrary units. Exact-EPMC theory at zero ~rn~ra~re is in error because it omits second-order (and higher-order) terms in the cumulam expansion. However, as the temperature increases, this error quickly goes to zero. The exact-
Leiters A 212 (1996) 161-166
EPMC theory is remarkably accurate at higher temperatures. On the other hand, the error in an LCAEPMC calculation initially grows, even though the size of the purely quantum fluctuations decreases. The error then slowly decreases, persisting well into the region where exact-EPMC calculations are virtually correct. It is also useful to know what systematic error the LCA introduces into LCA-EPMC calculations. To study this, we plot in Fig. 2 the difference between the LCA-EPMC free energy and the exact-EPMC free energy. The LCA-EPMC calculation gives the exact-EPMC result at zero degrees, as discussed earlier, and its error increases quickly to a maximum, and then decreases slowly. It is important to note that the size of this maximum error is of the same order as the exact-EPMC error at zero degrees (from Fig. 1). Although a single particle model cannot give any information on the volume dependence of the free energy, such as the pressure, we should expect similar behavior there. Before making a similar comparison for a threedimensional crystal, we discuss some computational aspects of an exact-EPMC simulation. It is. first and
0.006
I
1
I
I
‘;;; Ifi 0.005 5 .z e
0.004
z 5
0.003
z 0 $
0.002
z W
L= : 0.001
0.000
3 Technically, these arc not EPMC calculations since the classical integration is not done by Monte Carlo simulation. However, it is still the same effective-potential tbeory, and we call it EPMC theory to avoid confusion.
163
0
,
I
I
I
1
2
3
4
Temperature
(orb.
5
units)
Fig. 2. Difference between the LCA-EPMC and the exact-EPMC free energies.
164
D. Acocella et al./Physics
foremost, a classical Monte Carlo simulation with an effective potential. We use standard Metropolis sampling: one atom is moved; the change in the effective potential AV,, is calculated; if the change is negative then the move is accepted unconditionally and if it is positive then it is accepted with a probability e-a*“e~f. However, we note that V,,, is not easily calculated; Eqs. (1) are coupled and must be solved self-consistently. We start with an initial guess for U and w, solve for K,,, diagonalize it to get a new U and w,..., and iterate until convergence. Then we use these final lJ and w to calculate V,,,. This iterative calculation of V,,, must be done at every R, i.e., at every move. For each iteration step, the diagonalization procedure is of order N3 and, thus, is very time consuming. Optimization efforts, such as taking the largest convergence parameter possible and making the best initial guess possible, are necessary. In the end, each move requires over five minutes on a Sparcstation 2 for a lOO-atom crystal with periodic boundary conditions 4. To obtain acceptable accuracy for the energy and pressure, we required approximately 300000 moves, corresponding to 25000 hours of computation time for a single temperature and volume! We therefore used the 5 12-node CM-5 massively parallel supercomputer, which has a peak speed of 64 Gflops and which can theoretically reduce the total computation time of 300000 hours, for twelve different temperatures, to 21 hours. Of course, peak rates are never achieved and our total computation time on the CM-5 was equivalent to about 84 hours. However, we did not use the CM-5 as a single, very fast machine. Monte Carlo simulations are well suited to parallel architectures because each node can be used to simultaneously accumulate independent statistics. In fact, the CM-5 has four independent vector units per node, and each vector unit can accumulate independent statistics. Our program ran one temperature and volume on sixteen nodes, dividing the statistics into 64 independent blocks, one block per vector unit. Thus, twelve separate data points were run in parallel on 64 nodes, requiring approximately 600 hours.
4We use a Sparcstation 2 as a unit for comparison CM-5 computer is composed of Spare 2 processors.
because the
Letters A 212 (1996) 161-166 -4.20
-3.9
I
I
1
I
I ,a-
-4.0 -
-4.25
-4.1 -4.2 -
-4.45
’ 0.0
I
1
I
0.1
0.2
0.3
Temperature
(s/k)
Fig. 3. Internal energy at zero pressure versus temperature. Solid line: exact-EPMC; long-dashed line: LCA-EPMC; short-dashed line: IEP; solid circles: PIMC.
One problem with this procedure is that each of the 64 independent blocks must be thermalized separately, thus discarding 64 times more moves than would be necessary on a scalar machine. To avoid this problem, we started each block with a random thennalized configuration, which was obtained from an LCA-EPMC simulation. For a three-dimensional, nearest-neighbor Lenmud-Jones model of neon, we compare the results of a PIMC calculation 1111of the zero-pressure internal energy with those of an exact-EPMC calculation, an LCA-EPMC calculation [4], and with those of our recently published improved effective-potential Monte Carlo theory (IEP) [5]. The neon potential used was
where r = I x I, E = 72.09 X lo-l6 erg, and u = 2.7012 X 10m8 cm. We use this potential purely for the sake of comparison with earlier work. We will report on comparison with experiment using a realistic potential at a later date. We compare, in Fig. 3, the exact-EPMC internal energy with the LCA-EPMC internal energy. At zero
D. Acocella et al./Physics
temperature, as expected, they coincide since the LCA is exact at that temperature; and at T = 0.06, they are still rather close. However, at T = 0.14, differences arise, with the LCA-EPMC energy increasing, but the exact-EPMC energy decreasing and closely approaching the PIMC results. Even at the highest temperatures, near melting, the systematic error caused by the use of the LCA is significant. We note that the exact-EPMC energy varies unphysically with temperature, and the associated specific heat will be negative in that temperature range, In particular, it is impossible to avoid the negative specific heat behavior because EPMC theory is a variational theory, so the ground state energy must be higher than the exact energy. As the temperature increases, the exact-EPMC error decreases and the energy must, therefore, decrease. Thus, the negative specific heat is built into the formalism and is unavoidable. This illustrates the serious problems with the LCA-EPMC and exact-EPMC methods for our model of neon. Next, we compare our exact-EPMC results with the PIMC results. The size of the PIMC points corresponds to one standard deviation. It is gratifying that, at T = 0.25 and T = 0.41, the exact-EPMC results are indistinguishable from the PIMC results, as expected from the single particle case at higher temperatures. At lower temperatures, however, the exact-EPMC results are significantly in error, although, as the temperature increases, this error quickly goes to zero. Recall that second-order terms in the cumulant expansion are omitted in zero-temperature EPMC theory. In particular, the cubic term from the potential is of leading order in perturbation theory, which explains why the error in EPMC theory at zero temperature is rather large. According to the observations made when discussing the quartic oscillator, this means that the error in~oduced by the LCA will be of this order, and that it will persist to higher temperatures. This is exactly what happens. Including the cubic term in the EPMC theory would dramatically reduce the zero-temperature error and, thus, similarly reduce the error caused by the LCA at all temperatures. This is an important point and, in fact, is not specific to our particular illustration using Ne22. In EPMC theory in general, the odd terms in the potential do not contribute to the zero-point energy since the smearing is even. But, as the temperature in-
Letters A 212 (1996) 161-166
165
creases from zero, the atoms are moved away from their equilib~um sites in the classical Monte Carlo simulation, so that the odd terms do con~ibute at finite temperature [S]. In an exact-EPMC simulation, both the phonon terms (the frequency dependent terms in the effective potential) and the smeared potential 5 are configuration-dependent and, therefore, receive con~ibutions from all missing terms in the expansion of the potential. This is the reason why the exact-EPMC method becomes so accurate as the temperature increases. When using the LCA, however, the missing terms contribute to the smeared potential as the temperature increases, but neuer to the ~~un~n terms since these are constant. Thus, the terms which are missing at low temperature continue to be inadequately accounted for at higher temperatures. This explains why the error caused by the LCA is of the order of the exact-EPMC error at zero temperature and why it persists to higher temperatures. And since this is a general feature of the LCA, estimating the zero-temperature exact-EPMC error gives us an estimate, at all temperatures, of the systematic error introduced by the LCA. We now compare the exact-EPMC results with our IEP results. IEP theory adds the purely quantum cubic term as the leading correction to EPMC theory. This correction was modeled after the cubic correction introduced by the improved self-consistent phonon theory USC) to correct the deficiencies of the first-order self-consistent phonon theory [ 121. The form of the purely quantum cubic correction was chosen by using simple consistency arguments, but we could not be certain that these were legitimate. Keeping in mind our previous remarks, if the cubic contribution is included properly then the error introduced by the LCA must be very much reduced
at all tem~ra~res. The fact that this is what happens, i.e., that the IEP correction nearly completely cancels the errors due to the LCA, gives much
’ At zero temperature,the smeared potential (K in Eqs. (1)) is the zero-point potential energy and the phonon terms ((V,, - K) in Eqs. (1)) represent the zerepoint kinetic energy. But, at finite temperature, this is not the case and we refer to these quantities as the smeared potential and the phonon terms, rather than the potential energy and kinetic energy. The analogy is useful, however.
D. Acocella et ~l./Ph~si~s Letters A 212 (i996) 161-166
2
tem~ratures. This allows an estimation of the systematic error introduced by the LCA, for example by calculating the contribution of the cubic term. Furthermore, by using IEP theory, this error is not only estimated but also corrected for at all temperatures. An IEP simulation is as fast as an LCA-EPMC simulation and, for neon, reproduces PIMC results at all temperatures, as well as ISC results at the lowest temperatures [5] where no PIMC results are available. We conclude that, for solids in which exchange effects can be neglected, IEP theory is the best available numerically implementable, all-temperature lattice dynamical theory.
0.03 -
: E i n
0.02 -
zs ki
I=
0.01 0.00
. 0.0
0.2
0.3
Temperature
(e/k)
0.1
0.4
Fig. 4. Difference between the LCA-EPMC and the exact-EPMC internal energies at zero pressure versus temperature.
credibility to the IEP theory. The improvement in the ground-state energy shown by the IEP theory over the EPMC theory is approximately 0.7% and persists beyond half the melting ~rn~ra~re. This is a very si~i~c~t improvement. IEP theory was designed to be very reliable at zero temperature, and we now see that it is also very accurate at all temperatures, even for such a severe test as our model for Nez2. These results for the internal energy are typical and representative of other properties. Finally, we compare in Fig. 4 the LCA-EPMC internal energy at zero pressure with the exact-EPMC internal energy by plotting, in analogy with the single-particle case, their differences. Qualitatively, we see the same picture as in Fig. 2, although we realize that internal energies and free energies are different quantities, which justifies our use of a single-particle model for qualitative information. Again, the systematic error caused by the LCA rises from zero to a maximum and then falls. However, because neon melts at a relatively low temperature, we cannot follow the fall-off as we did in the single-particle model. Thus, we find that, although the LCA is necessary for practical reasons, LCA-EPMC theory is not satisfactory as a theory of solid neon. In general, the error the LCA introduces is of the size of the exactEPMC error at zero temperature and persists to high
This work was partially supported by the U.S. National Science Foundation under Grant No. DMR92-02907 and by the National Center for Supercomputing Applications under Grant No. DMR93-OOI5N and Grant No. DMR95-0005N.
References
[f] R.P. Feynman, Statistical mechanics (Addison-WesIey, Reading.. MA, 1988). [Z] R. Giachetti and V. Tognetti, Phys. Rev. Lett. 55 (1985) 912; Phys. Rev. B 33 (1986) 7647. [3] R.P. Feynman and H. Kleinert, Phys. Rev. A 34 (1986) 5080. [41 D. Acocella, G.K. Horton and E.R. Cowley, Phys. Rev. B 51 (1995) 11406. [51 D. Acocella, G.K. Horton and E.R. Cowley, Phys. Rev. Lett. 74 ( 1995) 4887. [6] K.E. Schmidt and D.M. Ceperly, in: The Monte Carlo method in condensed matter physics, ed. K. Binder (Springer, Berlin, 1992); A.R. McGum, in: Dynamical properties of solids, Vol. 7, eds. G.K. Horton and A.A. Maradudin (North-Holland, Amsterdam, 1995). [7] S. Srivastava and Vishwamittar, Phys. Rev. A 44 (1991) 8006. 181 S. Liu, G.K. Horton and E.R. Cowley, Phys. Rev. B 44 (1991) 11714. 191 A. Cuccoli, A. Macchi, V. Tognetti and R. Vaia, Phys. Rev. B 47 (1993) 14923. [lOI E.R. Cowley and G.K. Horton, in: Dynamical Properties of Solids, Vol. 7, eds. G.K. Horton and A.A. Maradudin (Noah-Hotl~d, Amsterdam, 1995). 1111 S. Liu, G.K. Horton, E.R. Cowley, A.R. McGum, A.A. M~adudin and R.F. Wallis, Phys. Rev. B 45 (1992) 9716. [I21 V.V. Gofdman, G.K. Horton and M.L. Klein, Phys. Rev. Lett. 21 (1%8) 1527.