Merge pull request #185 from gpetretto/workflow_refactoring

small updates and fixes

abipy/abio/factories.py

@@ -1369,7 +1369,7 @@ def dfpt_from_gsinput(gs_inp, ph_ngqpt=None, qpoints=None, do_ddk=True, do_dde=T
         do_strain: If True inputs for the strain perturbations will be included.
         do_dte: If True inputs for the non-linear perturbations will be included. The phonon non-linear perturbations
             will be included only if a phonon calculation at gamma is present. The caller is responsible for
-            adding it.
+            adding it. Automatically sets with_dde=True.
         ph_tol: a dictionary with a single key defining the type of tolerance used for the phonon calculations and
             its value. Default: {"tolvrs": 1.0e-10}.
         ddk_tol: a dictionary with a single key defining the type of tolerance used for the DDK calculations and
@@ -1400,6 +1400,9 @@ def dfpt_from_gsinput(gs_inp, ph_ngqpt=None, qpoints=None, do_ddk=True, do_dde=T
     if do_dde:
         do_ddk = True
 
+    if do_dte:
+        do_dde = True
+
     multi = MultiDataset.from_inputs([gs_inp])
     multi[0].add_tags(SCF)
 
@@ -1435,7 +1438,7 @@ def dfpt_from_gsinput(gs_inp, ph_ngqpt=None, qpoints=None, do_ddk=True, do_dde=T
         gs_inp_copy.set_vars(nband=nband)
         gs_inp_copy.pop('nbdbuf', None)
         multi_dte = gs_inp_copy.make_dte_inputs(phonon_pert=do_phonons and has_gamma,
-                                           skip_permutations=skip_dte_permutations, manager=manager)
+                                                skip_permutations=skip_dte_permutations, manager=manager)
         multi_dte.add_tags([DTE, DFPT])
         multi.extend(multi_dte)
 

abipy/abio/inputs.py

@@ -2758,10 +2758,23 @@ with the Abinit version you are using. Please contact the AbiPy developers.""" %
             anaddb_input["instrflag"] = 1
 
         if dte:
-            anaddb_input.set_vars(nlflag=1,
+            ramansr = 0
-                                  ramansr=1,
+            alphon = 0
-                                  alphon=1,
+            prtmbm = 0
-                                  prtmbm=1)
+
+            # if there are phonons at gamma
+            if ngqpt and (not q1shft or np.allclose(q1shft, [0, 0, 0])):
+                nlflag = 1
+                ramansr = 1
+                alphon = 1
+                prtmbm = 1
+            else:
+                nlflag = 3
+
+            anaddb_input.set_vars(nlflag=nlflag,
+                                  ramansr=ramansr,
+                                  alphon=alphon,
+                                  prtmbm=prtmbm)
 
         anaddb_args = [] if anaddb_args is None else anaddb_args
         anaddb_kwargs = {} if anaddb_kwargs is None else anaddb_kwargs

abipy/dfpt/ddb.py

@@ -1898,7 +1898,7 @@ class DielectricTensorGenerator(Has_Structure):
         Return: |matplotlib-Figure|
         """
         if w_max is None:
-            w_max = np.max(self.phfreqs) * phfactor_ev2units(units) + gamma_ev * 10
+            w_max = (np.max(self.phfreqs) + gamma_ev * 10) * phfactor_ev2units(units)
 
         wmesh = np.linspace(w_min, w_max, num, endpoint=True)
         t = np.zeros((num, 3, 3), dtype=complex)

abipy/dfpt/phtk.py

@@ -5,6 +5,7 @@ Phonon Toolkit: This module gathers low-level tools to operate on phonons.
 from __future__ import print_function, division, absolute_import  # unicode_literals,
 
 import warnings
+import sys
 import numpy as np
 import abipy.core.abinit_units as abu
 

abipy/dfpt/vsound.py

@@ -38,7 +38,7 @@ class SoundVelocity(Has_Structure, NotebookWriter):
             qpts: array with shape (len(directions), num_points, 3) with the coordinates of
                 the qpoints in fractional used to fit the phonon frequencies.
         """
-        self.directions = directions
+        self.directions = np.array(directions)
         self.sound_velocities = np.array(sound_velocities)
         self.mode_types = mode_types
         self.labels = labels
@@ -239,8 +239,11 @@ class SoundVelocity(Has_Structure, NotebookWriter):
 
             cart_versor = qpt_cart_coords[end -1] / np.linalg.norm(qpt_cart_coords[end -1])
             for k in range(3):
-                slope, se, _, _ = np.linalg.lstsq(qpt_cart_norms[start:end][:, np.newaxis],
+                start_fit = 0
-                                                  acoustic_freqs[:, k] * eV_to_Ha, rcond=None)
+                if ignore_neg_freqs and first_positive_freq_ind > 1:
+                    start_fit = first_positive_freq_ind
+                slope, se, _, _ = np.linalg.lstsq(qpt_cart_norms[start+start_fit:end][:, np.newaxis],
+                                                  acoustic_freqs[start_fit:, k] * eV_to_Ha, rcond=None)
                 sv.append(slope[0] * abu.velocity_at_to_si)
 
                 # identify the type of the mode (longitudinal/transversal) based on the