Merge branch 'release/v.1.0.0'

Author: ttadano

ChangeLog.md

@@ -1,3 +1,24 @@
+# Ver. 1.0.0 (2017-11-21)
+
+## New
+
+- Self-consistent phonon calculation (``MODE = SCPH``) 
+- Ewald summation for the non-analytic term of the dynamical matrix (``NONANALYTIC = 3``).
+- Support of the ``CLASSICAL`` option.
+- Python auxiliary script for LAMMPS
+
+## Changes
+
+- P+ and P- are printed seperately in PREFIX.sps when ``SPS = 1``
+- Use C++11 standard. From this version, the C++ compiler must support the C++11 standard.
+- The **anphon** code symmetrize the Born effective charges
+
+## Fix
+
+- Loosen the tolerance to detect the multiplicity of force constants. This is an important fix for low symmetry structures.
+- Fixed a problem of the restart mode of ``MODE = RTA``
+
+
 # Ver. 0.9.8 (2016-7-14)
 
 ## New

README.md

@@ -1,7 +1,12 @@
-# ALAMODE 
-### Version 0.9.8 (Beta)
+# ALAMODE
+ 
+[![License][license-image]][license-url]
+[![Doc status][docs-image]][docs-url]
+
+### Version 1.0.0
 ![alt ALAMODE](./docs/img/alamode.png)
 
+
 - - -
 
 ## Introduction 
@@ -20,7 +25,7 @@ linewidth, and other anharmonic phonon properties from first principles.
 * Extraction of harmonic and anharmonic force constants based on the supercell approach
 * Applicable to any crystal structures and low-dimensional systems
 * Accurate treatment of translational and rotational invariance
-* Interface to VASP, Quantum-ESPRESSO, and xTAPP codes
+* Interface to VASP, Quantum-ESPRESSO, xTAPP, and LAMMPS codes
 * Parallelization with MPI+OpenMP
 
 ### Harmonic properties
@@ -43,12 +48,15 @@ linewidth, and other anharmonic phonon properties from first principles.
 * Phonon linewidth due to 3-phonon interactions
 * Phonon frequency shift due to 3- and 4-phonon interactions
 * Temperature-dependent effective potential method
+* Self-consistent phonon calculation
 
 ## Prerequisite
 * C++ compiler
-* LAPACK libarary
+* LAPACK library
 * MPI library
 * Boost C++ library
+* FFT library
+* Eigen3 library
 
 ## Download
 
@@ -78,9 +86,23 @@ http://alamode.readthedocs.io
 
 
 ## License
-Copyright (c) 2014, 2015, 2016 Terumasa Tadano
+Copyright (c) 2014--2017 Terumasa Tadano
 This software is released under the MIT license. 
 For license rights and limitations, see LICENSE.txt file.
 
 ## Author
-Terumasa Tadano (The University of Tokyo, Japan)
+Terumasa Tadano (National Institute for Materials Science, Japan)
+
+## Contributors
+
+* Tatsuro Nishimoto (Univ. Tokyo)
+* Yusuke Oba (Univ. Tokyo)
+
+
+
+[license-image]: https://img.shields.io/github/license/ttadano/alamode.svg
+[license-url]:  https://github.com/ttadano/alamode/blob/develop/LICENSE.txt
+
+[docs-image]:  https://readthedocs.org/projects/alamode/badge/?version=latest
+[docs-url]: https://alamode.readthedocs.io/en/latest/?badge=latest
+

alm/Makefile.linux

@@ -6,7 +6,7 @@
 #-----------------------------------------------
 
 CXX = icpc 
-CXXFLAGS = -O2 -xHOST -openmp 
+CXXFLAGS = -O2 -xHOST -openmp -std=c++11
 INCLUDE = -I../include
 
 CXXL = ${CXX}

alm/Makefile.osx

@@ -7,8 +7,8 @@
 
 # Use gcc >= 4.8 to use OpenMP 
 # OpenMP-enabled gcc can be installed via homebrew
-CXX = g++-6
-CXXFLAGS = -O2 -fopenmp 
+CXX = g++-7
+CXXFLAGS = -O2 -fopenmp -std=c++11
 INCLUDE = -I../include
 
 CXXL = ${CXX}

alm/alamode.cpp

@@ -70,7 +70,7 @@ ALM::ALM(int narg, char **arg)
 
     finalize();
 
-    std::cout << std::endl << " Job finished at " 
+    std::cout << std::endl << " Job finished at "
         << timer->DateAndTime() << std::endl;
 }
 
@@ -116,4 +116,3 @@ void ALM::finalize()
     delete writes;
     delete memory;
 }
-

alm/constraint.cpp

@@ -28,7 +28,9 @@
 
 using namespace ALM_NS;
 
