Develop and Download Open Source Software

Browse Subversion Repository

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations) (download) (as text)
Fri Feb 22 04:51:31 2013 UTC (11 years, 1 month ago) by rtomiyasu
File MIME type: text/x-c++src
File size: 24988 byte(s)


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     if( LatticeFigureOfMerit::cmpFOMdeWolff( latFOM2, latFOM.putLatticeFigureOfMerit() ) )
255     {
256     latFOM.setLatticeFigureOfMerit(latFOM2);
257     }
258     }
259    
260     const SetOfFigureOfMerit& setFOM = latFOM.putLatticeFigureOfMerit().putFiguresOfMerit();
261    
262     if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
263     if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;
264    
265     #ifdef _OPENMP
266     #pragma omp critical
267     #endif
268     {
269     lattice_result_tri.push_back( latFOM );
270     }
271     }
272     }
273    
274     /* 2011.10.19 VIC Tamura */
275     CHECK_INTERRUPTION();
276    
277     sort( lattice_result_tri.begin(), lattice_result_tri.end() );
278     }
279    
280    
281     void SortingLattice::putLatticeCandidatesForEachBravaisTypes(
282     const Double& MIN_NormM,
283     const Double& MIN_RevM,
284     const eABCaxis& abc_axis,
285     const eRHaxis& rh_axis,
286     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
287     {
288     try{
289    
290     for(ArrayIndex i=1; i<NUM_LS; i++)
291     {
292     lattice_result[i].clear();
293     }
294     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(ArrayIndex)Triclinic];
295     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_P = lattice_result[(ArrayIndex)Monoclinic_P];
296     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_B = lattice_result[(ArrayIndex)Monoclinic_B];
297     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_P = lattice_result[(ArrayIndex)Orthorhombic_P];
298     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_B = lattice_result[(ArrayIndex)Orthorhombic_C];
299     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_I = lattice_result[(ArrayIndex)Orthorhombic_I];
300     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_F = lattice_result[(ArrayIndex)Orthorhombic_F];
301     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_P = lattice_result[(ArrayIndex)Tetragonal_P];
302     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_I = lattice_result[(ArrayIndex)Tetragonal_I];
303     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_rhom = lattice_result[(ArrayIndex)Rhombohedral];
304     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_hex = lattice_result[(ArrayIndex)Hexagonal];
305     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_P = lattice_result[(ArrayIndex)Cubic_P];
306     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_I = lattice_result[(ArrayIndex)Cubic_I];
307     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_F = lattice_result[(ArrayIndex)Cubic_F];
308    
309     const Int4 num_tri = lattice_result_tri.size();
310    
311     /* 2011.10.19 VIC Tamura */
312     Int4 LOOP_COUNTER = 0;
313    
314     #ifdef _OPENMP
315     #pragma omp parallel
316     #endif
317     {
318     vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
319    
320     #ifdef _OPENMP
321     #pragma omp for
322     #endif
323     for(Int4 n=0; n<num_tri; n++)
324     {
325     /* 2011.10.19 VIC Tamura */
326     SET_PROGRESS(99*(LOOP_COUNTER++)/num_tri, 66, 30); // critical, but works
327     if(IS_CANSELED()) continue;
328    
329     LatticeFigureOfMeritToCheckSymmetry& latFOM = lattice_result_tri[n];
330     latFOM.setLabel(n+1);
331     const SymMat43_VCData& S_red = latFOM.putInitialForm();
332    
333     // Calculate figures of merit as triclinic
334     const ReducedLatticeToCheckBravais RLCB(abc_axis, rh_axis, m_DoesPrudentSymSearch, m_cv2, S_red);
335    
336     if( JudgeSymmetry[Monoclinic_B] )
337     {
338     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_monoclinic_b_type(abc_axis) ), latFOM_tray);
339     #ifdef _OPENMP
340     #pragma omp critical(monoB)
341     #endif
342     {
343     lattice_result_mono_B.insert(lattice_result_mono_B.end(), latFOM_tray.begin(), latFOM_tray.end());
344     }
345     }
346     if( JudgeSymmetry[Orthorhombic_I] )
347     {
348     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Inner, D2h) ), latFOM_tray);
349     #ifdef _OPENMP
350     #pragma omp critical(orthoI)
351     #endif
352     {
353     lattice_result_ortho_I.insert(lattice_result_ortho_I.end(), latFOM_tray.begin(), latFOM_tray.end());
354     }
355     }
356     if( JudgeSymmetry[Orthorhombic_F] )
357     {
358     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Face, D2h) ), latFOM_tray);
359     #ifdef _OPENMP
360     #pragma omp critical(orthoF)
361     #endif
362     {
363     lattice_result_ortho_F.insert(lattice_result_ortho_F.end(), latFOM_tray.begin(), latFOM_tray.end());
364     }
365     }
366     if( JudgeSymmetry[Rhombohedral] )
367     {
368     putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_rhombohedral_type(rh_axis) ), latFOM_tray);
369     #ifdef _OPENMP
370     #pragma omp critical(rhom)
371     #endif
372     {
373     lattice_result_rhom.insert(lattice_result_rhom.end(), latFOM_tray.begin(), latFOM_tray.end());
374     }
375     }
376    
377     if( JudgeSymmetry[Monoclinic_P] )
378     {
379     putLatticeFigureOfMerit(latFOM, put_monoclinic_p_type(abc_axis), m_cv2, latFOM_tray);
380     #ifdef _OPENMP
381     #pragma omp critical(monoP)
382     #endif
383     {
384     lattice_result_mono_P.insert(lattice_result_mono_P.end(), latFOM_tray.begin(), latFOM_tray.end());
385     }
386     }
387     if( JudgeSymmetry[Orthorhombic_P] )
388     {
389     putLatticeFigureOfMerit(latFOM, D2h, m_cv2, latFOM_tray);
390     #ifdef _OPENMP
391     #pragma omp critical (ortho_P)
392     #endif
393     {
394     lattice_result_ortho_P.insert(lattice_result_ortho_P.end(), latFOM_tray.begin(), latFOM_tray.end());
395     }
396     }
397     }
398     }
399    
400     /* 2011.10.19 VIC Tamura */
401     CHECK_INTERRUPTION();
402    
403     ZLOG_INFO( "All the lattice constants are being optimized by linear least squares...\n" );
404    
405     /* 2011.10.19 VIC Tamura */
406     LOOP_COUNTER = 0;
407     Int4 SUM = 0;
408     for(Int4 i=0; i<NUM_LS; i++) { SUM += lattice_result[i].size(); }
409    
410     for(ArrayIndex i=1; i<NUM_LS; i++)
411     {
412     if( !JudgeSymmetry[i] ) continue;
413     sort( lattice_result[i].begin(), lattice_result[i].end() );
414    
415     const Int4 num_lattice = lattice_result[i].size();
416    
417     #ifdef _OPENMP
418     #pragma omp parallel
419     #endif
420     {
421     vector< VecDat3<Int4> > closest_hkl_tray;
422     vector<bool> is_cal_Q_observed_tray;
423     vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
424    
425     #ifdef _OPENMP
426     #pragma omp for
427     #endif
428     for(Int4 index=0; index<num_lattice; index++)
429     {
430     /* 2011.10.19 VIC Tamura */
431     SET_PROGRESS(99+1*(LOOP_COUNTER++)/SUM, 66, 30); // critical, but works
432     if(IS_CANSELED()) continue;
433    
434     LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
435     latFOM0.setLabel(index+1);
436    
437     latFOM0.setFigureOfMerit(m_num_ref_figure_of_merit,
438     VCData::putPeakQData(),
439     closest_hkl_tray, is_cal_Q_observed_tray);
440    
441     LatticeFigureOfMeritZeroShift latFOM2 = latFOM0.putLatticeFigureOfMerit();
442     pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
443     closest_hkl_tray, is_cal_Q_observed_tray, false);
444     if( ans.first )
445     {
446     assert( latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
447     latFOM2.setFigureOfMerit(latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
448     VCData::putPeakQData(),
449     closest_hkl_tray, is_cal_Q_observed_tray);
450     if( LatticeFigureOfMerit::cmpFOMdeWolff( latFOM2, latFOM0.putLatticeFigureOfMerit() ) )
451     {
452     latFOM0.setLatticeFigureOfMerit(latFOM2);
453     }
454     }
455    
456     const SetOfFigureOfMerit& setFOM = latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit();
457     if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
458     if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;
459    
460     if( eCrystalSystem(i) == Monoclinic_P )
461     {
462     if( JudgeSymmetry[Hexagonal] )
463     {
464     putLatticeFigureOfMerit(latFOM0, D6h, m_cv2, latFOM_tray);
465     #ifdef _OPENMP
466     #pragma omp critical (hex)
467     #endif
468     {
469     lattice_result_hex.insert(lattice_result_hex.end(), latFOM_tray.begin(), latFOM_tray.end());
470     }
471     }
472     }
473     else if( eCrystalSystem(i) == Monoclinic_B )
474     {
475     if( JudgeSymmetry[Orthorhombic_C] )
476     {
477     putLatticeFigureOfMerit(latFOM0, D2h, m_cv2, latFOM_tray);
478     #ifdef _OPENMP
479     #pragma omp critical (ortho_B)
480     #endif
481     {
482     lattice_result_ortho_B.insert(lattice_result_ortho_B.end(), latFOM_tray.begin(), latFOM_tray.end());
483     }
484     }
485     }
486     else if( eCrystalSystem(i) == Orthorhombic_P )
487     {
488     if( JudgeSymmetry[Tetragonal_P] )
489     {
490     putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
491     #ifdef _OPENMP
492     #pragma omp critical (tetra_P)
493     #endif
494     {
495     lattice_result_tetra_P.insert(lattice_result_tetra_P.end(), latFOM_tray.begin(), latFOM_tray.end());
496     }
497     }
498     if( JudgeSymmetry[Cubic_P] )
499     {
500     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
501     #ifdef _OPENMP
502     #pragma omp critical (cubic_P)
503     #endif
504     {
505     lattice_result_cubic_P.insert(lattice_result_cubic_P.end(), latFOM_tray.begin(), latFOM_tray.end());
506     }
507     }
508     }
509     else if( eCrystalSystem(i) == Orthorhombic_I )
510     {
511     if( JudgeSymmetry[Tetragonal_I] )
512     {
513     putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
514     #ifdef _OPENMP
515     #pragma omp critical (tetra_I)
516     #endif
517     {
518     lattice_result_tetra_I.insert(lattice_result_tetra_I.end(), latFOM_tray.begin(), latFOM_tray.end());
519     }
520     }
521     if( JudgeSymmetry[Cubic_I] )
522     {
523     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
524     #ifdef _OPENMP
525     #pragma omp critical (cubic_I)
526     #endif
527     {
528     lattice_result_cubic_I.insert(lattice_result_cubic_I.end(), latFOM_tray.begin(), latFOM_tray.end());
529     }
530     }
531     }
532     else if( eCrystalSystem(i) == Orthorhombic_F )
533     {
534     if( JudgeSymmetry[Cubic_F] )
535     {
536     putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
537     #ifdef _OPENMP
538     #pragma omp critical (cubic_F)
539     #endif
540     {
541     lattice_result_cubic_F.insert(lattice_result_cubic_F.end(), latFOM_tray.begin(), latFOM_tray.end());
542     }
543     }
544     }
545     }
546     }
547     /* 2011.10.19 VIC Tamura */
548     CHECK_INTERRUPTION();
549    
550     ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), abc_axis)
551     + " : " + num2str<Int4>( lattice_result[i].size() ) + "\n" );
552     }
553     ZLOG_INFO( "\n" );
554     }
555     catch(bad_alloc& ball)
556     {
557     throw nerror(ball, __FILE__, __LINE__, __FUNCTION__);
558     }
559     }
560    
561    
562     void SortingLattice::putLatticeCandidatesForEachBravaisTypes(const vector<SymMat43_VCData>& S_super,
563     const Double& MIN_NormM,
564     const Double& MIN_RevM,
565     const eABCaxis& abc_axis,
566     const eRHaxis& rh_axis,
567     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
568     {
569     vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(Int4)Triclinic];
570     putLatticeCandidatesForTriclinic(S_super, MIN_NormM, MIN_RevM, lattice_result_tri);
571    
572     ZLOG_INFO( "Determination of the Bravais type is being carried out...\n(Solutions with " + putLabel(SCWuM) + " less than " + num2str(MIN_NormM)
573     + " or " + putLabel(SCRevM) + " less than " + num2str(MIN_RevM)
574     + " are automatically removed).\n\n" );
575    
576     ZLOG_INFO( "The program has removed " + num2str<Int4>( S_super.size() - lattice_result_tri.size() )
577     + " solutions because their " + putLabel(SCWuM) + " is less than " + num2str(MIN_NormM)
578     + " or their " + putLabel(SCRevM) + " is less than " + num2str(MIN_RevM) + ".\n\n" );
579     ZLOG_INFO( "Determination of the Bravais type is being carried out...\n" );
580     putLatticeCandidatesForEachBravaisTypes(MIN_NormM, MIN_RevM, abc_axis, rh_axis, lattice_result);
581     }
582    
583    
584     void SortingLattice::setNumberOfNeighbors(const eABCaxis& baxis_flag,
585     vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
586     {
587    
588     #ifdef _OPENMP
589     #pragma omp for
590     #endif
591     for(ArrayIndex i=0; i<NUM_LS; i++)
592     {
593     if( !OutputSymmetry[i] ) continue;
594    
595     stable_sort( lattice_result[i].begin(), lattice_result[i].end() ); // Sort by the unit-cell volume.
596     for(vector<LatticeFigureOfMeritToCheckSymmetry>::iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
597     {
598     it->putNumberOfLatticesInNeighborhood() = 0;
599     }
600     }
601    
602     const Double coef_lower = 1.0 - sqrt(3.0)*m_resol2;
603     const Double coef_upper = 1.0 + sqrt(3.0)*m_resol2;
604    
605     Vec_INT index_tray(put_cs_num());
606    
607     /* 2011.10.19 VIC Tamura */
608     Int4 SUM=0, LOOP_COUNTER=0;
609     for(Int4 i=0; i<NUM_LS; i++ ) { SUM += lattice_result[i].size(); }
610    
611     #ifdef _OPENMP
612     #pragma omp for
613     #endif
614     for(ArrayIndex i=0; i<NUM_LS; i++)
615     {
616     if( !OutputSymmetry[i] ) continue;
617    
618     const Int4 num_lattice = lattice_result[i].size();
619    
620     for(Int4 index=0; index<num_lattice; index++)
621     {
622     /* 2011.10.19 VIC Tamura */
623     SET_PROGRESS(100*(LOOP_COUNTER++)/SUM, 97, 1); // critical, but works
624     if(IS_CANSELED()) continue;
625    
626     LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
627     const LatticeFigureOfMerit& latFOM0_prim = latFOM0.putLatticeFigureOfMerit();
628     if( latFOM0.putNumberOfLatticesInNeighborhood() < 0 ) continue;
629    
630     const Double& detS = latFOM0_prim.putDeterminantOfGramMatrix();
631     const Int4 ibegin = distance( lattice_result[i].begin(), lower_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_lower ) );
632     const Int4 iend = distance( lattice_result[i].begin(), upper_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_upper ) );
633    
634     Int4 count=0;
635     if( i == (Int4)Triclinic )
636     {
637     const ReducedLatticeToCheckEquiv RLCS(m_resol2, latFOM0_prim.putOptimizedForm());
638     for(Int4 index2=ibegin; index2<iend; index2++)
639     {
640     if( index2 == index ) continue;
641    
642     LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
643     const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();
644    
645     // lattice_result_tri[index2] equals trans_mat * RLCB.m_S_super_obtuse * Transpose(trans_mat)
646     if( RLCS.equiv( latFOM2_prim.putSellingReducedForm() ) )
647     {
648     // Compare the figures of merit.
649     if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
650     && LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
651     {
652     count = -1;
653     break;
654     }
655     else
656     {
657     count++;
658     latFOM2.putNumberOfLatticesInNeighborhood() = -1;
659     }
660     }
661     }
662     }
663     else
664     {
665     for(Int4 index2=ibegin; index2<iend; index2++)
666     {
667     if( index2 == index ) continue;
668    
669     LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
670     const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();
671    
672     // *it2 equals trans_mat * RLCS.m_S_super_obtuse * Transpose(trans_mat)
673     if( check_equiv_m( latFOM0_prim.putInverseOfMinkowskiReducedForm(), latFOM2_prim.putInverseOfMinkowskiReducedForm(), m_resol2 ) )
674     {
675     // Compare the figures of merit.
676     if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
677     && LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
678     {
679     count = -1;
680     break;
681     }
682     else
683     {
684     count++;
685     latFOM2.putNumberOfLatticesInNeighborhood() = -1;
686     }
687     }
688     }
689     }
690    
691     latFOM0.putNumberOfLatticesInNeighborhood() = count;
692     }
693    
694     Int4& index = index_tray[i];
695     for(vector<LatticeFigureOfMeritToCheckSymmetry>::const_iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
696     {
697     if( it->putNumberOfLatticesInNeighborhood() >= 0 ) index++;
698     }
699     }
700    
701     /* 2011.10.19 VIC Tamura */
702     CHECK_INTERRUPTION();
703    
704     for(ArrayIndex i=0; i<NUM_LS; i++)
705     {
706     if( !OutputSymmetry[i] ) continue;
707     ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), baxis_flag)
708     + " : " + num2str( lattice_result[i].size() ) + " -> " + num2str( index_tray[i] ) + "\n" );
709     }
710     ZLOG_INFO( "\n" );
711     }

Back to OSDN">Back to OSDN
ViewVC Help
Powered by ViewVC 1.1.26