Merge branch 'merge-dftb-interface' into develop

This commit is contained in:
Atsushi Togo 2019-01-12 07:40:45 +09:00
commit 7cc8ca802e
10 changed files with 1764 additions and 7 deletions

51
doc/dftb+.rst Normal file
View File

@ -0,0 +1,51 @@
.. _dftbp_interface:
DFTB+ & phonopy calculation
=========================================
How to run
-----------
DFTB+ phonon band structures are created as follows:
1) Create a DFTB+ input file dftb_in.hsd that is set up to perform a
single-point energy and force calculation for a structure which is named
``geo.gen``. The dftb_in.hsd file should turn on force evaluation by setting
``CalculateForces = Yes`` in its analysis block, and write the tagged results
by enabling ``WriteResultsTag = Yes`` in its options.
2) Generate the the required set of structures and the ``disp.yaml`` file by
issuing the command ::
% phonopy -d --dim="4 4 4" --dftb+
This example builds 2 x 2 x 2 supercell files. The undistorted supercell is
stored in ``geo.genS``, while the required displacements are stored in files
matching the pattern ``geo.genS-*``. Note that you have to increase the
supercell dimension until you reach convergence of the band structure.
2) For each each ``geo.genS-*`` structure perform a DFTB+ calculations,
retaining the resulting ``detailed.out`` file.
3) Create the ``FORCE_SETS`` file with the command ::
% phonopy -f disp-*/results.tag --dftb+ ...
Where the location of all of the ``results.tag`` files is given on the
command line. To run this command, the ``disp.yaml`` file has to be located
in the current directory, because the atomic displacements are written into
the FORCE_SETS file.
4) Create a ``band.conf`` file to specify the path in the Brillouin zone you are
interested in (see the phonopy documentation). Then post-process the phonopy
data, either in the settings file (DIM) or by providing the dimensions of the
the supercell repeat on the command line ::
% phonopy -p band.conf --dim="4 4 4" --dftb+
5) Create a band structure in gnuplot format ::
% phonopy-bandplot --gnuplot band.yaml > band.dat
All major phonopy options should be available for the DFTB+ interface.

View File

@ -0,0 +1,18 @@
This is an example of the use of the DFTB+ interface. An input is
given in the dftb_in.hsd file along with the conventional unit cell in
geo.gen. This example also requires the pbc-0-3 Slater-Koster
parameteters for carbon
(https://www.dftb.org/parameters/download/pbc/pbc-0-3-cc/).
1) Create the (one) displaced supercell structure(s):
phonopy -d --dim="4 4 4" --dftb+
2) Run the supercell input with dftb+. This has been pre-calculated
with DFTB+ version 18.2, with the required data stored in the file
results.tag.
3) Collect the forces:
phonopy -f results.tag --dftb+ ...
4) Calculate phonon dispersion data into band.yaml and save band.pdf:
phonopy --dftb+ --dim="4 4 4" -p -s band.conf

View File

@ -0,0 +1,3 @@
BAND = 0 0 0 1/2 0 1/2 1/2 1/4 3/4 3/8 3/8 3/4 0 0 0 1/2 1/2 1/2 5/8 1/4 5/8 1/2 1/4 3/4 1/2 1/2 1/2 3/8 3/8 3/4, 5/8 1/4 5/8 1/2 0 1/2
BAND_LABELS = $\Gamma$ X W K $\Gamma$ L U W L K:U X
BAND_POINTS = 101

View File

