Manual

Installation

PhonTS installation requires some software to be present on your system. PhonTS is written in Fortran90, hence Fortran90 compiler is required. It is most rigorously tested with Intel compiler, hence it is a recommended one. Currently only parallel version of the code is available, with the parallelization via the MPI set of libraries. Any version should work, but it was tested on OpenMPI and MPICH. Further, PhonTS relies on LAPACK set of libraries for the linear algebra operations. Any version of it should work, while supplied makefile links to Intel version (mkl libraries, version 10.0.011). Once all the needed pieces are available proceed with the PhonTS installation as follows:

  1. Download PhonTS zip archive here, if you have not done it yet.
  2. Unzip it by
    • unzip PhonTS-ver_num.zip (you will be prompted for a password: please contact developers for one)
  3. PhonTS directory appears: go to src subdirectory
    • cd PhonTS/src
  4. Two Makefiles are currently provided with the code: Makefile_ifort (for intel compiler, default) and Makefile_gnu (for gfortran).
  5. For intel:
    • cp Makefile_ifort Makefile
  6. For gnu fortran:
    • patch -p1 < Phonts.gnu.patch
    • cp Makefile_gnu Makefile
  7. Type:
    • make
    1. and PhonTS executable should be created. You are good to go!
  8. If that does not work, open up Makefile and adjust the location and/or version of the fortran compiler/mpi libraries (typically the latter will provide wrapper “mpif90” with all the paths taken care of), as well as desired LAPACK libraries. The following variables might need adjustments:
    • FF: desired fortran compiler/MPI wrapper
    • LDFLAGS: desired LAPACK set of libraries
    • FFLAGS: Fortran90 compilation options – adjust if compiler different from Intel is in use

If version for debugging needs to be created, type:

  • make clean; make debug

and PhonTS executable with many optional debugging features (timing, debugging information, run-time checks and no optimization) will be created

If you attempted to compile with other F90/MPI/LAPACK combination and compilation failed, please contact authors for support.

Getting started

PhonTS is driven by the single input file, called phonons_input.dat. This input file is keyword driven, freestyle type of file: it is does not matter what is the order of the keywords. Any line that starts from a “#” considered to be a comment line and is ignored, as well as empty lines. Most of the parameters have default values specified (for details see here), but a few are necessary for specification of the system under consideration. Directory “examples” in the PhonTS directory contains a few sample problems, one of the simplest one being “Argon”. phonons_input.dat found there contains absolute minimum amount of keywords, thus we consider this example here. Here is how it looks like:

cell 1. 1. 1. 
species 1
Ar 40.0 0.00
Lattice 3.9620
exp-6 1
Ar Ar 0.01051316 13. 3.85 5.625 9.625
natoms 4
fractional
Ar 0.0000 0.0000 0.0000
Ar 0.5000 0.5000 0.0000
Ar 0.5000 0.0000 0.5000
Ar 0.0000 0.5000 0.5000
end

Keyword “cell” determines geometry of the cell. Note that currently only the orthorhombic cells are supported, thus for the fcc crystal, one has to specify a conventional cubic cell (defined in this example by 1. 1. 1.). Keyword “species” determines a list of atomic species in the system, 1 being the number of different species. In this example, there is only one type of atoms, Ar, with mass 40 a.u. and charge of zero electrons. Keyword “lattice” determine the lattice constant at which simulations will be performed. “exp-6” is one of the available pair-wise potentials, so this keywords determines interactions between the atoms. “natoms” defines number of atoms in the cell, and “fractional” means that fractional coordinates of the atoms will be supplied. “end” specifies the end of the file.

Executing PhonTS in this directory should generate the following screen output:

Running on           16  processors...
 ================================================================
 This code does phonons calculations for the supplied structure.
 It is based on the HELL code, but requires phonon_input.dat file
 for keyword driven input, whose name is self-explanatory.
 See input.f or phonts.mse.ufl.edu for possible keywords and manual.
 ================================================================
     Generating a perfect crystal structure
     Currently satisfies minimum image
       Supercell size:           3           3           3
       Total number of ions/atoms=         108
 ===============================================
     Control parameters for the phonon calculations:
     Will do iterative solution, no symmetry implemented so far!!!
     Number of iterations:           1
     Mixing parameter: 0.50
     Number of temperature gradients:  1
     Temperature dependence will be done for   9 points, from   300.0 K to  1200.0 K
