Merge branch 'develop' of https://github.com/ttadano/alamode into dev…

…elop
ttadano · Feb 23, 2018 · 9637704 · 9637704
2 parents 026ba26 + 4b85010
commit 9637704
Show file tree

Hide file tree

Showing 52 changed files with 26,700 additions and 10,035 deletions.
diff --git a/ChangeLog.md b/ChangeLog.md
@@ -1,3 +1,26 @@
+# Ver. 1.0.2 (2018-1-29)
+
+## New
+
+- Phonon band connection by eigenvector similarity (``BCONNECT`` tag) 
+- New option to turn on/off the symmetrization of Born effective charge (``BORNSYM`` tag).
+
+## Changes
+
+- Improve the performance of the "suggest" mode for hexagonal systems
+- Use \<unorderd_set\> instead of \<set\> for better performance
+
+## Fix
+
+- Fix a bug in the symmetrization of the Born effective charge
+
+
+# Ver. 1.0.1 (2017-11-21)
+
+## Fix
+- Fixed a minor issue in the previous version
+
+
 # Ver. 1.0.0 (2017-11-21)
 
 ## New

diff --git a/README.md b/README.md
@@ -3,7 +3,7 @@
 [![License][license-image]][license-url]
 [![Doc status][docs-image]][docs-url]
 
-### Version 1.0.0
+### Version 1.0.2
 ![alt ALAMODE](./docs/img/alamode.png)
 
 

diff --git a/alm/constraint.cpp b/alm/constraint.cpp
@@ -846,7 +846,7 @@ void Constraint::translational_invariance()
 
         list_found.clear();
 
