Octopus
symmetries.F90
Go to the documentation of this file.
1!! Copyright (C) 2009 X. Andrade
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22 use debug_oct_m
23 use global_oct_m
24 use iso_c_binding
25 use, intrinsic :: iso_fortran_env
27 use math_oct_m
29 use mpi_oct_m
31 use parser_oct_m
33 use string_oct_m
34 use space_oct_m
35 use spglib_f08
37
38 implicit none
39
40 private
41
42 public :: &
52
53 type symmetries_t
54 private
55 type(symm_op_t), allocatable, public :: ops(:)
56 integer, public :: nops = 0
57 real(real64) :: breakdir(1:3)
58 integer, public :: periodic_dim
59 integer :: space_group
60 logical :: any_non_spherical
61 logical :: symmetries_compute = .false.
62 character(len=6) :: group_name = ""
63 character(len=30) :: group_elements = ""
64 character(len=11) :: symbol = ""
65 character(len=7) :: schoenflies = ""
66
67 type(symm_op_t), allocatable, public :: non_symmorphic_ops(:)
68 integer, public :: nops_nonsymmorphic = 0
69 contains
70 procedure :: copy => symmetries_copy
71 generic :: assignment(=) => copy
73 end type symmetries_t
74
75 real(real64), parameter, public :: default_symprec = 1.e-5_real64
76 real(real64), public :: symprec
77
79 interface
80 subroutine symmetries_finite_init(natoms, types, positions, verbosity, point_group)
81 use, intrinsic :: iso_fortran_env
82 integer, intent(in) :: natoms
83 integer, intent(in) :: types
84 real(real64), intent(in) :: positions
85 integer, intent(in) :: verbosity
86 integer, intent(out) :: point_group
87 end subroutine symmetries_finite_init
88
89 subroutine symmetries_finite_get_group_name(point_group, name) bind(c)
90 use iso_c_binding
91 integer, intent(in) :: point_group
92 character(kind=c_char), intent(out) :: name(*)
94
95 subroutine symmetries_finite_get_group_elements(point_group, elements) bind(c)
96 use iso_c_binding
97 integer, intent(in) :: point_group
98 character(kind=c_char), intent(out) :: elements(*)
100
101 subroutine symmetries_finite_end()
102 end subroutine symmetries_finite_end
103 end interface
104
105 interface symmetries_t
106 module procedure symmetries_constructor
107 end interface symmetries_t
108
109contains
110
111 function symmetries_constructor(namespace, space, latt, n_sites, site_pos, site_type, spherical_site) result(this)
112 type(namespace_t), intent(in) :: namespace
113 class(space_t), intent(in) :: space
114 type(lattice_vectors_t), intent(in) :: latt
115 integer, intent(in) :: n_sites
116 real(real64), intent(in) :: site_pos(:,:)
117 integer, intent(in) :: site_type(1:n_sites)
118 logical, intent(in) :: spherical_site(:)
119 type(symmetries_t) :: this
120
121 integer :: dim4syms
122 integer :: idir, isite, isite_symm, iop, verbosity, point_group
123 real(real64), allocatable :: position(:, :)
124 type(block_t) :: blk
125 type(symm_op_t) :: tmpop
126 integer :: identity(3,3)
127 logical :: found_identity, is_supercell
128 logical :: def_sym_comp
129 real(real64) :: rsite(space%dim)
130 character(kind=c_char) :: c_str(31)
131
132 type(SpglibDataset) :: spg_dataset
133 type(SpglibSpacegroupType) :: spg_spacegroup
134
135 push_sub(symmetries_constructor)
136
137 ! if someone cares, they could try to analyze the symmetry point group of the individual species too
138 this%any_non_spherical = any(.not. spherical_site)
139
140 this%periodic_dim = space%periodic_dim
141
142 dim4syms = min(3, space%dim)
143
144 def_sym_comp = (n_sites < 100) .or. space%is_periodic()
145 def_sym_comp = def_sym_comp .and. space%dim == 3
146
147 !%Variable SymmetriesCompute
148 !%Type logical
149 !%Section Execution::Symmetries
150 !%Description
151 !% If disabled, <tt>Octopus</tt> will not compute
152 !% nor print the symmetries.
153 !%
154 !% By default, symmetries are computed when running in 3
155 !% dimensions for systems with less than 100 atoms.
156 !% For periodic systems, the default is always true, irrespective of the number of atoms.
157 !%End
158 call parse_variable(namespace, 'SymmetriesCompute', def_sym_comp, this%symmetries_compute)
160 !%Variable SymmetriesTolerance
161 !%Type float
162 !%Section Execution::Symmetries
163 !%Description
164 !% For periodic systems, this variable controls the tolerance used by the symmetry finder
165 !% (spglib) to find the spacegroup and symmetries of the crystal.
166 !%End
167 call parse_variable(namespace, 'SymmetriesTolerance', default_symprec, symprec)
168
169
170 if (this%symmetries_compute .and. space%dim /= 3) then
171 call messages_experimental('symmetries for non 3D systems', namespace=namespace)
172 end if
173
174 if (this%any_non_spherical .or. .not. this%symmetries_compute) then
177 return
178 end if
179
180 ! In all cases, we must check that the grid respects the symmetries. --DAS
181
182 if (space%periodic_dim == 0) then
183
185
186 ! for the moment symmetries are only used for information, so we compute them only on one node.
187 if (mpi_world%is_root()) then
188 verbosity = -1
189
190 safe_allocate(position(1:3, 1:n_sites))
191
192 do isite = 1, n_sites
193 position(1:3, isite) = m_zero
194 position(1:dim4syms, isite) = site_pos(1:dim4syms, isite)
195 end do
197 if (this%symmetries_compute) then
198 call symmetries_finite_init(n_sites, site_type(1), position(1, 1), verbosity, point_group)
199 if (point_group > -1) then
200 c_str = c_null_char
201 call symmetries_finite_get_group_name(point_group, c_str)
202 call string_c_to_f(c_str, this%group_name)
203
204 c_str = c_null_char
205 call symmetries_finite_get_group_elements(point_group, c_str)
206 call string_c_to_f(c_str, this%group_elements)
207 else
208 this%group_name = ""
209 this%group_elements = ""
210 this%symbol = ""
211 this%schoenflies = ""
212 this%symmetries_compute = .false.
213 end if
215 end if
216 safe_deallocate_a(position)
217 end if
218
219 else
220
221 spg_dataset = symmetries_get_spg_dataset(namespace, space, latt, n_sites, site_pos, site_type)
222 if (spg_dataset%spglib_error /= 0) then
223 call init_identity()
225 return
226 end if
227
228 this%space_group = spg_dataset%spacegroup_number
229 this%symbol = spg_dataset%international_symbol
230 spg_spacegroup = spg_get_spacegroup_type(spg_dataset%hall_number)
231 this%schoenflies = spg_spacegroup%schoenflies
232
233 do iop = 1, spg_dataset%n_operations
234 ! transpose due to array order difference between C and fortran
235 spg_dataset%rotations(:,:,iop) = transpose(spg_dataset%rotations(:,:,iop))
236 ! sometimes spglib may return lattice vectors as 'fractional' translations
237 spg_dataset%translations(:, iop) = spg_dataset%translations(:, iop) - &
238 anint(spg_dataset%translations(:, iop) + m_half * tol_translation)
239 end do
240
241 ! we need to check that it is not a supercell, as in the QE routine (sgam_at)
242 ! they disable fractional translations if the identity has one, because the sym ops might not form a group.
243 ! spglib may return duplicate operations in this case!
244 is_supercell = (spg_dataset%n_operations > 48)
245 found_identity = .false.
246 identity = diagonal_matrix(3, 1)
247 do iop = 1, spg_dataset%n_operations
248 if (all(spg_dataset%rotations(1:3, 1:3, iop) == identity(1:3, 1:3))) then
249 found_identity = .true.
250 if (any(abs(spg_dataset%translations(1:3, iop)) > tol_translation )) then
251 is_supercell = .true.
252 write(message(1),'(a,3f12.6)') 'Identity has a fractional translation ', spg_dataset%translations(1:3, iop)
253 call messages_info(1, namespace=namespace)
254 end if
255 end if
256 end do
257 if (.not. found_identity) then
258 message(1) = "Symmetries internal error: Identity is missing from symmetry operations."
259 call messages_fatal(1, namespace=namespace)
260 end if
261
262 if (is_supercell) then
263 message(1) = "Disabling fractional translations. System appears to be a supercell."
264 call messages_info(1, namespace=namespace)
265 end if
266 ! actually, we do not use fractional translations regardless currently
267
268 ! this is a hack to get things working, this variable should be
269 ! eliminated and the direction calculated automatically from the
270 ! perturbations.
271
272 !%Variable SymmetryBreakDir
273 !%Type block
274 !%Section Mesh::Simulation Box
275 !%Description
276 !% This variable specifies a direction in which the symmetry of
277 !% the system will be broken. This is useful for generating <i>k</i>-point
278 !% grids when an external perturbation is applied.
279 !%End
280
281 this%breakdir(1:3) = m_zero
282
283 if (parse_block(namespace, 'SymmetryBreakDir', blk) == 0) then
284
285 do idir = 1, dim4syms
286 call parse_block_float(blk, 0, idir - 1, this%breakdir(idir))
287 end do
288
289 call parse_block_end(blk)
290
291 end if
292
293 safe_allocate(this%ops(1:spg_dataset%n_operations))
294 safe_allocate(this%non_symmorphic_ops(1:spg_dataset%n_operations))
295
296 ! check all operations and leave those that kept the symmetry-breaking
297 ! direction invariant and (for the moment) that do not have a translation
298 this%nops = 0
299 ! We also keep track of all non-symmorphic operations
300 this%nops_nonsymmorphic = 0
301 do iop = 1, spg_dataset%n_operations
302 ! If spglib returns tiny translations, we set them to zero to avoid numerical artifacts
303 if (all(abs(spg_dataset%translations(1:3, iop)) < tol_translation)) then
304 spg_dataset%translations(1:3, iop) = m_zero
305 end if
306
307 ! Reject all symmetries involving fractional translation in the non-periodic dimension
308 if (any(abs(spg_dataset%translations(space%periodic_dim+1:3, iop)) > tol_translation)) cycle
309
310 call symm_op_init(tmpop, spg_dataset%rotations(1:3, 1:3, iop), latt, dim4syms, spg_dataset%translations(1:3, iop))
311
312 ! If we do not specify SymmetryBreakDir, we want to use all symmetries
313 if( all(abs(this%breakdir) < m_epsilon) ) then
314 if (.not. symm_op_has_translation(tmpop, tol_translation)) then
315 this%nops = this%nops + 1
316 call symm_op_copy(tmpop, this%ops(this%nops))
317 else if(.not. is_supercell) then
318 this%nops_nonsymmorphic = this%nops_nonsymmorphic + 1
319 call symm_op_copy(tmpop, this%non_symmorphic_ops(this%nops_nonsymmorphic))
320 end if
321 else ! We only test symmorphic symmetries
322 if (symm_op_invariant_cart(tmpop, this%breakdir, symprec) .and. &
323 .not. symm_op_has_translation(tmpop, symprec)) then
324 this%nops = this%nops + 1
325 call symm_op_copy(tmpop, this%ops(this%nops))
326 end if
327 end if
328 end do
329
330
331 end if
332
333 call symmetries_write_info(this, space, namespace=namespace)
334
335 ! Checks if that the atomic coordinates are compatible with the symmetries
336 !
337 ! We want to use for instance that
338 !
339 ! \int dr f(Rr) V_isite(r) \nabla f(R(v)) = R\int dr f(r) V_isite(R*r) f(r)
340 !
341 ! and that the operator R should map the position of atom
342 ! isite to the position of some other atom isite_symm, so that
343 !
344 ! V_isite(R*r) = V_isite_symm(r)
345 !
346 do iop = 1, symmetries_number(this)
347 if (iop == symmetries_identity_index(this)) cycle
348
349 do isite = 1, n_sites
350 rsite = symm_op_apply_cart(this%ops(iop), site_pos(:, isite))
351
352 rsite = latt%fold_into_cell(rsite)
353
354 ! find isite_symm
355 do isite_symm = 1, n_sites
356 if (all(abs(rsite - site_pos(:, isite_symm)) < symprec)) exit
357 end do
358
359 if (isite_symm > n_sites) then
360 write(message(1),'(a,i6)') 'Internal error: could not find symmetric partner for atom number', isite
361 write(message(2),'(a,i3,a)') 'with symmetry operation number ', iop, '.'
362 call messages_fatal(2, namespace=namespace)
363 end if
364
365 end do
366 end do
367
368 ! Check non-symmorphic operations
369 do iop = 1, this%nops_nonsymmorphic
370 do isite = 1, n_sites
371 rsite = symm_op_apply_cart(this%non_symmorphic_ops(iop), site_pos(:, isite))
372
373 rsite = latt%fold_into_cell(rsite)
374
375 ! find isite_symm
376 do isite_symm = 1, n_sites
377 if (all(abs(rsite - site_pos(:, isite_symm)) < symprec)) exit
378 end do
379
380 if (isite_symm > n_sites) then
381 write(message(1),'(a,i6)') 'Internal error: could not find symmetric partner for atom number', isite
382 write(message(2),'(a,i3,a)') 'with non-symmorphic symmetry operation number ', iop, '.'
383 call messages_fatal(2, namespace=namespace)
384 end if
385
386 end do
387 end do
388
390
391 contains
392
393 subroutine init_identity()
394
395 push_sub(symmetries_init.init_identity)
396
397 safe_allocate(this%ops(1))
398 this%nops = 1
399 call symm_op_init(this%ops(1), diagonal_matrix(3, 1), latt, dim4syms)
400 this%nops_nonsymmorphic = 0
401 this%breakdir = m_zero
402 this%space_group = 1
403
404 pop_sub(symmetries_init.init_identity)
405
406 end subroutine init_identity
407
408 end function symmetries_constructor
409
410 ! -------------------------------------------------------------------------------
411 subroutine symmetries_copy(lhs, rhs)
412 class(symmetries_t), intent(out) :: lhs
413 class(symmetries_t), intent(in) :: rhs
414
415 integer :: iop
416
417 push_sub(symmetries_copy)
418
419 safe_allocate(lhs%ops(1:rhs%nops))
420 do iop = 1, rhs%nops
421 call symm_op_copy(rhs%ops(iop), lhs%ops(iop))
422 end do
423 lhs%nops = rhs%nops
424 safe_allocate(lhs%non_symmorphic_ops(1:rhs%nops_nonsymmorphic))
425 do iop = 1, rhs%nops_nonsymmorphic
426 call symm_op_copy(rhs%non_symmorphic_ops(iop), lhs%non_symmorphic_ops(iop))
427 end do
428 lhs%nops_nonsymmorphic = rhs%nops_nonsymmorphic
429 lhs%breakdir = rhs%breakdir
430 lhs%periodic_dim = rhs%periodic_dim
431 lhs%space_group = rhs%space_group
432 lhs%any_non_spherical = rhs%any_non_spherical
433 lhs%symmetries_compute = rhs%symmetries_compute
434 lhs%group_name = rhs%group_name
435 lhs%group_elements = rhs%group_elements
436 lhs%symbol = rhs%symbol
437 lhs%schoenflies = rhs%schoenflies
438
439 pop_sub(symmetries_copy)
440 end subroutine symmetries_copy
441
442 ! -------------------------------------------------------------------------------
443
444 subroutine symmetries_finalizer(this)
445 type(symmetries_t), intent(inout) :: this
446
447 push_sub(symmetries_finalizer)
448
449 safe_deallocate_a(this%ops)
450 safe_deallocate_a(this%non_symmorphic_ops)
451
452 pop_sub(symmetries_finalizer)
453 end subroutine symmetries_finalizer
454
455 ! -------------------------------------------------------------------------------
456
457 integer pure function symmetries_number(this) result(number)
458 type(symmetries_t), intent(in) :: this
459
460 number = this%nops
461 end function symmetries_number
462
463 ! -------------------------------------------------------------------------------
464
465 subroutine symmetries_apply_kpoint_red(this, iop, aa, bb)
466 type(symmetries_t), intent(in) :: this
467 integer, intent(in) :: iop
468 real(real64), intent(in) :: aa(1:3)
469 real(real64), intent(out) :: bb(1:3)
470
472
473 assert(iop /= 0 .and. abs(iop) <= this%nops)
474
475 !Time-reversal symmetry
476 if (iop < 0) then
477 bb(1:3) = symm_op_apply_transpose_red(this%ops(abs(iop)), -aa(1:3))
478 else
479 bb(1:3) = symm_op_apply_transpose_red(this%ops(abs(iop)), aa(1:3))
480 end if
481
483 end subroutine symmetries_apply_kpoint_red
484
485 ! -------------------------------------------------------------------------------
486
487 integer pure function symmetries_space_group_number(this) result(number)
488 type(symmetries_t), intent(in) :: this
489
490 number = this%space_group
492
493 ! -------------------------------------------------------------------------------
494
495 logical pure function symmetries_have_break_dir(this) result(have)
496 type(symmetries_t), intent(in) :: this
497
498 have = any(abs(this%breakdir(1:3)) > m_epsilon)
499 end function symmetries_have_break_dir
500
501 ! -------------------------------------------------------------------------------
502
503 integer pure function symmetries_identity_index(this) result(index)
504 type(symmetries_t), intent(in) :: this
505
506 integer :: iop
507
508 do iop = 1, this%nops
509 if (symm_op_is_identity(this%ops(iop))) then
510 index = iop
511 cycle
512 end if
513 end do
514
515 end function symmetries_identity_index
516
517 ! ---------------------------------------------------------
518 subroutine symmetries_write_info(this, space, iunit, namespace)
519 type(symmetries_t), intent(in) :: this
520 class(space_t), intent(in) :: space
521 integer, optional, intent(in) :: iunit
522 type(namespace_t), optional, intent(in) :: namespace
523
524 integer :: iop
525
526 push_sub(symmetries_write_info)
527
528 call messages_print_with_emphasis(msg='Symmetries', iunit=iunit, namespace=namespace)
529
530 if (this%any_non_spherical) then
531 message(1) = "Symmetries are disabled since non-spherically symmetric species may be present."
532 call messages_info(1,iunit = iunit, namespace=namespace)
533 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
534 pop_sub(symmetries_write_info)
535 return
536 end if
537
538 if (.not. this%symmetries_compute) then
539 message(1) = "Symmetries have been disabled by SymmetriesCompute = false."
540 call messages_info(1, iunit=iunit, namespace=namespace)
541 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
542 pop_sub(symmetries_write_info)
543 return
544 end if
545
546 if (space%periodic_dim == 0) then
547 ! At the moment only the root node has information about symetries of finite systems.
548 if (mpi_world%is_root()) then
549 if (this%symmetries_compute) then
550 call messages_write('Symmetry elements : '//trim(this%group_elements), new_line = .true.)
551 call messages_write('Symmetry group : '//trim(this%group_name))
552 call messages_info(iunit=iunit, namespace=namespace)
553 end if
554 end if
555 else
556 write(message(1),'(a, i4)') 'Space group No. ', this%space_group
557 write(message(2),'(2a)') 'International: ', trim(this%symbol)
558 write(message(3),'(2a)') 'Schoenflies: ', trim(this%schoenflies)
559 call messages_info(3, iunit=iunit, namespace=namespace)
561 write(message(1), '(a,i5,a)') 'Info: The system has ', this%nops, ' symmorphic symmetries that can be used.'
562 call messages_info(1, iunit=iunit, namespace=namespace)
563
564 write(message(1),'(a7,a31,12x,a33)') 'Index', 'Rotation matrix', 'Fractional translations'
565 call messages_info(1, iunit=iunit, namespace=namespace)
566 do iop = 1, this%nops
567 ! list all operations and leave those that kept the symmetry-breaking
568 ! direction invariant and (for the moment) that do not have a translation
569 if (space%dim == 1) then
570 write(message(1),'(i5,1x,a,2x,1(1i4,2x),1f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), &
571 symm_op_translation_vector_red(this%ops(iop))
572 end if
573 if (space%dim == 2) then
574 write(message(1),'(i5,1x,a,2x,2(2i4,2x),2f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), &
575 symm_op_translation_vector_red(this%ops(iop))
576 end if
577 if (space%dim == 3) then
578 write(message(1),'(i5,1x,a,2x,3(3i4,2x),3f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), &
579 symm_op_translation_vector_red(this%ops(iop))
580 end if
581 call messages_info(1, iunit=iunit, namespace=namespace)
582 end do
583
584 write(message(1), '(a,i5,a)') 'Info: The system also has ', this%nops_nonsymmorphic, &
585 ' nonsymmorphic symmetries (not used for electrons).'
586 call messages_info(1, iunit=iunit, namespace=namespace)
587 do iop = 1, this%nops_nonsymmorphic
588 ! list all operations and leave those that kept the symmetry-breaking
589 ! direction invariant and (for the moment) that do not have a translation
590 if (space%dim == 1) then
591 write(message(1),'(i5,1x,a,2x,1(1i4,2x),1f12.6)') iop, ':', symm_op_rotation_matrix_red(this%non_symmorphic_ops(iop)), &
592 symm_op_translation_vector_red(this%non_symmorphic_ops(iop))
593 end if
594 if (space%dim == 2) then
595 write(message(1),'(i5,1x,a,2x,2(2i4,2x),2f12.6)') iop, ':', symm_op_rotation_matrix_red(this%non_symmorphic_ops(iop)), &
596 symm_op_translation_vector_red(this%non_symmorphic_ops(iop))
597 end if
598 if (space%dim == 3) then
599 write(message(1),'(i5,1x,a,2x,3(3i4,2x),3f12.6)') iop, ':', symm_op_rotation_matrix_red(this%non_symmorphic_ops(iop)), &
600 symm_op_translation_vector_red(this%non_symmorphic_ops(iop))
601 end if
602 call messages_info(1, iunit=iunit, namespace=namespace)
603 end do
604 end if
605 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
606
607 pop_sub(symmetries_write_info)
608 end subroutine symmetries_write_info
609
610 ! ---------------------------------------------------------
612 subroutine symmetries_update_lattice_vectors(this, latt, dim)
613 type(symmetries_t), intent(inout) :: this
614 type(lattice_vectors_t), intent(in) :: latt
615 integer, intent(in) :: dim
616
617 integer :: iop
618
620
621 do iop = 1, this%nops
622 call symm_op_build_cartesian(this%ops(iop), latt, dim)
623 end do
624
626 end subroutine
627
629 type(spglibdataset) function symmetries_get_spg_dataset(namespace, space, latt, n_sites, site_pos, site_type) result(spg_dataset)
630 type(namespace_t), intent(in) :: namespace
631 class(space_t), intent(in) :: space
632 type(lattice_vectors_t), intent(in) :: latt
633 integer, intent(in) :: n_sites
634 real(real64), intent(in) :: site_pos(:,:)
635 integer, intent(in) :: site_type(1:n_sites)
636
637 integer :: isite, idir, dim4syms
638 real(real64), allocatable :: position(:, :)
639 real(real64) :: lattice(1:3, 1:3)
640
642
643 assert(space%is_periodic())
644
645 dim4syms = min(3, space%dim)
646
647 safe_allocate(position(1:3, 1:n_sites))
648
649 do isite = 1, n_sites
650 position(1:3,isite) = m_zero
651
652 ! Transform atomic positions to reduced coordinates
653 position(1:dim4syms, isite) = latt%cart_to_red(site_pos(1:dim4syms, isite))
654 position(1:dim4syms, isite) = position(1:dim4syms, isite) - m_half
655 do idir = 1, dim4syms
656 position(idir, isite) = position(idir, isite) - anint(position(idir, isite))
657 end do
658 position(1:dim4syms, isite) = position(1:dim4syms,isite) + m_half
659 end do
660
661 lattice = m_zero
662 !NOTE: Why "inverse matrix" ? (NTD)
663 ! get inverse matrix to extract reduced coordinates for spglib
664 lattice(1:space%periodic_dim, 1:space%periodic_dim) = latt%rlattice(1:space%periodic_dim, 1:space%periodic_dim)
665 ! transpose the lattice vectors for use in spglib as row-major matrix
666 lattice(:,:) = transpose(lattice(:,:))
667 ! spglib always assume 3D periodicity, so we need to set the lattice vectors along the non-periodic dimensions.
668 ! Here we choose a value that should guarantee that no incorrect symmetries are found along the non-periodic dimensions.
669 do idir = space%periodic_dim + 1, 3
670 if (idir <= space%dim) then
671 lattice(idir, idir) = m_two*(maxval(site_pos(idir, :)) - minval(site_pos(idir, :))) + m_one
672 else
673 lattice(idir, idir) = m_one
674 end if
675 end do
676
677 spg_dataset = spg_get_dataset(lattice, position, site_type, n_sites, symprec)
678 if (spg_dataset%spglib_error /= 0) then
679 message(1) = "Symmetry analysis failed in spglib. Disabling symmetries."
680 call messages_warning(1, namespace=namespace)
681
682 do isite = 1, n_sites
683 write(message(1),'(a,i6,a,3f12.6,a,3f12.6)') 'type ', site_type(isite), &
684 ' reduced coords ', position(:, isite), ' cartesian coords ', site_pos(:, isite)
685 call messages_info(1, namespace=namespace)
686 end do
687 end if
688 safe_deallocate_a(position)
689
691 end function symmetries_get_spg_dataset
692
693end module symmetries_oct_m
694
695!! Local Variables:
696!! mode: f90
697!! coding: utf-8
698!! End:
NOTE: unfortunately, these routines use global variables shared among them.
Definition: symmetries.F90:175
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public string_c_to_f(c_string, f_string)
convert a C string to a Fortran string
Definition: string.F90:292
subroutine, public symm_op_copy(inp, outp)
Definition: symm_op.F90:278
logical pure function, public symm_op_has_translation(this, prec)
Definition: symm_op.F90:295
real(real64), parameter, public tol_translation
Definition: symm_op.F90:185
subroutine, public symm_op_init(this, rot, latt, dim, trans)
Definition: symm_op.F90:192
subroutine, public symmetries_update_lattice_vectors(this, latt, dim)
Updates the symmetry operations when lattice vectors are updated.
Definition: symmetries.F90:708
subroutine, public symmetries_apply_kpoint_red(this, iop, aa, bb)
Definition: symmetries.F90:561
integer pure function, public symmetries_space_group_number(this)
Definition: symmetries.F90:583
integer pure function, public symmetries_identity_index(this)
Definition: symmetries.F90:599
type(symmetries_t) function symmetries_constructor(namespace, space, latt, n_sites, site_pos, site_type, spherical_site)
Definition: symmetries.F90:207
subroutine, public symmetries_write_info(this, space, iunit, namespace)
Definition: symmetries.F90:614
real(real64), public symprec
Definition: symmetries.F90:171
integer pure function, public symmetries_number(this)
Definition: symmetries.F90:553
subroutine symmetries_copy(lhs, rhs)
Definition: symmetries.F90:507
logical pure function, public symmetries_have_break_dir(this)
Definition: symmetries.F90:591
type(spglibdataset) function, public symmetries_get_spg_dataset(namespace, space, latt, n_sites, site_pos, site_type)
Returns the SpglibDataset for the system.
Definition: symmetries.F90:725
subroutine symmetries_finalizer(this)
Definition: symmetries.F90:540
subroutine init_identity()
Definition: symmetries.F90:489
int true(void)