----------------------------------
     Number of types of atoms:     1
     Type     Charge       Mass 
      Ar      0.0000      40.0000
----------------------------------
     Interatomic potentials found: 2body                                             
     Type   Atom  Atom A_lj         B_lj      R0_lj       Rcut
     Exp-6  Ar    Ar    0.0105    13.0000     3.8500     5.6250  
----------------------------------
     K-point grid:  9 x  9 x 81
     Z-axis refinement:  10
     Total number of k-points:   6561
     Density of states will be done in 200 bins, from .00 Thz to 30.00 THz.
     Gaussian smoothing factor:   10.00000
----------------------------------
     Hell unit (lattice constant):   3.96200 Ang
     Unit cell vectors in hell units:   1.00000   0.00000   0.00000
     Unit cell vectors in hell units:   0.00000   1.00000   0.00000
     Unit cell vectors in hell units:   0.00000   0.00000   1.00000
     Supercell dimensions are (in terms of unit cells):   3   3   3
----------------------------------
     Number of atoms in the unit cell:    4
     Coordinates in lattice basis    :
   c  Ar        0.0000000000        0.0000000000        0.0000000000
   c  Ar        0.5000000000        0.5000000000        0.0000000000
   c  Ar        0.5000000000        0.0000000000        0.5000000000
   c  Ar        0.0000000000        0.5000000000        0.5000000000
     Coordinates in cartesian basis  :
   c  Ar        0.0000000000        0.0000000000        0.0000000000
   c  Ar        0.5000000000        0.5000000000        0.0000000000
   c  Ar        0.5000000000        0.0000000000        0.5000000000
   c  Ar        0.0000000000        0.5000000000        0.5000000000
     Parallelization is over k-points...
     Size of the dynamical matrix:     12 x     12
     Third derivatives cutoff:   5.625 Angstr
     External pressure               :   0.00000 GPa
     No electrostatics!
     Average Gruneisen parameter is  1.00000
     Basic output level...
     All Done!
 ================================================
  
     Computing energy of the initial structure: 
       Energy of the system:  0.1055285080350E+03
       Energy per cell     :  0.3908463260555E+01
     Entering BTE_driver...
     Max number of neighbours for D3 is:    55
     Shortest interatomic distance is:  2.801557 A
     Running parallel job: partitioning:         410  k-points per node...
     ENTERING memory_estimator...
        Memory consumption (memory_estimator):
          Current allocations:  0.15419 Gb.
          Available memory   :  0.81730 Gb.
            Total allocations:   2.46698 Gb.
     End of initialization
     Running parallel job: partitioning:         410  k-points per node...
     Starting to fold and diagonalize dynamical matrix...
     Speed of sound estimates:
     Transverse:      4353.91 m/s, Longitudinal:      6993.95 m/s, Average:      4798.33 m/s
     Debye temperature:       572.46 K
     Diagonalization time:     1.455
     Do a density of states now...
     Cleaning up space and calculating perturbation matrices...
     Perturbation calculations time:     0.367
     Iterative solution requested: starting iterations...
     Running parallel job: partitioning:          45  k-points per node...
     ENTERING memory_estimator...
        Memory consumption (memory_estimator):
          Current allocations:  0.20617 Gb.
          Available memory   :  0.74693 Gb.
            Total allocations:   3.29877 Gb.
     Done with preparations...
 ==========================================================
     Thermal Conductivity:
     Temp(K)       k_xx        k_xy        k_xz        k_yx        k_yy        k_yz        k_zx        k_zy        k_zz
    300.0000     25.3455      0.0072     -0.0028      0.0072     25.3307     -0.0005     -0.0028     -0.0005     24.4998
    412.5000     18.5007      0.0056     -0.0021      0.0056     18.4895     -0.0003     -0.0021     -0.0003     17.8975
    525.0000     14.5657      0.0046     -0.0017      0.0046     14.5567     -0.0003     -0.0017     -0.0003     14.0958
    637.5000     12.0092      0.0038     -0.0014      0.0038     12.0017     -0.0002     -0.0014     -0.0002     11.6241
    750.0000     10.2151      0.0033     -0.0012      0.0033     10.2088     -0.0002     -0.0012     -0.0002      9.8887
    862.5000       8.8869      0.0029     -0.0010      0.0029       8.8813     -0.0002     -0.0010     -0.0002      8.6035
    975.0000       7.8640      0.0025     -0.0009      0.0025       7.8591     -0.0001     -0.0009     -0.0001      7.6136
   1087.5000      7.0521      0.0023     -0.0008      0.0023      7.0477     -0.0001     -0.0008     -0.0001      6.8278
   1200.0000      6.3921      0.0021     -0.0007      0.0021      6.3881     -0.0001     -0.0007     -0.0001      6.1889
 
     Iteration    1 completed, Time:  0.5557855E+02 sec.
     All Done
 ===========================================================

