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