Müller-Plathe method

From VASP Wiki

The thermal conductivity can be obtained by Fourier's law

,

where is the heat-flux vector and is the temperature gradient. In reverse nonequilibrium molecular dynamics proposed by Müller-Plathe[1], a temperature gradient () along selected lattice axis is created by splitting the simulation box into 2N slabs of identical thickness orthogonal to and swapping the velocity of the coolest particle of type from the slab 1 () with the velocity of the hottest particle of the same type from the slab N+1 (). Assuming that is orthogonal to the remaining two lattice vectors and , is computed as

where the first summation is taken over all swapping events, is the mass of a particle of the type , is the total simulation time, and represents an ensemble average of the quantity enclosed. The method is invoked by defining the parameter FMP_ACTIVE (see below). The simulation has to be started from a well-equilibrated POSCAR-file and one of the following thermostats can be utilized in the NVE ensemble:

thermostats tags
Langevin thermostat MDALGO = 3 and LANGEVIN_GAMMA=NTYP×0
Andersen thermostat MDALGO = 1 and ANDERSEN_PROB = 0.0
Mind:
  • For measurements of thermal conductivity , the Langevin thermostat is preferred because it utilizes the velocity Verlet algorithm. The velocity Verlet algorithm computes positions and velocities in a synchronized manner. Nevertheless, setting MDALGO = 1 can be useful when studying particle redistributions under a temperature gradient.
  • The Müller-Plathe method can be utilized in the canonical ensemble. In the canonical ensemble, the thermal conductivity can not be measured because the energy is not conserved. Still, it can be useful to study a material's properties with a temperature gradient under canonical conditions.

The behavior of the simulation is controlled via the following parameters defined in the INCAR file:

  • FMP_ACTIVE=[logical array] invokes the method and specifies whether or not the atomic type (as defined in POSCAR) is allowed for swapping.
  • FMP_DIRECTION=[integer] index of the lattice vector along which the gradient is created.
  • FMP_SNUMBER=[integer] number of slabs.
  • FMP_SWAPNUM=[integer] number pairs of particles exchanged in a single swapping event.
  • FMP_PERIOD=[integer] number of time steps between two swapping events.

The output needed to evaluate the expression for is written out to the REPORT file. In particular, lines introduced by the string "tsl>" contain the ID number of the slab (item 2), the number of atoms in the slab (item 3), and the instantaneous temperature of the slab (item 4). This information, is needed to evaluate and is written for each MD step: 

grep "tsl>" REPORT 
tsl>        1     128   348.740
tsl>        2     138   387.874
tsl>        3     136   391.949
tsl>        4     127   380.342
tsl>        5     113   432.304
tsl>        6     122   409.074
tsl>        7     121   406.230
tsl>        8     120   370.238
tsl>        9     118   377.384
tsl>       10     134   383.762
tsl>       11     109   376.807
tsl>       12     146   377.061
tsl>        1     131   374.098
 ...

After a swapping event, FMP_SWAPNUM lines starting with the string "fmp>" are written (only written if NSW>FMP_PERIOD): 

grep "fmp>" REPORT
fmp> swapping atoms:     1133    1080     1158.178       15.604
fmp> swapping atoms:     1109    1485     1160.495       56.120
fmp> swapping atoms:     1059    1281     1400.528       18.375
fmp> swapping atoms:     1054    1357     1162.322       21.041
...

These lines contain the ID numbers (items 4 and 5) and instantaneous temperatures of hot (item 6) and cold atoms (item 7). Since the instantaneous temperature of a single particle is defined via , the instantaneous temperatures of hot and cold atoms can be used, in a straightforward manner, to evaluate the numerator of the equation for given above: 

 grep "fmp>" REPORT | awk 'BEGIN {bolkEV=8.6173857e-5; Q=0.} {Q+= $6-$7} END {print 1.5*Q*bolkEV}'


Mind:
  • The swapping process defined above leaves the total energy unchanged. Thus, for instance, the total energy is a conserved quantity if the simulation with the Müller-Plathe method is realized in the NVE ensemble.
  • The lattice vector along which the gradient of temperature is evaluated must be orthogonal to the remaining two lattice vectors.
  • To obtain meaningful results, a large supercell must be used in this type of simulation. Additionally averaging over several MD runs is recommended.

Example

INCAR: 

ML_LMLFF       = .TRUE. # switch on machine learning
ML_MODE        = RUN    # execute in production mode
ML_OUTPUT_MODE = 0      # reduce written output
ML_OUTBLOCK    = 100    # write every 100th step

TEBEG          = 320    # simulation temperature
ISIF           = 2      # compute stress tensor without box update
IBRION         = 0      # run molecular dynamics
MDALGO         = 3      # use Langevin thermostat
LANGEVIN_GAMMA = 2*0    # define microcanonical ensemble
POTIM          = 0.50   # time step 0.5fs
NSW            = 100000 # number of time steps
NBLOCK         = 100    # reduce output writing for speed
NWRITE         = 0      # write only little output
KBLOCK         = 1000   # frequency of pair correlation computation

FMP_SNUMBER   = 12    # number of slabs
FMP_DIRECTION = 3     # heat will be transfered along 3rd lattice vector
FMP_PERIOD    = 400   # swapping will be made every 400 steps
FMP_SWAPNUM   = 1     # one pair will be swapped during every swapping event
FMP_ACTIVE    = .FALSE. .TRUE. # swap only second species defined in POSCAR

POSCAR:

Click to show POSCAR

To extract the data from a collection of REPORT files the following analysis script can be used. AnalysisScript.sh:

Click to show AnalysisScript.sh

Convergence Analysis for thermal conductivity

In general, a lot of molecular dynamics steps and a large supercell are needed to get converged results for thermal conductivity. Therefore, convergence analysis is a crucial step when measuring thermal conductivities. To measure the convergence of it is recommended to do several MD runs. The thermal conductivity can then be computed with the help of the provided analysis script. It should be noted that the analysis script assumes a linear temperature gradient. Hence, it is essential to check the linearity of the temperature gradient when using the script.

Convergence of . As the simulation time progresses, the averages of the plot incorporate an increasing amount of data until convergence is achieved

Related tags and articles

FMP DIRECTION, FMP_ACTIVE, FMP_SNUMBER, FMP_SWAPNUM, FMP_PERIOD

References

  1. ↑ F. Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, J. Chem. Phys. 106, 6082 (1997).
Retrieved from "https://www.vasp.at/wiki/index.php?title=Müller-Plathe_method&oldid=27500"