25 use,
intrinsic :: iso_fortran_env
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 =
""
67 type(symm_op_t),
allocatable,
public :: non_symmorphic_ops(:)
68 integer,
public :: nops_nonsymmorphic = 0
71 generic ::
assignment(=) => copy
75 real(real64),
parameter,
public :: default_symprec = 1.e-5_real64
76 real(real64),
public :: symprec
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
91 integer,
intent(in) :: point_group
92 character(kind=c_char),
intent(out) :: name(*)
97 integer,
intent(in) :: point_group
98 character(kind=c_char),
intent(out) :: elements(*)
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
122 integer :: idir, isite, isite_symm, iop, verbosity, point_group
123 real(real64),
allocatable :: position(:, :)
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)
132 type(SpglibDataset) :: spg_dataset
133 type(SpglibSpacegroupType) :: spg_spacegroup
138 this%any_non_spherical = any(.not. spherical_site)
140 this%periodic_dim = space%periodic_dim
142 dim4syms = min(3, space%dim)
144 def_sym_comp = (n_sites < 100) .or. space%is_periodic()
145 def_sym_comp = def_sym_comp .and. space%dim == 3
158 call parse_variable(namespace,
'SymmetriesCompute', def_sym_comp, this%symmetries_compute)
170 if (this%symmetries_compute .and. space%dim /= 3)
then
174 if (this%any_non_spherical .or. .not. this%symmetries_compute)
then
182 if (space%periodic_dim == 0)
then
190 safe_allocate(position(1:3, 1:n_sites))
192 do isite = 1, n_sites
193 position(1:3, isite) =
m_zero
194 position(1:dim4syms, isite) = site_pos(1:dim4syms, isite)
197 if (this%symmetries_compute)
then
199 if (point_group > -1)
then
209 this%group_elements =
""
211 this%schoenflies =
""
212 this%symmetries_compute = .false.
216 safe_deallocate_a(position)
222 if (spg_dataset%spglib_error /= 0)
then
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
233 do iop = 1, spg_dataset%n_operations
235 spg_dataset%rotations(:,:,iop) = transpose(spg_dataset%rotations(:,:,iop))
237 spg_dataset%translations(:, iop) = spg_dataset%translations(:, iop) - &
244 is_supercell = (spg_dataset%n_operations > 48)
245 found_identity = .false.
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)
257 if (.not. found_identity)
then
258 message(1) =
"Symmetries internal error: Identity is missing from symmetry operations."
262 if (is_supercell)
then
263 message(1) =
"Disabling fractional translations. System appears to be a supercell."
281 this%breakdir(1:3) =
m_zero
283 if (
parse_block(namespace,
'SymmetryBreakDir', blk) == 0)
then
285 do idir = 1, dim4syms
293 safe_allocate(this%ops(1:spg_dataset%n_operations))
294 safe_allocate(this%non_symmorphic_ops(1:spg_dataset%n_operations))
300 this%nops_nonsymmorphic = 0
301 do iop = 1, spg_dataset%n_operations
303 if (all(abs(spg_dataset%translations(1:3, iop)) <
tol_translation))
then
304 spg_dataset%translations(1:3, iop) =
m_zero
308 if (any(abs(spg_dataset%translations(space%periodic_dim+1:3, iop)) >
tol_translation)) cycle
310 call symm_op_init(tmpop, spg_dataset%rotations(1:3, 1:3, iop), latt, dim4syms, spg_dataset%translations(1:3, iop))
313 if( all(abs(this%breakdir) <
m_epsilon) )
then
315 this%nops = this%nops + 1
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))
324 this%nops = this%nops + 1
349 do isite = 1, n_sites
352 rsite = latt%fold_into_cell(rsite)
355 do isite_symm = 1, n_sites
356 if (all(abs(rsite - site_pos(:, isite_symm)) < symprec))
exit
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,
'.'
369 do iop = 1, this%nops_nonsymmorphic
370 do isite = 1, n_sites
373 rsite = latt%fold_into_cell(rsite)
376 do isite_symm = 1, n_sites
377 if (all(abs(rsite - site_pos(:, isite_symm)) < symprec))
exit
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,
'.'
397 safe_allocate(this%ops(1))
400 this%nops_nonsymmorphic = 0
419 safe_allocate(lhs%ops(1: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))
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
449 safe_deallocate_a(this%ops)
450 safe_deallocate_a(this%non_symmorphic_ops)
467 integer,
intent(in) :: iop
468 real(real64),
intent(in) :: aa(1:3)
469 real(real64),
intent(out) :: bb(1:3)
473 assert(iop /= 0 .and. abs(iop) <= this%nops)
477 bb(1:3) = symm_op_apply_transpose_red(this%ops(abs(iop)), -aa(1:3))
479 bb(1:3) = symm_op_apply_transpose_red(this%ops(abs(iop)), aa(1:3))
490 number = this%space_group
495 logical pure function symmetries_have_break_dir(this) result(have)
498 have = any(abs(this%breakdir(1:3)) > m_epsilon)
503 integer pure function symmetries_identity_index(this) result(index)
508 do iop = 1, this%nops
509 if (symm_op_is_identity(this%ops(iop)))
then
520 class(space_t),
intent(in) :: space
521 integer,
optional,
intent(in) :: iunit
522 type(namespace_t),
optional,
intent(in) :: namespace
528 call messages_print_with_emphasis(msg=
'Symmetries', iunit=iunit, namespace=namespace)
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)
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)
546 if (space%periodic_dim == 0)
then
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)
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)
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
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))
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))
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))
581 call messages_info(1, iunit=iunit, namespace=namespace)
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
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))
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))
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))
602 call messages_info(1, iunit=iunit, namespace=namespace)
605 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
614 type(lattice_vectors_t),
intent(in) :: latt
615 integer,
intent(in) :: dim
621 do iop = 1, this%nops
622 call symm_op_build_cartesian(this%ops(iop), latt, dim)
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)
637 integer :: isite, idir, dim4syms
638 real(real64),
allocatable :: position(:, :)
639 real(real64) :: lattice(1:3, 1:3)
643 assert(space%is_periodic())
645 dim4syms = min(3, space%dim)
647 safe_allocate(position(1:3, 1:n_sites))
649 do isite = 1, n_sites
650 position(1:3,isite) = m_zero
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))
658 position(1:dim4syms, isite) = position(1:dim4syms,isite) + m_half
664 lattice(1:space%periodic_dim, 1:space%periodic_dim) = latt%rlattice(1:space%periodic_dim, 1:space%periodic_dim)
666 lattice(:,:) = transpose(lattice(:,:))
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
673 lattice(idir, idir) = m_one
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)
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)
688 safe_deallocate_a(position)
NOTE: unfortunately, these routines use global variables shared among them.
real(real64), parameter, public m_zero
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
This module is intended to contain "only mathematical" functions and procedures.
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
type(mpi_grp_t), public mpi_world
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public string_c_to_f(c_string, f_string)
convert a C string to a Fortran string
subroutine, public symm_op_copy(inp, outp)
logical pure function, public symm_op_has_translation(this, prec)
real(real64), parameter, public tol_translation
subroutine, public symm_op_init(this, rot, latt, dim, trans)
subroutine, public symmetries_update_lattice_vectors(this, latt, dim)
Updates the symmetry operations when lattice vectors are updated.
subroutine, public symmetries_apply_kpoint_red(this, iop, aa, bb)
integer pure function, public symmetries_space_group_number(this)
integer pure function, public symmetries_identity_index(this)
type(symmetries_t) function symmetries_constructor(namespace, space, latt, n_sites, site_pos, site_type, spherical_site)
subroutine, public symmetries_write_info(this, space, iunit, namespace)
real(real64), public symprec
integer pure function, public symmetries_number(this)
subroutine symmetries_copy(lhs, rhs)
logical pure function, public symmetries_have_break_dir(this)
type(spglibdataset) function, public symmetries_get_spg_dataset(namespace, space, latt, n_sites, site_pos, site_type)
Returns the SpglibDataset for the system.
subroutine symmetries_finalizer(this)
subroutine init_identity()