Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

/[conograph]/Conograph/trunk/src/lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.cc

Diff of /Conograph/trunk/src/lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 17 by rtomiyasu, Tue Apr 30 05:04:41 2013 UTC revision 25 by rtomiyasu, Mon Jul 7 02:35:51 2014 UTC
# Line 1  Line 1 
1  /*  /*
2   * The MIT License   * The MIT License
3    
4     Conograph (powder auto-indexing program)     BLDConograph (Bravais lattice determination module in Conograph)
5    
6  Copyright (c) <2012> <Ryoko Oishi-Tomiyasu, KEK>  Copyright (c) <2012> <Ryoko Oishi-Tomiyasu, KEK>
7    
# Line 24  OUT OF OR IN CONNECTION WITH THE SOFTWAR Line 24  OUT OF OR IN CONNECTION WITH THE SOFTWAR
24  THE SOFTWARE.  THE SOFTWARE.
25   *   *
26   */   */
27  #include "../utility_func/chToDouble.hh"  #include <limits>
28  #include "../utility_lattice_reduction/put_Minkowski_reduced_lattice.hh"  #include "../utility_lattice_reduction/put_Buerger_reduced_lattice.hh"
29  #include "OutputInfo.hh"  #include "check_equiv.hh"
 #include "ReducedLatticeToCheckBravais.hh"  
30  #include "LatticeFigureOfMeritToCheckSymmetry.hh"  #include "LatticeFigureOfMeritToCheckSymmetry.hh"
31    
 const string LatticeFigureOfMeritToCheckSymmetry::CS_LABEL[NUM_LS] =  
                 {       "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14" };  
   
 // default c'tor for GUI  
32  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry()  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry()
33          : m_label(-1),          : m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
           m_S_red( SymMat43_VCData( SymMat<VCData>(3), NRMat<Int4>(4,3) ) )  