@ -0,0 +1,48 @@
# Input for DFTB+ diamond primitive cell example, see
# https://www.dftbplus.org/ for details of code useage
Geometry = GenFormat {
# geometry input file
<<< geo.genS-001
}
Hamiltonian = DFTB {
SCC = No # non-self-consistent, as diamond
MaxAngularMomentum = {
C = "p" # s and p functions, read parameter file for details
}
# local location of Slater-Koster parameter data files.
# These are available at https://www.dftb.org/parameters
SlaterKosterFiles = Type2FileNames {
Separator = "-"
Suffix = ".skf"
}
# Example uses a fairly large supercell, so not many k-points are included
KPointsAndWeights = SupercellFolding {
2 0 0
0 2 0
0 0 2
0.5 0.5 0.5
}
}
Options = {
# Required options for storing data needed by phonopy
WriteResultsTag = Yes
}
ParserOptions = {
# input parser version for DFTB+ 18.2, but this is backward
# compatible
ParserVersion = 6
}
Analysis = {
# required option for phonopy
CalculateForces = Yes
}

View File

@ -0,0 +1,14 @@
8 F
C
1 1 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
2 1 0.2500000000E+00 0.2500000000E+00 0.2500000000E+00
3 1 0.5000000000E+00 0.5000000000E+00 0.0000000000E+00
4 1 0.7500000000E+00 0.7500000000E+00 0.2500000000E+00
5 1 0.5000000000E+00 0.0000000000E+00 0.5000000000E+00
6 1 0.7500000000E+00 0.2500000000E+00 0.7500000000E+00
7 1 0.0000000000E+00 0.5000000000E+00 0.5000000000E+00
8 1 0.2500000000E+00 0.7500000000E+00 0.7500000000E+00
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.3561772125E+01 0.0000000000E+00 0.0000000000E+00
0.0000000000E+00 0.3561772125E+01 0.0000000000E+00
0.0000000000E+00 0.0000000000E+00 0.3561772125E+01

File diff suppressed because it is too large Load Diff

View File

@ -83,7 +83,8 @@ def get_parser():
elk_mode=False,
siesta_mode=False,
cp2k_mode=False,
fc_symmetry=False,
dftbp_mode=False,
fc_symmetry=None,
fc_format=None,
fc_spg_symmetry=False,
fits_debye_model=False,
@ -220,6 +221,9 @@ def get_parser():
parser.add_argument(
"-d", "--displacement", dest="is_displacement", action="store_true",
help="Create supercells with displacements")
parser.add_argument(
"--dftb+", dest="dftbp_mode", action="store_true",
help="Invoke dftb+ mode")
parser.add_argument(
"--dim", nargs='+', dest="supercell_dimension",
help="Same behavior as DIM tag")

View File

