Added an option to enforce the ASR for harmonic terms.

To use this function, please set force_asr = true in the examine_translational_invariance().
ttadano · Dec 28, 2015 · f14d546 · f14d546
1 parent 285a347
commit f14d546
Show file tree

Hide file tree

Showing 3 changed files with 47 additions and 2 deletions.
diff --git a/anphon/dynamical.h b/anphon/dynamical.h
@@ -96,3 +96,4 @@ namespace PHON_NS {
             std::complex<double> *beta, std::complex<double> *c, int *ldc);
     }
 }
+
diff --git a/anphon/fcs_phonon.cpp b/anphon/fcs_phonon.cpp
@@ -389,7 +389,7 @@ void Fcs_phonon::MPI_Bcast_fc2_ext()
 
 
 void Fcs_phonon::examine_translational_invariance(const int n, const unsigned int nat, const unsigned int natmin, 
-                                                  double *ret, std::vector<FcsClassExtent> fc2,
+                                                  double *ret, std::vector<FcsClassExtent> &fc2,
                                                   std::vector<FcsArrayWithCell> *fcs)
 {
     int i, j, k, l, m;
@@ -400,6 +400,11 @@ void Fcs_phonon::examine_translational_invariance(const int n, const unsigned in
     double ***sum3;
     double ****sum4;
 
+    bool force_asr = false;
+    FcsClassExtent fc2_tmp;
+    std::vector<FcsClassExtent>::iterator it_target;
+
+
     for (i = 0; i < n; ++i) ret[i] = 0.0;
 
     for (i = 0; i < n; ++i) {
@@ -419,6 +424,39 @@ void Fcs_phonon::examine_translational_invariance(const int n, const unsigned in
                   sum2[3 * (*it).atm1 + (*it).xyz1][(*it).xyz2] += (*it).fcs_val;
             }
 
+            if (force_asr) {
+                std::cout << "  force_asr = true: Modify harmonic force constans so that the ASR is satisfied." << std::endl;
+                for (j = 0; j < natmin; ++j) {
+                    for (k = 0; k < 3; ++k) {
+                        for (m = 0; m < 3; ++m) {
+                            fc2_tmp.atm1 = j;
+                            fc2_tmp.xyz1 = k;
+                            fc2_tmp.atm2 = system->map_p2s[j][0];
+                            fc2_tmp.xyz2 = m;
+                            fc2_tmp.cell_s = 0;
+                            fc2_tmp.fcs_val = sum2[3*j+k][m];
+                            it_target = std::find(fc2.begin(),fc2.end(),fc2_tmp);
+                            if (std::abs(fc2_tmp.fcs_val) > eps12) {
+                                if (it_target != fc2.end()) {
+                                    fc2[it_target-fc2.begin()].fcs_val -= fc2_tmp.fcs_val;
+                                } else {
+                                    error->exit("examine_translational_invariance", "Corresponding IFC not found.");
+                                }
+                            }
+                        }
+                    }
+                }
+                for (j = 0; j < 3 * natmin; ++j) {
+                    for (k = 0; k < 3; ++k) {
+                        sum2[j][k] = 0.0;
+                    }
+                }
+
+                for (std::vector<FcsClassExtent>::const_iterator it = fc2.begin(); it != fc2.end(); ++it) {
+                    sum2[3 * (*it).atm1 + (*it).xyz1][(*it).xyz2] += (*it).fcs_val;
+                }
+            }
+
             for (j = 0; j < 3 * natmin; ++j) {
                 for (k = 0; k < 3; ++k) {
                     dev = std::abs(sum2[j][k]);

diff --git a/anphon/fcs_phonon.h b/anphon/fcs_phonon.h
@@ -79,6 +79,12 @@ namespace PHON_NS {
             cell_s = obj.cell_s;
             fcs_val = obj.fcs_val;
         }
+
+        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;
+             }
     };
 
     struct AtomCellSuper {
@@ -141,7 +147,7 @@ namespace PHON_NS {
         void load_fcs_xml();
 
         void examine_translational_invariance(const int, const unsigned int, const unsigned int,
-            double *, std::vector<FcsClassExtent>, std::vector<FcsArrayWithCell> *);
+            double *, std::vector<FcsClassExtent> &, std::vector<FcsArrayWithCell> *);
 
 
         void MPI_Bcast_fc_class(const unsigned int);