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

…elop

ChangeLog.md

@@ -1,3 +1,22 @@
+# Ver. 0.9.5 (2015-06-28)
+
+## New
+
+- PERIODIC tag for low dimensional systems
+
+## Changes
+
+- Print distances for anharmonic terms in PREFIX.fcs
+
+- Improved the way to define the multiplicity of force constants
+
+- OpenMP parallelization for generating constraints between force constants
+
+## Fix
+
+- Fixed a minor bug in interaction.cpp
+
+- Fixed a bug in displace.py for QE
 
 # Ver. 0.9.4 (2015-02-16)
 

README.md

@@ -1,5 +1,5 @@
 # ALAMODE
-* Version 0.9.4 (Beta)
+* Version 0.9.5 (Beta)
 
 - - -
 
@@ -38,4 +38,4 @@ This software is released under the MIT license.
 For license rights and limitations, see LICENSE.txt file.
 
 ## Author
-Terumasa Tadano (Department of Physics, The University of Tokyo, Japan)
+Terumasa Tadano (Department of Applied Physics, The University of Tokyo, Japan)

alm/fcs.cpp

@@ -149,6 +149,12 @@ void Fcs::generate_fclists(int maxorder)
                     if (!symmetry->sym_available[isym]) continue;
 
                     for (i = 0; i < order + 2; ++i) atmn_mapped[i] = symmetry->map_sym[atmn[i]][isym];
