Merge branch 'feature/without_nnp' into develop

alm/symmetry.h

@@ -22,47 +22,15 @@
 namespace ALM_NS {
 
     class SymmetryOperation {
-    public:
-        std::vector<int> symop;
-
-        SymmetryOperation();
-
-        // Declaration construction
-
-        SymmetryOperation(const SymmetryOperation &a)
-        {
-            for(std::vector<int>::const_iterator p = a.symop.begin(); p != a.symop.end(); ++p){
-                symop.push_back(*p);
-            }
-        }
-
-        SymmetryOperation(const int rot[3][3], const int trans[3])
-        {
-            for (int i = 0; i < 3; ++i){
-                for (int j = 0; j < 3; ++j){
-                    symop.push_back(rot[i][j]);
-                }
-            }
-            for (int i = 0; i < 3; ++i){
-                symop.push_back(trans[i]);
-            }
-        }
-    };
-
-    inline bool operator<(const SymmetryOperation a, const SymmetryOperation b){
-        return std::lexicographical_compare(a.symop.begin(), a.symop.end(), b.symop.begin(), b.symop.end());
-    }
-
-    class SymmetryOperationTransFloat {
     public:
         int rot[3][3];
         double tran[3];
 
-        SymmetryOperationTransFloat();
+        SymmetryOperation();
 
         // Declaration construction
 
-        SymmetryOperationTransFloat(const int rot_in[3][3], const double tran_in[3])
+        SymmetryOperation(const int rot_in[3][3], const double tran_in[3])
         {
             for (int i = 0; i < 3; ++i){
                 for (int j = 0; j < 3; ++j){
@@ -95,62 +63,46 @@ namespace ALM_NS {
         ~Symmetry();
 
         void init();
-        void setup_symmetry_operation(int, unsigned int&, unsigned int&, double[3][3], double[3][3], 
-            double **, int *);
-
-        unsigned int nsym, nnp;
-        int ntran, natmin;
-        int nsym_s, ntran_s, natmin_s; // for reference system (supercell?)
-        int ntran_ref;
 
+        unsigned int nsym, ntran, natmin;
         int is_printsymmetry;
         int multiply_data;
-
-        int ***symrel_int;
         int *symnum_tran;
+
         double tolerance;
         double ***symrel;
         double **tnons;
 
         int **map_sym;
         int **map_p2s;
-        int **map_p2s_s;
+
         class Maps {
         public:
             int atom_num;
             int tran_num;
         };
-        Maps *map_s2p, *map_s2p_s;
-
-        void genmaps(int, double **, int **, int **, class Symmetry::Maps *);
-
+        Maps *map_s2p;
         bool *sym_available;
 
-        std::string file_sym;
-        std::ofstream ofs_sym;
-        std::ifstream ifs_sym;
-
-        void gensym_notran(std::vector<SymmetryOperation> &);
-        int numsymop(int, double **x, double);
-
     private:
 
+        void setup_symmetry_operation(int, unsigned int&, double[3][3], double[3][3], 
+				      double **, int *);
+        void genmaps(int, double **, int **, int **, class Symmetry::Maps *);
         void findsym(int, double [3][3], double **, std::vector<SymmetryOperation> &);
         bool is_translation(int **);
-
         void symop_in_cart(double [3][3], double[3][3]);
         void pure_translations();
         void print_symmetrized_coordinate(double **);
         void symop_availability_check(double ***, bool *, const int, int &);
-
         void find_lattice_symmetry(double [3][3], std::vector<RotationMatrix> &);
         void find_crystal_symmetry(int, int, std::vector<unsigned int> *, double **x,
-            std::vector<RotationMatrix>, std::vector<SymmetryOperationTransFloat> &);
-        void find_nnp_for_translation(unsigned int &, std::vector<SymmetryOperationTransFloat>);
+            std::vector<RotationMatrix>, std::vector<SymmetryOperation> &);
+	//        void find_nnp_for_translation(unsigned int &, std::vector<SymmetryOperationTransFloat>);
 
+        std::string file_sym;
+        int ***symrel_int;
         std::vector<SymmetryOperation> SymmList;
     };
-
-
 }
 

alm/system.cpp

@@ -290,12 +290,14 @@ void System::load_reference_system()
     int iat, jat;
     int icrd;
 
-    int nat_s, nkd_s;
+    unsigned int nat_s, nkd_s;
+    unsigned int natmin_ref, ntran_ref;
     double lavec_s[3][3];
     int *kd_s;
     double **xcoord_s;
     int *map_ref;
-
+    int **map_p2s_s;
+    Symmetry::Maps *map_s2p_s;
     std::ifstream ifs_fc2;
 
     ifs_fc2.open(constraint->fc2_file.c_str(), std::ios::in);
@@ -326,9 +328,9 @@ void System::load_reference_system()
             std::getline(ifs_fc2, str_tmp);
             std::getline(ifs_fc2, str_tmp);
 
-            ifs_fc2 >> nat_s >> symmetry->natmin_s >> symmetry->ntran_s;
+            ifs_fc2 >> nat_s >> natmin_ref >> ntran_ref;
 
-            if (symmetry->natmin_s != symmetry->natmin) {
+            if (natmin_ref != symmetry->natmin) {
                 error->exit("load_reference_system", "The number of atoms in the primitive cell is not consistent");
             }
 
@@ -339,18 +341,18 @@ void System::load_reference_system()
 
             memory->allocate(xcoord_s, nat_s, 3);
             memory->allocate(kd_s, nat_s);
-            memory->allocate(symmetry->map_p2s_s, symmetry->natmin_s, symmetry->ntran_s);
-            memory->allocate(symmetry->map_s2p_s, nat_s);
+            memory->allocate(map_p2s_s, natmin_ref, ntran_ref);
+            memory->allocate(map_s2p_s, nat_s);
 
             unsigned int ikd, itran, icell;
             std::getline(ifs_fc2, str_tmp);
             std::getline(ifs_fc2, str_tmp);
             for (i = 0; i < nat_s; ++i) {
                 ifs_fc2 >> str_tmp >> ikd >> xcoord_s[i][0] >> xcoord_s[i][1] >> xcoord_s[i][2] >> itran >> icell;
                 kd_s[i] = ikd;
-                symmetry->map_p2s_s[icell - 1][itran - 1] = i;
-                symmetry->map_s2p_s[i].atom_num = icell - 1;
-                symmetry->map_s2p_s[i].tran_num = itran - 1;
+                map_p2s_s[icell - 1][itran - 1] = i;
+                map_s2p_s[i].atom_num = icell - 1;
+                map_s2p_s[i].tran_num = itran - 1;
             }
         }
     }

anphon/fcs_phonon.h

@@ -81,9 +81,9 @@ namespace PHON_NS {
         }
 
         bool operator==(const FcsClassExtent &a) const {
-            return this->atm1 == a.atm1 & this->atm2 == a.atm2 
-                 & this->xyz1 == a.xyz1 & this->xyz2 == a.xyz2
-                 & this->cell_s == a.cell_s;
+            return (this->atm1 == a.atm1) & (this->atm2 == a.atm2) 
+                 & (this->xyz1 == a.xyz1) & (this->xyz2 == a.xyz2)
+                 & (this->cell_s == a.cell_s);
              }
     };
 

anphon/kpoint.cpp

@@ -602,7 +602,7 @@ void Kpoint::reduce_kpoints(const unsigned int nsym, double **xkr, const unsigne
 
         for (i = 0; i < 3; ++i) {
             for (j = 0; j < 3; ++j) {
-                srot[i][j] = static_cast<double>(symmetry->SymmList[isym].symop[3 * i + j]);
+                srot[i][j] = static_cast<double>(symmetry->SymmList[isym].rot[i][j]);
             }
         }
 

anphon/phonon_dos.cpp

@@ -107,7 +107,7 @@ void Dos::setup()
         for (i = 0; i < symmetry->nsym; ++i) {
             for (int j = 0; j < 3; ++j) {
                 for (int k = 0; k < 3; ++k) {
-                    symmetry_tmp[i][j][k] = symmetry->SymmList[i].symop[3 * j + k];
+                    symmetry_tmp[i][j][k] = symmetry->SymmList[i].rot[j][k];
                 }
             }
         }

anphon/relaxation.cpp

@@ -2062,7 +2062,7 @@ void Relaxation::print_momentum_resolved_final_state(const unsigned int NT, doub
         for (isym = 0; isym < symmetry->nsym; ++isym) {
             for (j = 0; j < 3; ++j) {
                 for (k = 0; k < 3; ++k) {
-                    srot[j][k] = static_cast<double>(symmetry->SymmList[isym].symop[3 * j + k]);
+                    srot[j][k] = static_cast<double>(symmetry->SymmList[isym].rot[j][k]);
                 }
             }
 
@@ -2374,7 +2374,7 @@ int Relaxation::knum_sym(const int nk_in, const int symop_num) {
 
     for (i = 0; i < 3; ++i) {
         for (j = 0; j < 3; ++j) {
-            srot[i][j] = static_cast<double>(symmetry->SymmList[symop_num].symop[3 * i + j]);
+            srot[i][j] = static_cast<double>(symmetry->SymmList[symop_num].rot[i][j]);
         }
     }
 
@@ -2402,7 +2402,7 @@ bool Relaxation::is_proper(const int isym)
 
     for (i = 0; i < 3; ++i) {
         for (j = 0; j < 3; ++j) {
-            S[i][j] = static_cast<double>(symmetry->SymmList[isym].symop[3 * i + j]);
+            S[i][j] = static_cast<double>(symmetry->SymmList[isym].rot[i][j]);
         }
     }
 
@@ -2427,9 +2427,9 @@ bool Relaxation::is_symmorphic(const int isym)
     int tran[3];
     bool ret;
 
-    for (i = 0; i < 3; ++i) tran[i] = symmetry->SymmList[isym].symop[9 + i];
+    for (i = 0; i < 3; ++i) tran[i] = symmetry->SymmList[isym].tran[i];
 
-    if (tran[0] == 0 && tran[1] == 0 && tran[2] == 0) {
+    if (std::abs(tran[0]) < eps && std::abs(tran[1]) < eps && std::abs(tran[2]) < eps) {
         ret = true;
     } else {
         ret = false;