Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 13 - (show annotations) (download) (as text)
Sat Mar 2 10:44:08 2013 UTC (11 years ago) by rtomiyasu
File MIME type: text/x-c++src
File size: 24984 byte(s) 
デバッグモードでコンパイルした際に発生するエラーの修正。
1 /*
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 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 sort( lattice_result_tri.begin(), lattice_result_tri.end() );
276 }
277
278
279 void SortingLattice::putLatticeCandidatesForEachBravaisTypes(
280 const Double& MIN_NormM,
281 const Double& MIN_RevM,
282 const eABCaxis& abc_axis,
283 const eRHaxis& rh_axis,
284 vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
285 {
286 try{
287
288 for(ArrayIndex i=1; i<NUM_LS; i++)
289 {
290 lattice_result[i].clear();
291 }
292 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(ArrayIndex)Triclinic];
293 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_P = lattice_result[(ArrayIndex)Monoclinic_P];
294 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_B = lattice_result[(ArrayIndex)Monoclinic_B];
295 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_P = lattice_result[(ArrayIndex)Orthorhombic_P];
296 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_B = lattice_result[(ArrayIndex)Orthorhombic_C];
297 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_I = lattice_result[(ArrayIndex)Orthorhombic_I];
298 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_F = lattice_result[(ArrayIndex)Orthorhombic_F];
299 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_P = lattice_result[(ArrayIndex)Tetragonal_P];
300 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_I = lattice_result[(ArrayIndex)Tetragonal_I];
301 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_rhom = lattice_result[(ArrayIndex)Rhombohedral];
302 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_hex = lattice_result[(ArrayIndex)Hexagonal];
303 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_P = lattice_result[(ArrayIndex)Cubic_P];
304 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_I = lattice_result[(ArrayIndex)Cubic_I];
305 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_F = lattice_result[(ArrayIndex)Cubic_F];
306
307 const Int4 num_tri = lattice_result_tri.size();
308
309 /* 2011.10.19 VIC Tamura */
310 Int4 LOOP_COUNTER = 0;
311
312 #ifdef _OPENMP
313 #pragma omp parallel
314 #endif
315 {
316 vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
317
318 #ifdef _OPENMP
319 #pragma omp for
320 #endif
321 for(Int4 n=0; n<num_tri; n++)
322 {
323 /* 2011.10.19 VIC Tamura */
324 SET_PROGRESS(99*(LOOP_COUNTER++)/num_tri, 66, 30); // critical, but works
325 if(IS_CANSELED()) continue;
326
327 LatticeFigureOfMeritToCheckSymmetry& latFOM = lattice_result_tri[n];
328 latFOM.setLabel(n+1);
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 /* 2011.10.19 VIC Tamura */
399 CHECK_INTERRUPTION();
400
401 ZLOG_INFO( "All the lattice constants are being optimized by linear least squares...\n" );
402
403 /* 2011.10.19 VIC Tamura */
404 LOOP_COUNTER = 0;
405 Int4 SUM = 0;
406 for(Int4 i=0; i<NUM_LS; i++) { SUM += lattice_result[i].size(); }
407
408 for(ArrayIndex i=1; i<NUM_LS; i++)
409 {
410 if( !JudgeSymmetry[i] ) continue;
411 sort( lattice_result[i].begin(), lattice_result[i].end() );
412
413 const Int4 num_lattice = lattice_result[i].size();
414
415 #ifdef _OPENMP
416 #pragma omp parallel
417 #endif
418 {
419 vector< VecDat3<Int4> > closest_hkl_tray;
420 vector<bool> is_cal_Q_observed_tray;
421 vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;
422
423 #ifdef _OPENMP
424 #pragma omp for
425 #endif
426 for(Int4 index=0; index<num_lattice; index++)
427 {
428 /* 2011.10.19 VIC Tamura */
429 SET_PROGRESS(99+1*(LOOP_COUNTER++)/SUM, 66, 30); // critical, but works
430 if(IS_CANSELED()) continue;
431
432 LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
433 latFOM0.setLabel(index+1);
434
435 latFOM0.setFigureOfMerit(m_num_ref_figure_of_merit,
436 VCData::putPeakQData(),
437 closest_hkl_tray, is_cal_Q_observed_tray);
438
439 LatticeFigureOfMeritZeroShift latFOM2 = latFOM0.putLatticeFigureOfMerit();
440 pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
441 closest_hkl_tray, is_cal_Q_observed_tray, false);
442 if( ans.first )
443 {
444 assert( latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
445 latFOM2.setFigureOfMerit(latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
446 VCData::putPeakQData(),
447 closest_hkl_tray, is_cal_Q_observed_tray);
448 if( LatticeFigureOfMerit::cmpFOMdeWolff( latFOM2, latFOM0.putLatticeFigureOfMerit() ) )
449 {
450 latFOM0.setLatticeFigureOfMerit(latFOM2);
451 }
452 }
453
454 const SetOfFigureOfMerit& setFOM = latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit();
455 if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
456 if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;
457
458 if( eCrystalSystem(i) == Monoclinic_P )
459 {
460 if( JudgeSymmetry[Hexagonal] )
461 {
462 putLatticeFigureOfMerit(latFOM0, D6h, m_cv2, latFOM_tray);
463 #ifdef _OPENMP
464 #pragma omp critical (hex)
465 #endif
466 {
467 lattice_result_hex.insert(lattice_result_hex.end(), latFOM_tray.begin(), latFOM_tray.end());
468 }
469 }
470 }
471 else if( eCrystalSystem(i) == Monoclinic_B )
472 {
473 if( JudgeSymmetry[Orthorhombic_C] )
474 {
475 putLatticeFigureOfMerit(latFOM0, D2h, m_cv2, latFOM_tray);
476 #ifdef _OPENMP
477 #pragma omp critical (ortho_B)
478 #endif
479 {
480 lattice_result_ortho_B.insert(lattice_result_ortho_B.end(), latFOM_tray.begin(), latFOM_tray.end());
481 }
482 }
483 }
484 else if( eCrystalSystem(i) == Orthorhombic_P )
485 {
486 if( JudgeSymmetry[Tetragonal_P] )
487 {
488 putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
489 #ifdef _OPENMP
490 #pragma omp critical (tetra_P)
491 #endif
492 {
493 lattice_result_tetra_P.insert(lattice_result_tetra_P.end(), latFOM_tray.begin(), latFOM_tray.end());
494 }
495 }
496 if( JudgeSymmetry[Cubic_P] )
497 {
498 putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
499 #ifdef _OPENMP
500 #pragma omp critical (cubic_P)
501 #endif
502 {
503 lattice_result_cubic_P.insert(lattice_result_cubic_P.end(), latFOM_tray.begin(), latFOM_tray.end());
504 }
505 }
506 }
507 else if( eCrystalSystem(i) == Orthorhombic_I )
508 {
509 if( JudgeSymmetry[Tetragonal_I] )
510 {
511 putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
512 #ifdef _OPENMP
513 #pragma omp critical (tetra_I)
514 #endif
515 {
516 lattice_result_tetra_I.insert(lattice_result_tetra_I.end(), latFOM_tray.begin(), latFOM_tray.end());
517 }
518 }
519 if( JudgeSymmetry[Cubic_I] )
520 {
521 putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
522 #ifdef _OPENMP
523 #pragma omp critical (cubic_I)
524 #endif
525 {
526 lattice_result_cubic_I.insert(lattice_result_cubic_I.end(), latFOM_tray.begin(), latFOM_tray.end());
527 }
528 }
529 }
530 else if( eCrystalSystem(i) == Orthorhombic_F )
531 {
532 if( JudgeSymmetry[Cubic_F] )
533 {
534 putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
535 #ifdef _OPENMP
536 #pragma omp critical (cubic_F)
537 #endif
538 {
539 lattice_result_cubic_F.insert(lattice_result_cubic_F.end(), latFOM_tray.begin(), latFOM_tray.end());
540 }
541 }
542 }
543 }
544 }
545 /* 2011.10.19 VIC Tamura */
546 CHECK_INTERRUPTION();
547
548 ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), abc_axis)
549 + " : " + num2str<Int4>( lattice_result[i].size() ) + "\n" );
550 }
551 ZLOG_INFO( "\n" );
552 }
553 catch(bad_alloc& ball)
554 {
555 throw nerror(ball, __FILE__, __LINE__, __FUNCTION__);
556 }
557 }
558
559
560 void SortingLattice::putLatticeCandidatesForEachBravaisTypes(const vector<SymMat43_VCData>& S_super,
561 const Double& MIN_NormM,
562 const Double& MIN_RevM,
563 const eABCaxis& abc_axis,
564 const eRHaxis& rh_axis,
565 vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
566 {
567 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(Int4)Triclinic];
568 putLatticeCandidatesForTriclinic(S_super, MIN_NormM, MIN_RevM, lattice_result_tri);
569
570 ZLOG_INFO( "Determination of the Bravais type is being carried out...\n(Solutions with " + putLabel(SCWuM) + " less than " + num2str(MIN_NormM)
571 + " or " + putLabel(SCRevM) + " less than " + num2str(MIN_RevM)
572 + " are automatically removed).\n\n" );
573
574 ZLOG_INFO( "The program has removed " + num2str<Int4>( S_super.size() - lattice_result_tri.size() )
575 + " solutions because their " + putLabel(SCWuM) + " is less than " + num2str(MIN_NormM)
576 + " or their " + putLabel(SCRevM) + " is less than " + num2str(MIN_RevM) + ".\n\n" );
577 ZLOG_INFO( "Determination of the Bravais type is being carried out...\n" );
578 putLatticeCandidatesForEachBravaisTypes(MIN_NormM, MIN_RevM, abc_axis, rh_axis, lattice_result);
579 }
580
581
582 void SortingLattice::setNumberOfNeighbors(const eABCaxis& baxis_flag,
583 vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
584 {
585
586 #ifdef _OPENMP
587 #pragma omp for
588 #endif
589 for(ArrayIndex i=0; i<NUM_LS; i++)
590 {
591 if( !OutputSymmetry[i] ) continue;
592
593 stable_sort( lattice_result[i].begin(), lattice_result[i].end() ); // Sort by the unit-cell volume.
594 for(vector<LatticeFigureOfMeritToCheckSymmetry>::iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
595 {
596 it->putNumberOfLatticesInNeighborhood() = 0;
597 }
598 }
599
600 const Double coef_lower = 1.0 - sqrt(3.0)*m_resol2;
601 const Double coef_upper = 1.0 + sqrt(3.0)*m_resol2;
602
603 Vec_INT index_tray(put_cs_num());
604
605 /* 2011.10.19 VIC Tamura */
606 Int4 SUM=0, LOOP_COUNTER=0;
607 for(Int4 i=0; i<NUM_LS; i++ ) { SUM += lattice_result[i].size(); }
608
609 #ifdef _OPENMP
610 #pragma omp for
611 #endif
612 for(ArrayIndex i=0; i<NUM_LS; i++)
613 {
614 if( !OutputSymmetry[i] ) continue;
615
616 const Int4 num_lattice = lattice_result[i].size();
617
618 for(Int4 index=0; index<num_lattice; index++)
619 {
620 /* 2011.10.19 VIC Tamura */
621 SET_PROGRESS(100*(LOOP_COUNTER++)/SUM, 97, 1); // critical, but works
622 if(IS_CANSELED()) continue;
623
624 LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
625 const LatticeFigureOfMerit& latFOM0_prim = latFOM0.putLatticeFigureOfMerit();
626 if( latFOM0.putNumberOfLatticesInNeighborhood() < 0 ) continue;
627
628 const Double& detS = latFOM0_prim.putDeterminantOfGramMatrix();
629 const Int4 ibegin = distance( lattice_result[i].begin(), lower_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_lower ) );
630 const Int4 iend = distance( lattice_result[i].begin(), upper_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_upper ) );
631
632 Int4 count=0;
633 if( i == (Int4)Triclinic )
634 {
635 const ReducedLatticeToCheckEquiv RLCS(m_resol2, latFOM0_prim.putOptimizedForm());
636 for(Int4 index2=ibegin; index2<iend; index2++)
637 {
638 if( index2 == index ) continue;
639
640 LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
641 const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();
642
643 // lattice_result_tri[index2] equals trans_mat * RLCB.m_S_super_obtuse * Transpose(trans_mat)
644 if( RLCS.equiv( latFOM2_prim.putSellingReducedForm() ) )
645 {
646 // Compare the figures of merit.
647 if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
648 && LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
649 {
650 count = -1;
651 break;
652 }
653 else
654 {
655 count++;
656 latFOM2.putNumberOfLatticesInNeighborhood() = -1;
657 }
658 }
659 }
660 }
661 else
662 {
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 // *it2 equals trans_mat * RLCS.m_S_super_obtuse * Transpose(trans_mat)
671 if( check_equiv_m( latFOM0_prim.putInverseOfMinkowskiReducedForm(), latFOM2_prim.putInverseOfMinkowskiReducedForm(), m_resol2 ) )
672 {
673 // Compare the figures of merit.
674 if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
675 && LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
676 {
677 count = -1;
678 break;
679 }
680 else
681 {
682 count++;
683 latFOM2.putNumberOfLatticesInNeighborhood() = -1;
684 }
685 }
686 }
687 }
688
689 latFOM0.putNumberOfLatticesInNeighborhood() = count;
690 }
691
692 Int4& index = index_tray[i];
693 for(vector<LatticeFigureOfMeritToCheckSymmetry>::const_iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
694 {
695 if( it->putNumberOfLatticesInNeighborhood() >= 0 ) index++;
696 }
697 }
698
699 /* 2011.10.19 VIC Tamura */
700 CHECK_INTERRUPTION();
701
702 for(ArrayIndex i=0; i<NUM_LS; i++)
703 {
704 if( !OutputSymmetry[i] ) continue;
705 ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), baxis_flag)
706 + " : " + num2str( lattice_result[i].size() ) + " -> " + num2str( index_tray[i] ) + "\n" );
707 }
708 ZLOG_INFO( "\n" );
709 }
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help