@ -48,14 +48,13 @@ def get_interface_mode(args):
"""
calculator_list = ['wien2k', 'abinit', 'qe', 'elk', 'siesta', 'cp2k',
'crystal', 'vasp']
'crystal', 'vasp','dftbp']
for calculator in calculator_list:
mode = "%s_mode" % calculator
if mode in args and args.__dict__[mode]:
return calculator
return None
def write_supercells_with_displacements(interface_mode,
supercell,
cells_with_disps,
@ -103,6 +102,9 @@ def write_supercells_with_displacements(interface_mode,
optional_structure_info[1],
num_unitcells_in_supercell,
template_file="TEMPLATE")
elif interface_mode == 'dftbp':
from phonopy.interface.dftbp import write_supercells_with_displacements
write_supercells_with_displacements(supercell, cells_with_disps)
def read_crystal_structure(filename=None,
@ -181,7 +183,10 @@ def read_crystal_structure(filename=None,
from phonopy.interface.crystal import read_crystal
unitcell, conv_numbers = read_crystal(cell_filename)
return unitcell, (cell_filename, conv_numbers)
elif interface_mode == 'dftbp':
from phonopy.interface.dftbp import read_dftbp
unitcell = read_dftbp(cell_filename)
return unitcell, (cell_filename,)
def get_default_cell_filename(interface_mode):
if interface_mode == 'phonopy_yaml':
@ -200,10 +205,11 @@ def get_default_cell_filename(interface_mode):
return "unitcell.inp"
elif interface_mode == 'crystal':
return "crystal.o"
elif interface_mode == 'dftbp':
return "geo.gen"
else:
return None
def get_default_supercell_filename(interface_mode):
if interface_mode == 'phonopy_yaml':
return "phonopy_disp.yaml"
@ -219,6 +225,8 @@ def get_default_supercell_filename(interface_mode):
return "supercell.inp"
elif interface_mode == 'crystal':
return None # supercell.ext can not be parsed by crystal interface.
elif interface_mode == 'dftbp':
return "geo.gen"
else:
return None
@ -244,12 +252,13 @@ def get_default_physical_units(interface_mode):
qe : Ry, au, AMU, Ry/au
siesta : eV, au, AMU, eV/Angstroem
CRYSTAL : eV, Angstrom, AMU, eV/Angstroem
DFTB+ : hartree, au, AMU hartree/au
"""
from phonopy.units import (Wien2kToTHz, AbinitToTHz, PwscfToTHz, ElkToTHz,
SiestaToTHz, VaspToTHz, CP2KToTHz, CrystalToTHz,
Hartree, Bohr)
DftbpToTHz, Hartree, Bohr)
units = {'factor': None,
'nac_factor': None,
@ -305,6 +314,12 @@ def get_default_physical_units(interface_mode):
units['distance_to_A'] = 1.0
units['force_constants_unit'] = 'eV/Angstrom^2'
units['length_unit'] = 'Angstrom'
elif interface_mode == 'dftbp':
units['factor'] = DftbpToTHz
units['nac_factor'] = Hartree * Bohr
units['distance_to_A'] = Bohr
units['force_constants_unit'] = 'hartree/au^2'
units['length_unit'] = 'Angstrom'
return units
@ -334,7 +349,7 @@ def create_FORCE_SETS(interface_mode,
"other files." % force_filenames[0])
if interface_mode in (None, 'vasp', 'abinit', 'elk', 'qe', 'siesta',
'cp2k', 'crystal'):
'cp2k', 'crystal', 'dftbp'):
disp_dataset = parse_disp_yaml(filename=disp_filename)
num_atoms = disp_dataset['natom']
num_displacements = len(disp_dataset['first_atoms'])
@ -415,6 +430,8 @@ def get_force_sets(interface_mode,
from phonopy.interface.cp2k import parse_set_of_forces
elif interface_mode == 'crystal':
from phonopy.interface.crystal import parse_set_of_forces
elif interface_mode == 'dftbp':
from phonopy.interface.dftbp import parse_set_of_forces
else:
return []

219
phonopy/interface/dftbp.py Normal file
View File

@ -0,0 +1,219 @@
# Copyright (C) 2015 Atsushi Togo
# All rights reserved.
#
# This file is part of phonopy.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import sys
import numpy as np
from phonopy.structure.atoms import Atoms
from phonopy.units import dftbpToBohr
from phonopy.file_IO import collect_forces
from phonopy.interface.vasp import (get_scaled_positions_lines, check_forces,
get_drift_forces)
def parse_set_of_forces(num_atoms, forces_filenames, verbose=True):
hook = 'forces :real:2:'
is_parsed = True
force_sets = []
for i, filename in enumerate(forces_filenames):
if verbose:
sys.stdout.write("%d. " % (i + 1))
f = open(filename)
dftbp_forces = collect_forces(f, num_atoms, hook, [0, 1, 2])
if check_forces(dftbp_forces, num_atoms, filename, verbose=verbose):
drift_force = get_drift_forces(dftbp_forces,
filename=filename,
verbose=verbose)
force_sets.append(np.array(dftbp_forces) - drift_force)
else:
is_parsed = False
if is_parsed:
return force_sets
else:
return []
#
# read dftbp-files
#
def read_dftbp(filename):
""" Reads DFTB+ structure files in gen format.
Args:
filename: name of the gen-file to be read
Returns:
atoms: an object of the phonopy.Atoms class, representing the structure
found in filename
"""
infile = open(filename, 'r')
lines = infile.readlines()
# remove any comments
for ss in lines:
if ss.strip().startswith('#'):
lines.remove(ss)
natoms = int(lines[0].split()[0])
symbols = lines[1].split()
if (lines[0].split()[1].lower() == 'f'):
is_scaled = True
scale_pos = 1
scale_latvecs = dftbpToBohr
else:
is_scaled = False
scale_pos = dftbpToBohr
scale_latvecs = dftbpToBohr
# assign positions and expanded symbols
positions = []
expaned_symbols = []
for ii in range(2, natoms+2):
lsplit = lines[ii].split()
expaned_symbols.append(symbols[int(lsplit[1]) - 1])
positions.append([float(ss)*scale_pos for ss in lsplit[2:5]])
# origin is ignored, may be used in future
origin = [float(ss) for ss in lines[natoms+2].split()]
# assign coords of unitcell
cell = []
for ii in range(natoms+3, natoms+6):
lsplit = lines[ii].split()
cell.append([float(ss)*scale_latvecs for ss in lsplit[:3]])
cell = np.array(cell)
if is_scaled:
atoms = Atoms(symbols=expaned_symbols,
cell=cell,
scaled_positions=positions)
else:
atoms = Atoms(symbols=expaned_symbols,
cell=cell,
positions=positions)
return atoms
#
# write dftb+ .gen-file
#
def get_reduced_symbols(symbols):
"""Reduces expanded list of symbols.
Args:
symbols: list containing any chemical symbols as often as
the atom appears in the structure
Returns:
reduced_symbols: any symbols appears only once
"""
reduced_symbols = []
for ss in symbols:
if not (ss in reduced_symbols):
reduced_symbols.append(ss)
return reduced_symbols
def write_dftbp(filename, atoms):
"""Writes DFTB+ readable, gen-formatted structure files
Args:
filename: name of the gen-file to be written
atoms: object containing information about structure
"""
scale_pos = dftbpToBohr
lines = ""
# 1. line, use absolute positions
natoms = atoms.get_number_of_atoms()
lines += str(natoms)
lines += ' S \n'
# 2. line
expaned_symbols = atoms.get_chemical_symbols()
symbols = get_reduced_symbols(expaned_symbols)
lines += ' '.join(symbols) + '\n'
atom_numbers = []
for ss in expaned_symbols:
atom_numbers.append(symbols.index(ss) + 1)
positions = atoms.get_positions()/scale_pos
for ii in range(natoms):
pos = positions[ii]
pos_str = "{:3d} {:3d} {:20.15f} {:20.15f} {:20.15f}\n".format(
ii + 1, atom_numbers[ii], pos[0], pos[1], pos[2])
lines += pos_str
# origin arbitrary
lines +='0.0 0.0 0.0\n'
cell = atoms.get_cell()/scale_pos
for ii in range(3):
cell_str = "{:20.15f} {:20.15f} {:20.15f}\n".format(
cell[ii][0], cell[ii][1], cell[ii][2])
lines += cell_str
outfile = open(filename, 'w')
outfile.write(lines)
def write_supercells_with_displacements(supercell, cells_with_disps, filename="geo.gen"):
"""Writes perfect supercell and supercells with displacements
Args:
supercell: perfect supercell
cells_with_disps: supercells with displaced atoms
filename: root-filename
"""
# original cell
write_dftbp(filename + "S", supercell)
# displaced cells
for ii in range(len(cells_with_disps)):
write_dftbp(filename + "S-{:03d}".format(ii+1), cells_with_disps[ii])

View File

@ -68,4 +68,6 @@ ElkToTHz = sqrt(Hartree*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 154.10794
SiestaToTHz = sqrt(EV/(AMU*Bohr))/Angstrom/(2*pi)/1e12 # [THz] 21.49068
CP2KToTHz = ElkToTHz # CP2K uses a.u. units (Hartree/Bohr)
CrystalToTHz = VaspToTHz
DftbpToTHz = sqrt(Hartree*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 154.10794344
dftbpToBohr = 0.188972598857892E+01
EVAngstromToGPa = EV * 1e21