PhonTS first reports how many processors it is running on and prints the greeting message. It then proceeds to construct a simulation supercell from the cell that is supplied in the input file that satisfy the minimum image requirement with respect to longest real-space cutoff in the system and reports its dimensions. The next block of data is a printout of the input parameters, including some that are set by default. For example, thermal conductivity is set to be computed by the iterative technique, but only a single iteration is performed by default, thus producing relaxation time approximation result. Particularly important is the k-mesh used for the calculations, which is set to default of 9x9x9, but typically is a part of the input file. Further, the energy of the ideal system is computed for the verification purposes, and various other quantities, such as speed of sound and Debye temperature, and parallelization parameters are reported. Memory estimator attempts to report the memory footprint of the code, since all the arrays are allocated withing the code, and memory is the main bottleneck for the scalability. This routiune accesses system information about PhonTS process (this is one of the places that unlickely to trasfer easily to other compilers) and reads how much memory it takes. Estimates of the memory leftovers are often misleading, however, and should not be taken seriously. Finally, once PhonTS reposts “Done with preparations…” All the preliminary things are done, and PhonTS enters the main loop that cycles through all possible 3-phonons processes and adds up their contribution to thermal conductivity. This is typically the most time consuming part of the calculations. Finally, thermal conductivity tensor is reported, together with the differences between current iteration and the previous one (in the case of the iterative solution). In the simulations directory, there is also a number of other files created. Among the important ones are “pdos.dat” that contains phonon density of states and “tc_dos.dat” that contains spectral thermal conductivty.

Recommended procedure

Based on our expirience performing thermal conductivity calculations and comparing with the results of the molecular dynamics simulations (which employs different set of approximations, thus comparison have to be done carefully) we recommend the following setup for the calculations.

  1. Identify the classical potential suitable for the simulations for the simulations of the system of interest. These might be available from the literature, online depositories, or even developed for this specific purpose. If the form of the potential is not implemented in PhonTS then one can attempt to implement it (see separate sections on this), or use interface to another code (LAMMPS).
  2. Perform a quench of the system of interest in order to identify 0K lattice constant. This will also verify that potential is correctly transfered to PhonTS.
  3. Perform QHA simulations in order to determine temperature dependent volume of the system. Since we employ very simple SZIZA approximation (simply scale the volume in QHA), if atomic relaxation is essential, then another code has to be used to determine equilibrium structure at given temperature.
  4. For the temperature of interest, perform thermal conductivity simulations.

This setup will reqiuire separate calculations for the thermal conductivity at each temperature, which might be too time consuming. An obvious shortcut is to use a single volume (for exampl T=0 volume) for all temperatures of interest, however, in our expirience this will slightly overestimate thermal conductivity, especially at high temperatures.

Structural optimization

PhonTS can perform a structural optimization of the supplied structure. Optimizsation performe via very simple steepest decent algorithm, namely, the atoms are moved in the direction of the forces with some step which can be controlled by the user. There is also a variation of the algorithm, which allow for a change of the minimization step on the fly, which may improve the performane of this step. Either constant volume, or simulatlations under hydrostatic pressure can be performed, which is controlled by the “cell_control” keyword. User also have control over the precision of the optimization procedure.