-Constraint::Constraint(ALM *alm) : Pointers(alm) {}
+Constraint::Constraint(ALM *alm) : Pointers(alm)
+{
+}
 
 Constraint::~Constraint()
 {
@@ -556,25 +558,26 @@ void Constraint::constraint_from_symmetry(std::vector<ConstraintClass> *const_ou
     int **xyzcomponent;
     int nparams;
 
+    double *arr_constraint;
     bool has_constraint_from_symm = false;
-
     std::set<FcProperty> list_found;
+    std::vector<std::vector<double>> const_mat;
 
     for (isym = 0; isym < symmetry->nsym; ++isym) {
         if (symmetry->sym_available[isym]) continue;
         has_constraint_from_symm = true;
     }
 
-    for (order = 0; order < maxorder; ++order) {
-        const_out[order].clear();
-    }
+    for (order = 0; order < maxorder; ++order) const_out[order].clear();
 
     if (has_constraint_from_symm) {
         std::cout << "  Generating constraints from crystal symmetry ..." << std::endl;
     }
 
     memory->allocate(index_tmp, maxorder + 1);
 
+    const_mat.clear();
+
     for (order = 0; order < maxorder; ++order) {
 
         nparams = fcs->ndup[order].size();
@@ -589,8 +592,7 @@ void Constraint::constraint_from_symmetry(std::vector<ConstraintClass> *const_ou
 
         // Generate temporary list of parameters
         list_found.clear();
-        for (std::vector<FcProperty>::iterator p = fcs->fc_set[order].begin();
-             p != fcs->fc_set[order].end(); ++p) {
+        for (auto p = fcs->fc_set[order].begin(); p != fcs->fc_set[order].end(); ++p) {
             for (i = 0; i < order + 2; ++i) index_tmp[i] = (*p).elems[i];
             list_found.insert(FcProperty(order + 2, (*p).coef,
                                          index_tmp, (*p).mother));
@@ -602,25 +604,33 @@ void Constraint::constraint_from_symmetry(std::vector<ConstraintClass> *const_ou
 
         int nfcs = fcs->fc_set[order].size();
 
-#pragma omp parallel 
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
         {
+            int j;
             int i_prim;
-
+            int loc_nonzero;
             int *ind;
             int *atm_index, *atm_index_symm;
             int *xyz_index;
             double c_tmp;
-            double *arr_constraint;
 
             std::set<FcProperty>::iterator iter_found;
+            std::vector<double> const_now_omp;
+            std::vector<std::vector<double>> const_omp;
 
-            memory->allocate(arr_constraint, nparams);
             memory->allocate(ind, order + 2);
             memory->allocate(atm_index, order + 2);
             memory->allocate(atm_index_symm, order + 2);
             memory->allocate(xyz_index, order + 2);
 
-#pragma omp for private(i, isym, ixyz) 
+            const_omp.clear();
+            const_now_omp.resize(nparams);
+
+#ifdef _OPENMP
+#pragma omp for private(i, isym, ixyz)
+#endif
             for (int ii = 0; ii < nfcs; ++ii) {
                 FcProperty list_tmp = fcs->fc_set[order][ii];
 
@@ -637,9 +647,9 @@ void Constraint::constraint_from_symmetry(std::vector<ConstraintClass> *const_ou
                         atm_index_symm[i] = symmetry->map_sym[atm_index[i]][isym];
                     if (!fcs->is_inprim(order + 2, atm_index_symm)) continue;
 
-                    for (i = 0; i < nparams; ++i) arr_constraint[i] = 0.0;
+                    for (i = 0; i < nparams; ++i) const_now_omp[i] = 0.0;
 
-                    arr_constraint[list_tmp.mother] = -list_tmp.coef;
+                    const_now_omp[list_tmp.mother] = -list_tmp.coef;
 
                     for (ixyz = 0; ixyz < nxyz; ++ixyz) {
                         for (i = 0; i < order + 2; ++i)
@@ -652,32 +662,64 @@ void Constraint::constraint_from_symmetry(std::vector<ConstraintClass> *const_ou
                         iter_found = list_found.find(FcProperty(order + 2, 1.0, ind, 1));
                         if (iter_found != list_found.end()) {
                             c_tmp = fcs->coef_sym(order + 2, isym, xyz_index, xyzcomponent[ixyz]);
-                            arr_constraint[(*iter_found).mother] += (*iter_found).coef * c_tmp;
+                            const_now_omp[(*iter_found).mother] += (*iter_found).coef * c_tmp;
                         }
                     }
 
-                    if (!is_allzero(nparams, arr_constraint)) {
+                    if (!is_allzero(const_now_omp, loc_nonzero)) {
+                        if (const_now_omp[loc_nonzero] < 0.0) {
+                            for (j = 0; j < nparams; ++j) const_now_omp[j] *= -1.0;
+                        }
+                        const_omp.push_back(const_now_omp);
+                    }
+
+                } // close isym loop
+
+
+                // sort-->uniq the array
+                std::sort(const_omp.begin(), const_omp.end());
+                const_omp.erase(std::unique(const_omp.begin(), const_omp.end()),
+                                const_omp.end());
+
+                // Merge vectors
 #pragma omp critical
-                        const_out[order].push_back(ConstraintClass(nparams,
-                                                                   arr_constraint));
+                {
+                    for (std::vector<std::vector<double>>::iterator it = const_omp.begin();
+                         it != const_omp.end(); ++it) {
+                        const_mat.push_back(*it);
                     }
                 }
-            }
+                const_omp.clear();
+
+            } // close ii loop
 
-            memory->deallocate(arr_constraint);
             memory->deallocate(ind);
             memory->deallocate(atm_index);
             memory->deallocate(atm_index_symm);
             memory->deallocate(xyz_index);
+
+        } // close openmp region
+
+        memory->allocate(arr_constraint, nparams);
+        for (std::vector<std::vector<double>>::reverse_iterator it = const_mat.rbegin();
+             it != const_mat.rend(); ++it) {
+            for (i = 0; i < (*it).size(); ++i) {
+                arr_constraint[i] = (*it)[i];
+            }
+            const_out[order].push_back(ConstraintClass(nparams,
+                                                       arr_constraint));
         }
+        const_mat.clear();
 
         memory->deallocate(xyzcomponent);
+        memory->deallocate(arr_constraint);
+
         remove_redundant_rows(nparams, const_out[order], eps8);
 
         if (has_constraint_from_symm) {
             std::cout << " done." << std::endl;
         }
-    }
+    } // close loop order
 
     memory->deallocate(index_tmp);
 
@@ -712,11 +754,11 @@ void Constraint::translational_invariance()
     std::vector<int> intlist, data;
     std::set<FcProperty> list_found;
     std::set<FcProperty>::iterator iter_found;
-    std::vector<std::vector<int> > data_vec;
+    std::vector<std::vector<int>> data_vec;
     std::vector<FcProperty> list_vec;
     std::vector<FcProperty>::iterator iter_vec;
     std::vector<int> const_now;
-    std::vector<std::vector<int> > const_mat;
+    std::vector<std::vector<int>> const_mat;
 
     std::cout << "  Generating constraints for translational invariance ..." << std::endl;
 
@@ -848,7 +890,7 @@ void Constraint::translational_invariance()
                     memory->allocate(intarr_omp, order + 2);
                     memory->allocate(intarr_copy_omp, order + 2);
 
-                    std::vector<std::vector<int> > const_omp;
+                    std::vector<std::vector<int>> const_omp;
                     std::vector<int> data_omp;
                     std::vector<int> const_now_omp;
 
@@ -913,7 +955,7 @@ void Constraint::translational_invariance()
                         // Merge vectors
 #pragma omp critical
                         {
-                            for (std::vector<std::vector<int> >::iterator it = const_omp.begin();
+                            for (std::vector<std::vector<int>>::iterator it = const_omp.begin();
                                  it != const_omp.end(); ++it) {
                                 const_mat.push_back(*it);
                             }
@@ -943,7 +985,7 @@ void Constraint::translational_invariance()
         // Copy to constraint class 
 
         const_translation[order].clear();
-        for (std::vector<std::vector<int> >::reverse_iterator it = const_mat.rbegin();
+        for (std::vector<std::vector<int>>::reverse_iterator it = const_mat.rbegin();
              it != const_mat.rend(); ++it) {
             for (i = 0; i < (*it).size(); ++i) {
                 arr_constraint[i] = static_cast<double>((*it)[i]);
@@ -1007,7 +1049,7 @@ void Constraint::rotational_invariance()
     CombinationWithRepetition<int> g;
 
     std::vector<int> atom_tmp;
-    std::vector<std::vector<int> > cell_dummy;
+    std::vector<std::vector<int>> cell_dummy;
     std::set<MinimumDistanceCluster>::iterator iter_cluster;
 
     setup_rotation_axis(valid_rotation_axis);
@@ -1614,7 +1656,6 @@ bool Constraint::is_allzero(const std::vector<int> vec, int &loc)
 {
     loc = -1;
     for (int i = 0; i < vec.size(); ++i) {
-        //    for(std::vector<int>::const_iterator it = vec.begin(); it != vec.end(); ++it){
         if (std::abs(vec[i]) > 0) {
             loc = i;
             return false;
@@ -1623,6 +1664,18 @@ bool Constraint::is_allzero(const std::vector<int> vec, int &loc)
     return true;
 }
 
+bool Constraint::is_allzero(const std::vector<double> vec, int &loc)
+{
+    loc = -1;
+    for (int i = 0; i < vec.size(); ++i) {
+        if (std::abs(vec[i]) > eps) {
+            loc = i;
+            return false;
+        }
+    }
+    return true;
+}
+
 void Constraint::setup_rotation_axis(bool flag[3][3])
 {
     unsigned int mu, nu;

alm/constraint.h

@@ -131,6 +131,7 @@ namespace ALM_NS
         void setup_rotation_axis(bool [3][3]);
         bool is_allzero(const int, const double *, const int nshift = 0);
         bool is_allzero(const std::vector<int>, int &);
+        bool is_allzero(const std::vector<double>, int &);
 
         void remove_redundant_rows(const int, std::vector<ConstraintClass> &,
                                    const double tolerance = eps12);