Merge branch 'feature/scph' into develop

README.md

@@ -1,5 +1,5 @@
 # ALAMODE 
-### Version 0.9.8 (Beta)
+### Version 1.0.0
 ![alt ALAMODE](./docs/img/alamode.png)
 
 - - -
@@ -43,12 +43,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 library
 * MPI library
 * Boost C++ library
+* FFT library
+* Eigen3 library
 
 ## Download
 

anphon/Makefile.linux

@@ -28,7 +28,8 @@ CXXSRC= phonons.cpp error.cpp fcs_phonon.cpp parsephon.cpp dynamical.cpp  \
 	main.cpp memory.cpp system.cpp timer.cpp write_phonons.cpp kpoint.cpp \
 	phonon_dos.cpp phonon_velocity.cpp integration.cpp relaxation.cpp \
 	thermodynamics.cpp conductivity.cpp symmetry_core.cpp \
-        mpi_common.cpp gruneisen.cpp isotope.cpp selfenergy.cpp
+        mpi_common.cpp gruneisen.cpp isotope.cpp selfenergy.cpp \
+	scph.cpp 
 
 OBJS= ${CXXSRC:.cpp=.o}
 

anphon/Makefile.osx

@@ -6,8 +6,8 @@
 #-----------------------------------------------
 
 MPICXX = mpic++
-CXXFLAGS = -O2 
-INCLUDE = -I../include 
+CXXFLAGS = -O3 -std=c++11
+INCLUDE = -I../include  -I$(HOME)/include
 
 LINKER = ${MPICXX}
 
@@ -17,7 +17,7 @@ LINKER = ${MPICXX}
 LDFLAGS = -framework Accelerate
 
 LAPACK = -llapack
-LIBS = ${LAPACK} 
+LIBS = ${LAPACK}  -lfftw3
 
 #----------------------------------------------
 # General rules
@@ -32,7 +32,8 @@ CXXSRC= phonons.cpp error.cpp fcs_phonon.cpp parsephon.cpp dynamical.cpp  \
 	main.cpp memory.cpp system.cpp timer.cpp write_phonons.cpp kpoint.cpp \
 	phonon_dos.cpp phonon_velocity.cpp integration.cpp relaxation.cpp \
 	thermodynamics.cpp conductivity.cpp symmetry_core.cpp \
-        mpi_common.cpp gruneisen.cpp isotope.cpp selfenergy.cpp
+        mpi_common.cpp gruneisen.cpp isotope.cpp selfenergy.cpp \
+	scph.cpp 
 
 OBJS= ${CXXSRC:.cpp=.o}
 

anphon/dynamical.cpp