First Principles Calculations

Key advantage of the BTE approach over molecular dynamics is the possibility of the first principles simulations of the thermal conductivity. In particular Density Functional Therory (DFT) caluclations strike the balance between the accuracy and computational resources required. As a result they are extremely popular in the community, with the number of the well-tested codes available to perform such simulations (VASP, QuantumEspresso, ABINIT to name just a few). PhonTS relies on these codes (interfaces to VASP and QuantumEspresso only are currently implemented, but authors will be happy to work work with you to produce an interface to Your favorite code) to perform first principles calculations. The quantities that are needed from the first principles calculations are the second and third derivatives of the total energy with respect to atomic positions. Those can be obtained in two different ways. Firstly, analytical calculations in principle is possible via the Density Functional Perturbation Theory (DFPT). Since all the codes mentioned above are the plain-wave codes, calclulations essentially proceeds in the recipical space, and one obtains Fourier transformed force constants on requested grid of k-points. Second derivatives are the function of a single k-point, while third one is a function of the pair of k-points. However, currently third derivatives calculations for all pairs of k-points is not included in any of the standard codes to the best of our knowledge. Moreover, if one is interested in using any of the formalism beyond standard DFT (such DFT+U, or hybrid functionals), then second derivatives are not implemented for this case either. One therefore has to resort to the second method of calculating derivatives, namely via the numerical differentiation. In this approach one starts with the supercell, and “manually” displace atoms from the equilibrium positions and compute derivative by numberically differentiating the forces, which are typically computed analytically. Currently, PhonTS is coupled to QuantumEspresso to obtain either second or first derivatives, and to VASP to obtain the first derivatives in real space, and then internally computes the rest via the fitnite differences. The key issue in this calculations is accuracy: numerical differentiation performed twice significantly reduces the accuracy of the third derivatives. It is thus recommended to set convergence parameters (especially for the self-consistency and internal Fourier grids) as tight as one can afford. In the example below, we consider calculations performed with VASP via the supercell approach.

In light of the discussion above, first principles calculations of thermal conductivity with PhonTS is performed in 3 stages:

  1. Preparation run, where phonon_input.dat file tells PhonTS the system of interest and PhonTS prepares the number of structure files for the ab-initio calculations.
  2. These files are then used by ab initio code to perform first principles calculations. All other input files required by those is user responsibility.
  3. PhonTS is run again, this time it reads the output of first principles calculations together with phonons_input.dat and performes actual thermal conductivity calculations from first principles.

All the intermediate steps, such as file copying can be performed by the shell scripts, and example included with distribution provides scripts used by the authors. Users can either adopt those, or make up their own, more suited for their needs.

