Merge branch 'hotfix/v.1.4.2'

ChangeLog.md

@@ -1,3 +1,17 @@
+# Ver. 1.4.2 (2022-09-28)
+
+## Fixes
+
+- Fix a bug in group velocity when NONANALYTIC=3
+
+
+# Ver. 1.4.1 (2022-06-22)
+
+## Fixes
+
+- Fix a bug in xml writer in the previous version
+
+
 # Ver. 1.4.0 (2022-06-21)
 
 ## New

README.md

@@ -3,7 +3,7 @@
 [![License][license-image]][license-url]
 [![Doc status][docs-image]][docs-url]
 
-### Version 1.4.0
+### Version 1.4.2
 ![alt ALAMODE](./docs/img/alamode.png)
 
 

anphon/conductivity.cpp

@@ -124,7 +124,6 @@ void Conductivity::setup_kappa()
 
     phonon_velocity->get_phonon_group_velocity_mesh_mpi(*dos->kmesh_dos,
                                                         system->lavec_p,
-                                                        fcs_phonon->fc2_ext,
                                                         vel);
     if (mympi->my_rank == 0) {
         for (i = 0; i < nk; ++i) {

anphon/phonon_velocity.cpp

@@ -13,6 +13,7 @@ or http://opensource.org/licenses/mit-license.php for information.
 #include "constants.h"
 #include "dynamical.h"
 #include "error.h"
+#include "ewald.h"
 #include "fcs_phonon.h"
 #include "kpoint.h"
 #include "mathfunctions.h"
@@ -77,6 +78,7 @@ void PhononVelocity::get_phonon_group_velocity_bandstructure(const KpointBandStr
                                                              const double lavec_p[3][3],
                                                              const double rlavec_p[3][3],
                                                              const std::vector<FcsClassExtent> &fc2_ext_in,
+                                                             const std::vector<FcsClassExtent> &fc2_without_dipole,
                                                              double **phvel_out) const
 {
     unsigned int i;
@@ -130,12 +132,19 @@ void PhononVelocity::get_phonon_group_velocity_bandstructure(const KpointBandStr
             rotvec(xk_shift[idiff], xk_shift[idiff], lavec_p, 'T');
             for (i = 0; i < 3; ++i) xk_shift[idiff][i] /= 2.0 * pi;
 
-            dynamical->eval_k(xk_shift[idiff],
-                              kpoint_bs_in->kvec_na[ik],
-                              fc2_ext_in,
-                              omega_shift[idiff],
-                              evec_tmp, false);
-
+            if (dynamical->nonanalytic == 3) {
+                dynamical->eval_k_ewald(xk_shift[idiff],
+                                        kpoint_bs_in->kvec_na[ik],
+                                        fc2_without_dipole,
+                                        omega_shift[idiff],
+                                        evec_tmp, false);
+            } else {
+                dynamical->eval_k(xk_shift[idiff],
+                                  kpoint_bs_in->kvec_na[ik],
+                                  fc2_ext_in,
+                                  omega_shift[idiff],
+                                  evec_tmp, false);
+            }
         }
 
         for (i = 0; i < n; ++i) {
@@ -159,7 +168,6 @@ void PhononVelocity::get_phonon_group_velocity_bandstructure(const KpointBandStr
 
 void PhononVelocity::get_phonon_group_velocity_mesh(const KpointMeshUniform &kmesh_in,
                                                     const double lavec_p[3][3],
-                                                    const std::vector<FcsClassExtent> &fc2_ext_in,
                                                     const bool irreducible_only,
                                                     double ***phvel3_out) const
 {
@@ -202,7 +210,6 @@ void PhononVelocity::get_phonon_group_velocity_mesh(const KpointMeshUniform &kme
 
 void PhononVelocity::get_phonon_group_velocity_mesh_mpi(const KpointMeshUniform &kmesh_in,
                                                         const double lavec_p[3][3],
-                                                        const std::vector<FcsClassExtent> &fc2_ext_in,
                                                         double ***phvel3_out) const
 {
     // This routine computes the group velocities for the given uniform k mesh
@@ -469,13 +476,21 @@ void PhononVelocity::phonon_vel_k(const double *xk_in,
 
         for (idiff = 0; idiff < ndiff; ++idiff) {
 
-            dynamical->eval_k(xk_shift[idiff],
-                              kvec_na_tmp[0],
-                              fcs_phonon->fc2_ext,
-                              omega_shift[idiff],
-                              evec_tmp,
-                              false);
-
+            if (dynamical->nonanalytic == 3) {
+                dynamical->eval_k_ewald(xk_shift[idiff],
+                                        kvec_na_tmp[idiff],
+                                        ewald->fc2_without_dipole,
+                                        omega_shift[idiff],
+                                        evec_tmp,
+                                        false);
+            } else {
+                dynamical->eval_k(xk_shift[idiff],
+                                  kvec_na_tmp[idiff],
+                                  fcs_phonon->fc2_ext,
+                                  omega_shift[idiff],
+                                  evec_tmp,
+                                  false);
+            }
         }
 
         for (j = 0; j < n; ++j) {

anphon/phonon_velocity.h

@@ -35,13 +35,11 @@ class PhononVelocity : protected Pointers {
 
     void get_phonon_group_velocity_mesh(const KpointMeshUniform &kmesh_in,
                                         const double lavec_p[3][3],
-                                        const std::vector<FcsClassExtent> &fc2_ext_in,
                                         const bool irreducible_only,
                                         double ***phvel3_out) const;
 
     void get_phonon_group_velocity_mesh_mpi(const KpointMeshUniform &kmesh_in,
                                             const double lavec_p[3][3],
-                                            const std::vector<FcsClassExtent> &fc2_ext_in,
                                             double ***phvel3_out) const;
 
     void calc_phonon_velmat_mesh(std::complex<double> ****velmat_out) const;
@@ -50,6 +48,7 @@ class PhononVelocity : protected Pointers {
                                                  const double lavec_p[3][3],
                                                  const double rlavec_p[3][3],
                                                  const std::vector<FcsClassExtent> &fc2_ext_in,
+                                                 const std::vector<FcsClassExtent> &fc2_without_dipole,
                                                  double **phvel_out) const;
 
     void velocity_matrix_analytic(const double *xk_in,

anphon/write_phonons.cpp

@@ -15,6 +15,7 @@ or http://opensource.org/licenses/mit-license.php for information.
 #include "dielec.h"
 #include "dynamical.h"
 #include "error.h"
+#include "ewald.h"
 #include "fcs_phonon.h"
 #include "gruneisen.h"
 #include "kpoint.h"
@@ -779,6 +780,7 @@ void Writes::write_phonon_vel() const
                                                              system->lavec_p,
                                                              system->rlavec_p,
                                                              fcs_phonon->fc2_ext,
+                                                             ewald->fc2_without_dipole,
                                                              phvel_bs);
 
     ofs_vel << "# k-axis, |Velocity| [m / sec]" << std::endl;
@@ -834,7 +836,6 @@ void Writes::write_phonon_vel_all() const
 
     phonon_velocity->get_phonon_group_velocity_mesh(*dos->kmesh_dos,
                                                     system->lavec_p,
-                                                    fcs_phonon->fc2_ext,
                                                     false,
                                                     phvel_xyz);
     unsigned int ik, is;

docs/source/conf.py

@@ -59,7 +59,7 @@
 # The short X.Y version.
 version = '1.4'
 # The full version, including alpha/beta/rc tags.
-release = '1.4.0'
+release = '1.4.2'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

include/version.h

@@ -12,4 +12,4 @@
 
 #include <string>
 
-static const std::string ALAMODE_VERSION = "1.4.0";
+static const std::string ALAMODE_VERSION = "1.4.2";

tools/interface/LAMMPS.py

@@ -380,24 +380,36 @@ def x_fractional(self):
     @staticmethod
     def _get_coordinate_and_force_lammps(lammps_dump_file):
 
-        add_flag = False
+        add_flag = None
         ret = []
 
         with open(lammps_dump_file) as f:
             for line in f:
-                if "ITEM:" in line and "ITEM: ATOMS id xu yu zu fx fy fz" not in line:
-                    add_flag = False
+                #if "ITEM:" in line and "ITEM: ATOMS id xu yu zu fx fy fz" not in line:
+                #    add_flag = False
+                #    continue
+                if "ITEM: ATOMS id xu yu zu fx fy fz" in line:
+                    add_flag = "id.xu"
                     continue
-                elif "ITEM: ATOMS id xu yu zu fx fy fz" in line:
-                    add_flag = True
+                elif "ITEM: ATOMS element xu yu zu fx fy fz" in line:
+                    add_flag = "element.xu"
+                    id_ = 0
                     continue
 
-                if add_flag:
-                    if line.strip():
-                        entries = line.strip().split()
-                        data_atom = [int(entries[0]),
-                                     [float(t) for t in entries[1:4]],
-                                     [float(t) for t in entries[4:]]]
+                if add_flag is not None:
+                    if add_flag == "id.xu":
+                        if line.strip():
+                            entries = line.strip().split()
+                            data_atom = [int(entries[0]),
+                                        [float(t) for t in entries[1:4]],
+                                        [float(t) for t in entries[4:]]]
+                    elif add_flag == "element.xu":
+                        if line.strip():
+                            entries = line.strip().split()
+                            id_ += 1
+                            data_atom = [int(id_), 
+                                        [float(t) for t in entries[1:4]],
+                                        [float(t) for t in entries[4:]]]
 
                         ret.append(data_atom)
 

tools/interface/QE.py

@@ -1024,13 +1024,13 @@ def _get_borninfo_phout(phout_file):
         borncharge = []
 
         search_tag1 = "Dielectric constant in cartesian axis"
-
+        search_tags = ["Px", "Py", "Pz", "Ex", "Ey", "Ez"]
         f = open(phout_file, 'r')
         line = f.readline()
 
         found_tag1 = False
         found_tag2 = False
-
+        
         while line:
             if search_tag1 in line:
                 found_tag1 = True
@@ -1039,7 +1039,7 @@ def _get_borninfo_phout(phout_file):
                     line = f.readline()
                     dielec.extend([float(t) for t in line.strip().split()[1:4]])
 
-            if "Px" in line or "Py" in line or "Pz" in line:
+            if any(s in line for s in search_tags):
                 found_tag2 = True
                 borncharge.extend(float(t) for t in line.strip().split()[2:5])
 
@@ -1051,9 +1051,9 @@ def _get_borninfo_phout(phout_file):
                   "in %s" % phout_file, file=sys.stderr)
             return None
 
-        nat = len(borncharge) // 9
+        nat2 = len(borncharge) // 9
         dielec = np.reshape(np.array(dielec[9:]), (3, 3))
-        borncharge = np.reshape(np.array(borncharge), (nat, 3, 3))
-        return dielec, borncharge
+        borncharge = np.reshape(np.array(borncharge), (nat2, 3, 3))
+        return dielec, borncharge[:nat2 // 2, :, :]