Order N Fourier Interpolation

This is where lone-q and LID happens.

The core part of the code, Fourier interpolation, needs to be abstracted from this worker class so that the codes are better structured and the FI can be more easily accessible for other uses.

Moreover, we now have: 1. A YAML based data storage system used in this code and the frozen phonon code 2. A tar-archive based system used in BID/HSBID These needs to be replaced by a HDF5 based storage system to improve flexibility, safety and IO speed.

Currently, the code does 2 jobs. The first is finding irreducible Q for a given system size and order based on symmetry, creating lone-q or LID runs using frozen_phonon_n or frozen_phonon_sym module, it also collects the dynamic tensors from those runs and rotated them to the full BZ. Second, it performs Fourier transformation and pack the force constants within the WS cell. These two can be completed 2 seperate features and keeping them seperate also helps modules to utilize the features without redundant methods dangling around.

The proposed feature consists of a core Fourier interpolation class, from which a class can be derived to add in symmetry operations allowing dynamic tensors from the irreducible Q to be rotated to the full BZ, which can not only be used in the LID approach but also the BID approaches. Then a class can be derived to create and collect LID calculations.

import numpy as np
from fourier_interp_Dqn import fourier_interp_Dqn
from coordinate_tools import coord_to_string as c2s

poscar = """\
graphene
   1.00000000000000
     2.1217113285639000    1.2249706066897641    0.0000000000000000
     2.1217113285639000   -1.2249706066897641    0.0000000000000000
     0.0000000000000000    0.0000000000000000  -15.0000000000000000
   2
Direct
  0.00000000  0.00000000  0.0000000000000000 C:p
  0.33333333  0.33333333  0.0000000000000000 C:p
""".strip()

f = fourier_interp_Dqn(
    pos=poscar,
    nq="2 -1 0  -1 2 0  0 0 1",
    order=3,
    pgn="D6h",
    pg_loc="1/3 1/3 0"
)

f.generate_qmesh_runs(delta=0.01, dry_run=True)

print "A dictionary containing the LID directories and the corresponding Q"
print "-" * 10
print "{:20s} {:20s}".format("Directory", "Q")
for k, v in f.fi_info_dict["qmesh"].items():
    print "{:20s} {:20s}".format(k, c2s(v))

A dictionary containing the LID directories and the corresponding Q
----------
Directory            Q
q1__000__000         ((0, 0, 0), (0, 0, 0))
q3__210__210         ((2/3, 1/3, 0), (2/3, 1/3, 0))
q3__210__120         ((2/3, 1/3, 0), (1/3, 2/3, 0))

In order to collect the LID results and perform Fourier interpolation:

f.gather_qmesh_runs(
    dft_force_io=vasp.get_forces,
    # Provide a callable interface for LID to extract forces from first principle calculation
)
f.fourier_transform()

Then the force constants before the WS packing can be found in f.phin and the one after the packing f.phin_WS.

If one desire to merely use the FI feature, then after instantiated the object, one can provide the dynamic tensor for either the irreducible Q or the full BZ to the attribute f.Dqn in the form of a dictionary and run f.fourier_transform().