Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Wien2k support have been added #323

Closed
wants to merge 16 commits into from
Closed
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7,436 changes: 7,436 additions & 0 deletions example/Si-wien2k/FORCES_FC3

Large diffs are not rendered by default.

32 changes: 32 additions & 0 deletions example/Si-wien2k/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
This is the example of silicon calculation. The Wien2k calculation was made to
obtain forces with 2x2x2 k-point mesh for the supercell, PBE, and the lattice
parameters in Si.struct.

`phono3py_disp.yaml` is generated by

```
% phono3py --wien2k -d --dim 2 2 2 -c Si.struct
```

To create `fc3.hdf5` and `fc2.hdf5`,

```
% phono3py-load
```

Using 11x11x11 sampling mesh, lattice thermal conductivity is calculated by

```
% phono3py-load --mesh 11 11 11 --br --ts 300
```

`kappa-m111111.hdf5` is written as the result. The lattice thermal conductivity
is calculated as 119.1 W/m-K at 300 K. This becomes, with 19x19x19 sampling
mesh, 128.1 W/m-K.

The .scf files for supercells are found in `supercell_out.zip`. If phono3py
is properly installed, the following command should work.

```
% phono3py --cf3 supercell_out/Si.structS-{00001..00111}/Si.structS-{00001..00111}.scf
```
53 changes: 53 additions & 0 deletions example/Si-wien2k/Si.struct
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
Si
P LATTICE,NONEQUIV.ATOMS: 8
MODE OF CALC=RELA unit=ang
10.329744 10.329744 10.329744 90.000000 90.000000 90.000000
ATOM -1: X=0.87500000 Y=0.87500000 Z=0.87500000
MULT= 1 ISPLIT= 8
SiSi1 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000
0.0000000 1.0000000 0.0000000
0.0000000 0.0000000 1.0000000
ATOM 2: X=0.87500000 Y=0.37500000 Z=0.37500000
MULT= 1 ISPLIT= 8
SiSi2 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 3: X=0.37500000 Y=0.87500000 Z=0.37500000
MULT= 1 ISPLIT= 8
SiSi3 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 4: X=0.37500000 Y=0.37500000 Z=0.87500000
MULT= 1 ISPLIT= 8
SiSi4 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 5: X=0.12500000 Y=0.12500000 Z=0.12500000
MULT= 1 ISPLIT= 8
SiSi5 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 6: X=0.12500000 Y=0.62500000 Z=0.62500000
MULT= 1 ISPLIT= 8
SiSi6 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 7: X=0.62500000 Y=0.12500000 Z=0.62500000
MULT= 1 ISPLIT= 8
SiSi7 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM 8: X=0.62500000 Y=0.62500000 Z=0.12500000
MULT= 1 ISPLIT= 8
SiSi8 NPT= 781 R0=0.00010000 RMT= 2.1100 Z: 14.000
LOCAL ROT MATRIX: 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0 NUMBER OF SYMMETRY OPERATIONS
587 changes: 587 additions & 0 deletions example/Si-wien2k/phono3py_disp.yaml

Large diffs are not rendered by default.

39 changes: 25 additions & 14 deletions phono3py/conductivity/direct_solution.py
Original file line number Diff line number Diff line change
@@ -1881,6 +1881,7 @@ def diagonalize_collision_matrix(
assert size == shape[1]

solver = select_colmat_solver(pinv_solver)
f_norm = np.linalg.norm(collision_matrices[i_sigma, i_temp].ravel())
trace = np.trace(collision_matrices[i_sigma, i_temp].reshape(size, size))

# [1] dsyev: safer and slower than dsyevd and smallest memory usage
@@ -1906,31 +1907,41 @@ def diagonalize_collision_matrix(
) # only diagonalization
elif solver == 3: # np.linalg.eigh depends on dsyevd.
if log_level:
sys.stdout.write("Diagonalize by np.linalg.eigh ")
sys.stdout.flush()
print("Diagonalize by np.linalg.eigh ", end="")
col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, col_mat[:] = np.linalg.eigh(col_mat)

info_str = ""
elif solver == 4: # fully scipy dsyev
if log_level:
sys.stdout.write("Diagonalize by scipy.linalg.lapack.dsyev ")
sys.stdout.flush()
print("Diagonalize by scipy.linalg.lapack.dsyev ", end="")
import scipy.linalg

col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, _, info = scipy.linalg.lapack.dsyev(col_mat.T, overwrite_a=1)
_col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size).copy()
print(_col_mat.shape)
print(_col_mat.flags, _col_mat.dtype)
sym_str = f"sym={scipy.linalg.issymmetric(_col_mat)} "
w, _, info = scipy.linalg.lapack.dsyev(_col_mat.T, overwrite_a=1)
collision_matrices[i_sigma, i_temp].reshape(size, size)[:, :] = _col_mat[:, :]
info_str = f"info={info} {sym_str} "
elif solver == 5: # fully scipy dsyevd
if log_level:
sys.stdout.write("Diagnalize by scipy.linalg.lapack.dsyevd ")
sys.stdout.flush()
print("Diagnalize by scipy.linalg.lapack.dsyevd ", end="")
import scipy.linalg

col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, _, info = scipy.linalg.lapack.dsyevd(col_mat.T, overwrite_a=1)

