1# coding: utf-8
2"""
3Work for computing Grüneisen parameters with finite differences on DFPT phonons
4
5WARNING: This code must be tested more carefully.
6"""
7import numpy as np
8
9from abipy.core.structure import Structure
10from .works import Work, PhononWork
11
12
13class GruneisenWork(Work):
14    """
15    This work computes the Grüneisen parameters (derivative of frequencies wrt Volume)
16    using finite differences and the phonons obtained with the DFPT part of Abinit.
17    The Anaddb input file needed to compute Grüneisen parameters will be generated
18    in the outdata directory of the flow.
19
20    It is necessary to run three DFPT phonon calculations.
21    One is calculated at the equilibrium volume and the remaining two are calculated
22    at the slightly larger volume and smaller volume than the equilibrium volume.
23    The unitcells at these volumes have to be fully relaxed under the constraint of each volume.
24    This Work automates the entire procedure starting from an input for GS calculations.
25    """
26
27    @classmethod
28    def from_gs_input(cls, gs_inp, voldelta, ngqpt, tolerance=None, with_becs=False,
29                      ddk_tolerance=None, workdir=None, manager=None):
30        """
31        Build the work from an |AbinitInput| representing a GS calculations.
32
33        Args:
34            gs_inp: |AbinitInput| representing a GS calculation in the initial unit cell.
35            voldelta: Absolute increment for unit cell volume. The three volumes are:
36                [v0 - voldelta, v0, v0 + voldelta] where v0 is taken from gs_inp.structure.
37            ngqpt: three integers defining the q-mesh for phonon calculations.
38            tolerance: dict {"varname": value} with the tolerance to be used in the phonon run.
39                None to use AbiPy default.
40            with_becs: Activate calculation of Electric field and Born effective charges.
41            ddk_tolerance: dict {"varname": value} with the tolerance used in the DDK run if with_becs.
42                None to use AbiPy default.
43        """
44        new = cls(workdir=workdir, manager=manager)
45        new.ngqpt = np.reshape(ngqpt, (3,))
46        new.with_becs = with_becs
47        new.ddk_tolerance = ddk_tolerance
48        new.tolerance = tolerance
49
50        if any(gs_inp["ngkpt"] % new.ngqpt != 0):
51            raise ValueError("Kmesh and Qmesh must be commensurate.\nGot ngkpt: `%s`\nand ngqpt: `%s`" % (
52                             str(gs_inp["ngkpt"]), str(new.ngqpt)))
53
54        # Build three tasks for structural optimization at constant volume.
55        v0 = gs_inp.structure.volume
56        if voldelta <= 0:
57            raise ValueError("voldelta must be > 0 but got %s" % voldelta)
58        volumes = [v0 - voldelta, v0, v0 + voldelta]
59        if any(v <= 0 for v in volumes):
60            raise ValueError("volumes must be > 0 but got %s" % str(volumes))
61
62        # Keep a copy of the GS input that will be used to generate the Phonon Works
63        new.gs_inp = gs_inp.deepcopy()
64
65        new.relax_tasks = []
66        for vol in volumes:
67            # Build new structure
68            new_lattice = gs_inp.structure.lattice.scale(vol)
69            new_structure = Structure(new_lattice, gs_inp.structure.species, gs_inp.structure.frac_coords)
70            new_input = gs_inp.new_with_structure(new_structure)
71            # Set variables for structural optimization at constant volume.
72            new_input.pop_tolerances()
73            new_input.set_vars(optcell=3, ionmov=3, tolvrs=1e-10, toldff=1.e-6)
74            new_input.set_vars_ifnotin(ecutsm=0.5, dilatmx=1.05)
75            t = new.register_relax_task(new_input)
76            new.relax_tasks.append(t)
77
78        return new
79
80    def on_all_ok(self):
81        """
82        This method is called once the `Work` is completed i.e. when all the tasks have reached status S_OK.
83        """
84        self.add_phonon_works_and_build()
85        return super().on_all_ok()
86
87    def add_phonon_works_and_build(self):
88        """
89        Get relaxed structures from the tasks, build Phonons works with new structures.
90        Add works to the flow and build new directories.
91        """
92        ddb_paths, struct_middle, middle_idx = [], None, None
93        relaxed_structs = []
94
95        for i, task in enumerate(self.relax_tasks):
96            relaxed_structure = task.get_final_structure()
97            relaxed_structs.append(relaxed_structure)
98            if i == len(self.relax_tasks) // 2:
99                middle_idx, struct_middle = i, relaxed_structure
100
101            # work to compute phonons with new structure.
102            gsinp_vol = self.gs_inp.new_with_structure(relaxed_structure)
103            work = PhononWork.from_scf_input(gsinp_vol, self.ngqpt, is_ngqpt=True, tolerance=self.tolerance,
104                                             with_becs=self.with_becs, ddk_tolerance=self.ddk_tolerance)
105            # Add it to the flow.
106            self.flow.register_work(work)
107
108            ddb_paths.append(work.outdir.path_in("out_DDB"))
109
110        # Write Anaddb input file to compute Gruneisen parameters in flow.outdata.
111        from abipy.abio.inputs import AnaddbInput
112        anaddb_inp = AnaddbInput.phbands_and_dos(struct_middle, self.ngqpt, nqsmall=20, ndivsm=20,
113                                                 chneut=1 if self.with_becs else 0,
114                                                 dipdip=1 if self.with_becs else 0,
115                                                 lo_to_splitting=self.with_becs,
116                                                 comment="Anaddb input file for Grunesein parameters")
117        # Add DDB files for Gruns
118        anaddb_inp["gruns_nddbs"] = len(ddb_paths)
119        anaddb_inp["gruns_ddbs"] = "\n" + "\n".join('"%s"' % p for p in ddb_paths)
120        in_path = self.flow.outdir.path_in("anaddb_gruns.in")
121        anaddb_inp.write(in_path)
122
123        files_file = []
124        app = files_file.append
125        app(in_path)                                        # 1) Path of the input file
126        app(self.flow.outdir.path_in("anaddb_gruns.out"))   # 2) Path of the output file
127        app(ddb_paths[middle_idx])                          # 3) Input derivative database (not used if Gruns)
128        for i in range(4):
129            app("FOOBAR")
130
131        with open(self.flow.outdir.path_in("anaddb_gruns.files"), "wt") as fh:
132            fh.write("\n".join(files_file))
133
134        #with_ebands = False
135        #if with_ebands:
136        #    bands_work = Work(manager=self.manager)
137        #    for i, structure in enumerate(relaxed_structs):
138        #        nscf_inp = self.gs_inp.new_with_structure(structure)
139        #        nscf_inp.pop_tolerances()
140        #        nscf_inp.set_kpath(ndivsm, kptbounds=None, iscf=-2)
141        #        nscf_inp["tolwfr"] = 1e-18
142        #        work.register_nscf_task(nscf_inp, deps={self.relax_tasks[i]: "DEN"})
143
144        #    # Add it to the flow.
145        #    self.flow.register_work(work)
146
147        # Allocate new works and update the pickle database.
148        self.flow.allocate()
149        self.flow.build_and_pickle_dump()
150