34  {  {
         m_num_lattice_found = 0;  
35  }  }
36    
37    
38  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const Double& rhs)  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const Double& rhs)
39          : m_label(-1), m_latfom(rhs),          : m_latfom(rhs),
40            m_S_red( SymMat43_VCData( SymMat<VCData>(3), NRMat<Int4>(4,3) ) )            m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
41  {  {
42  }  }
43    
44    
45  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const BravaisType& brat,  LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const BravaisType& brat,
46                  const SymMat43_VCData& S,                  const SymMat43_Double& S)
47                  const ePeakShiftFunctionType& type,          : m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
48                  const Double& wave_length,  {
49                  const vector<ZParawError>& peak_shift_param_rad)          setLatticeConstants43(brat, S);
         : m_label(-1),  
           m_S_red( SymMat43_VCData( SymMat<VCData>(3), NRMat<Int4>(4,3) ) )  
 {  
         this->setLatticeConstants43(brat, S);  
         m_latfom.setPeakShiftParamRadian(type, wave_length, peak_shift_param_rad);  
         m_num_lattice_found = 0;  
50  }  }
51    
52    
53    
54    
55  #ifdef DEBUG  #ifdef DEBUG
56  static bool checkInitialLatticeParameters(  static bool checkInitialLatticeParameters(
57                  const BravaisType& brat,                  const BravaisType& brat,
58                  const SymMat43_VCData& S_red)                  const SymMat43_Double& S_red)
59  {  {
60          const SymMat<Double> dbl_S_red( chToDouble(S_red.first) );          const SymMat<Double> dbl_S_red( S_red.first );
61    
62          if( brat.enumLaueGroup() == Ci && brat.enumBravaisLattice() == Prim )          if( brat.enumLaueGroup() == Ci && brat.enumCentringType() == Prim )
63          {          {
64                  assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)                  assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
65                                  && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1)                                  && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1)
66                                  && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2)                                  && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2)
67                                  && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );                                  && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
68          }          }
69          else if( brat.enumLaueGroup() == C2h_Y && brat.enumBravaisLattice() == Prim )          else if( brat.enumLaueGroup() == C2h_Y && brat.enumCentringType() == Prim )
70          {          {
71                  assert( 0.0 <= dbl_S_red(0,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(0,0)                  assert( 0.0 <= dbl_S_red(0,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(0,0)
72                                  && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(0,0) );                                  && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(0,0) );
73          }          }
74          else if( brat.enumLaueGroup() == C2h_Z && brat.enumBravaisLattice() == Prim )          else if( brat.enumLaueGroup() == C2h_Z && brat.enumCentringType() == Prim )
75          {          {
76                  assert( 0.0 <= dbl_S_red(0,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)                  assert( 0.0 <= dbl_S_red(0,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
77                                  && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1) );                                  && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1) );
78          }          }
79          else if( brat.enumLaueGroup() == C2h_X && brat.enumBravaisLattice() == Prim )          else if( brat.enumLaueGroup() == C2h_X && brat.enumCentringType() == Prim )
80          {          {
81                  assert( 0.0 <= dbl_S_red(1,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1)                  assert( 0.0 <= dbl_S_red(1,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1)
82                                  && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );                                  && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
83          }          }
84          else if( brat.enumLaueGroup() == C2h_Y && brat.enumBravaisLattice() == BaseZ )          else if( brat.enumLaueGroup() == C2h_Y && brat.enumCentringType() == BaseZ )
85          {          {
86                  assert( 0.0 <= dbl_S_red(0,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 0.9999 < dbl_S_red(0,0) );                  assert( 0.0 <= dbl_S_red(0,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 0.9999 < dbl_S_red(0,0) );
87          }          }
88          else if( brat.enumLaueGroup() == C2h_Z && brat.enumBravaisLattice() == BaseX )          else if( brat.enumLaueGroup() == C2h_Z && brat.enumCentringType() == BaseX )
89          {          {
90                  assert( 0.0 <= dbl_S_red(0,1) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 0.9999 < dbl_S_red(1,1) );                  assert( 0.0 <= dbl_S_red(0,1) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 0.9999 < dbl_S_red(1,1) );
91          }          }
92          else if( brat.enumLaueGroup() == C2h_X && brat.enumBravaisLattice() == BaseY )          else if( brat.enumLaueGroup() == C2h_X && brat.enumCentringType() == BaseY )
93          {          {
94                  assert( 0.0 <= dbl_S_red(1,2) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 0.9999  < dbl_S_red(2,2) );                  assert( 0.0 <= dbl_S_red(1,2) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 0.9999  < dbl_S_red(2,2) );
95          }          }
96          else if( brat.enumLaueGroup() == D2h          else if( brat.enumLaueGroup() == D2h
97                          && brat.enumBravaisLattice() != BaseX                          && brat.enumCentringType() != BaseX
98                          && brat.enumBravaisLattice() != BaseY                          && brat.enumCentringType() != BaseY
99                          && brat.enumBravaisLattice() != BaseZ )                          && brat.enumCentringType() != BaseZ )
100          {          {
101                  assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0) );                  assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0) );
102          }          }
103    
104          const SymMat<VCData> S_super = transform_sym_matrix(S_red.second, S_red.first);          const SymMat<Double> S_super = transform_sym_matrix(S_red.second, S_red.first);
105          assert( S_super(0,1) <= VCData()          assert( S_super(0,1) <= 0.0
106                          && S_super(0,2) <= VCData()                          && S_super(0,2) <= 0.0
107                          && S_super(0,3) <= VCData()                          && S_super(0,3) <= 0.0
108                          && S_super(1,2) <= VCData()                          && S_super(1,2) <= 0.0
109                          && S_super(1,3) <= VCData()                          && S_super(1,3) <= 0.0
110                          && S_super(2,3) <= VCData() );                          && S_super(2,3) <= 0.0 );
   
         SymMat<VCData> S_red_cp = S_red.first;  
         cal_average_crystal_system(brat.enumLaueGroup(), S_red_cp);  
         assert( S_red.first(0,0).putVecCoef() == S_red_cp(0,0).putVecCoef() );  
         assert( S_red.first(1,1).putVecCoef() == S_red_cp(1,1).putVecCoef() );  
         assert( S_red.first(2,2).putVecCoef() == S_red_cp(2,2).putVecCoef() );  
         assert( S_red.first(0,1).putVecCoef() == S_red_cp(0,1).putVecCoef() );  
         assert( S_red.first(0,2).putVecCoef() == S_red_cp(0,2).putVecCoef() );  
         assert( S_red.first(1,2).putVecCoef() == S_red_cp(1,2).putVecCoef() );  