-        for (auto p = fcs->fc_table[order].begin(); p != fcs->fc_table[order].end(); ++p) {
+        for (auto p = fcs->fc_table[order].cbegin(); p != fcs->fc_table[order].cend(); ++p) {
             for (i = 0; i < order + 2; ++i) {
                 ind[i] = (*p).elems[i];
             }
@@ -1539,6 +1539,7 @@ void Constraint::rotational_invariance()
     } // order
 
     for (order = 0; order < maxorder; ++order) {
+        nparam_sub = nparams[order] + nparams[order - 1];
         remove_redundant_rows(nparam_sub, const_rotation_cross[order], eps6);
         remove_redundant_rows(nparams[order], const_rotation_self[order], eps6);
     }

diff --git a/alm/interaction.cpp b/alm/interaction.cpp
@@ -412,11 +412,8 @@ void Interaction::search_interactions(std::vector<int> **interaction_list_out,
     int ikd, jkd;
 
     double cutoff_tmp;
-    std::set<IntList> *interacting_atom_pairs;
     std::vector<int> intlist;
 
-    memory->allocate(interacting_atom_pairs, maxorder);
-
     for (order = 0; order < maxorder; ++order) {
         for (i = 0; i < natmin; ++i) {
             interaction_list_out[order][i].clear();
@@ -451,8 +448,6 @@ void Interaction::search_interactions(std::vector<int> **interaction_list_out,
 
     for (order = 0; order < maxorder; ++order) {
 
-        interacting_atom_pairs[order].clear();
-
         std::cout << std::endl << "   ***" << str_order[order] << "***" << std::endl;
 
         for (i = 0; i < natmin; ++i) {
@@ -465,11 +460,11 @@ void Interaction::search_interactions(std::vector<int> **interaction_list_out,
             iat = symmetry->map_p2s[i][0];
 
             intlist.clear();
-            for (std::vector<int>::const_iterator it = interaction_list_out[order][i].begin();
+            for (auto it = interaction_list_out[order][i].begin();
                  it != interaction_list_out[order][i].end(); ++it) {
                 intlist.push_back((*it));
             }
-            std::sort(intlist.begin(), intlist.end()); // Necessarily to sort here
+            std::sort(intlist.begin(), intlist.end()); 
 
             // write atoms inside the cutoff radius
             int id = 0;
@@ -493,44 +488,10 @@ void Interaction::search_interactions(std::vector<int> **interaction_list_out,
             std::cout << std::endl << std::endl;
             std::cout << "    Number of total interaction pairs = "
                 << interaction_list_out[order][i].size() << std::endl << std::endl;
-
-            int *intarr;
-            memory->allocate(intarr, order + 2);
-
-            if (intlist.size() > 0) {
-                if (order == 0) {
-                    for (unsigned int ielem = 0; ielem < intlist.size(); ++ielem) {
-                        intarr[0] = iat;
-                        intarr[1] = intlist[ielem];
-                        insort(2, intarr);
-
-                        interacting_atom_pairs[0].insert(IntList(2, intarr));
-                    }
-                } else if (order > 0) {
-                    CombinationWithRepetition<int> g(intlist.begin(), intlist.end(), order + 1);
-                    do {
-                        std::vector<int> data = g.now();
-                        intarr[0] = iat;
-                        intarr[1] = data[0];
-                        for (unsigned int isize = 1; isize < data.size(); ++isize) {
-                            intarr[isize + 1] = data[isize];
-                        }
-                        if (!is_incutoff(order + 2, intarr, order)) continue;
-
-                        insort(order + 2, intarr);
-                        interacting_atom_pairs[order].insert(IntList(order + 2, intarr));
-
-                    } while (g.next());
-                }
-            }
-            intlist.clear();
-            memory->deallocate(intarr);
         }
     }
 
     std::cout << std::endl;
-
-    memory->deallocate(interacting_atom_pairs);
 }
 
 bool Interaction::is_incutoff(const int n,
@@ -752,8 +713,8 @@ void Interaction::calc_mindist_clusters(std::vector<int> **interaction_pair_in,
 
             // List of 2-body interaction pairs
             intlist.clear();
-            for (std::vector<int>::const_iterator it = interaction_pair_in[order][i].begin();
-                 it != interaction_pair_in[order][i].end(); ++it) {
+            for (auto it = interaction_pair_in[order][i].cbegin();
+                 it != interaction_pair_in[order][i].cend(); ++it) {
                 intlist.push_back((*it));
             }
             std::sort(intlist.begin(), intlist.end()); // Need to sort here
@@ -840,8 +801,8 @@ void Interaction::calc_mindist_clusters(std::vector<int> **interaction_pair_in,
 
                         // Loop over the cell images of atom 'jat' and add to the list 
                         // as a candidate for the minimum distance cluster
-                        for (std::vector<DistInfo>::const_iterator it = distance_image[iat][jat].begin();
-                             it != distance_image[iat][jat].end(); ++it) {
+                        for (auto it = distance_image[iat][jat].cbegin();
+                             it != distance_image[iat][jat].cend(); ++it) {
                             if (exist[(*it).cell]) {
                                 if (rc_tmp < 0.0 || (*it).dist <= rc_tmp) {
                                     cell_vector.push_back((*it).cell);
@@ -930,42 +891,9 @@ void Interaction::calc_mindist_clusters(std::vector<int> **interaction_pair_in,
                         mindist_cluster_out[order][i].insert(MinimumDistanceCluster(data_now,
                                                                                     comb_cell_atom_center,
                                                                                     distmax));
-
                     }
                 }
                 memory->deallocate(list_now);
-
-                /*
-                  std::sort(distance_list.begin(), distance_list.end(), MinDistList::compare_sum_distance);
-                  comb_cell_min.clear();
-                  
-                  sum_dist_min = 0.0;
-                  for (j = 0; j < distance_list[0].dist.size(); ++j) {
-                  sum_dist_min += distance_list[0].dist[j];
-                  }
-                  
-                  for (j = 0; j < distance_list.size(); ++j) {
-                  sum_dist = 0.0;
-                  
-                  for (k = 0; k < distance_list[j].dist.size(); ++k) {
-                  sum_dist += distance_list[j].dist[k];
-                  }
-                  
-                  // In the following, only pairs having minimum sum of distances
-                  // are stored. However, we found that this treatment didn't
-                  // return a reliable value of phonon linewidth.
-                  if (std::abs(sum_dist - sum_dist_min) < eps6) {
-                  comb_cell_min.push_back(distance_list[j].cell);
-                  } else {
-                  break;
-                  }
-                  // Therefore, we consider all duplicate pairs
-                  //comb_cell_min.push_back(distance_list[j].cell);
-                  }
-                  // mindist_cluster_out[order][i].insert(MinimumDistanceCluster(data, comb_cell_min));
-                  */
-
-
             }
         }
     }

diff --git a/docs/img/Si_phband_DFT.png b/docs/img/Si_phband_DFT.png
diff --git a/docs/img/Si_phdos_DFT.png b/docs/img/Si_phdos_DFT.png
diff --git a/docs/img/si_cumulative.png b/docs/img/si_cumulative.png
diff --git a/docs/img/si_kappa.png b/docs/img/si_kappa.png
diff --git a/docs/img/si_kappa_spec.png b/docs/img/si_kappa_spec.png
diff --git a/docs/img/si_tau.png b/docs/img/si_tau.png
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -62,7 +62,7 @@
 # The short X.Y version.
 version = '1.0'
 # The full version, including alpha/beta/rc tags.
-release = '1.0.0'
+release = '1.0.2'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

diff --git a/docs/source/input/inputanphon.rst b/docs/source/input/inputanphon.rst
@@ -114,7 +114,6 @@ List of input variables
 
  === =======================================================
   0   Symmetry operations won’t be saved in “SYMM_INFO_PRIM”
-
   1   Symmetry operations will be saved in “SYMM_INFO_PRIM”
  === =======================================================
 
@@ -157,6 +156,18 @@ List of input variables
 
 ````
 
+* BORNSYM-tag = 0 | 1
+
+ === =================================================================
+  0   Do not symmetrize Born effective charges
+  1   Symmetrize Born effective charges by using point group symmetry
+ === =================================================================
+
+ :Default: 0
+ :Type: Integer
+
+````
+
 * TMIN, TMAX, DT-tags : Temperature range and its stride in units of Kelvin
 
  :Default: ``TMIN = 0``, ``TMAX = 1000``, ``DT = 10``
@@ -193,6 +204,25 @@ List of input variables
 
 ````
 
+* BCONNECT-tag = 0 | 1 | 2 
+
+ === ===================================================================================
+  0   | Phonon band is saved without change (sorted in order of energy)
+
+  1   | Phonon band is connected by using the similarity of eigenvectors.
+
+  2   | Same as ``BCONNECT=1``. In addition, information of the connectivity is 
+      | saved as ``PREFIX.connection``.
+ === ===================================================================================
+
+ :Default: 0
+ :Type: Integer
+ :Description: The algorithm for connecting a band structure is described here_.
+
+ .. _here : https://www.slideshare.net/TakeshiNishimatsu/two-efficient-algorithms-for-drawing-accurate-and-beautiful-phonon-dispersion
+
+````
+
 * CLASSICAL-tag = 0 | 1
 
  === =======================================================

diff --git a/docs/source/tutorial_pages/01_silicon.rst b/docs/source/tutorial_pages/01_silicon.rst
@@ -69,20 +69,25 @@ Using the script :red:`displace.py` in the tools/ directory, you can generate th
     **QE**
     ::
 
-        $ python displace.py --QE=si222.pw.in --mag=0.02 si222.pattern_HARMONIC
+        $ python displace.py --QE=si222.pw.in --mag=0.01 si222.pattern_HARMONIC
 
     **VASP**
     ::
 
-        $ python displace.py --VASP=POSCAR.orig --mag=0.02 si222.pattern_HARMONIC
+        $ python displace.py --VASP=POSCAR.orig --mag=0.01 si222.pattern_HARMONIC
 
     **xTAPP**
     ::
 
-        $ python displace.py --xTAPP=si222.cg --mag=0.02 si222.pattern_HARMONIC
+        $ python displace.py --xTAPP=si222.cg --mag=0.01 si222.pattern_HARMONIC
+
+    **OpenMX**
+    ::
+
+        $ python displace.py --OpenMX=si222.dat --mag=0.01 si222.pattern_HARMONIC
 
 The ``--mag`` option specifies the displacement length in units of Angstrom. 
-You need to specify an input file with equilibrium atomic positions either by the ``--QE``, ``--VASP``, ``--xTAPP``, or ``--LAMMPS``.
+You need to specify an input file with equilibrium atomic positions either by the ``--QE``, ``--VASP``, ``--xTAPP``, ``--OpenMX`` or ``--LAMMPS``.
 
 Then, calculate atomic forces for all the configurations. This can be done with a simple shell script 
 as follows::
@@ -124,6 +129,12 @@ The next step is to collect the displacement data and force data by the Python s
     $ python extract.py --xTAPP=si222.cg --get=disp *.str > disp.dat
     $ python extract.py --xTAPP=si222.cg --get=force *.str > force.dat
 
+    **OpenMX**
+    ::
+
+    $ python extract.py --OpenMX=si222.dat --get=disp *.out > disp.dat
+    $ python extract.py --OpenMX=si222.dat --get=force *.out > force.dat
+
 In the above examples, atomic displacements of all the configurations are merged as *disp.dat*, and the corresponding atomic forces are saved in the file *force.dat*. These files will be used in the following fitting procedure as ``DFILE`` and ``FFILE``. (See :ref:`Format of DFILE and FFILE<label_format_DFILE>`).
 
 
@@ -178,7 +189,7 @@ Try
 ::
 
   $ grep "Fitting error" si_alm2.log
-  Fitting error (%) : 1.47766
+  Fitting error (%) : 0.567187
 
 The other file :red:`si222.xml` contains crystal structure, symmetry, IFCs, and all other information necessary for subsequent phonon calculations.
 
@@ -418,13 +429,13 @@ Phonon lifetime
 The file :red:`si222.result` contains phonon linewidths at irreducible :math:`k` points. 
 You can extract phonon lifetime from this file as follows::
 
-    $ analyze_phonons.py --calc tau --temp 300 si222.result > tau300K.dat
+    $ analyze_phonons.py --calc tau --temp 300 si222.result > tau300K_10.dat
     $ gnuplot
     gnuplot> set xrange [1:]
     gnuplot> set logscale y
     gnuplot> set xlabel "Phonon frequency (cm^{-1})"
     gnuplot> set ylabel "Phonon lifetime (ps)"
-    gnuplot> plot "tau300K.dat" using 3:4 w p
+    gnuplot> plot "tau300K_10.dat" using 3:4 w p
 
 .. figure:: ../img/si_tau.png
    :scale: 40%
@@ -440,12 +451,12 @@ Cumulative thermal conductivity
 
 Following the procedure below, you can obtain the :ref:`cumulative thermal conductivity <cumulative_kappa>`::
 
-    $ analyze_phonons.py --calc cumulative --temp 300 --length 10000:5 si222.result > cumulative_300K.dat
+    $ analyze_phonons.py --calc cumulative --temp 300 --length 10000:5 si222.result > cumulative_300K_10.dat
     $ gnuplot
     gnuplot> set logscale x
     gnuplot> set xlabel "L (nm)"
     gnuplot> set ylabel "Cumulative kappa (W/mK)"
-    gnuplot> plot "cumulative_300K.dat" using 1:2 w lp
+    gnuplot> plot "cumulative_300K_10.dat" using 1:2 w lp
 
 .. figure:: ../img/si_cumulative.png
    :scale: 40%
@@ -495,11 +506,11 @@ The frequency range is specified with the ``EMIN``, ``EMAX``, and ``DELTA_E`` ta
 
 After the calculation finishes, you can find the file :red:`si222.kl_spec` which contains the spectra of thermal conductivity at various temperatures. You can plot the data at room temperature as follows::
     
-    $ awk '{if ($1 == 300.0) print $0}' si222.kl_spec > si222_300K.kl_spec
+    $ awk '{if ($1 == 300.0) print $0}' si222.kl_spec > si222_300K_10.kl_spec
     $ gnuplot
     gnuplot> set xlabel "Frequency (cm^{-1})"
     gnuplot> set ylabel "Spectrum of kappa (W/mK/cm^{-1})"
-    gnuplot> plot "si222_300K.kl_spec" using 2:3 w l lt 2 lw 2
+    gnuplot> plot "si222_300K_10.kl_spec" using 2:3 w l lt 2 lw 2
 
 .. figure:: ../img/si_kappa_spec.png
    :scale: 40%

diff --git a/example/Si/reference/disp.dat b/example/Si/reference/disp.dat
@@ -1,4 +1,4 @@
-      0.0377940       0.0000000       0.0000000
+      0.0188980       0.0000000       0.0000000
       0.0000000       0.0000000       0.0000000
       0.0000000       0.0000000       0.0000000
       0.0000000       0.0000000       0.0000000