Below we discuss these steps in turn on the example of the thermal conductivity of solid argon (at the pressure of 50 GPa) from the recent paper. Note, that scripts used here need the directory structure in PhonTS/example/AbInitio folder. AbInitio_reference contains all the data after the calculations for comparison. We start in PhonTS/example/AbInitio/PhonTS folder.

  1. In order to tell PhonTS that ab initio calculations will be performed, one uses the keyword AbInitio followed by two logical variables. First one specifys that preparation run will be done, while second one signifies production run. Therefore only “AbInitio T F” or “AbInitio F T” is allowed (somewhat redundant, but…). One also have to specify the interface to be used, or in other words, what ab initio code is intended to be used since they have different input files formats. This is done by “FP_interface” keyword, with possible values of VASP or QE. We intend to use VASP here hence one has to have in phonons_input.dat for the interaction keywords:
    1. AbInitio T F
    2. FP_interface VASP
    3. delta 0.005
    4. The last key word specifies the magnitude of finite displacement in the units of lattice constant. Default values here, 0.0005 is typically too small for first principles calculations. After that, PhonTS needs to be run. This run is very simple and fast (no need to run it in parallel, or submit to a queue), since PhonTS will only analyze supplied structure, create supercell and will make all independent displacements of atoms in order to determine second and third derivatives. Note, that this will work in exactly the same way if “phonons_only T” is specified (much fewer number of files is produced), or if QHA is specifyed. Thus these two types of calculations from first principles are also supported by themselves. The number of files produced might be rather large, as our summetry analyser is not very sophisticated, and there is some (small) redundancy. In this example, 182 input structures is computed, which covers third derivatives out to 5th nearest neighbours. All the files produced have an extension of “.str” and have a basename according to a naming convention which labels what kind of atoms are displaced, and which direction (central differences are used). Also, the lattice constant is included in the file name in order to distinguish files in the QHA type of calculations.
  2. Now, that all input structures are created, we move to ../VASP folder. There one can find some scripts (prep_script, submit_script, post_script and submission script submit.sh) and files that are needed for the VASP calculations: INCAR, POTCAR and KPOINTS. VASP executable also need to be placed in this folder (not provided for copyright reasons). Note that INCAR instructs VASP to perform a single point calculations: our goal is to obtain forces for these configuration of atoms. One starts with running preparation script “prep_script” which sets everything up for the VASP calculations. Namely, it copies all the structure files into VASP directory from PhonTS directory, for each structure file it creates a separate directory (named basename.dir), copies INCAR, POTCAR, KPOINTS and corresponding structure file as a POSCAR there and also adds a submission script (need to be modified for your computer). Each job gets its own name with the specified prefix for convinience.
    1. Now we are in position to run VASP calculations. Submit script conviniently dives into every directory created and submits job there: one ends up with 182 jobs. It is time for user to have a break – first principles calculations take some time.
  3. Coming back after the break, we find all VASP calculations done and can continue now. One execute “./post_script” that handles postprocessing. All it does is dives into each directory and copy OUTCAR file back to ../PhonTS directory and renaming it to the “basename.out”. phonon_inputs.dat needs to be adjusted for a production run by changing the logical variables in “AbInitio” key word:
    1. AbInitio F T
    2. FP_interface VASP
    3. delta 0.005
    4. And PhonTS needs to be executed again. This time it will read in the results of the first principles calculations and perform calculations of the thermal conductivity, the latter step being quite lenghly, so one might want to execute it in parallel. Note, that we request full iterative solution here – this makes about 10% difference in Argon.

Description of the output files

screen outputMirrors the input parameters specified in phonons_input.dat, as well as default parameters of the simulations. Reports thermal conductivity.
pdos.datReports phonons density of states. First three coloumns is the energy in different units (THz, cm -1 and meV), 4th coloumn is the density of states as calculated, and the 5th one is density of states smoothed out by the gaussian lineshape.
tc_dos.datThis file contains spectral thermal conductivity (contribution to thermal conductivity from phonons in the frequency interval) for xx, yy and zz components. Its calculation is similar top phonon density of states calculations, thus raw data typically contains unphisical sharp spikes due to relatively coarse grid. PhonTS therefore spoothes it out using the same gaussian smoothing as for the PDOS. First coloumn is the frequency in THz, folloed by xx component raw data, then smoothed out, then the same for yy and zz components. Each temperature has its own block.
phon_lifetimes.datHere we have phonons lifetimes in picoseconds as a function of the k-point and temperature. For each k-point (specified by a number and k x ,k y ,k z ) there are 3N lifeitimes which correspond to 3N phonons at this k-point, where N is the number of atoms.
eigenvectors.datThis file contained phonons eigenvectors. The format of the file is first column specifies the eigenvector at given k-point, the second specifies the component of the eigenvetor, followed by real and imaginary part of the value. The k-point number is the last number in the line.
scatt.datThis file contains the scattering amplitudes map, which essentially shows with which phonons a phonon at particular frequency interact the most. It is a 3D data, with x and y (first two coloumns) being the energy of the original phonon and phonons it scatters into, while the following 3 coloumns samples the amplitude of different scattering processes. The last 3 coloumns have the same data but smoothed out by gaussians. There are as many blocks of data as there are temperature points.
freq.dat
group_vel.dat