_col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size).copy()
print(_col_mat.shape)
print(_col_mat.flags, _col_mat.dtype)
sym_str = f"sym={scipy.linalg.issymmetric(_col_mat)} "
w, _, info = scipy.linalg.lapack.dsyevd(_col_mat.T, overwrite_a=1)
collision_matrices[i_sigma, i_temp].reshape(size, size)[:, :] = _col_mat[:, :]
info_str = f"info={info} {sym_str}"
if log_level:
print(f"sum={w.sum():<.1e} d={trace - w.sum():<.1e} ", end="")
print("[%.3fs]" % (time.time() - start))
print(f"[{time.time() - start:.3f}s]")
print(
" "
f"{info_str}"
f"f_norm={f_norm:<.1e} tr={w.sum():<.1e} d={trace - w.sum():<.1e} ",
)
sys.stdout.flush()

return w
23 changes: 17 additions & 6 deletions phono3py/cui/create_force_sets.py
Original file line number Diff line number Diff line change
@@ -63,6 +63,7 @@
Phono3pyYaml,
displacements_yaml_lines_type1,
)
from phono3py.interface.wien2k import get_fc3_calc_dataset_wien2k


def create_FORCES_FC3_and_FORCES_FC2(
@@ -385,12 +386,22 @@ def _get_force_sets_fc3(
): # type-1
calc_dataset = {"forces": []}
else: # type-2
calc_dataset = get_calc_dataset(
interface_mode,
num_atoms,
force_filenames,
verbose=(log_level > 0),
)
if interface_mode == "wien2k":
calc_dataset = get_fc3_calc_dataset_wien2k(
force_filenames,
supercell,
disp_dataset,
# wien2k_P1_mode=wien2k_P1_mode,
# symmetry_tolerance=symmetry_tolerance,
verbose=(log_level > 0),
)
else:
calc_dataset = get_calc_dataset(
interface_mode,
num_atoms,
force_filenames,
verbose=(log_level > 0),
)
force_sets = calc_dataset["forces"]
if "points" in calc_dataset:
if filename := check_agreements_of_displacements(
3 changes: 3 additions & 0 deletions phono3py/cui/create_supercells.py
Original file line number Diff line number Diff line change
@@ -116,6 +116,9 @@ def create_phono3py_supercells(
additional_info = get_additional_info_to_write_supercells(
interface_mode, phono3py.supercell_matrix
)

additional_info["supercell_matrix"] = phono3py.supercell_matrix

write_supercells_with_displacements(
interface_mode,
phono3py.supercell,
3 changes: 1 addition & 2 deletions phono3py/interface/calculator.py
Original file line number Diff line number Diff line change
@@ -50,8 +50,7 @@
# 'help': "Invoke Siesta mode"}},
"turbomole": {"option": {"name": "--turbomole", "help": "Invoke TURBOMOLE mode"}},
"vasp": {"option": {"name": "--vasp", "help": "Invoke Vasp mode"}},
# 'wien2k': {'option': {'name': "--wien2k",
# 'help': "Invoke Wien2k mode"}},
"wien2k": {"option": {"name": "--wien2k", "help": "Invoke Wien2k mode"}},
}


61 changes: 61 additions & 0 deletions phono3py/interface/wien2k.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
"""ALM interface for force constants calculation."""

# Copyright (C) 2016 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# 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.

from phono3py.phonon3.dataset import get_displacements_fc3


def get_fc3_calc_dataset_wien2k(
force_filenames,
supercell,
disp_dataset,
wien2k_P1_mode=False,
symmetry_tolerance=None,
verbose=False,
):
"""Read Wien2k output files and parse force sets."""
from phonopy.interface.wien2k import parse_set_of_forces

# disps, _ = get_displacements_and_forces_fc3(disp_dataset)
disps = get_displacements_fc3(disp_dataset)
force_sets = parse_set_of_forces(
disps,
force_filenames,
supercell,
wien2k_P1_mode=wien2k_P1_mode,
symmetry_tolerance=symmetry_tolerance,
verbose=verbose,
)
return {"forces": force_sets}
55 changes: 55 additions & 0 deletions phono3py/phonon3/dataset.py
Original file line number Diff line number Diff line change
@@ -98,6 +98,61 @@ def get_displacements_and_forces_fc3(disp_dataset):
raise RuntimeError("disp_dataset doesn't contain correct information.")


def get_displacements_fc3(disp_dataset):
"""Return displacements and forces from disp_dataset.
Note
----
Dipslacements and forces of all atoms in supercells are returned.
Parameters
----------
disp_dataset : dict
Displacement dataset.
Returns
-------
displacements : ndarray
Displacements of all atoms in all supercells.
shape=(snapshots, supercell atoms, 3), dtype='double', order='C'
forces : ndarray or None
Forces of all atoms in all supercells.
shape=(snapshots, supercell atoms, 3), dtype='double', order='C'
None is returned when forces don't exist.
"""
if "first_atoms" in disp_dataset:
natom = disp_dataset["natom"]
ndisp = len(disp_dataset["first_atoms"])
for disp1 in disp_dataset["first_atoms"]:
ndisp += len(disp1["second_atoms"])
displacements = np.zeros((ndisp, natom, 3), dtype="double", order="C")
indices = []
count = 0
for disp1 in disp_dataset["first_atoms"]:
indices.append(count)
displacements[count, disp1["number"]] = disp1["displacement"]
count += 1

for disp1 in disp_dataset["first_atoms"]:
for disp2 in disp1["second_atoms"]:
if "included" in disp2:
if disp2["included"]:
indices.append(count)
else:
indices.append(count)
displacements[count, disp1["number"]] = disp1["displacement"]
displacements[count, disp2["number"]] = disp2["displacement"]
count += 1

return np.array(displacements[indices], dtype="double", order="C")

elif "forces" in disp_dataset and "displacements" in disp_dataset:
return disp_dataset["displacements"]
else:
raise RuntimeError("disp_dataset doesn't contain correct information.")


def forces_in_dataset(dataset: Optional[dict]) -> bool:
"""Return whether forces in dataset or not."""
if dataset is None: