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 21module symmetries_oct_m 22 use iso_c_binding 23 use geometry_oct_m 24 use global_oct_m 25 use messages_oct_m 26 use mpi_oct_m 27 use namespace_oct_m 28 use parser_oct_m 29 use profiling_oct_m 30 use species_oct_m 31 use spglib_f08 32 use symm_op_oct_m 33 34 implicit none 35 36 private 37 38 public :: & 39 symmetries_t, & 40 symmetries_init, & 41 symmetries_copy, & 42 symmetries_end, & 43 symmetries_number, & 44 symmetries_apply_kpoint_red, & 45 symmetries_space_group_number, & 46 symmetries_have_break_dir, & 47 symmetries_identity_index, & 48 symmetries_write_info 49 50 type symmetries_t 51 private 52 type(symm_op_t), allocatable, public :: ops(:) 53 integer, public :: nops 54 FLOAT :: breakdir(1:3) 55 integer :: space_group 56 logical :: any_non_spherical 57 logical :: symmetries_compute 58 character(len=6) :: group_name 59 character(len=30) :: group_elements 60 character(len=10) :: symbol 61 character(len=6) :: schoenflies 62 end type symmetries_t 63 64 FLOAT, parameter, public :: SYMPREC = CNST(1e-5) 65 66 !> NOTE: unfortunately, these routines use global variables shared among them 67 interface 68 subroutine symmetries_finite_init(natoms, types, positions, verbosity, point_group) 69 integer, intent(in) :: natoms 70 integer, intent(in) :: types !< (natoms) 71 FLOAT, intent(in) :: positions !< (3, natoms) 72 integer, intent(in) :: verbosity 73 integer, intent(out) :: point_group 74 end subroutine symmetries_finite_init 75 76 subroutine symmetries_finite_get_group_name(point_group, name) 77 integer, intent(in) :: point_group 78 character(len=*), intent(out) :: name 79 end subroutine symmetries_finite_get_group_name 80 81 subroutine symmetries_finite_get_group_elements(point_group, elements) 82 integer, intent(in) :: point_group 83 character(len=*), intent(out) :: elements 84 end subroutine symmetries_finite_get_group_elements 85 86 subroutine symmetries_finite_end() 87 end subroutine symmetries_finite_end 88 end interface 89 90contains 91 92 subroutine symmetries_init(this, namespace, geo, dim, periodic_dim, rlattice, klattice) 93 type(symmetries_t), intent(out) :: this 94 type(namespace_t), intent(in) :: namespace 95 type(geometry_t), intent(in) :: geo 96 integer, intent(in) :: dim 97 integer, intent(in) :: periodic_dim 98 FLOAT, intent(in) :: rlattice(:, :) 99 FLOAT, intent(in) :: klattice(:, :) 100 101 integer :: max_size, dim4syms 102 integer :: idir, iatom, iop, verbosity, point_group 103 FLOAT :: lattice(1:3, 1:3) 104 FLOAT, allocatable :: position(:, :) 105 integer, allocatable :: typs(:) 106 type(block_t) :: blk 107 type(symm_op_t) :: tmpop 108 integer :: natoms, identity(3,3) 109 logical :: found_identity, is_supercell 110 integer :: fullnops 111 integer, allocatable :: rotation(:, :, :) 112 FLOAT, allocatable :: translation(:, :) 113 character(kind=c_char) :: c_symbol(11), c_schoenflies(7) 114 logical :: def_sym_comp 115 116 PUSH_SUB(symmetries_init) 117 118 ! if someone cares, they could try to analyze the symmetry point group of the individual species too 119 this%any_non_spherical = .false. 120 do iatom = 1, geo%natoms 121 this%any_non_spherical = this%any_non_spherical .or. & 122 species_type(geo%atom(iatom)%species) == SPECIES_USDEF .or. & 123 species_type(geo%atom(iatom)%species) == SPECIES_JELLIUM_SLAB .or. & 124 species_type(geo%atom(iatom)%species) == SPECIES_CHARGE_DENSITY .or. & 125 species_type(geo%atom(iatom)%species) == SPECIES_FROM_FILE .or. & 126 species_type(geo%atom(iatom)%species) == SPECIES_FROZEN 127 if(this%any_non_spherical)exit 128 end do 129 130 dim4syms = min(3,dim) 131 132 def_sym_comp = (geo%natoms < 100) .or. periodic_dim > 0 133 def_sym_comp = def_sym_comp .and. dim == 3 134 135 !%Variable SymmetriesCompute 136 !%Type logical 137 !%Section Execution::Symmetries 138 !%Description 139 !% If disabled, <tt>Octopus</tt> will not compute 140 !% nor print the symmetries. 141 !% 142 !% By default, symmetries are computed when running in 3 143 !% dimensions for systems with less than 100 atoms. 144 !% For periodic systems, the default is always true, irrespective of the number of atoms. 145 !%End 146 call parse_variable(namespace, 'SymmetriesCompute', def_sym_comp, this%symmetries_compute) 147 148 if(this%symmetries_compute .and. dim /= 3) then 149 call messages_experimental('symmetries for non 3D systems') 150 end if 151 152 if(this%any_non_spherical .or. .not. this%symmetries_compute) then 153 call init_identity() 154 155 POP_SUB(symmetries_init) 156 return 157 end if 158 159 ! In all cases, we must check that the grid respects the symmetries. --DAS 160 161 if (periodic_dim == 0) then 162 163 call init_identity() 164 165 ! for the moment symmetries are only used for information, so we compute them only on one node. 166 if(mpi_grp_is_root(mpi_world)) then 167 natoms = max(1,geo%natoms) 168 verbosity = -1 169 170 SAFE_ALLOCATE(position(1:3, natoms)) 171 SAFE_ALLOCATE(typs(1:natoms)) 172 173 do iatom = 1, geo%natoms 174 position(1:3, iatom) = M_ZERO 175 position(1:dim4syms, iatom) = geo%atom(iatom)%x(1:dim4syms) 176 typs(iatom) = species_index(geo%atom(iatom)%species) 177 end do 178 179 if (this%symmetries_compute) then 180 call symmetries_finite_init(geo%natoms, typs(1), position(1, 1), verbosity, point_group) 181 if(point_group > -1) then 182 call symmetries_finite_get_group_name(point_group, this%group_name) 183 call symmetries_finite_get_group_elements(point_group, this%group_elements) 184 end if 185 call symmetries_finite_end() 186 end if 187 SAFE_DEALLOCATE_A(position) 188 SAFE_DEALLOCATE_A(typs) 189 end if 190 191 else 192 193 SAFE_ALLOCATE(position(1:3,1:geo%natoms)) ! transpose!! 194 SAFE_ALLOCATE(typs(1:geo%natoms)) 195 196 do iatom = 1, geo%natoms 197 position(1:3,iatom) = M_ZERO 198 199 if(.not. geo%reduced_coordinates) then 200 ! Transform atomic positions to reduced coordinates 201 position(1:dim4syms,iatom) = matmul(geo%atom(iatom)%x(1:dim4syms),klattice(1:dim4syms,1:dim4syms))/(M_TWO*M_PI) 202 else 203 position(1:dim4syms,iatom) = geo%atom(iatom)%x(1:dim4syms) 204 end if 205 position(1:dim4syms,iatom) = position(1:dim4syms,iatom)- M_HALF 206 do idir = 1, dim4syms 207 position(idir,iatom) = position(idir,iatom) - anint(position(idir,iatom)) 208 end do 209 position(1:dim4syms,iatom) = position(1:dim4syms,iatom) + M_HALF 210 211 typs(iatom) = species_index(geo%atom(iatom)%species) 212 end do 213 214 lattice = M_ZERO 215 !NOTE: Why "inverse matrix" ? (NTD) 216 ! get inverse matrix to extract reduced coordinates for spglib 217 lattice(1:dim, 1:dim) = rlattice(1:dim, 1:dim) 218 ! transpose the lattice vectors for use in spglib as row-major matrix 219 lattice(:,:) = transpose(lattice(:,:)) 220 ! fix things for low-dimensional systems: higher dimension lattice constants set to 1 221 do idir = dim + 1, 3 222 lattice(idir, idir) = M_ONE 223 end do 224 225 this%space_group = spg_get_international(c_symbol, lattice(1,1), & 226 position(1,1), typs(1), geo%natoms, SYMPREC) 227 this%space_group = spg_get_schoenflies(c_schoenflies, lattice(1, 1), & 228 position(1, 1), typs(1), geo%natoms, SYMPREC) 229 230 if(this%space_group == 0) then 231 message(1) = "Symmetry analysis failed in spglib. Disabling symmetries." 232 call messages_warning(1, namespace=namespace) 233 234 do iatom = 1, geo%natoms 235 write(message(1),'(a,i6,a,3f12.6,a,3f12.6)') 'type ', typs(iatom), & 236 ' reduced coords ', position(:, iatom), ' cartesian coords ', geo%atom(iatom)%x(:) 237 call messages_info(1) 238 end do 239 240 call init_identity() 241 POP_SUB(symmetries_init) 242 return 243 end if 244 245 call c_to_f_string(c_symbol, this%symbol) 246 call c_to_f_string(c_schoenflies, this%schoenflies) 247 248 max_size = spg_get_multiplicity(lattice(1, 1), position(1, 1), & 249 typs(1), geo%natoms, SYMPREC) 250 251 ! spglib returns row-major not column-major matrices!!! --DAS 252 SAFE_ALLOCATE(rotation(1:3, 1:3, 1:max_size)) 253 SAFE_ALLOCATE(translation(1:3, 1:max_size)) 254 255 fullnops = spg_get_symmetry(rotation(1, 1, 1), translation(1, 1), & 256 max_size, lattice(1, 1), position(1, 1), typs(1), geo%natoms, SYMPREC) 257 258 do iop = 1, fullnops 259 ! transpose due to array order difference between C and fortran 260 rotation(:,:,iop) = transpose(rotation(:,:,iop)) 261 ! sometimes spglib may return lattice vectors as 'fractional' translations 262 translation(:, iop) = translation(:, iop) - anint(translation(:, iop) + M_HALF*SYMPREC) 263 end do 264 265 ! we need to check that it is not a supercell, as in the QE routine (sgam_at) 266 ! they disable fractional translations if the identity has one, because the sym ops might not form a group. 267 ! spglib may return duplicate operations in this case! 268 is_supercell = (fullnops > 48) 269 found_identity = .false. 270 identity = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), shape(identity)) 271 do iop = 1, fullnops 272 if(all(rotation(1:3, 1:3, iop) == identity(1:3, 1:3))) then 273 found_identity = .true. 274 if(any(abs(translation(1:3, iop)) > TOFLOAT(SYMPREC))) then 275 is_supercell = .true. 276 write(message(1),'(a,3f12.6)') 'Identity has a fractional translation ', translation(1:3, iop) 277 call messages_info(1) 278 end if 279 end if 280 end do 281 if(.not. found_identity) then 282 message(1) = "Symmetries internal error: Identity is missing from symmetry operations." 283 call messages_fatal(1, namespace=namespace) 284 end if 285 286 if(is_supercell) then 287 message(1) = "Disabling fractional translations. System appears to be a supercell." 288 call messages_info(1) 289 end if 290 ! actually, we do not use fractional translations regardless currently 291 292 ! this is a hack to get things working, this variable should be 293 ! eliminated and the direction calculated automatically from the 294 ! perturbations. 295 296 !%Variable SymmetryBreakDir 297 !%Type block 298 !%Section Mesh::Simulation Box 299 !%Description 300 !% This variable specifies a direction in which the symmetry of 301 !% the system will be broken. This is useful for generating <i>k</i>-point 302 !% grids when an external perturbation is applied. 303 !%End 304 305 this%breakdir(1:3) = M_ZERO 306 307 if(parse_block(namespace, 'SymmetryBreakDir', blk) == 0) then 308 309 do idir = 1, dim4syms 310 call parse_block_float(blk, 0, idir - 1, this%breakdir(idir)) 311 end do 312 313 call parse_block_end(blk) 314 315 end if 316 317 SAFE_ALLOCATE(this%ops(1:fullnops)) 318 319 320 ! check all operations and leave those that kept the symmetry-breaking 321 ! direction invariant and (for the moment) that do not have a translation 322 this%nops = 0 323 do iop = 1, fullnops 324 call symm_op_init(tmpop, rotation(1:3, 1:3, iop), rlattice(1:dim4syms,1:dim4syms), & 325 klattice(1:dim4syms,1:dim4syms), dim4syms, & 326 TOFLOAT(translation(1:3, iop))) 327 328 if(symm_op_invariant_cart(tmpop, this%breakdir, TOFLOAT(SYMPREC)) & 329 .and. .not. symm_op_has_translation(tmpop, TOFLOAT(SYMPREC))) then 330 this%nops = this%nops + 1 331 call symm_op_copy(tmpop, this%ops(this%nops)) 332 end if 333 end do 334 335 SAFE_DEALLOCATE_A(position) 336 SAFE_DEALLOCATE_A(typs) 337 SAFE_DEALLOCATE_A(rotation) 338 SAFE_DEALLOCATE_A(translation) 339 340 end if 341 342 call symmetries_write_info(this, namespace, dim, periodic_dim, stdout) 343 344 POP_SUB(symmetries_init) 345 346 contains 347 348 subroutine init_identity() 349 350 PUSH_SUB(symmetries_init.init_identity) 351 352 SAFE_ALLOCATE(this%ops(1:1)) 353 this%nops = 1 354 call symm_op_init(this%ops(1), reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/)), & 355 rlattice, klattice, dim4syms) 356 this%breakdir = M_ZERO 357 this%space_group = 1 358 359 POP_SUB(symmetries_init.init_identity) 360 361 end subroutine init_identity 362 363 subroutine c_to_f_string(c_string, f_string) 364 character(kind=c_char,len=1), intent(in) :: c_string(*) 365 character(len=*), intent(out) :: f_string 366 367 integer :: i 368 369 i = 1 370 do while(c_string(i) /= C_NULL_CHAR .and. i <= len(f_string)) 371 f_string(i:i) = c_string(i) 372 i = i + 1 373 end do 374 375 if (i <= len(f_string)) f_string(i:) = ' ' 376 377 end subroutine c_to_f_string 378 379 end subroutine symmetries_init 380 381 382 383 ! ------------------------------------------------------------------------------- 384 subroutine symmetries_copy(inp, outp) 385 type(symmetries_t), intent(in) :: inp 386 type(symmetries_t), intent(out) :: outp 387 388 integer :: iop 389 390 PUSH_SUB(symmetries_copy) 391 392 outp%nops = inp%nops 393 outp%breakdir = inp%breakdir 394 395 SAFE_ALLOCATE(outp%ops(1:outp%nops)) 396 397 do iop = 1, outp%nops 398 call symm_op_copy(inp%ops(iop), outp%ops(iop)) 399 end do 400 401 outp%group_name = inp%group_name 402 outp%group_elements = inp%group_elements 403 outp%symbol = inp%symbol 404 outp%schoenflies = inp%schoenflies 405 406 POP_SUB(symmetries_copy) 407 end subroutine symmetries_copy 408 409 ! ------------------------------------------------------------------------------- 410 411 subroutine symmetries_end(this) 412 type(symmetries_t), intent(inout) :: this 413 414 PUSH_SUB(symmetries_end) 415 416 SAFE_DEALLOCATE_A(this%ops) 417 418 POP_SUB(symmetries_end) 419 end subroutine symmetries_end 420 421 ! ------------------------------------------------------------------------------- 422 423 integer pure function symmetries_number(this) result(number) 424 type(symmetries_t), intent(in) :: this 425 426 number = this%nops 427 end function symmetries_number 428 429 ! ------------------------------------------------------------------------------- 430 431 subroutine symmetries_apply_kpoint_red(this, iop, aa, bb) 432 type(symmetries_t), intent(in) :: this 433 integer, intent(in) :: iop 434 FLOAT, intent(in) :: aa(1:3) 435 FLOAT, intent(out) :: bb(1:3) 436 437 PUSH_SUB(symmetries_apply_kpoint_red) 438 439 ASSERT(0 < iop .and. iop <= this%nops) 440 441 bb(1:3) = symm_op_apply_transpose_red(this%ops(iop), aa(1:3)) 442 443 POP_SUB(symmetries_apply_kpoint_red) 444 end subroutine symmetries_apply_kpoint_red 445 446 ! ------------------------------------------------------------------------------- 447 448 integer pure function symmetries_space_group_number(this) result(number) 449 type(symmetries_t), intent(in) :: this 450 451 number = this%space_group 452 end function symmetries_space_group_number 453 454 ! ------------------------------------------------------------------------------- 455 456 logical pure function symmetries_have_break_dir(this) result(have) 457 type(symmetries_t), intent(in) :: this 458 459 have = any(abs(this%breakdir(1:3)) > M_EPSILON) 460 end function symmetries_have_break_dir 461 462 ! ------------------------------------------------------------------------------- 463 464 integer pure function symmetries_identity_index(this) result(index) 465 type(symmetries_t), intent(in) :: this 466 467 integer :: iop 468 469 do iop = 1, this%nops 470 if(symm_op_is_identity(this%ops(iop))) then 471 index = iop 472 cycle 473 end if 474 end do 475 476 end function symmetries_identity_index 477 478 ! --------------------------------------------------------- 479 subroutine symmetries_write_info(this, namespace, dim, periodic_dim, iunit) 480 type(symmetries_t), intent(in) :: this 481 type(namespace_t), intent(in) :: namespace 482 integer, intent(in) :: dim, periodic_dim 483 integer, intent(in) :: iunit 484 485 integer :: iop 486 487 PUSH_SUB(symmetries_write_info) 488 489 call messages_print_stress(iunit, 'Symmetries', namespace=namespace) 490 491 if(this%any_non_spherical) then 492 message(1) = "Symmetries are disabled since non-spherically symmetric species may be present." 493 call messages_info(1,iunit = iunit) 494 call messages_print_stress(iunit, namespace=namespace) 495 end if 496 497 if(.not. this%symmetries_compute) then 498 message(1) = "Symmetries have been disabled by SymmetriesCompute = false." 499 call messages_info(1,iunit = iunit) 500 call messages_print_stress(iunit, namespace=namespace) 501 POP_SUB(symmetries_write_info) 502 return 503 end if 504 505 if (periodic_dim == 0) then 506 ! At the moment only the root node has information about symetries of finite systems. 507 if(mpi_grp_is_root(mpi_world)) then 508 if (this%symmetries_compute) then 509 call messages_write('Symmetry elements : '//trim(this%group_elements), new_line = .true.) 510 call messages_write('Symmetry group : '//trim(this%group_name)) 511 call messages_info(iunit = iunit) 512 end if 513 end if 514 else 515 write(message(1),'(a, i4)') 'Space group No. ', this%space_group 516 write(message(2),'(2a)') 'International: ', trim(this%symbol) 517 write(message(3),'(2a)') 'Schoenflies: ', trim(this%schoenflies) 518 call messages_info(3,iunit = iunit) 519 520 write(message(1),'(a7,a31,12x,a33)') 'Index', 'Rotation matrix', 'Fractional translations' 521 call messages_info(1,iunit = iunit) 522 do iop = 1, this%nops 523 ! list all operations and leave those that kept the symmetry-breaking 524 ! direction invariant and (for the moment) that do not have a translation 525 if(dim == 1) & 526 write(message(1),'(i5,1x,a,2x,1(1i4,2x),1f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), & 527 symm_op_translation_vector_red(this%ops(iop)) 528 if(dim == 2) & 529 write(message(1),'(i5,1x,a,2x,2(2i4,2x),2f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), & 530 symm_op_translation_vector_red(this%ops(iop)) 531 if(dim == 3) & 532 write(message(1),'(i5,1x,a,2x,3(3i4,2x),3f12.6)') iop, ':', symm_op_rotation_matrix_red(this%ops(iop)), & 533 symm_op_translation_vector_red(this%ops(iop)) 534 call messages_info(1,iunit = iunit) 535 end do 536 write(message(1), '(a,i5,a)') 'Info: The system has ', this%nops, ' symmetries that can be used.' 537 call messages_info(iunit = iunit) 538 end if 539 call messages_print_stress(iunit, namespace=namespace) 540 541 POP_SUB(symmetries_write_info) 542 end subroutine symmetries_write_info 543 544 545end module symmetries_oct_m 546 547!! Local Variables: 548!! mode: f90 549!! coding: utf-8 550!! End: 551