111    
112          return true;          return true;
113  }  }
114  #endif  #endif
115    
116    
117  void LatticeFigureOfMeritToCheckSymmetry::setLatticeConstants43(const BravaisType& brat, const SymMat43_VCData& S)  void LatticeFigureOfMeritToCheckSymmetry::setLatticeConstants43(const BravaisType& brat, const SymMat43_Double& S)
118  {  {
119          m_S_red = S;          m_S_red = S;
120          assert( checkInitialLatticeParameters(brat, m_S_red) );          assert( checkInitialLatticeParameters(brat, m_S_red) );
121    
122          m_latfom.setLatticeConstants43(brat, SymMat43_Double(chToDouble(m_S_red.first), m_S_red.second));          m_latfom.setLatticeConstants43(brat, S);
         m_num_lattice_found = 0;  
 //      for(Int4 i=0; i<NUM_LS; i++) m_lattice_equiv[i].clear();  
123  }  }
124    
125    
126    
127    static bool operator<(const SymMat<Double>& lhs, const SymMat<Double>& rhs)
128    {
129            static const Double EPS_1 = 1.0+sqrt( numeric_limits<double>::epsilon() );
130            assert( lhs.size() == 3 );
131            assert( rhs.size() == 3 );
132    
133            static const Int4 ISIZE = 3;
134            for(Int4 i=0; i<ISIZE; i++)
135            {
136                    if( lhs(i,i)*EPS_1 < rhs(i,i) ) return true;
137                    if( rhs(i,i)*EPS_1 < lhs(i,i) ) return false;
138    
139                    for(Int4 j=0; j<i; j++)
140                    {
141                            const Double lhs_ij = lhs(i,i)+lhs(j,j)+lhs(i,j)*2.0;
142                            const Double rhs_ij = rhs(i,i)+rhs(j,j)+rhs(i,j)*2.0;
143    
144                            if( lhs_ij*EPS_1 < rhs_ij ) return true;
145                            if( rhs_ij*EPS_1 < lhs_ij ) return false;
146                    }
147            }
148    
149            return false;
150    }
151    
152    
153    
154  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsMonoclinic(const ePointGroup& epg_new,  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsMonoclinic(const ePointGroup& epg_new,
155                  const Double& cv2,                  const Double& resol,
156                  map< SymMat<VCData>, NRMat<Int4> >& ans) const                  map< SymMat<Double>, NRMat<Int4> >& ans) const
157  {  {
158          ans.clear();          ans.clear();
159    
160          SymMat<VCData> ans0 = m_S_red.first;          SymMat<Double> ans0 = m_S_red.first;
161          cal_average_crystal_system(C2h_X, ans0);          cal_average_crystal_system(C2h_X, ans0);
162    
163          SymMat<VCData> S_red(3);          SymMat<Double> S_red(3);
164          NRMat<Int4> trans_mat2;          NRMat<Int4> trans_mat2;
165          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
166          {          {
167                  if( epg_new == C2h_X )                  if( epg_new == C2h_X )
168                  {                  {
169                          S_red = ans0;                          S_red = ans0;
170                          trans_mat2 = m_S_red.second;                          trans_mat2 = m_S_red.second;
171                          putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(1, 2, S_red, trans_mat2);
172                  }                  }
173                  else if( epg_new == C2h_Y )                  else if( epg_new == C2h_Y )
174                  {                  {
175                          S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);                          S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
176                          trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());                          trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
177                          putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 2, S_red, trans_mat2);
178                  }                  }
179                  else // if( epg_new == C2h_Z )                  else // if( epg_new == C2h_Z )
180                  {                  {
181                          S_red = transform_sym_matrix(put_matrix_YZX(), ans0);                          S_red = transform_sym_matrix(put_matrix_YZX(), ans0);
182                          trans_mat2 = mprod(m_S_red.second, put_matrix_ZXY());                          trans_mat2 = mprod(m_S_red.second, put_matrix_ZXY());
183                          putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 1, S_red, trans_mat2);
184                  }                  }
185                  ans.insert( SymMat43_VCData( S_red, trans_mat2) );                  ans.insert( SymMat43_Double( S_red, trans_mat2) );
186          }          }
187    
188          ans0 = m_S_red.first;          ans0 = m_S_red.first;
189          cal_average_crystal_system(C2h_Y, ans0);          cal_average_crystal_system(C2h_Y, ans0);
190          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
191          {          {
192                  if( epg_new == C2h_X )                  if( epg_new == C2h_X )
193                  {                  {
194                          S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);                          S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
195                          trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());                          trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
196                          putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(1, 2, S_red, trans_mat2);
197                  }                  }
198                  else if( epg_new == C2h_Y )                  else if( epg_new == C2h_Y )
199                  {                  {
200                          S_red = ans0;                          S_red = ans0;
201                          trans_mat2 = m_S_red.second;                          trans_mat2 = m_S_red.second;
202                          putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 2, S_red, trans_mat2);
203                  }                  }
204                  else // if( epg_new == C2h_Z )                  else // if( epg_new == C2h_Z )
205                  {                  {
206                          S_red = transform_sym_matrix(put_matrix_XZY(), ans0);                          S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
207                          trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());                          trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
208                          putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 1, S_red, trans_mat2);
209                  }                  }
210                  ans.insert( SymMat43_VCData( S_red, trans_mat2) );                  ans.insert( SymMat43_Double( S_red, trans_mat2) );
211          }          }
212    
213          ans0 = m_S_red.first;          ans0 = m_S_red.first;
214          cal_average_crystal_system(C2h_Z, ans0);          cal_average_crystal_system(C2h_Z, ans0);
215          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
216          {          {
217                  if( epg_new == C2h_X )                  if( epg_new == C2h_X )
218                  {                  {
219                          S_red = transform_sym_matrix(put_matrix_ZXY(), ans0);                          S_red = transform_sym_matrix(put_matrix_ZXY(), ans0);
220                          trans_mat2 = mprod(m_S_red.second, put_matrix_YZX());                          trans_mat2 = mprod(m_S_red.second, put_matrix_YZX());
221                          putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(1, 2, S_red, trans_mat2);
222                  }                  }
223                  else if( epg_new == C2h_Y )                  else if( epg_new == C2h_Y )
224                  {                  {
225                          S_red = transform_sym_matrix(put_matrix_XZY(), ans0);                          S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
226                          trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());                          trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
227                          putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 2, S_red, trans_mat2);
228                  }                  }
229                  else // if( epg_new == C2h_Z )                  else // if( epg_new == C2h_Z )
230                  {                  {
231                          S_red = ans0;                          S_red = ans0;
232                          trans_mat2 = m_S_red.second;                          trans_mat2 = m_S_red.second;
233                          putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);                          putBuergerReducedMonoclinicP(0, 1, S_red, trans_mat2);
234                  }                  }
235                  ans.insert( SymMat43_VCData( S_red, trans_mat2) );                  ans.insert( SymMat43_Double( S_red, trans_mat2) );
236          }          }
237    
238          return !( ans.empty() );          return !( ans.empty() );
239  }  }
240    
241    
242  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsOrthorhombic(const Double& cv2,  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsOrthorhombic(const Double& resol,
243                                                          map< SymMat<VCData>, NRMat<Int4> >& ans) const                                                          map< SymMat<Double>, NRMat<Int4> >& ans) const
244  {  {
245          ans.clear();          ans.clear();
246    
247          const BravaisType& brat = m_latfom.putBravaisType();          const BravaisType& brat = m_latfom.putBravaisType();
248    
249          SymMat<VCData> ans0 = m_S_red.first;          SymMat<Double> ans0 = m_S_red.first;
250          cal_average_crystal_system(D2h, ans0);          cal_average_crystal_system(D2h, ans0);
251          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
252          {          {
253                  if( brat.enumBravaisLattice() == BaseX )                  if( brat.enumCentringType() == BaseX )
254                  {                  {
255                          if( ans0(1,1) < ans0(2,2) )                          if( ans0(1,1) < ans0(2,2) )
256                          {                          {
257                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_ZYX(), ans0), mprod( m_S_red.second, put_matrix_ZYX() ) ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_ZYX(), ans0), mprod( m_S_red.second, put_matrix_ZYX() ) ) );
258                          }                          }
259                          else                          else
260                          {                          {
261                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
262                          }                          }
263                  }                  }
264                  else if( brat.enumBravaisLattice() == BaseY )                  else if( brat.enumCentringType() == BaseY )
265                  {                  {
266                          if( ans0(0,0) < ans0(2,2) )                          if( ans0(0,0) < ans0(2,2) )
267                          {                          {
268                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_ZXY(), ans0), mprod( m_S_red.second, put_matrix_YZX() ) ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_ZXY(), ans0), mprod( m_S_red.second, put_matrix_YZX() ) ) );
269                          }                          }
270                          else                          else
271                          {                          {
272                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
273                          }                          }
274                  }                  }
275                  else if( brat.enumBravaisLattice() == BaseZ )                  else if( brat.enumCentringType() == BaseZ )
276                  {                  {
277                          if( ans0(0,0) < ans0(1,1) )                          if( ans0(0,0) < ans0(1,1) )
278                          {                          {
279                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_YXZ(), ans0), mprod( m_S_red.second, put_matrix_YXZ() ) ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_YXZ(), ans0), mprod( m_S_red.second, put_matrix_YXZ() ) ) );
280                          }                          }
281                          else                          else
282                          {                          {
283                                  ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_XYZ(), ans0), m_S_red.second ) );                                  ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_XYZ(), ans0), m_S_red.second ) );
284                          }                          }
285                  }                  }
286                  else                  else
287                  {                  {
288                          NRMat<Int4> trans_mat = m_S_red.second;                          NRMat<Int4> trans_mat = m_S_red.second;
289                          putMinkowskiReducedOrthorhombic(brat.enumBravaisLattice(), ans0, trans_mat);                          putBuergerReducedOrthorhombic(brat.enumCentringType(), ans0, trans_mat);
290                          ans.insert( SymMat43_VCData(ans0, trans_mat ) );                          ans.insert( SymMat43_Double(ans0, trans_mat ) );
291                  }                  }
292                  return true;                  return true;
293          }          }
# Line 291  bool LatticeFigureOfMeritToCheckSymmetry Line 295  bool LatticeFigureOfMeritToCheckSymmetry
295  }  }
296    
297    
298  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsTetragonal(const Double& cv2,  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsTetragonal(const Double& resol,
299                  map< SymMat<VCData>, NRMat<Int4> >& ans) const                  map< SymMat<Double>, NRMat<Int4> >& ans) const
300  {  {
301          ans.clear();          ans.clear();
302    
303          SymMat<VCData> ans0 = m_S_red.first;          SymMat<Double> ans0 = m_S_red.first;
304          cal_average_crystal_system(D4h_X, ans0);          cal_average_crystal_system(D4h_X, ans0);
305          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
306          {          {
307                  ans.insert( SymMat43_VCData(                  ans.insert( SymMat43_Double(
308                                  transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );                                  transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
309          }          }
310    
311          ans0 = m_S_red.first;          ans0 = m_S_red.first;
312          cal_average_crystal_system(D4h_Y, ans0);          cal_average_crystal_system(D4h_Y, ans0);
313          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
314          {          {
315                  ans.insert( SymMat43_VCData(                  ans.insert( SymMat43_Double(
316                                  transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );                                  transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
317          }          }
318    
319          ans0 = m_S_red.first;          ans0 = m_S_red.first;
320          cal_average_crystal_system(D4h_Z, ans0);          cal_average_crystal_system(D4h_Z, ans0);
321          if( check_equiv_m(ans0, m_S_red.first, cv2 ) )          if( check_equiv_m(ans0, m_S_red.first, resol ) )
322          {          {
323                  ans.insert( SymMat43_VCData(ans0, m_S_red.second ) );                  ans.insert( SymMat43_Double(ans0, m_S_red.second ) );
324          }          }
325    
326          return !( ans.empty() );          return !( ans.empty() );
# Line 325  bool LatticeFigureOfMeritToCheckSymmetry Line 329  bool LatticeFigureOfMeritToCheckSymmetry
329    
330    
331    
332  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsHexagonal(const ePointGroup& epg_new, const Double& cv2,  bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsHexagonal(const ePointGroup& epg_new, const Double& resol,
333                  map< SymMat<VCData>, NRMat<Int4> >& ans) const                  map< SymMat<Double>, NRMat<Int4> >& ans) const
334  {  {
335          ans.clear();          ans.clear();
336          const BravaisType& brat = m_latfom.putBravaisType();          const BravaisType& brat = m_latfom.putBravaisType();
337    
338          SymMat43_VCData ans2(SymMat<VCData>(3), NRMat<Int4>(3,3));          SymMat43_Double ans2(SymMat<Double>(3), NRMat<Int4>(3,3));
339    
340          if( brat.enumLaueGroup() == C2h_X )          if( brat.enumLaueGroup() == C2h_X )
341          {          {
# Line 349  bool LatticeFigureOfMeritToCheckSymmetry Line 353  bool LatticeFigureOfMeritToCheckSymmetry
353                  ans2.second = m_S_red.second;                  ans2.second = m_S_red.second;
354          }          }
355    
356          if( ans2.first(0,1) < VCData() )          if( ans2.first(0,1) < 0.0 )
357          {          {
358                  ans2.first(0,1) *= -1;                  ans2.first(0,1) *= -1;
359                  ans2.second[0][0] *= -1;                  ans2.second[0][0] *= -1;
# Line 357  bool LatticeFigureOfMeritToCheckSymmetry Line 361  bool LatticeFigureOfMeritToCheckSymmetry
361                  ans2.second[2][0] *= -1;                  ans2.second[2][0] *= -1;
362          }          }
363    
364          SymMat<VCData> ans0 = ans2.first;          SymMat<Double> ans0 = ans2.first;
365          cal_average_crystal_system(epg_new, ans2.first);          cal_average_crystal_system(epg_new, ans2.first);
366          if( check_equiv_m(ans2.first, ans0, cv2 ) )          if( check_equiv_m(ans2.first, ans0, resol ) )
367          {          {
368                  ans.insert( ans2 );                  ans.insert( ans2 );
369                  return true;                  return true;
# Line 368  bool LatticeFigureOfMeritToCheckSymmetry Line 372  bool LatticeFigureOfMeritToCheckSymmetry
372  }  }
373    
374    
375  bool LatticeFigureOfMeritToCheckSymmetry::checkLatticeSymmetry(const ePointGroup& epg_new, const Double& cv2,  bool LatticeFigureOfMeritToCheckSymmetry::checkLatticeSymmetry(const ePointGroup& epg_new, const Double& resol,
376                  map< SymMat<VCData>, NRMat<Int4> >& ans) const                  map< SymMat<Double>, NRMat<Int4> >& ans) const
377  {  {
378          ans.clear();          ans.clear();
379          const BravaisType& brat = m_latfom.putBravaisType();          const BravaisType& brat = m_latfom.putBravaisType();
# Line 382  bool LatticeFigureOfMeritToCheckSymmetry Line 386  bool LatticeFigureOfMeritToCheckSymmetry
386          if( epg_new == C2h_X || epg_new == C2h_Y ||  epg_new == C2h_Z )          if( epg_new == C2h_X || epg_new == C2h_Y ||  epg_new == C2h_Z )
387          {          {
388                  assert( brat.enumLaueGroup() == Ci );                  assert( brat.enumLaueGroup() == Ci );
389                  assert( brat.enumBravaisLattice() == Prim );                  assert( brat.enumCentringType() == Prim );
390    
391                  return checkIfLatticeIsMonoclinic(epg_new, cv2, ans);                  return checkIfLatticeIsMonoclinic(epg_new, resol, ans);
392          }          }
393          else if( epg_new == D4h_Z )          else if( epg_new == D4h_Z )
394          {          {
395                  assert( brat.enumLaueGroup() == D2h );                  assert( brat.enumLaueGroup() == D2h );
396                  assert( brat.enumBravaisLattice() == Prim                  assert( brat.enumCentringType() == Prim
397                                  || brat.enumBravaisLattice() == Inner );                                  || brat.enumCentringType() == Inner );
398    
399                  return checkIfLatticeIsTetragonal(cv2, ans);                  return checkIfLatticeIsTetragonal(resol, ans);
400          }          }
401          else if( epg_new == D2h )          else if( epg_new == D2h )
402          {          {
403                  assert( brat.enumLaueGroup() != Ci || brat.enumBravaisLattice() == Prim );                  assert( brat.enumLaueGroup() != Ci || brat.enumCentringType() == Prim );
404                  assert( brat.enumLaueGroup() != C2h_Z || brat.enumBravaisLattice() == BaseX );                  assert( brat.enumLaueGroup() != C2h_Z || brat.enumCentringType() == BaseX );
405                  assert( brat.enumLaueGroup() != C2h_X || brat.enumBravaisLattice() == BaseY );                  assert( brat.enumLaueGroup() != C2h_X || brat.enumCentringType() == BaseY );
406                  assert( brat.enumLaueGroup() != C2h_Y || brat.enumBravaisLattice() == BaseZ );                  assert( brat.enumLaueGroup() != C2h_Y || brat.enumCentringType() == BaseZ );
407                  assert( brat.enumBravaisLattice() != Rhom_hex );                  assert( brat.enumCentringType() != Rhom_hex );
408    
409                  return checkIfLatticeIsOrthorhombic(cv2, ans);                  return checkIfLatticeIsOrthorhombic(resol, ans);
410          }          }
411          else if( epg_new == D6h )          else if( epg_new == D6h )
412          {          {
413                  assert( brat.enumBravaisLattice() == Prim );                  assert( brat.enumCentringType() == Prim );
414                  assert( brat.enumLaueGroup() == C2h_X                  assert( brat.enumLaueGroup() == C2h_X
415                                  || brat.enumLaueGroup() == C2h_Y                                  || brat.enumLaueGroup() == C2h_Y
416                                  || brat.enumLaueGroup() == C2h_Z );                                  || brat.enumLaueGroup() == C2h_Z );
417                  return checkIfLatticeIsHexagonal(epg_new, cv2, ans);                  return checkIfLatticeIsHexagonal(epg_new, resol, ans);
418          }          }
419          else          else
420          {          {
421                  assert( epg_new == Oh );                  assert( epg_new == Oh );
422                  assert( brat.enumBravaisLattice() == Prim                  assert( brat.enumCentringType() == Prim
423                                  || brat.enumBravaisLattice() == Inner                                  || brat.enumCentringType() == Inner
424                                  || brat.enumBravaisLattice() == Face );                                  || brat.enumCentringType() == Face );
425    
426                  SymMat43_VCData ans2 = m_S_red;                  SymMat43_Double ans2 = m_S_red;
427                  cal_average_crystal_system(epg_new, ans2.first);                  cal_average_crystal_system(epg_new, ans2.first);
428                  if( check_equiv_m(ans2.first, m_S_red.first, cv2 ) )                  if( check_equiv_m(ans2.first, m_S_red.first, resol ) )
429                  {                  {
430                          ans.insert( ans2 );                          ans.insert( ans2 );
431                          return true;                          return true;
# Line 431  bool LatticeFigureOfMeritToCheckSymmetry Line 435  bool LatticeFigureOfMeritToCheckSymmetry
435  }  }
436    
437    
438  void LatticeFigureOfMeritToCheckSymmetry::putEquivalentLatticeConstantsDegreeWithOtherCentring(  void LatticeFigureOfMeritToCheckSymmetry::putLatticesOfHigherSymmetry(
439                  const eABCaxis& abc_axis, const eRHaxis& rh_axis,                  const ePointGroup& epg, const Double& resol,
440                  vector< pair< eCrystalSystem, SymMat<Double> > >& ans) const                  vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
441  {  {
442          ans.clear();          lattice_result.clear();
443            map< SymMat<Double>, NRMat<Int4> > S_red_tray;
444            if( !this->checkLatticeSymmetry(epg, resol, S_red_tray) ) return;
445    
446          static const Double cv2 = 0.04;          const BravaisType& ebrat_original = this->putLatticeFigureOfMerit().putBravaisType();
447          static const Double resol2 = 0.06;          const eCentringType eblat = (ebrat_original.enumBravaisType()==Monoclinic_B?
448                                                                            (epg==D31d_rho?Prim:(epg==D3d_1_hex?Rhom_hex:BaseZ)):ebrat_original.enumCentringType());
449    
450          // Calculate figures of merit as triclinic          const NRMat<Int4> matrix_min_to_sell = this->putInitialForm().second;
         const ReducedLatticeToCheckBravais RLCB(abc_axis, rh_axis, false, cv2, this->putInitialForm());  
         const SymMat<Double> S_obtuse = this->putInitialSellingReducedForm();  
451    
452          if( m_latfom.enumCrystalSystem() != Rhombohedral )          SymMat<Double> S_super(4);
453          {          NRMat<Int4> trans_mat(4,3);
                 const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Rhombohedral, abc_axis, rh_axis));  
                 for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)  
                 {  
                         const SymMat<Double> S_red = chToDouble(it->first);  
                         if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;  
                         ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Rhombohedral, S_red) );  
                 }  
         }  
         if( m_latfom.enumBravaisLattice() != Face )  
         {  
                 const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Orthorhombic_F, abc_axis, rh_axis));  
                 for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)  
                 {  
                         const SymMat<Double> S_red = chToDouble(it->first);  
                         if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;  
                         ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Orthorhombic_F, S_red ) );  
                 }  
         }  
         if( m_latfom.enumBravaisLattice() != Inner )  
         {  
                 const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Orthorhombic_I, abc_axis, rh_axis));  
                 for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)  
                 {  
                         const SymMat<Double> S_red = chToDouble(it->first);  
                         if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;  
                         ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Orthorhombic_I, S_red ) );  
                 }  
         }  
         if( m_latfom.enumBravaisLattice() != BaseX && m_latfom.enumBravaisLattice() != BaseY && m_latfom.enumBravaisLattice() != BaseZ )  
         {  
                 const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Monoclinic_B, abc_axis, rh_axis));  
                 for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)  
                 {  
                         const SymMat<Double> S_red = chToDouble(it->first);  
                         if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;  
                         ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Monoclinic_B, S_red ) );  
                 }  
         }  