@@ -91,7 +91,7 @@ void Dynamical::setup_dynamical(std::string mode)
         eigenvectors = true;
 
         /*
-        if (phon->mode == "RTA") {
+        if (phon->mode == "RTA" || phon->mode == "SCPH") {
             eigenvectors = true;
         } else {
             if (print_eigenvectors || writes->print_msd || writes->print_xsf || writes->print_anime

anphon/dynamical.h

@@ -67,25 +67,30 @@ namespace PHON_NS
                     double *, std::complex<double> **, bool);
         void modify_eigenvectors();
 
+
+
         double fold(double);
         double freq(const double);
 
         void calc_participation_ratio_all(std::complex<double> ***,
                                           double **,
                                           double ***);
 
-    private:
-
-        void load_born();
         void calc_analytic_k(double *,
-                             std::vector<FcsClassExtent>,
-                             std::complex<double> **);
+            std::vector<FcsClassExtent>,
+            std::complex<double> **);
         void calc_nonanalytic_k(double *, double *,
-                                std::complex<double> **);
+                std::complex<double> **);
         void calc_nonanalytic_k2(double *, double *,
-                                 std::vector<FcsClassExtent>,
-                                 std::complex<double> **);
+                 std::vector<FcsClassExtent>,
+                 std::complex<double> **);
 
+
+    private:
+
+        void load_born();
+
+
         void prepare_mindist_list(std::vector<int> **);
         void calc_atomic_participation_ratio(std::complex<double> *, double *);
         double distance(double *, double *);

anphon/fcs_phonon.cpp

@@ -80,6 +80,11 @@ void Fcs_phonon::setup(std::string mode)
             maxorder = 2;
             require_quartic = false;
         }
+    } else if (mode == "SCPH") {
+        require_cubic = true;
+        require_quartic = true;
+        maxorder = 3;
+        relaxation->quartic_mode = 1;
     }
 
     memory->allocate(force_constant_with_cell, maxorder);

anphon/fcs_phonon.h

@@ -134,6 +134,12 @@ namespace PHON_NS
             index_a.push_back(a.pairs[i].index);
             index_b.push_back(b.pairs[i].index);
         }
+        for (int i = 0; i < a.pairs.size(); ++i) {
+            index_a.push_back(a.pairs[i].tran);
+            index_a.push_back(a.pairs[i].cell_s);
+            index_b.push_back(b.pairs[i].tran);
+            index_b.push_back(b.pairs[i].cell_s);
+        }
         return lexicographical_compare(index_a.begin(), index_a.end(), index_b.begin(), index_b.end());
     }
 

anphon/kpoint.cpp

@@ -62,7 +62,7 @@ Kpoint::~Kpoint()
 
 void Kpoint::kpoint_setups(std::string mode)
 {
-    small_group_of_k = NULL;
+    small_group_of_k = nullptr;
     symmetry->symmetry_flag = true;
 
     unsigned int i, j;

anphon/parsephon.cpp

@@ -34,6 +34,9 @@
 #include "isotope.h"
 #include "phonon_velocity.h"
 #include "integration.h"
+#include "scph.h"
+#include <boost/lexical_cast.hpp>
+
 #include <boost/algorithm/string.hpp>
 #include <boost/lexical_cast.hpp>
 
@@ -93,27 +96,36 @@ void Input::parse_general_vars()
 
     int i;
     int nsym, nbands, ismear, nkd;
+    unsigned int maxiter;
     unsigned int nonanalytic;
+    unsigned int ialgo_scph;
     double *masskd;
     double Tmin, Tmax, dT, na_sigma, epsilon;
     double emin, emax, delta_e;
-    double tolerance;
+    double tolerance, tolerance_scph;
+    double mixalpha;
     bool printsymmetry;
-    bool restart;
+    bool restart, restart_scph;
     bool sym_time_reversal, use_triplet_symmetry;
+    bool selenergy_offdiagonal;
     bool update_fc2;
+    bool lower_temp, warm_start;
 
     struct stat st;
     std::string prefix, mode, fcsinfo, fc2info;
     std::string borninfo, file_result;
+    std::string file_dymat;
     std::string *kdname;
     std::string str_tmp;
     std::string str_allowed_list = "PREFIX MODE NSYM TOLERANCE PRINTSYM FCSXML FC2XML TMIN TMAX DT \
                                    NBANDS NONANALYTIC BORNINFO NA_SIGMA ISMEAR EPSILON EMIN EMAX DELTA_E \
-                                   RESTART TREVSYM NKD KD MASS TRISYM";
+                                   RESTART TREVSYM NKD KD MASS TRISYM KMESH_SCPH KMESH_INTERPOLATE \
+                                   MIXALPHA MAXITER RESTART_SCPH IALGO SELF_OFFDIAG PREC_EWALD TOL_SCPH \
+                                   LOWER_TEMP WARMSTART";
     std::string str_no_defaults = "PREFIX MODE FCSXML NKD KD MASS";
     std::vector<std::string> no_defaults, celldim_v;
     std::vector<std::string> kdname_v, masskd_v;
+    std::vector<int> kmesh_v, kmesh_interpolate_v;
     std::map<std::string, std::string> general_var_dict;
 
     if (from_stdin) {
@@ -141,6 +153,7 @@ void Input::parse_general_vars()
     mode = general_var_dict["MODE"];
 
     file_result = prefix + ".result";
+    file_dymat = prefix + ".scph_dymat";
 
     std::transform(mode.begin(), mode.end(), mode.begin(), toupper);
     assign_val(nsym, "NSYM", general_var_dict);
@@ -185,6 +198,7 @@ void Input::parse_general_vars()
     nonanalytic = 0;
     nsym = 0;
     tolerance = 1.0e-6;
+    tolerance_scph = 1.0e-10;
     printsymmetry = false;
     sym_time_reversal = false;
     use_triplet_symmetry = true;
@@ -198,6 +212,15 @@ void Input::parse_general_vars()
         restart = false;
     }
 
+    // if file_dymat exists in the current directory, 
+    // restart mode will be automatically turned on for SCPH calculations.
+
+    if (stat(file_dymat.c_str(), &st) == 0) {
+        restart_scph = true;
+    } else {
+        restart_scph = false;
+    }
+
     nbands = -1;
     borninfo = "";
     fc2info = "";
@@ -206,6 +229,12 @@ void Input::parse_general_vars()
     epsilon = 10.0;
     na_sigma = 0.1;
 
+    maxiter = 1000;
+    mixalpha = 0.1;
+    selenergy_offdiagonal = true;
+    ialgo_scph = 0;
+    lower_temp = true;
+    warm_start = true;
 
     // Assign given values
 
@@ -224,6 +253,7 @@ void Input::parse_general_vars()
 
     assign_val(nonanalytic, "NONANALYTIC", general_var_dict);
     assign_val(restart, "RESTART", general_var_dict);
+    assign_val(restart_scph, "RESTART_SCPH", general_var_dict);
 
     assign_val(nbands, "NBANDS", general_var_dict);
     assign_val(borninfo, "BORNINFO", general_var_dict);
@@ -236,6 +266,62 @@ void Input::parse_general_vars()
 
     assign_val(use_triplet_symmetry, "TRISYM", general_var_dict);
 
+    assign_val(maxiter, "MAXITER", general_var_dict);
+    assign_val(mixalpha, "MIXALPHA", general_var_dict);
+    assign_val(selenergy_offdiagonal, "SELF_OFFDIAG", general_var_dict);
+    assign_val(ialgo_scph, "IALGO", general_var_dict);
+    assign_val(tolerance_scph, "TOL_SCPH", general_var_dict);
+    assign_val(lower_temp, "LOWER_TEMP", general_var_dict);
+    assign_val(warm_start, "WARMSTART", general_var_dict);
+
+
+    str_tmp = general_var_dict["KMESH_SCPH"];
+
+        if (!str_tmp.empty()) {
+
+            std::istringstream is(str_tmp);
+
+            while (1) {
+                str_tmp.clear();
+                is >> str_tmp;
+                if (str_tmp.empty()) {
+                    break;
+                }
+                kmesh_v.push_back(my_cast<unsigned int>(str_tmp));
+            }
+
+            if (kmesh_v.size() != 3) {
+                error->exit("parse_general_vars", "The number of entries for KMESH_SCPH has to be 3.");
+            }
+        } else {
+            if (mode == "SCPH") {
+                error->exit("parse_general_vars", "Please specify KMESH_SCPH for mode = SCPH");
+            }
+        }
+
+        str_tmp = general_var_dict["KMESH_INTERPOLATE"];
+        if (!str_tmp.empty()) {
+
+            std::istringstream is(str_tmp);
+
+            while (1) {
+                str_tmp.clear();
+                is >> str_tmp;
+                if (str_tmp.empty()) {
+                    break;
+                }
+                kmesh_interpolate_v.push_back(my_cast<unsigned int>(str_tmp));
+            }
+
+            if (kmesh_interpolate_v.size() != 3) {
+                error->exit("parse_general_vars", "The number of entries for KMESH_INTERPOLATE has to be 3.");
+            }
+        } else {
+            if (mode == "SCPH") {
+                error->exit("parse_general_vars", "Please specify KMESH_INTERPOLATE for mode = SCPH");
+            }
+        }
+
     if (nonanalytic > 2) {
         error->exit("parse_general_vars",
                     "NONANALYTIC should be 0, 1, or 2.");
@@ -287,6 +373,23 @@ void Input::parse_general_vars()
     integration->ismear = ismear;
     relaxation->use_triplet_symmetry = use_triplet_symmetry;
 
+    if (mode == "SCPH") {
+        for (i = 0; i < 3; ++i) {
+            scph->kmesh_scph[i] = kmesh_v[i];
+            scph->kmesh_interpolate[i] = kmesh_interpolate_v[i];
+        }
+        scph->mixalpha = mixalpha;
+        scph->maxiter = maxiter;
+        scph->restart_scph = restart_scph;
+        scph->selfenergy_offdiagonal = selenergy_offdiagonal;
+        scph->ialgo = ialgo_scph;
+        scph->tolerance_scph = tolerance_scph;
+        scph->lower_temp = lower_temp;
+        scph->warmstart_scph = warm_start;
+    }
+    kmesh_v.clear();
+    kmesh_interpolate_v.clear();
+
     general_var_dict.clear();
 }
 
@@ -307,6 +410,7 @@ void Input::parse_analysis_vars(const bool use_default_values)
     bool two_phonon_dos;
     bool print_xsf, print_anime;
     bool print_V3, participation_ratio;
+    bool print_self_consistent_fc2;
     bool bubble_omega;
 
     int quartic_mode;
@@ -351,6 +455,8 @@ void Input::parse_analysis_vars(const bool use_default_values)
 
     calculate_kappa_spec = 0;
 
+    print_self_consistent_fc2 = false;
+
 
     // Assign values to variables
 
@@ -489,6 +595,10 @@ void Input::parse_analysis_vars(const bool use_default_values)
         memory->deallocate(isotope_factor);
     }
 
+    if (phon->mode == "SCPH") {
+        scph->print_self_consistent_fc2 = print_self_consistent_fc2;
+    }
+
     analysis_var_dict.clear();
 }
 

anphon/phonon_dos.cpp

@@ -525,6 +525,35 @@ void Dos::calc_total_scattering_phase_space(double **omega,
     }
 }
 
+void Dos::calc_dos_scph(double ***eval_anharm, double **dos_scph)
+{
+    int i;
+    unsigned int j, k;
+    unsigned int iT;
+    unsigned int nk = kpoint->nk;
+    unsigned int neval = dynamical->neval;
+    double **eval;
+
+    double Tmin = system->Tmin;
+    double Tmax = system->Tmax;
+    double dT = system->dT;
+    unsigned int NT = static_cast<unsigned int>((Tmax - Tmin) / dT) + 1;
+
+    memory->allocate(eval, neval, nk);
+
+    for (iT = 0; iT < NT; ++iT) {
+
+        for (j = 0; j < nk; ++j) {
+            for (k = 0; k < neval; ++k) {
+                eval[k][j] = writes->in_kayser(eval_anharm[iT][j][k]);
+                //                std::cout << eval[k][j] << std::endl;
+            }
+        }
+
+        calc_dos(nk_irreducible, kmap_irreducible, eval, n_energy, energy_dos,
+                 dos_scph[iT], neval, integration->ismear, kpoint->kpoint_irred_all);
+    }
+}
 
 void Dos::calc_scattering_phase_space_with_Bose(double **eval,
                                                 const int smearing_method,