Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

/[conograph]/Conograph/trunk/src/SortingLattice.cc

Annotation of /Conograph/trunk/src/SortingLattice.cc

Parent Directory Parent Directory | Revision Log Revision Log

Revision 16 - (hide annotations) (download) (as text)
Fri Apr 26 08:50:24 2013 UTC (10 years, 10 months ago) by rtomiyasu
File MIME type: text/x-c++src
File size: 26075 byte(s) 
Modifications for GUI development.
1 rtomiyasu 3 /*
2     * The MIT License
3    
4     Conograph (powder auto-indexing program)
5    
6     Copyright (c) <2012> <Ryoko Oishi-Tomiyasu, KEK>
7    
8     Permission is hereby granted, free of charge, to any person obtaining a copy
9     of this software and associated documentation files (the "Software"), to deal
10     in the Software without restriction, including without limitation the rights
11     to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12     copies of the Software, and to permit persons to whom the Software is
13     furnished to do so, subject to the following conditions:
14    
15     The above copyright notice and this permission notice shall be included in
16     all copies or substantial portions of the Software.
17    
18     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21     AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22     LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23     OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24     THE SOFTWARE.
25     *
26     */
27     #ifdef _OPENMP
28     # include <omp.h>
29     #endif
30     #include "utility_data_structure/index_set.hh"
31     #include "utility_func/chToDouble.hh"
32     #include "utility_func/transform_sym_matrix.hh"
33     #include "utility_lattice_reduction/super_basis3.hh"
34     #include "lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.hh"
35     #include "lattice_symmetry/ReducedLatticeToCheckBravais.hh"
36     #include "lattice_symmetry/ReducedLatticeToCheckEquiv.hh"
37     #include "p_out_indexing.hh"
38     #include "zerror_type/error_out.hh"
39     #include "zlog/zlog.hh"
40     #include "ControlParam.hh"
41     #include "SortingLattice.hh"
42     #include "utility_func/stopx.hh"
43    
44     const bool SortingLattice::m_DoesPrudentSymSearch = false;
45     const Double SortingLattice::m_cv2 = 0.5;
46    
47     SortingLattice::SortingLattice()
48     {
49     for(ArrayIndex i=0; i<NUM_LS; i++)
50     {
51     OutputSymmetry[i] = false;
52     JudgeSymmetry[i] = false;
53     }
54    
55     m_resol2 = 0.0;
56     m_num_ref_figure_of_merit = 20;
57     m_etype_peak_shift = kPeakShiftFunction_Type0;
58     m_WlengthX = 1.54056;
59     }
60    
61    
62     SortingLattice::~SortingLattice()
63     {
64     }
65    
66    
67     // Set the member variables.
68     void SortingLattice::setParam(const ControlParam& cont)
69     {
70     OutputSymmetry[(Int4)Triclinic] = cont.putOutputSymmetry(Triclinic);
71     JudgeSymmetry[(Int4)Triclinic] = false;
72     for(ArrayIndex i=1; i<NUM_LS; i++)
73     {
74     OutputSymmetry[i] = cont.putOutputSymmetry(eCrystalSystem(i));
75     JudgeSymmetry[i] = cont.putOutputSymmetry(eCrystalSystem(i));
76     }
77    
78     if( JudgeSymmetry[(Int4)Cubic_P] )
79     {
80     JudgeSymmetry[(Int4)Tetragonal_P] = true;
81     }
82     if( JudgeSymmetry[(Int4)Hexagonal] )
83     {
84     JudgeSymmetry[(Int4)Monoclinic_P] = true;
85     }
86     if( JudgeSymmetry[(Int4)Tetragonal_P] )
87     {
88     JudgeSymmetry[(Int4)Orthorhombic_P] = true;
89     }
90     if( JudgeSymmetry[(Int4)Orthorhombic_P] )
91     {
92     JudgeSymmetry[(Int4)Monoclinic_P] = true;
93     }
94    
95     if( JudgeSymmetry[(Int4)Orthorhombic_C] )
96     {
97     JudgeSymmetry[(Int4)Monoclinic_B] = true;
98     }
99    
100     if( JudgeSymmetry[(Int4)Cubic_I] )
101     {
102     JudgeSymmetry[(Int4)Tetragonal_I] = true;
103     }
104     if( JudgeSymmetry[(Int4)Tetragonal_I] )
105     {
106     JudgeSymmetry[(Int4)Orthorhombic_I] = true;
107     }
108    
109     if( JudgeSymmetry[(Int4)Cubic_F] )
110     {
111     JudgeSymmetry[(Int4)Orthorhombic_F] = true;
112     }
113    
114     m_resol2 = cont.putResolution() * 2.0;
115     m_num_ref_figure_of_merit = cont.putNumberOfReflectionsForFigureOfMerit();
116     m_etype_peak_shift = cont.putPeakShiftFunctionType();
117     m_WlengthX = cont.putWaveLength();
118    
119     const vector<Double>& peak_shift_param_rad = cont.putPeakShiftParamRadian();
120     const Int4 param_num = peak_shift_param_rad.size();
121     assert( m_etype_peak_shift != kPeakShiftFunction_Type0 || param_num == 0 );
122     assert( m_etype_peak_shift != kPeakShiftFunction_Type1 || param_num == 1 );
123    
124     m_peak_shift_param_rad.resize(param_num);
125     for(Int4 i=0; i<param_num; i++)
126     {
127     m_peak_shift_param_rad[i] = peak_shift_param_rad[i];
128     }
129     }
130    
131    
132     void SortingLattice::putLatticeFigureOfMerit(const LatticeFigureOfMeritToCheckSymmetry& lattice_original,
133     const ePointGroup& epg, const Double& cv2,
134     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
135     {
136     lattice_result.clear();
137     map< SymMat<VCData>, NRMat<Int4> > S_red_tray;
138     if( !lattice_original.checkLatticeSymmetry(epg, cv2, S_red_tray) ) return;
139    
140     const BravaisType& ebrat_original = lattice_original.putLatticeFigureOfMerit().putBravaisType();
141     const eBravaisLattice eblat = (ebrat_original.enumCrystalSystem()==Monoclinic_B?
142     (epg==D31d_rho?Prim:(epg==D3d_1_hex?Rhom_hex:BaseZ)):ebrat_original.enumBravaisLattice());
143    
144     const NRMat<Int4> matrix_min_to_sell = lattice_original.putInitialForm().second;
145    
146     SymMat<Double> S_super(4);
147     Int4 itnum;
148     NRMat<Int4> trans_mat(4,3);
149    
150     for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
151     {
152     // S_super = it->second * it->first * Transpose(it->second) is close to
153     // Delone-reduced form of the original lattice.
154     S_super = transform_sym_matrix(it->second, chToDouble(it->first) );
155    
156     trans_mat = identity_matrix<Int4>(4);
157    
158     // S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second).
159     put_super_Gram_matrix_obtuse_angle< Double, SymMat<Double> >(trans_mat, S_super, itnum);
160     moveSmallerDiagonalLeftUpper< Double, SymMat<Double> >(S_super, trans_mat);
161    
162     lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( BravaisType( pair<eBravaisLattice, ePointGroup>(eblat, epg) ),
163     SymMat43_VCData(it->first, mprod(trans_mat, it->second) ),
164     lattice_original.putLatticeFigureOfMerit().putPeakShiftFunctionType(),
165     lattice_original.putLatticeFigureOfMerit().putWaveLength(),
166     lattice_original.putLatticeFigureOfMerit().putPeakShiftParamRadian() ) );
167     }
168     }
169    
170    
171     void SortingLattice::putBravaisLatticeFigureOfMerit(const ReducedLatticeToCheckBravais& RLCB,
172     const LatticeFigureOfMeritToCheckSymmetry& lattice_original,
173     const BravaisType& brat,
174     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
175     {
176     lattice_result.clear();
177    
178     const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(brat);
179     if( S_red_tray.empty() ) return;
180    
181     // The lattice of RLCB has at least the symmetry given by eblat.
182     SymMat<VCData> S_super_obtuse(4);
183     Int4 itnum;
184     NRMat<Int4> trans_mat(4,3);
185    
186     for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
187     {
188     S_super_obtuse = transform_sym_matrix(it->second, it->first);
189     trans_mat = identity_matrix<Int4>(4);
190    
191     // S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second) is Delone reduced.
192     if( !put_super_Gram_matrix_obtuse_angle< VCData, SymMat<VCData> >(trans_mat, S_super_obtuse, itnum) )
193     {
194     assert( false );
195     }
196     moveSmallerDiagonalLeftUpper< VCData, SymMat<VCData> >(S_super_obtuse, trans_mat);
197     lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( brat, SymMat43_VCData(it->first, mprod(trans_mat, it->second) ),
198     lattice_original.putLatticeFigureOfMerit().putPeakShiftFunctionType(),
199     lattice_original.putLatticeFigureOfMerit().putWaveLength(),
200     lattice_original.putLatticeFigureOfMerit().putPeakShiftParamRadian() ) );
201     }
202     }
203    
204    
205    
206    
207    
208    
209     void SortingLattice::putLatticeCandidatesForTriclinic(const vector<SymMat43_VCData>& S_super,
210     const Double& MIN_NormM,
211     const Double& MIN_RevM,
212     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri) const
213     {
214     const Int4 num_topo = S_super.size();
215     lattice_result_tri.clear();
216    
217     /* 2011.10.19 VIC Tamura */
218     Int4 LOOP_COUNTER = 0;
219    
220     #ifdef _OPENMP
221     #pragma omp parallel
222     #endif
223     {
224     vector< VecDat3<Int4> > closest_hkl_tray;
225     vector<bool> is_cal_Q_observed_tray;
226     vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
227    
228     #ifdef _OPENMP
229     #pragma omp for
230     #endif
231     for(Int4 n=0; n<num_topo; n++)
232     {
233     /* 2011.10.19 VIC Tamura */
234     SET_PROGRESS(100*(LOOP_COUNTER++)/num_topo, 65, 1); // critical, but works
235     if(IS_CANSELED()) continue;
236    
237     LatticeFigureOfMeritToCheckSymmetry latFOM(BravaisType( pair<eBravaisLattice, ePointGroup>(Prim, Ci) ), S_super[n],
238     m_etype_peak_shift, m_WlengthX, m_peak_shift_param_rad);
239    
240     latFOM.setFigureOfMerit(m_num_ref_figure_of_merit,
241     VCData::putPeakQData(),
242     closest_hkl_tray, is_cal_Q_observed_tray);
243    
244     LatticeFigureOfMeritZeroShift latFOM2 = latFOM.putLatticeFigureOfMerit();
245     pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
246     closest_hkl_tray, is_cal_Q_observed_tray, false);
247    
248     if( ans.first )
249     {
250     assert( latFOM.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
251     latFOM2.setFigureOfMerit(latFOM.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
252     VCData::putPeakQData(),
253     closest_hkl_tray, is_cal_Q_observed_tray);
254 rtomiyasu 16 if( LatticeFigureOfMerit::cmpFOMWu( latFOM2, latFOM.putLatticeFigureOfMerit() ) )
255 rtomiyasu 3 {
256     latFOM.setLatticeFigureOfMerit(latFOM2);
257     }
258     }
259     const SetOfFigureOfMerit& setFOM = latFOM.putLatticeFigureOfMerit().putFiguresOfMerit();
260     if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
261     if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;
262    
263     #ifdef _OPENMP
264     #pragma omp critical
265     #endif
266     {
267     lattice_result_tri.push_back( latFOM );
268     }
269     }
270     }
271    
272     /* 2011.10.19 VIC Tamura */
273     CHECK_INTERRUPTION();
274    
275 rtomiyasu 16 // sort( lattice_result_tri.begin(), lattice_result_tri.end() );
276 rtomiyasu 3 }
277    
278    
279     void SortingLattice::putLatticeCandidatesForEachBravaisTypes(
280     const Double& MIN_NormM,
281     const Double& MIN_RevM,
282 rtomiyasu 16 const size_t& MAX_SIZE,
283 rtomiyasu 3 const eABCaxis& abc_axis,
284     const eRHaxis& rh_axis,
285     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
286     {
287     try{
288    
289 rtomiyasu 16 for(ArrayIndex i=1; i<NUM_LS; i++)
290 rtomiyasu 3 {
291     lattice_result[i].clear();
292     }
293     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(ArrayIndex)Triclinic];
294     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_P = lattice_result[(ArrayIndex)Monoclinic_P];
295     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_B = lattice_result[(ArrayIndex)Monoclinic_B];
296     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_P = lattice_result[(ArrayIndex)Orthorhombic_P];
297     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_B = lattice_result[(ArrayIndex)Orthorhombic_C];
298     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_I = lattice_result[(ArrayIndex)Orthorhombic_I];
299     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_F = lattice_result[(ArrayIndex)Orthorhombic_F];
300     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_P = lattice_result[(ArrayIndex)Tetragonal_P];
301     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_I = lattice_result[(ArrayIndex)Tetragonal_I];
302     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_rhom = lattice_result[(ArrayIndex)Rhombohedral];
303     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_hex = lattice_result[(ArrayIndex)Hexagonal];
304     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_P = lattice_result[(ArrayIndex)Cubic_P];
305     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_I = lattice_result[(ArrayIndex)Cubic_I];
306     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_F = lattice_result[(ArrayIndex)Cubic_F];
307    
308     const Int4 num_tri = lattice_result_tri.size();
309    
310     /* 2011.10.19 VIC Tamura */
311     Int4 LOOP_COUNTER = 0;
312    
313     #ifdef _OPENMP
314     #pragma omp parallel
315     #endif
316     {
317     vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
318    
319     #ifdef _OPENMP
320     #pragma omp for
321     #endif
322     for(Int4 n=0; n<num_tri; n++)
323     {
324     /* 2011.10.19 VIC Tamura */
325     SET_PROGRESS(99*(LOOP_COUNTER++)/num_tri, 66, 30); // critical, but works
326     if(IS_CANSELED()) continue;
327    
328     LatticeFigureOfMeritToCheckSymmetry& latFOM = lattice_result_tri[n];
329     const SymMat43_VCData& S_red = latFOM.putInitialForm();
330    
331     // Calculate figures of merit as triclinic
332     const ReducedLatticeToCheckBravais RLCB(abc_axis, rh_axis, m_DoesPrudentSymSearch, m_cv2, S_red);
333    
334     if( JudgeSymmetry[Monoclinic_B] )
335     {
336     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_monoclinic_b_type(abc_axis) ), latFOM_tray);
337     #ifdef _OPENMP
338     #pragma omp critical(monoB)
339     #endif
340     {
341     lattice_result_mono_B.insert(lattice_result_mono_B.end(), latFOM_tray.begin(), latFOM_tray.end());
342     }
343     }
344     if( JudgeSymmetry[Orthorhombic_I] )
345     {
346     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Inner, D2h) ), latFOM_tray);
347     #ifdef _OPENMP
348     #pragma omp critical(orthoI)
349     #endif
350     {
351     lattice_result_ortho_I.insert(lattice_result_ortho_I.end(), latFOM_tray.begin(), latFOM_tray.end());
352     }
353     }
354     if( JudgeSymmetry[Orthorhombic_F] )
355     {
356     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Face, D2h) ), latFOM_tray);
357     #ifdef _OPENMP
358     #pragma omp critical(orthoF)
359     #endif
360     {
361     lattice_result_ortho_F.insert(lattice_result_ortho_F.end(), latFOM_tray.begin(), latFOM_tray.end());
362     }
363     }
364     if( JudgeSymmetry[Rhombohedral] )
365     {
366     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_rhombohedral_type(rh_axis) ), latFOM_tray);
367     #ifdef _OPENMP
368     #pragma omp critical(rhom)
369     #endif
370     {
371     lattice_result_rhom.insert(lattice_result_rhom.end(), latFOM_tray.begin(), latFOM_tray.end());
372     }
373     }
374    
375     if( JudgeSymmetry[Monoclinic_P] )
376     {
377     putLatticeFigureOfMerit(latFOM, put_monoclinic_p_type(abc_axis), m_cv2, latFOM_tray);
378     #ifdef _OPENMP
379     #pragma omp critical(monoP)
380     #endif
381     {
382     lattice_result_mono_P.insert(lattice_result_mono_P.end(), latFOM_tray.begin(), latFOM_tray.end());
383     }
384     }
385     if( JudgeSymmetry[Orthorhombic_P] )
386     {
387     putLatticeFigureOfMerit(latFOM, D2h, m_cv2, latFOM_tray);
388     #ifdef _OPENMP
389     #pragma omp critical (ortho_P)
390     #endif
391     {
392     lattice_result_ortho_P.insert(lattice_result_ortho_P.end(), latFOM_tray.begin(), latFOM_tray.end());
393     }
394     }
395     }
396     }
397    
398 rtomiyasu 16 sort( lattice_result_tri.begin(), lattice_result_tri.end(), LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff );
399     if( MAX_SIZE < lattice_result_tri.size() )
400     {
401     lattice_result_tri.erase( lattice_result_tri.begin() + MAX_SIZE, lattice_result_tri.end() );
402     }
403    
404     const size_t num_tri2 = lattice_result_tri.size();
405     for(size_t n=0; n<num_tri2; n++)
406     {
407     lattice_result_tri[n].setLabel(n+1);
408     }
409    
410     ZLOG_INFO( "The program has selected " + num2str<Int4>( lattice_result_tri.size() )
411     + " triclinic solutions using the Wu figure of merit.\n\n" );
412    
413    
414 rtomiyasu 3 /* 2011.10.19 VIC Tamura */
415     CHECK_INTERRUPTION();
416    
417     /* 2011.10.19 VIC Tamura */
418     LOOP_COUNTER = 0;
419     Int4 SUM = 0;
420     for(Int4 i=0; i<NUM_LS; i++) { SUM += lattice_result[i].size(); }
421    
422     for(ArrayIndex i=1; i<NUM_LS; i++)
423     {
424     if( !JudgeSymmetry[i] ) continue;
425 rtomiyasu 16 // sort( lattice_result[i].begin(), lattice_result[i].end() );
426 rtomiyasu 3
427     const Int4 num_lattice = lattice_result[i].size();
428    
429     #ifdef _OPENMP
430     #pragma omp parallel
431     #endif
432     {
433     vector< VecDat3<Int4> > closest_hkl_tray;
434     vector<bool> is_cal_Q_observed_tray;
435     vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
436    
437     #ifdef _OPENMP
438     #pragma omp for
439     #endif
440     for(Int4 index=0; index<num_lattice; index++)
441     {
442     /* 2011.10.19 VIC Tamura */
443     SET_PROGRESS(99+1*(LOOP_COUNTER++)/SUM, 66, 30); // critical, but works
444     if(IS_CANSELED()) continue;
445    
446     LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
447     latFOM0.setFigureOfMerit(m_num_ref_figure_of_merit,
448     VCData::putPeakQData(),
449     closest_hkl_tray, is_cal_Q_observed_tray);
450    
451     LatticeFigureOfMeritZeroShift latFOM2 = latFOM0.putLatticeFigureOfMerit();
452     pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
453     closest_hkl_tray, is_cal_Q_observed_tray, false);
454     if( ans.first )
455     {
456     assert( latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
457     latFOM2.setFigureOfMerit(latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
458     VCData::putPeakQData(),
459     closest_hkl_tray, is_cal_Q_observed_tray);
460 rtomiyasu 16 if( LatticeFigureOfMerit::cmpFOMWu( latFOM2, latFOM0.putLatticeFigureOfMerit() ) )
461 rtomiyasu 3 {
462     latFOM0.setLatticeFigureOfMerit(latFOM2);
463     }
464     }
465    
466     const SetOfFigureOfMerit& setFOM = latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit();
467     if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
468     if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;
469    
470     if( eCrystalSystem(i) == Monoclinic_P )
471     {
472     if( JudgeSymmetry[Hexagonal] )
473     {
474     putLatticeFigureOfMerit(latFOM0, D6h, m_cv2, latFOM_tray);
475     #ifdef _OPENMP
476     #pragma omp critical (hex)
477     #endif
478     {
479     lattice_result_hex.insert(lattice_result_hex.end(), latFOM_tray.begin(), latFOM_tray.end());
480     }
481     }
482     }
483     else if( eCrystalSystem(i) == Monoclinic_B )
484     {
485     if( JudgeSymmetry[Orthorhombic_C] )
486     {
487     putLatticeFigureOfMerit(latFOM0, D2h, m_cv2, latFOM_tray);
488     #ifdef _OPENMP
489     #pragma omp critical (ortho_B)
490     #endif
491     {
492     lattice_result_ortho_B.insert(lattice_result_ortho_B.end(), latFOM_tray.begin(), latFOM_tray.end());
493     }
494     }
495     }
496     else if( eCrystalSystem(i) == Orthorhombic_P )
497     {
498     if( JudgeSymmetry[Tetragonal_P] )
499     {
500     putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
501     #ifdef _OPENMP
502     #pragma omp critical (tetra_P)
503     #endif
504     {
505     lattice_result_tetra_P.insert(lattice_result_tetra_P.end(), latFOM_tray.begin(), latFOM_tray.end());
506     }
507     }
508     if( JudgeSymmetry[Cubic_P] )
509     {
510     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
511     #ifdef _OPENMP
512     #pragma omp critical (cubic_P)
513     #endif
514     {
515     lattice_result_cubic_P.insert(lattice_result_cubic_P.end(), latFOM_tray.begin(), latFOM_tray.end());
516     }
517     }
518     }
519     else if( eCrystalSystem(i) == Orthorhombic_I )
520     {
521     if( JudgeSymmetry[Tetragonal_I] )
522     {
523     putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
524     #ifdef _OPENMP
525     #pragma omp critical (tetra_I)
526     #endif
527     {
528     lattice_result_tetra_I.insert(lattice_result_tetra_I.end(), latFOM_tray.begin(), latFOM_tray.end());
529     }
530     }
531     if( JudgeSymmetry[Cubic_I] )
532     {
533     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
534     #ifdef _OPENMP
535     #pragma omp critical (cubic_I)
536     #endif
537     {
538     lattice_result_cubic_I.insert(lattice_result_cubic_I.end(), latFOM_tray.begin(), latFOM_tray.end());
539     }
540     }
541     }
542     else if( eCrystalSystem(i) == Orthorhombic_F )
543     {
544     if( JudgeSymmetry[Cubic_F] )
545     {
546     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
547     #ifdef _OPENMP
548     #pragma omp critical (cubic_F)
549     #endif
550     {
551     lattice_result_cubic_F.insert(lattice_result_cubic_F.end(), latFOM_tray.begin(), latFOM_tray.end());
552     }
553     }
554     }
555     }
556     }
557     /* 2011.10.19 VIC Tamura */
558     CHECK_INTERRUPTION();
559    
560 rtomiyasu 16 sort( lattice_result[i].begin(), lattice_result[i].end(), LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff );
561     if( MAX_SIZE < lattice_result[i].size() )
562     {
563     lattice_result[i].erase( lattice_result[i].begin() + MAX_SIZE, lattice_result[i].end() );
564     }
565    
566     const size_t num_lattice2 = lattice_result[i].size();
567     for(size_t n=0; n<num_lattice2; n++)
568     {
569     lattice_result[i][n].setLabel(n+1);
570     }
571    
572 rtomiyasu 3 ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), abc_axis)
573     + " : " + num2str<Int4>( lattice_result[i].size() ) + "\n" );
574     }
575     ZLOG_INFO( "\n" );
576     }
577     catch(bad_alloc& ball)
578     {
579     throw nerror(ball, __FILE__, __LINE__, __FUNCTION__);
580     }
581     }
582    
583    
584     void SortingLattice::putLatticeCandidatesForEachBravaisTypes(const vector<SymMat43_VCData>& S_super,
585     const Double& MIN_NormM,
586     const Double& MIN_RevM,
587 rtomiyasu 16 const size_t& MAX_SIZE,
588 rtomiyasu 3 const eABCaxis& abc_axis,
589     const eRHaxis& rh_axis,
590     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
591     {
592     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(Int4)Triclinic];
593     putLatticeCandidatesForTriclinic(S_super, MIN_NormM, MIN_RevM, lattice_result_tri);
594    
595     ZLOG_INFO( "Determination of the Bravais type is being carried out...\n(Solutions with " + putLabel(SCWuM) + " less than " + num2str(MIN_NormM)
596     + " or " + putLabel(SCRevM) + " less than " + num2str(MIN_RevM)
597 rtomiyasu 16 + " are automatically removed).\n" );
598     ZLOG_INFO( "All the unit-cell parameters are being optimized by linear least squares...\n" );
599 rtomiyasu 3
600 rtomiyasu 16 //ZLOG_INFO( "The program has removed " + num2str<Int4>( S_super.size() - lattice_result_tri.size() )
601     // + " triclinic solutions because their " + putLabel(SCWuM) + " is less than " + num2str(MIN_NormM)
602     // + " or their " + putLabel(SCRevM) + " is less than " + num2str(MIN_RevM) + ".\n\n" );
603     //ZLOG_INFO( "Determination of the Bravais type is being carried out...\n" );
604     putLatticeCandidatesForEachBravaisTypes(MIN_NormM, MIN_RevM, MAX_SIZE, abc_axis, rh_axis, lattice_result);
605 rtomiyasu 3 }
606    
607    
608     void SortingLattice::setNumberOfNeighbors(const eABCaxis& baxis_flag,
609 rtomiyasu 16 bool (*CmpFunc)(const LatticeFigureOfMeritToCheckSymmetry&, const LatticeFigureOfMeritToCheckSymmetry&),
610     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
611 rtomiyasu 3 {
612    
613     #ifdef _OPENMP
614     #pragma omp for
615     #endif
616     for(ArrayIndex i=0; i<NUM_LS; i++)
617     {
618     if( !OutputSymmetry[i] ) continue;
619    
620     stable_sort( lattice_result[i].begin(), lattice_result[i].end() ); // Sort by the unit-cell volume.
621     for(vector<LatticeFigureOfMeritToCheckSymmetry>::iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
622     {
623     it->putNumberOfLatticesInNeighborhood() = 0;
624     }
625     }
626    
627     const Double coef_lower = 1.0 - sqrt(3.0)*m_resol2;
628     const Double coef_upper = 1.0 + sqrt(3.0)*m_resol2;
629    
630     Vec_INT index_tray(put_cs_num());
631    
632     /* 2011.10.19 VIC Tamura */
633     Int4 SUM=0, LOOP_COUNTER=0;
634     for(Int4 i=0; i<NUM_LS; i++ ) { SUM += lattice_result[i].size(); }
635    
636     #ifdef _OPENMP
637     #pragma omp for
638     #endif
639     for(ArrayIndex i=0; i<NUM_LS; i++)
640     {
641     if( !OutputSymmetry[i] ) continue;
642    
643     const Int4 num_lattice = lattice_result[i].size();
644    
645     for(Int4 index=0; index<num_lattice; index++)
646     {
647     /* 2011.10.19 VIC Tamura */
648     SET_PROGRESS(100*(LOOP_COUNTER++)/SUM, 97, 1); // critical, but works
649     if(IS_CANSELED()) continue;
650    
651     LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
652     const LatticeFigureOfMerit& latFOM0_prim = latFOM0.putLatticeFigureOfMerit();
653     if( latFOM0.putNumberOfLatticesInNeighborhood() < 0 ) continue;
654    
655     const Double& detS = latFOM0_prim.putDeterminantOfGramMatrix();
656     const Int4 ibegin = distance( lattice_result[i].begin(), lower_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_lower ) );
657     const Int4 iend = distance( lattice_result[i].begin(), upper_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_upper ) );
658    
659     Int4 count=0;
660     if( i == (Int4)Triclinic )
661     {
662     const ReducedLatticeToCheckEquiv RLCS(m_resol2, latFOM0_prim.putOptimizedForm());
663     for(Int4 index2=ibegin; index2<iend; index2++)
664     {
665     if( index2 == index ) continue;
666    
667     LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
668     const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();
669    
670     // lattice_result_tri[index2] equals trans_mat * RLCB.m_S_super_obtuse * Transpose(trans_mat)
671     if( RLCS.equiv( latFOM2_prim.putSellingReducedForm() ) )
672     {
673     // Compare the figures of merit.
674 rtomiyasu 16 if( CmpFunc( latFOM2, latFOM0 ) )
675 rtomiyasu 3 {
676 rtomiyasu 16 if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0 )
677     {
678     count = -1;
679     break;
680     }
681 rtomiyasu 3 }
682     else
683     {
684     count++;
685     latFOM2.putNumberOfLatticesInNeighborhood() = -1;
686     }
687     }
688     }
689     }
690     else
691     {
692     for(Int4 index2=ibegin; index2<iend; index2++)
693     {
694     if( index2 == index ) continue;
695    
696     LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
697     const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();
698    
699     // *it2 equals trans_mat * RLCS.m_S_super_obtuse * Transpose(trans_mat)
700     if( check_equiv_m( latFOM0_prim.putInverseOfMinkowskiReducedForm(), latFOM2_prim.putInverseOfMinkowskiReducedForm(), m_resol2 ) )
701     {
702     // Compare the figures of merit.
703 rtomiyasu 16 if( CmpFunc( latFOM2, latFOM0 ) )
704 rtomiyasu 3 {
705 rtomiyasu 16 if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0 )
706     {
707     count = -1;
708     break;
709     }
710 rtomiyasu 3 }
711     else
712     {
713     count++;
714     latFOM2.putNumberOfLatticesInNeighborhood() = -1;
715     }
716     }
717     }
718     }
719    
720     latFOM0.putNumberOfLatticesInNeighborhood() = count;
721     }
722    
723     Int4& index = index_tray[i];
724 rtomiyasu 16 index = 0;
725 rtomiyasu 3 for(vector<LatticeFigureOfMeritToCheckSymmetry>::const_iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
726     {
727     if( it->putNumberOfLatticesInNeighborhood() >= 0 ) index++;
728     }
729     }
730    
731     /* 2011.10.19 VIC Tamura */
732     CHECK_INTERRUPTION();
733    
734     for(ArrayIndex i=0; i<NUM_LS; i++)
735     {
736     if( !OutputSymmetry[i] ) continue;
737     ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), baxis_flag)
738     + " : " + num2str( lattice_result[i].size() ) + " -> " + num2str( index_tray[i] ) + "\n" );
739     }
740     ZLOG_INFO( "\n" );
741     }
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help