454    
455          NRMat<Int4> trans_mat;          for(map< SymMat<Double>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
         SymMat<Double> S_red(3);  
         if( m_latfom.enumCrystalSystem() == Rhombohedral || m_latfom.enumBravaisLattice() != Prim )  
456          {          {
457                  const SymMat<Double> S_super = put_sym_matrix_size4to3(this->putInitialSellingReducedForm());                  // S_super = it->second * it->first * Transpose(it->second) is close to
458                  putTransformMatrixToMinkowskiReduced(Inverse3(S_super), trans_mat);                  // Delone-reduced form of the original lattice.
459                  transpose_square_matrix(trans_mat);                  S_super = transform_sym_matrix(it->second, it->first );
                 ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Triclinic, transform_sym_matrix(Inverse3(trans_mat), S_super) ) );  
         }  
 }  
460    
461                    trans_mat = identity_matrix<Int4>(4);
462    
463  void LatticeFigureOfMeritToCheckSymmetry::putEquivalentLatticeConstantsDegreeWithOtherCentring(                  // S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second).
464                  const eABCaxis& abc_axis, const eRHaxis& rh_axis,                  put_Selling_reduced_dim_less_than_4(S_super, trans_mat);
465                  vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >& ans) const                  moveSmallerDiagonalLeftUpper(S_super, trans_mat);
 {  
         vector< pair< eCrystalSystem, SymMat<Double> > > ans0;  
         putEquivalentLatticeConstantsDegreeWithOtherCentring(abc_axis, rh_axis, ans0);  
466    
467          ans.clear();                  lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( BravaisType( pair<eCentringType, ePointGroup>(eblat, epg) ),
468          ans.resize( ans0.size() );                                                                                  SymMat43_Double(it->first, mprod(trans_mat, it->second) ) ) );
469          vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >::iterator it2 = ans.begin();          }
         for(vector< pair< eCrystalSystem, SymMat<Double> > >::const_iterator it=ans0.begin(); it<ans0.end(); it++, it2++)  
         {  
                 it2->first = it->first;  
                 LatticeFigureOfMerit::putLatticeConstantsDegree( BravaisType(it->first, abc_axis, rh_axis), it->second, abc_axis, rh_axis, it2->second.first, it2->second.second );  
         }  
 }  
   
   
 void LatticeFigureOfMeritToCheckSymmetry::printLatticeInformation(  
                                         const vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[],  
                                         const OutputInfo outinfo[],  
                                         const eABCaxis& abc_axis,  
                                         const eRHaxis& rh_axis,  
                                         const Int4& label_start0,  
                                         ostream* os) const  
 {  
         m_latfom.printLatticeInformation(abc_axis, rh_axis, label_start0, os);  
   
 //      const FracMat mat_sell_to_min = FInverse3( put_transform_matrix_row4to3( putLatticeFigureOfMerit().putOptimizedFormWithTransformMatrixToSellingReduced().second) );  
 //      const SymMat<Double> dbl_S_init = transform_sym_matrix(mat_sell_to_min.mat,  
 //                                                                              put_sym_matrix_size4to3( this->putInitialSellingReducedForm() ) ) / (mat_sell_to_min.denom*mat_sell_to_min.denom);  
 //      const SymMat<Double>& dbl_S_init = this->putInitialSellingReducedForm();  
   
 //      os->width(label_start); *os << "";  
 //      *os << "<!-- A*, B*, C*, D*, E*, F*(angstrom^(-2)) first given by peak-positions.-->\n";  
 //      os->width(label_start); *os << "";  
 //      *os << "<InitialReciprocalLatticeParameters>";  
 //      os->width(14);  
 //      *os << dbl_S_init(0,0);  
 //      os->width(14);  
 //      *os << dbl_S_init(1,1);  
 //      os->width(14);  
 //      *os << dbl_S_init(2,2);  
 //      os->width(14);  
 //      *os << dbl_S_init(1,2);  
 //      os->width(14);  
 //      *os << dbl_S_init(0,2);  
 //      os->width(14);  
 //      *os << dbl_S_init(0,1);  
 //      *os << " </InitialReciprocalLatticeParameters>\n";  
   
         Int4 label_start = label_start0;  
         os->width(label_start);  
         *os << "" << "<NumberOfLatticesInNeighborhood>";  
         os->width(14);  
         *os << this->putNumberOfLatticesInNeighborhood();  
         *os << " </NumberOfLatticesInNeighborhood>\n\n";  
   
         os->width(label_start);  
         *os << "" << "<EquivalentLatticeCandidates>\n";  
         label_start++;  
   
         vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > > lattice_equiv;  
         this->putEquivalentLatticeConstantsDegreeWithOtherCentring(abc_axis, rh_axis, lattice_equiv);  
   
                 for(vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >::const_iterator it=lattice_equiv.begin(); it<lattice_equiv.end(); it++)  
                 {  
                         os->width(label_start); *os << "";  
                 *os << "<LatticeCandidate>\n";  
                         label_start++;  
   
                         os->width(label_start); *os << "";  
                 *os << "<CrystalSystem>";  
                         os->width(17);  
                         *os << put_cs_name(it->first, abc_axis);  
                 *os << " </CrystalSystem>\n";  
   
                         os->width(label_start);  
                 *os << "" << "<LatticeParameters>";  
                 os->width(14);  
                 *os << it->second.first[0];  
                 os->width(14);  
                         *os << it->second.first[1];  
                 os->width(14);  
                         *os << it->second.first[2];  
                 os->width(14);  
                         *os << it->second.second[0];  
                 os->width(14);  
                         *os << it->second.second[1];  
                 os->width(14);  
                         *os << it->second.second[2];  
                 *os << " </LatticeParameters>\n";  
   
                         label_start--;  
                         os->width(label_start); *os << "";  
                 *os << "</LatticeCandidate>\n";  
                 }  
   
         label_start--;  
         os->width(label_start); *os << "";  
         *os << "</EquivalentLatticeCandidates>\n\n";  
470  }  }
Legend:
Removed from v.17  
changed lines
  Added in v.25
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help