+                    // Special treatment for non-periodic cases
+                    if (!(interaction->is_periodic[0] & 
+                          interaction->is_periodic[1] & 
+                          interaction->is_periodic[2])) {
+                        if (!interaction->is_incutoff(order+2, atmn_mapped)) continue;
+                    }
                     if (!is_inprim(order + 2, atmn_mapped)) continue;
 
                     for (i2 = 0; i2 < nxyz; ++i2) {

alm/input.cpp

@@ -90,8 +90,8 @@ void Input::parse_general_vars()
     int i;
     std::string prefix, mode, str_tmp, str_disp_basis;
     int nat, nkd, nsym;
-    int is_printsymmetry;
-    bool is_periodic[3], trim_dispsign_for_evenfunc;
+    int is_printsymmetry, is_periodic[3];
+    bool trim_dispsign_for_evenfunc;
     std::string *kdname;
     double tolerance;
 

alm/interaction.cpp

@@ -62,7 +62,7 @@ void Interaction::init()
         std::cout << "  " <<  std::setw(9) << str_order[i] << std::endl; 
         for (j = 0; j < nkd; ++j) {
             for (k = 0; k < nkd; ++k) {
-                if (rcs[i][j][k] < 0) {
+                if (rcs[i][j][k] < 0.0) {
                     std::cout << std::setw(9) << "None";
                 } else {
                     std::cout << std::setw(9) << rcs[i][j][k];
@@ -438,7 +438,7 @@ bool Interaction::is_incutoff(const int n, int *atomnumlist)
 {
     int i, j;
     int iat, jat, kat;
-    int jkd, kkd;
+    int ikd, jkd, kkd;
     int ncheck = n - 1;
     int order = n - 2;   
     bool in_cutoff_tmp;
@@ -447,17 +447,24 @@ bool Interaction::is_incutoff(const int n, int *atomnumlist)
     std::vector<DistInfo>::const_iterator it, it2;
 
     iat = atomnumlist[0];
+    ikd = system->kd[iat] - 1;
 
     for (i = 0; i < ncheck; ++i) {
 
         jat = atomnumlist[i + 1];
         jkd = system->kd[jat] - 1;
 
+        if (rcs[order][ikd][jkd] >= 0.0 && 
+            (mindist_pairs[iat][jat][0].dist > rcs[order][ikd][jkd])) return false;
+
         for (j = i + 1; j < ncheck; ++j) {
 
             kat = atomnumlist[j + 1];
             kkd = system->kd[kat] - 1;
 
+            if (rcs[order][ikd][kkd] >= 0.0 && 
+                (mindist_pairs[iat][kat][0].dist > rcs[order][ikd][kkd])) return false;
+
             cutoff_tmp = rcs[order][jkd][kkd];
 
             if (cutoff_tmp >= 0.0) {
@@ -644,7 +651,8 @@ void Interaction::calc_mindist_clusters(std::vector<int> **interaction_pair_in,
 }
 
 
-void Interaction::cell_combination(std::vector<std::vector<int> > array, int i, std::vector<int> accum, std::vector<std::vector<int> > &comb)
+void Interaction::cell_combination(std::vector<std::vector<int> > array, int i, 
+                                   std::vector<int> accum, std::vector<std::vector<int> > &comb)
 {
     if (i == array.size())  {
         comb.push_back(accum); 

alm/interaction.h

@@ -133,10 +133,9 @@ namespace ALM_NS {
         Interaction(class ALM *);
         ~Interaction();
 
-        bool is_periodic[3];
+        int is_periodic[3];
         int nneib;
         int maxorder;
-  //      int interaction_type;
 
         int *nbody_include;
 

alm/patterndisp.cpp

@@ -119,7 +119,7 @@ void Displace::gen_displacement_pattern()
 
                 group_tmp.clear();
 
-                // Store first order + 1 indexes as a necessarily displacement pattern.
+                // Store first order + 1 indexes as a necessary displacement pattern.
                 // Here, duplicate entries will be removed. 
                 // For example, (iij) will be reduced to (ij).
                 for (j = 0; j < order + 1; ++j) {

alm/symmetry.cpp

@@ -18,6 +18,7 @@
 #include "constants.h"
 #include "timer.h"
 #include "error.h"
+#include "interaction.h"
 #include "files.h"
 #include <cmath>
 #include <vector>
@@ -55,7 +56,8 @@ void Symmetry::init()
     std::cout << " SYMMETRY" << std::endl;
     std::cout << " ========" << std::endl << std::endl;
 
-    setup_symmetry_operation(nat, nsym, nnp, system->lavec, system->rlavec, system->xcoord, system->kd);
+    setup_symmetry_operation(nat, nsym, nnp, system->lavec, system->rlavec, 
+                             system->xcoord, system->kd);
 
     memory->allocate(tnons, nsym, 3);
     memory->allocate(symrel_int, nsym, 3, 3);
@@ -119,7 +121,8 @@ void Symmetry::init()
     std::cout << std::endl;
 }
 
-void Symmetry::setup_symmetry_operation(int nat, unsigned int &nsym, unsigned int &nnp, double aa[3][3], double bb[3][3], double **x, int *kd)
+void Symmetry::setup_symmetry_operation(int nat, unsigned int &nsym, unsigned int &nnp, 
+                                        double aa[3][3], double bb[3][3], double **x, int *kd)
 {
     int i, j;
 
@@ -220,7 +223,8 @@ void Symmetry::findsym(int nat, double aa[3][3], double **x, std::vector<Symmetr
     LatticeSymmList.clear();
     find_lattice_symmetry(aa, LatticeSymmList);
     CrystalSymmList.clear();
-    find_crystal_symmetry(nat, system->nclassatom, system->atomlist_class, x, LatticeSymmList, CrystalSymmList);
+    find_crystal_symmetry(nat, system->nclassatom, system->atomlist_class, x,
+                          LatticeSymmList, CrystalSymmList);
 
     //	findsym(nat, nnp, kd, aa, bb, x);
     find_nnp_for_translation(nnp, CrystalSymmList);
@@ -361,7 +365,9 @@ void Symmetry::find_lattice_symmetry(double aa[3][3], std::vector<RotationMatrix
     }
 }
 
-void Symmetry::find_crystal_symmetry(int nat, int nclass, std::vector<unsigned int> *atomclass, double **x, std::vector<RotationMatrix> LatticeSymmList, std::vector<SymmetryOperationTransFloat> &CrystalSymmList)
+void Symmetry::find_crystal_symmetry(int nat, int nclass, std::vector<unsigned int> *atomclass, double **x, 
+                                     std::vector<RotationMatrix> LatticeSymmList, 
+                                     std::vector<SymmetryOperationTransFloat> &CrystalSymmList)
 {
     unsigned int i, j;
     unsigned int iat, jat, kat, lat;
@@ -398,7 +404,8 @@ void Symmetry::find_crystal_symmetry(int nat, int nclass, std::vector<unsigned i
     CrystalSymmList.push_back(SymmetryOperationTransFloat(rot_int, tran));
 
 
-    for (std::vector<RotationMatrix>::iterator it_latsym = LatticeSymmList.begin(); it_latsym != LatticeSymmList.end(); ++it_latsym) {
+    for (std::vector<RotationMatrix>::iterator it_latsym = LatticeSymmList.begin(); 
+                                               it_latsym != LatticeSymmList.end(); ++it_latsym) {
 
         iat = atomclass[0][0];
 
@@ -421,16 +428,19 @@ void Symmetry::find_crystal_symmetry(int nat, int nclass, std::vector<unsigned i
                 tran[i] = tran[i] - nint(tran[i]);
             }
 
-            isok = true;
+            if ((std::abs(tran[0]) > eps12 && !interaction->is_periodic[0]) ||
+                (std::abs(tran[1]) > eps12 && !interaction->is_periodic[1]) ||
+                (std::abs(tran[2]) > eps12 && !interaction->is_periodic[2])) continue;
 
             is_identity_matrix = 
                 ( std::pow(rot[0][0] - 1.0, 2) + std::pow(rot[0][1], 2) + std::pow(rot[0][2], 2) 
                 + std::pow(rot[1][0], 2) + std::pow(rot[1][1] - 1.0, 2) + std::pow(rot[1][2], 2)
                 + std::pow(rot[2][0], 2) + std::pow(rot[2][1], 2) + std::pow(rot[2][2] - 1.0, 2)
                 + std::pow(tran[0], 2) + std::pow(tran[1], 2) + std::pow(tran[2], 2) ) < eps12;
-
             if (is_identity_matrix) continue;
 
+            isok = true;
+
             for (itype = 0; itype < nclass; ++itype) {
 
                 for (jj = 0; jj < atomclass[itype].size(); ++jj) {

alm/symmetry.h

@@ -97,7 +97,6 @@ namespace ALM_NS {
         void init();
         void setup_symmetry_operation(int, unsigned int&, unsigned int&, double[3][3], double[3][3], 
             double **, int *);
-        void findsym(int, double [3][3], double **, std::vector<SymmetryOperation> &);
 
         unsigned int nsym, nnp;
         int ntran, natmin;
@@ -136,6 +135,7 @@ namespace ALM_NS {
 
     private:
 
+        void findsym(int, double [3][3], double **, std::vector<SymmetryOperation> &);
         bool is_translation(int **);
 
         void symop_in_cart(double [3][3], double[3][3]);
@@ -144,7 +144,7 @@ namespace ALM_NS {
         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, 
+        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>);