1! 2! Copyright (C) 2010-2017 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!--------------------------------------------------------------------------------- 9! 10MODULE symm_base 11 ! 12 !! This module contains the variables needed to describe the symmetry properties 13 !! and the routines to find crystal symmetries. 14 ! 15 USE kinds, ONLY : DP 16 USE io_global, ONLY : stdout 17 USE cell_base, ONLY : at, bg 18 ! 19 ! ... these are acceptance criteria 20 ! 21 REAL(DP), PARAMETER :: eps1 = 1.0d-6, eps2 = 1.0d-5 22 ! 23 SAVE 24 ! 25 PRIVATE 26 ! 27 ! ... Exported variables 28 ! 29 PUBLIC :: s, sr, sname, ft, nrot, nsym, nsym_ns, nsym_na, t_rev, & 30 no_t_rev, time_reversal, irt, invs, invsym, d1, d2, d3, & 31 allfrac, nofrac, nosym, nosym_evc, fft_fact, spacegroup 32 ! 33 INTEGER :: s(3,3,48) 34 !! symmetry matrices, in crystal axis 35 INTEGER :: invs(48) 36 !! index of inverse operation: S^{-1}_i=S(invs(i)) 37 INTEGER :: fft_fact(3) = 1 38 !! FFT dimensions must be multiple of fft_fact 39 INTEGER :: nrot 40 !! number of bravais lattice symmetries 41 INTEGER :: spacegroup = 0 42 !! space group index, as read from input 43 INTEGER :: nsym = 1 44 !! total number of crystal symmetries 45 INTEGER :: nsym_ns = 0 46 !! nonsymmorphic (fractional translation) symms 47 INTEGER :: nsym_na = 0 48 !! excluded nonsymmorphic symmetries because 49 !! fract. transl. is noncommensurate with FFT grid 50 REAL(DP) :: ft (3,48) 51 !! fractional translations, in crystal axis 52 REAL(DP) :: sr (3,3,48) 53 !! symmetry matrices, in cartesian axis 54 REAL(DP) :: accep = 1.0d-5 55 !! initial value of the acceptance threshold 56 !! for position comparison by eqvect in checksym 57 ! 58 CHARACTER(LEN=45) :: sname(48) 59 !! name of the symmetries 60 INTEGER :: t_rev(48) = 0 61 !! time reversal flag, for noncolinear magnetism 62 INTEGER, ALLOCATABLE :: irt(:,:) 63 !! symmetric atom for each atom and sym.op. 64 LOGICAL :: time_reversal = .TRUE. 65 !! if .TRUE. the system has time reversal symmetry 66 LOGICAL :: invsym 67 !! if .TRUE. the system has inversion symmetry 68 LOGICAL :: nofrac = .FALSE. 69 !! if .TRUE. fract. translations are not allowed 70 LOGICAL :: allfrac = .FALSE. 71 !! if .TRUE. all fractionary translations allowed, 72 !! even those not commensurate with FFT grid 73 LOGICAL :: nosym = .FALSE. 74 !! if .TRUE. no symmetry is used 75 LOGICAL :: nosym_evc = .FALSE. 76 !! if .TRUE. symmetry is used only to symmetrize 77 !! k points 78 LOGICAL :: no_t_rev = .FALSE. 79 !! if .TRUE. remove the symmetries that 80 !! require time reversal 81 REAL(DP), TARGET :: d1(3,3,48) 82 !! matrix d1 for rotation of spherical harmonics (d1 for l=1, ...) 83 REAL(DP), TARGET :: d2(5,5,48) 84 !! matrix d2 for rotation of spherical harmonics 85 REAL(DP), TARGET :: d3(7,7,48) 86 !! matrix d3 for rotation of spherical harmonics 87 ! 88 ! ... Exported routines 89 ! 90 PUBLIC :: find_sym, inverse_s, copy_sym, checkallsym, & 91 s_axis_to_cart, set_sym, set_sym_bl, check_grid_sym 92 PUBLIC :: find_sym_ifc ! FIXME: should be merged with find_sym 93 PUBLIC :: remove_sym ! FIXME: is this still useful? 94 ! 95CONTAINS 96 ! 97 !----------------------------------------------------------------------- 98 SUBROUTINE inverse_s() 99 !----------------------------------------------------------------------- 100 !! Locate index of \(S^{-1}\). 101 ! 102 IMPLICIT NONE 103 ! 104 INTEGER :: isym, jsym, ss(3,3) 105 LOGICAL :: found 106 ! 107 DO isym = 1, nsym 108 found = .FALSE. 109 DO jsym = 1, nsym 110 ! 111 ss = MATMUL( s(:,:,jsym), s(:,:,isym) ) 112 ! s(:,:,1) is the identity 113 IF (ALL( s(:,:,1) == ss(:,:) )) THEN 114 invs(isym) = jsym 115 found = .TRUE. 116 ENDIF 117 ! 118 ENDDO 119 IF (.NOT.found) CALL errore( 'inverse_s', ' Not a group', 1 ) 120 ENDDO 121 ! 122 END SUBROUTINE inverse_s 123 ! 124 ! 125 !----------------------------------------------------------------------- 126 SUBROUTINE set_sym_bl() 127 !--------------------------------------------------------------------- 128 !! Provides symmetry operations for all bravais lattices. 129 !! Tests the 24 proper rotations for the cubic lattice first, then 130 !! the 8 rotations specific for the hexagonal axis (special axis c), 131 !! then inversion is added. 132 ! 133 USE matrix_inversion 134 ! 135 IMPLICIT NONE 136 ! 137 CHARACTER(LEN=6), EXTERNAL :: int_to_char 138 ! 139 ! sin3 = sin(pi/3), cos3 = cos(pi/3), msin3 = -sin(pi/3), mcos3 = -cos(pi/3) 140 ! 141 REAL(DP), PARAMETER :: sin3 = 0.866025403784438597d0, cos3 = 0.5d0, & 142 msin3 =-0.866025403784438597d0, mcos3 = -0.5d0 143 ! 144 REAL(DP) :: s0(3,3,32), overlap(3,3), rat(3), rot(3,3), value 145 ! s0: the s matrices in cartesian axis 146 ! overlap: inverse overlap matrix between direct lattice 147 ! rat: the rotated of a direct vector ( cartesian ) 148 ! rot: the rotated of a direct vector ( crystal axis ) 149 ! value: component of the s matrix in axis basis 150 INTEGER :: jpol, kpol, mpol, irot, imat(32) 151 ! counters over the polarizations and the rotations 152 ! 153 CHARACTER(LEN=45) :: s0name(64) 154 ! full name of the rotational part of each symmetry operation 155 ! 156 DATA s0/ 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, & 157 -1.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, 1.d0, & 158 -1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, -1.d0, & 159 1.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, -1.d0, & 160 0.d0, 1.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, -1.d0, & 161 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, -1.d0, & 162 0.d0, -1.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0, & 163 0.d0, 1.d0, 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0, & 164 0.d0, 0.d0, 1.d0, 0.d0, -1.d0, 0.d0, 1.d0, 0.d0, 0.d0, & 165 0.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, 0.d0, & 166 0.d0, 0.d0, -1.d0, 0.d0, 1.d0, 0.d0, 1.d0, 0.d0, 0.d0, & 167 0.d0, 0.d0, 1.d0, 0.d0, 1.d0, 0.d0, -1.d0, 0.d0, 0.d0, & 168 -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 1.d0, 0.d0, & 169 -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, & 170 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, 1.d0, 0.d0, & 171 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, -1.d0, 0.d0, & 172 0.d0, 0.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, & 173 0.d0, 0.d0, -1.d0, -1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, & 174 0.d0, 0.d0, -1.d0, 1.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, & 175 0.d0, 0.d0, 1.d0, -1.d0, 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, & 176 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 1.d0, 0.d0, 0.d0, & 177 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, -1.d0, 1.d0, 0.d0, 0.d0, & 178 0.d0, -1.d0, 0.d0, 0.d0, 0.d0, 1.d0, -1.d0, 0.d0, 0.d0, & 179 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, -1.d0, -1.d0, 0.d0, 0.d0, & 180 cos3, sin3, 0.d0, msin3, cos3, 0.d0, 0.d0, 0.d0, 1.d0, & 181 cos3, msin3, 0.d0, sin3, cos3, 0.d0, 0.d0, 0.d0, 1.d0, & 182 mcos3, sin3, 0.d0, msin3, mcos3, 0.d0, 0.d0, 0.d0, 1.d0, & 183 mcos3, msin3, 0.d0, sin3, mcos3, 0.d0, 0.d0, 0.d0, 1.d0, & 184 cos3, msin3, 0.d0, msin3, mcos3, 0.d0, 0.d0, 0.d0, -1.d0, & 185 cos3, sin3, 0.d0, sin3, mcos3, 0.d0, 0.d0, 0.d0, -1.d0, & 186 mcos3, msin3, 0.d0, msin3, cos3, 0.d0, 0.d0, 0.d0, -1.d0, & 187 mcos3, sin3, 0.d0, sin3, cos3, 0.d0, 0.d0, 0.d0, -1.d0 / 188 ! 189 DATA s0name/ 'identity ',& 190 '180 deg rotation - cart. axis [0,0,1] ',& 191 '180 deg rotation - cart. axis [0,1,0] ',& 192 '180 deg rotation - cart. axis [1,0,0] ',& 193 '180 deg rotation - cart. axis [1,1,0] ',& 194 '180 deg rotation - cart. axis [1,-1,0] ',& 195 ' 90 deg rotation - cart. axis [0,0,-1] ',& 196 ' 90 deg rotation - cart. axis [0,0,1] ',& 197 '180 deg rotation - cart. axis [1,0,1] ',& 198 '180 deg rotation - cart. axis [-1,0,1] ',& 199 ' 90 deg rotation - cart. axis [0,1,0] ',& 200 ' 90 deg rotation - cart. axis [0,-1,0] ',& 201 '180 deg rotation - cart. axis [0,1,1] ',& 202 '180 deg rotation - cart. axis [0,1,-1] ',& 203 ' 90 deg rotation - cart. axis [-1,0,0] ',& 204 ' 90 deg rotation - cart. axis [1,0,0] ',& 205 '120 deg rotation - cart. axis [-1,-1,-1] ',& 206 '120 deg rotation - cart. axis [-1,1,1] ',& 207 '120 deg rotation - cart. axis [1,1,-1] ',& 208 '120 deg rotation - cart. axis [1,-1,1] ',& 209 '120 deg rotation - cart. axis [1,1,1] ',& 210 '120 deg rotation - cart. axis [-1,1,-1] ',& 211 '120 deg rotation - cart. axis [1,-1,-1] ',& 212 '120 deg rotation - cart. axis [-1,-1,1] ',& 213 ' 60 deg rotation - cryst. axis [0,0,1] ',& 214 ' 60 deg rotation - cryst. axis [0,0,-1] ',& 215 '120 deg rotation - cryst. axis [0,0,1] ',& 216 '120 deg rotation - cryst. axis [0,0,-1] ',& 217 '180 deg rotation - cryst. axis [1,-1,0] ',& 218 '180 deg rotation - cryst. axis [2,1,0] ',& 219 '180 deg rotation - cryst. axis [0,1,0] ',& 220 '180 deg rotation - cryst. axis [1,1,0] ',& 221 'inversion ',& 222 'inv. 180 deg rotation - cart. axis [0,0,1] ',& 223 'inv. 180 deg rotation - cart. axis [0,1,0] ',& 224 'inv. 180 deg rotation - cart. axis [1,0,0] ',& 225 'inv. 180 deg rotation - cart. axis [1,1,0] ',& 226 'inv. 180 deg rotation - cart. axis [1,-1,0] ',& 227 'inv. 90 deg rotation - cart. axis [0,0,-1] ',& 228 'inv. 90 deg rotation - cart. axis [0,0,1] ',& 229 'inv. 180 deg rotation - cart. axis [1,0,1] ',& 230 'inv. 180 deg rotation - cart. axis [-1,0,1] ',& 231 'inv. 90 deg rotation - cart. axis [0,1,0] ',& 232 'inv. 90 deg rotation - cart. axis [0,-1,0] ',& 233 'inv. 180 deg rotation - cart. axis [0,1,1] ',& 234 'inv. 180 deg rotation - cart. axis [0,1,-1] ',& 235 'inv. 90 deg rotation - cart. axis [-1,0,0] ',& 236 'inv. 90 deg rotation - cart. axis [1,0,0] ',& 237 'inv. 120 deg rotation - cart. axis [-1,-1,-1]',& 238 'inv. 120 deg rotation - cart. axis [-1,1,1] ',& 239 'inv. 120 deg rotation - cart. axis [1,1,-1] ',& 240 'inv. 120 deg rotation - cart. axis [1,-1,1] ',& 241 'inv. 120 deg rotation - cart. axis [1,1,1] ',& 242 'inv. 120 deg rotation - cart. axis [-1,1,-1] ',& 243 'inv. 120 deg rotation - cart. axis [1,-1,-1] ',& 244 'inv. 120 deg rotation - cart. axis [-1,-1,1] ',& 245 'inv. 60 deg rotation - cryst. axis [0,0,1] ',& 246 'inv. 60 deg rotation - cryst. axis [0,0,-1] ',& 247 'inv. 120 deg rotation - cryst. axis [0,0,1] ',& 248 'inv. 120 deg rotation - cryst. axis [0,0,-1] ',& 249 'inv. 180 deg rotation - cryst. axis [1,-1,0] ',& 250 'inv. 180 deg rotation - cryst. axis [2,1,0] ',& 251 'inv. 180 deg rotation - cryst. axis [0,1,0] ',& 252 'inv. 180 deg rotation - cryst. axis [1,1,0] ' / 253 ! 254 ! ... compute the overlap matrix for crystal axis 255 DO jpol = 1, 3 256 DO kpol = 1, 3 257 rot(kpol,jpol) = at(1,kpol)*at(1,jpol) + & 258 at(2,kpol)*at(2,jpol) + & 259 at(3,kpol)*at(3,jpol) 260 ENDDO 261 ENDDO 262 ! 263 ! ... then its inverse (rot is used as work space) 264 CALL invmat( 3, rot, overlap ) 265 ! 266 nrot = 1 267 ! 268 DO irot = 1, 32 269 ! 270 ! ... for each possible symmetry 271 DO jpol = 1, 3 272 DO mpol = 1, 3 273 ! 274 ! ... compute, in cartesian coordinates the rotated vector 275 rat(mpol) = s0(mpol,1,irot)*at(1,jpol) + & 276 s0(mpol,2,irot)*at(2,jpol) + & 277 s0(mpol,3,irot)*at(3,jpol) 278 ENDDO 279 280 DO kpol = 1, 3 281 ! 282 ! ... the rotated vector is projected on the direct lattice 283 rot(kpol,jpol) = at(1,kpol)*rat(1) + & 284 at(2,kpol)*rat(2) + & 285 at(3,kpol)*rat(3) 286 ENDDO 287 ENDDO 288 ! 289 ! ... and the inverse of the overlap matrix is applied 290 DO jpol = 1,3 291 DO kpol = 1,3 292 value = overlap(jpol,1)*rot(1,kpol) + & 293 overlap(jpol,2)*rot(2,kpol) + & 294 overlap(jpol,3)*rot(3,kpol) 295 ! 296 IF ( ABS(DBLE(NINT(value))-value) > eps1 ) THEN 297 ! 298 ! ... if a noninteger is obtained, this implies that this operation 299 ! is not a symmetry operation for the given lattice 300 ! 301 GOTO 10 302 ENDIF 303 ! 304 s(kpol,jpol,nrot) = NINT(value) 305 ENDDO 306 ENDDO 307 ! 308 sname(nrot) = s0name(irot) 309 imat(nrot) = irot 310 nrot = nrot+1 311 ! 312 10 CONTINUE 313 ! 314 ENDDO 315 ! 316 nrot = nrot-1 317 ! 318 IF ( nrot /= 1 .AND. nrot /= 2 .AND. nrot /= 4 .AND. nrot /= 6 .AND. & 319 nrot /= 8 .AND. nrot /=12 .AND. nrot /=24 ) THEN 320 WRITE (stdout, '(80("-"),/,"NOTICE: Bravais lattice has wrong number (",& 321 & i2,") of symmetries - symmetries are disabled",/,80("-"))' ) nrot 322 nrot = 1 323 ENDIF 324 ! 325 ! ... set the inversion symmetry (Bravais lattices have always inversion symmetry) 326 DO irot = 1, nrot 327 sname(irot+nrot) = s0name(imat(irot)+32) 328 DO kpol = 1, 3 329 DO jpol = 1, 3 330 s(kpol,jpol,irot+nrot) = -s(kpol,jpol,irot) 331 ENDDO 332 ENDDO 333 ENDDO 334 ! 335 nrot = 2*nrot 336 ! 337 ! ... reset fractional translations to zero before checking the group 338 ft(:,:) = 0.0_dp 339 IF ( .NOT. is_group(nrot) ) THEN 340 ! ... This happens for instance for an hexagonal lattice with one axis 341 ! oriented at 15 degrees from the x axis, the other along (-1,1,0) 342 CALL infomsg( 'set_sym_bl', 'NOTICE: Symmetry group for Bravais lattice & 343 &is not a group (' // TRIM(int_to_char(nrot)) // & 344 ') - symmetries are disabled' ) 345 nrot = 1 346 ENDIF 347 ! 348 RETURN 349 ! 350 END SUBROUTINE set_sym_bl 351 ! 352 ! 353 !----------------------------------------------------------------------- 354 SUBROUTINE find_sym( nat, tau, ityp, magnetic_sym, m_loc, no_z_inv ) 355 !----------------------------------------------------------------------- 356 !! This routine finds the point group of the crystal, by eliminating 357 !! the symmetries of the Bravais lattice which are not allowed 358 !! by the atomic positions (or by the magnetization if present). 359 ! 360 IMPLICIT NONE 361 ! 362 INTEGER, INTENT(IN) :: nat 363 !! total number of atoms of all species 364 INTEGER, INTENT(IN) :: ityp(nat) 365 !! the type of each i-th atom in stdin 366 REAL(DP), INTENT(IN) :: tau(3,nat) 367 !! atomic positions 368 REAL(DP), INTENT(IN) :: m_loc(3,nat) 369 !! local integrated magnetization 370 LOGICAL, INTENT(IN) :: magnetic_sym 371 !! magnetic_sym = noncolin .AND. domag 372 LOGICAL, INTENT(IN), OPTIONAL :: no_z_inv 373 !! if .TRUE., disable symmetries sending z into -z. 374 !! Some calculations (e.g. gate fields) require this. 375 ! 376 ! ... local variables 377 ! 378 INTEGER :: i 379 LOGICAL :: sym(48) 380 ! if true the corresponding operation is a symmetry operation 381 ! 382 IF ( .NOT. ALLOCATED(irt) ) ALLOCATE( irt(48,nat) ) 383 irt(:,:) = 0 384 ! 385 ! Here we find the true symmetries of the crystal 386 ! 387 symm: DO i = 1, 3 !emine: if it is not resolved in 3 steps it is sth else? 388 IF ( PRESENT(no_z_inv) ) THEN 389 CALL sgam_at( nat, tau, ityp, sym, no_z_inv ) 390 ELSE 391 CALL sgam_at( nat, tau, ityp, sym ) 392 ENDIF 393 ! 394 ! ... Here we check for magnetic symmetries 395 IF ( magnetic_sym ) CALL sgam_at_mag( nat, m_loc, sym ) 396 ! 397 ! ... If nosym_evc is true from now on we do not use the symmetry any more 398 IF (nosym_evc) THEN 399 sym = .FALSE. 400 sym(1) = .TRUE. 401 ENDIF 402 ! 403 ! ... Here we re-order all rotations in such a way that true sym.ops 404 ! are the first nsym; rotations that are not sym.ops. follow 405 nsym = copy_sym( nrot, sym ) 406 ! 407 IF ( .NOT. is_group( nsym ) ) THEN 408 IF (i == 1) CALL infomsg( 'find_sym', & 409 'Not a group! Trying with lower acceptance parameter...' ) 410 accep = accep * 0.5d0 411 IF (i == 3) THEN 412 CALL infomsg( 'find_sym', 'Still not a group! symmetry disabled' ) 413 nsym = 1 414 ENDIF 415 CYCLE symm 416 ELSE 417 IF (i > 1) CALL infomsg( 'find_sym', 'Symmetry operations form a group' ) 418 EXIT symm 419 ENDIF 420 ENDDO symm 421 ! 422 ! ... check if inversion (I) is a symmetry. 423 ! If so, it should be the (nsym/2+1)-th operation of the group 424 ! 425 invsym = ALL( s(:,:,nsym/2+1) == -s(:,:,1) ) 426 ! 427 CALL inverse_s() 428 ! 429 CALL s_axis_to_cart() 430 ! 431 RETURN 432 ! 433 END SUBROUTINE find_sym 434 ! 435 ! 436 !----------------------------------------------------------------------- 437 SUBROUTINE sgam_at( nat, tau, ityp, sym, no_z_inv ) 438 !----------------------------------------------------------------------- 439 !! Given the point group of the Bravais lattice, this routine finds 440 !! the subgroup which is the point group of the considered crystal. 441 !! Non symmorphic groups are allowed, provided that fractional 442 !! translations are allowed (nofrac=.false) and that the unit cell 443 !! is not a supercell. 444 ! 445 !! On output, the array sym is set to .TRUE.. for each operation 446 !! of the original point group that is also a symmetry operation 447 !! of the crystal symmetry point group. 448 ! 449 IMPLICIT NONE 450 ! 451 INTEGER, INTENT(IN) :: nat 452 !! number of atoms in the unit cell 453 INTEGER, INTENT(IN) :: ityp(nat) 454 !! species of each atom in the unit cell 455 REAL(DP), INTENT(IN) :: tau(3,nat) 456 !! cartesian coordinates of the atoms 457 LOGICAL, INTENT(IN), OPTIONAL :: no_z_inv 458 !! if .TRUE., disable symmetry operations sending z into -z. 459 !! Some calculations (e.g. gate fields) require this 460 LOGICAL, INTENT(OUT) :: sym(48) 461 !! flag indicating if sym.op. isym in the parent group 462 !! is a true symmetry operation of the crystal. 463 ! 464 ! ... local variables 465 ! 466 INTEGER :: na, kpol, nb, irot, i, j 467 ! counters 468 REAL(DP) , ALLOCATABLE :: xau(:,:), rau(:,:) 469 ! atomic coordinates in crystal axis 470 LOGICAL :: fractional_translations, no_z 471 INTEGER :: nfrac 472 REAL(DP) :: ft_(3), ftaux(3) 473 ! 474 ALLOCATE( xau(3,nat) ) 475 ALLOCATE( rau(3,nat) ) 476 ! 477 ! ... Compute the coordinates of each atom in the basis of 478 ! the direct lattice vectors 479 DO na = 1, nat 480 xau(:,na) = bg(1,:)*tau(1,na) + bg(2,:)*tau(2,na) + bg(3,:)*tau(3,na) 481 ENDDO 482 ! 483 ! ... check if the identity has fractional translations (this means 484 ! that the cell is actually a supercell). When this happens, fractional 485 ! translations are disabled, because there is no guarantee that the 486 ! generated sym.ops. form a group. 487 ! 488 nb = 1 489 irot = 1 490 ! 491 fractional_translations = .NOT. nofrac 492 ! 493 IF ( fractional_translations ) THEN 494 DO na = 2, nat 495 IF (ityp(nb) == ityp(na)) THEN 496 ! 497 ft_(:) = xau(:,na) - xau(:,nb) - NINT( xau(:,na) - xau(:,nb) ) 498 sym(irot) = checksym ( irot, nat, ityp, xau, xau, ft_ ) 499 IF (sym(irot)) THEN 500 fractional_translations = .FALSE. 501 WRITE( stdout, '(5x,"Found symmetry operation: I + (",& 502 & 3f8.4, ")",/,5x,"This is a supercell,", & 503 & " fractional translations are disabled")') ft_ 504 GOTO 10 505 ENDIF 506 ! 507 ENDIF 508 ENDDO 509 ENDIF 510 ! 511 10 CONTINUE 512 ! 513 nsym_ns = 0 514 fft_fact(:) = 1 515 ! 516 DO irot = 1, nrot 517 ! 518 DO na = 1, nat 519 ! rau = rotated atom coordinates 520 rau(:,na) = s(1,:,irot) * xau(1,na) + & 521 s(2,:,irot) * xau(2,na) + & 522 s(3,:,irot) * xau(3,na) 523 ENDDO 524 ! 525 ! ... first attempt: no fractional translation 526 ft(:,irot) = 0 527 ft_(:) = 0.d0 528 ! 529 sym(irot) = checksym( irot, nat, ityp, xau, rau, ft_ ) 530 ! 531 IF (.NOT.sym(irot) .AND. fractional_translations) THEN 532 nb = 1 533 DO na = 1, nat 534 IF (ityp(nb) == ityp(na)) THEN 535 ! 536 ! ... second attempt: check all possible fractional translations 537 ft_(:) = rau(:,na) - xau(:,nb) - NINT( rau(:,na) - xau(:,nb) ) 538 ! 539 ! ... ft_ is in crystal axis and is a valid fractional translation 540 ! only if ft_(i)=0 or ft_(i)=1/n, with n=2,3,4, 541 ! 542 DO i = 1, 3 543 IF ( ABS (ft_(i)) > eps2 ) THEN 544 ftaux(i) = ABS (1.0_dp/ft_(i) - NINT(1.0_dp/ft_(i)) ) 545 nfrac = NINT(1.0_dp/ABS(ft_(i))) 546 IF ( ftaux(i) < eps2 .AND. nfrac /= 2 .AND. & 547 nfrac /= 3 .AND. nfrac /= 4 .AND. nfrac /= 6 ) & 548 ftaux(i) = 2*eps2 549 ELSE 550 ftaux(i) = 0.0_dp 551 ENDIF 552 ENDDO 553 ! 554 IF ( ANY( ftaux(:) > eps2 ) ) CYCLE 555 ! 556 sym(irot) = checksym( irot, nat, ityp, xau, rau, ft_ ) 557 ! 558 IF ( sym(irot) ) THEN 559 nsym_ns = nsym_ns + 1 560 ft(:,irot) = ft_(:) 561 ! 562 ! ... Find factors that must be present in FFT grid dimensions 563 ! in order to ensure that fractional translations are 564 ! commensurate with FFT grids. 565 DO i = 1, 3 566 IF ( ABS (ft_(i)) > eps2 ) THEN 567 nfrac = NINT(1.0_dp/ABS(ft_(i))) 568 ELSE 569 nfrac = 0 570 END IF 571 fft_fact(i) = mcm(fft_fact(i),nfrac) 572 ENDDO 573 ! 574 GOTO 20 575 ENDIF 576 ENDIF 577 ENDDO 578 ! 579 ENDIF 580 ! 581 20 CONTINUE 582 ! 583 ENDDO 584 ! 585 ! ... disable all symmetries z -> -z 586 IF ( PRESENT(no_z_inv) ) THEN 587 IF ( no_z_inv ) THEN 588 DO irot = 1, nrot 589 IF (s(3,3,irot) == -1) sym(irot) = .FALSE. 590 ENDDO 591 ENDIF 592 ENDIF 593 ! 594 ! ... deallocate work space 595 DEALLOCATE( rau ) 596 DEALLOCATE( xau ) 597 ! 598 RETURN 599 ! 600 END SUBROUTINE sgam_at 601 ! 602 ! 603 !----------------------------------------------------------------------- 604 SUBROUTINE sgam_at_mag( nat, m_loc, sym ) 605 !----------------------------------------------------------------------- 606 !! Find magnetic symmetries, i.e. point-group symmetries that are 607 !! also symmetries of the local magnetization - including 608 !! rotation + time reversal operations. 609 ! 610 IMPLICIT NONE 611 ! 612 INTEGER, INTENT(IN) :: nat 613 !! numbero of atoms of all species 614 REAL(DP), INTENT(IN) :: m_loc(3,nat) 615 !! local magnetization, must be invariant under the sym.op. 616 LOGICAL, INTENT(INOUT) :: sym(48) 617 !! .TRUE. if rotation isym is a sym.op. of the crystal 618 !! (i.e. not of the bravais lattice only) 619 ! 620 ! ... local variables 621 ! 622 INTEGER :: na, nb, irot 623 LOGICAL :: t1, t2 624 REAL(DP) , ALLOCATABLE :: mxau(:,:), mrau(:,:) 625 ! magnetization and rotated magnetization in crystal axis 626 ! 627 ALLOCATE( mxau(3,nat), mrau(3,nat) ) 628 ! 629 ! ... Compute the local magnetization of each atom in the basis of 630 ! the direct lattice vectors 631 DO na = 1, nat 632 mxau(:,na) = bg(1,:) * m_loc(1,na) + & 633 bg(2,:) * m_loc(2,na) + & 634 bg(3,:) * m_loc(3,na) 635 ENDDO 636 ! 637 DO irot = 1, nrot 638 ! 639 t_rev(irot) = 0 640 ! 641 IF ( sym(irot) ) THEN 642 ! 643 ! ... mrau = rotated local magnetization 644 DO na = 1, nat 645 mrau(:,na) = s(1,:,irot) * mxau(1,na) + & 646 s(2,:,irot) * mxau(2,na) + & 647 s(3,:,irot) * mxau(3,na) 648 ENDDO 649 ! 650 IF (sname(irot)(1:3) == 'inv') mrau = -mrau 651 ! 652 ! ... check if this a magnetic symmetry 653 t1 = .TRUE. 654 t2 = .TRUE. 655 ! 656 DO na = 1, nat 657 ! 658 nb = irt(irot,na) 659 IF ( nb < 1 .OR. nb > nat ) CALL errore( 'check_mag_sym', & 660 'internal error: out-of-bound atomic index', na ) 661 ! 662 t1 = ( ABS(mrau(1,na) - mxau(1,nb)) + & 663 ABS(mrau(2,na) - mxau(2,nb)) + & 664 ABS(mrau(3,na) - mxau(3,nb)) < eps2 ) .AND. t1 665 t2 = ( ABS(mrau(1,na) + mxau(1,nb)) + & 666 ABS(mrau(2,na) + mxau(2,nb)) + & 667 ABS(mrau(3,na) + mxau(3,nb)) < eps2 ) .AND. t2 668 ! 669 ENDDO 670 ! 671 IF ( .NOT. t1 .AND. .NOT. t2 ) THEN 672 ! not a magnetic symmetry 673 sym(irot) = .FALSE. 674 ELSEIF( t2 .AND. .NOT. t1 ) THEN 675 ! magnetic symmetry with time reversal, if allowed 676 IF (no_t_rev) THEN 677 sym(irot) = .FALSE. 678 ELSE 679 t_rev(irot) = 1 680 ENDIF 681 ENDIF 682 IF ( (.NOT. sym(irot) ) .AND. & 683 ( ABS(ft(1,irot)) > eps2 .OR. & 684 ABS(ft(2,irot)) > eps2 .OR. & 685 ABS(ft(3,irot)) > eps2 ) ) nsym_ns = nsym_ns-1 686 ! 687 ENDIF 688 ! 689 ENDDO 690 ! 691 ! ... deallocate work space 692 DEALLOCATE( mrau, mxau ) 693 ! 694 RETURN 695 ! 696 END SUBROUTINE sgam_at_mag 697 ! 698 ! 699 !------------------------------------------------------------------------- 700 SUBROUTINE set_sym( nat, tau, ityp, nspin_mag, m_loc ) 701 !----------------------------------------------------------------------- 702 !! This routine receives as input atomic types and positions, if there 703 !! is noncollinear magnetism and the initial magnetic moments 704 !! it sets the symmetry elements of this module. 705 !! Note that \(at\) and \(bg\) are those in \(\textrm{cell_base}\). It sets nrot, nsym, s, 706 !! sname, sr, invs, ft, irt, t_rev, time_reversal, and invsym. 707 ! 708 IMPLICIT NONE 709 ! 710 INTEGER, INTENT(IN) :: nat 711 !! number of atoms in the unit cell 712 INTEGER, INTENT(IN) :: ityp(nat) 713 !! species of each atom in the unit cell 714 INTEGER, INTENT(IN) :: nspin_mag 715 !! =1 when nspin=1,4 (domag=.false.), =2 when 716 !! nspin=2, =4 nspin=4 (domag=.true.) 717 REAL(DP), INTENT(IN) :: tau(3,nat) 718 !! cartesian coordinates of the atoms 719 REAL(DP), INTENT(IN) :: m_loc(3,nat) 720 !! local magnetization, must be invariant under the sym.op. 721 ! 722 time_reversal = (nspin_mag /= 4) 723 t_rev(:) = 0 724 ! 725 CALL set_sym_bl() 726 CALL find_sym( nat, tau, ityp, .NOT.time_reversal, m_loc ) 727 ! 728 RETURN 729 ! 730 END SUBROUTINE set_sym 731 ! 732 ! 733 !----------------------------------------------------------------------- 734 INTEGER FUNCTION copy_sym( nrot_, sym ) 735 !----------------------------------------------------------------------- 736 !! Copy symmetry operations in sequential order so that: 737 ! 738 !! * \(s(i,j,\text{irot})\), with \(\text{irot} \leq \text{nsym}\) are the symmetry 739 !! operations of the crystal; 740 !! * \(s(i,j,\text{irot})\), with \(\text{nsym}+1<\text{irot}\leq \text{nrot}\) are 741 !! the symmetry operations of the lattice. 742 ! 743 !! On exit \(\textrm{copy_sym}\) returns nsym. 744 ! 745 IMPLICIT NONE 746 ! 747 INTEGER, INTENT(IN) :: nrot_ 748 !! number of rotations 749 LOGICAL, INTENT(INOUT) :: sym(48) 750 !! .TRUE. if rotation isym is a sym.op. of the crystal 751 !! (i.e. not of the bravais lattice only) 752 ! 753 ! ... local variables 754 ! 755 INTEGER :: stemp(3,3), ftemp(3), ttemp, irot, jrot 756 REAL(DP) :: ft_(3) 757 INTEGER, ALLOCATABLE :: irtemp(:) 758 CHARACTER(LEN=45) :: nametemp 759 ! 760 ! 761 ALLOCATE ( irtemp(SIZE(irt,2)) ) 762 ! 763 jrot = 0 764 ! 765 DO irot = 1, nrot_ 766 IF ( sym(irot) ) THEN 767 jrot = jrot + 1 768 IF (irot > jrot) THEN 769 stemp = s(:,:,jrot) 770 s(:,:,jrot) = s(:,:,irot) 771 s(:,:,irot) = stemp 772 ft_(:) = ft(:,jrot) 773 ft(:,jrot) = ft(:,irot) 774 ft(:,irot) = ft_(:) 775 irtemp(:) = irt(jrot,:) 776 irt(jrot,:) = irt(irot,:) 777 irt(irot,:) = irtemp (:) 778 nametemp = sname(jrot) 779 sname(jrot) = sname(irot) 780 sname(irot) = nametemp 781 ttemp = t_rev(jrot) 782 t_rev(jrot) = t_rev(irot) 783 t_rev(irot) = ttemp 784 ENDIF 785 ENDIF 786 ENDDO 787 ! 788 sym(1:jrot) = .TRUE. 789 sym(jrot+1:nrot_) = .FALSE. 790 ! 791 DEALLOCATE( irtemp ) 792 ! 793 copy_sym = jrot 794 ! 795 RETURN 796 ! 797 END FUNCTION copy_sym 798 ! 799 ! 800 !----------------------------------------------------------------------- 801 LOGICAL FUNCTION is_group( nsym_ ) 802 !----------------------------------------------------------------------- 803 !! Checks that {S} is a group. 804 ! 805 IMPLICIT NONE 806 ! 807 INTEGER, INTENT(IN) :: nsym_ 808 INTEGER :: isym, jsym, ksym, ss (3,3) 809 REAL(DP) :: st(3), dt(3) 810 LOGICAL :: found 811 ! 812 DO isym = 1, nsym_ 813 DO jsym = 1, nsym_ 814 ! 815 ss = MATMUL(s(:,:,isym), s(:,:,jsym)) 816 st(:) = ft(:,jsym) + s(1,:,jsym)*ft(1,isym) + & 817 s(2,:,jsym)*ft(2,isym) + & 818 s(3,:,jsym)*ft(3,isym) 819 ! 820 ! ... here we check that the input matrices really form a group: 821 ! S(k) = S(i)*S(j) 822 ! ftau_k = S(j)*ftau_i+ftau_j (modulo a lattice vector) 823 ! 824 found = .FALSE. 825 ! 826 DO ksym = 1, nsym_ 827 dt(:) = ft(:,ksym) - st(:) - NINT(ft(:,ksym) - st(:)) 828 IF ( ALL( s(:,:,ksym) == ss(:,:) ) .AND. & 829 ( ABS( dt(1) ) < eps2 ) .AND. & 830 ( ABS( dt(2) ) < eps2 ) .AND. & 831 ( ABS( dt(3) ) < eps2 ) ) THEN 832 IF (found) THEN 833 is_group = .FALSE. 834 RETURN 835 ENDIF 836 found = .TRUE. 837 ENDIF 838 ENDDO 839 ! 840 IF ( .NOT. found ) THEN 841 is_group = .FALSE. 842 RETURN 843 ENDIF 844 ! 845 ENDDO 846 ENDDO 847 ! 848 is_group=.TRUE. 849 ! 850 RETURN 851 ! 852 END FUNCTION is_group 853 ! 854 ! 855 !----------------------------------------------------------------------- 856 LOGICAL FUNCTION checksym( irot, nat, ityp, xau, rau, ft_ ) 857 !----------------------------------------------------------------------- 858 !! This function receives as input all the atomic positions xau, 859 !! and the rotated rau by the symmetry operation ir. It returns 860 !! .TRUE. if, for each atom na, it is possible to find an atom nb 861 !! which is of the same type of na, and coincides with it after the 862 !! symmetry operation. Fractional translations are allowed. 863 ! 864 IMPLICIT NONE 865 ! 866 INTEGER, INTENT(IN) :: nat 867 !! number of atoms 868 INTEGER, INTENT(IN) :: ityp(nat) 869 !! the type of each atom 870 INTEGER, INTENT(IN) :: irot 871 !! rotation index 872 REAL(DP), INTENT(IN) :: xau(3,nat) 873 !! the initial vectors (in crystal coordinates) 874 REAL(DP), INTENT(IN) :: rau(3,nat) 875 !! the rotated vectors (as above) 876 REAL(DP), INTENT(IN) :: ft_(3) 877 !! fractionary translation (as above) 878 ! 879 ! ... local variables 880 ! 881 INTEGER :: na, nb 882 LOGICAL, EXTERNAL :: eqvect 883 ! the testing function 884 ! 885 DO na = 1, nat 886 DO nb = 1, nat 887 ! 888 IF( ityp (nb) == ityp (na) ) THEN 889 checksym = eqvect( rau(1,na), xau(1,nb), ft_ , accep ) 890 IF ( checksym ) THEN 891 ! 892 ! ... the rotated atom does coincide with one of the like atoms 893 ! keep track of which atom the rotated atom coincides with 894 irt (irot, na) = nb 895 GOTO 10 896 ! 897 ENDIF 898 ENDIF 899 ! 900 ENDDO 901 ! 902 ! ... the rotated atom does not coincide with any of the like atoms 903 ! s(ir) + ft is not a symmetry operation 904 checksym = .FALSE. 905 RETURN 906 ! 907 10 CONTINUE 908 ENDDO 909 ! 910 ! ... s(ir) + ft is a symmetry operation 911 ! 912 RETURN 913 ! 914 END FUNCTION checksym 915 ! 916 ! 917 !----------------------------------------------------------------------- 918 SUBROUTINE checkallsym( nat, tau, ityp ) 919 !----------------------------------------------------------------------- 920 !! Given a crystal group this routine checks that the actual atomic 921 !! positions and bravais lattice vectors are compatible with it. 922 !! Used in relaxation/MD runs to check that atomic motion is 923 !! consistent with assumed symmetry. 924 ! 925 IMPLICIT NONE 926 ! 927 INTEGER, INTENT(IN) :: nat 928 !! number of atoms 929 INTEGER, INTENT(IN) :: ityp(nat) 930 !! the type of each atom 931 REAL(DP), INTENT(IN) :: tau(3,nat) 932 !! postion of each atom 933 ! 934 ! ... local variables 935 ! 936 INTEGER :: na, kpol, isym, i, j, k, l 937 LOGICAL :: loksym(48) 938 REAL(DP) :: sx(3,3), sy(3,3) 939 REAL(DP), ALLOCATABLE :: xau(:,:), rau(:,:) 940 ! 941 ALLOCATE( xau(3,nat) ) 942 ALLOCATE( rau(3,nat) ) 943 ! 944 ! ... check that s(i,j, isym) is an orthogonal operation 945 DO isym = 1, nsym 946 sx = DBLE( s(:,:,isym) ) 947 sy = MATMUL( bg, sx ) 948 sx = MATMUL( sy, TRANSPOSE(at) ) 949 ! sx is s in cartesian axis 950 sy = MATMUL( TRANSPOSE(sx), sx ) 951 ! sy = s*TRANSPOSE(s) = I 952 DO i = 1, 3 953 sy(i,i) = sy(i,i) - 1.0_dp 954 ENDDO 955 IF ( ANY(ABS(sy) > eps1) ) & 956 CALL errore( 'checkallsym', 'not orthogonal operation', isym ) 957 ENDDO 958 ! 959 ! ... Compute the coordinates of each atom in the basis of the lattice 960 DO na = 1, nat 961 DO kpol = 1, 3 962 xau(kpol,na) = bg(1,kpol) * tau(1,na) + & 963 bg(2,kpol) * tau(2,na) + & 964 bg(3,kpol) * tau(3,na) 965 ENDDO 966 ENDDO 967 ! 968 ! ... Generate the coordinates of the rotated atoms 969 DO isym = 1, nsym 970 DO na = 1, nat 971 DO kpol = 1, 3 972 rau(kpol,na) = s(1,kpol,isym) * xau(1,na) + & 973 s(2,kpol,isym) * xau(2,na) + & 974 s(3,kpol,isym) * xau(3,na) 975 ENDDO 976 ENDDO 977 ! 978 loksym(isym) = checksym( isym, nat, ityp, xau, rau, ft(1,isym) ) 979 ENDDO 980 ! 981 ! ... deallocate work space 982 ! 983 DEALLOCATE( rau ) 984 DEALLOCATE( xau ) 985 ! 986 DO isym = 1,nsym 987 IF ( .NOT.loksym(isym) ) CALL errore( 'checkallsym', & 988 'the following symmetry operation is not satisfied ', -isym ) 989 ENDDO 990 ! 991 IF (ANY(.NOT.loksym(1:nsym) ) ) THEN 992 !call symmetrize_at (nsym, s, invs, ft, irt, nat, tau, at, bg, & 993 ! alat, omega) 994 CALL errore( 'checkallsym', & 995 'some of the original symmetry operations not satisfied ',1 ) 996 ENDIF 997 ! 998 RETURN 999 ! 1000 END SUBROUTINE checkallsym 1001 ! 1002 ! 1003 !---------------------------------------------------------------------- 1004 SUBROUTINE s_axis_to_cart() 1005 !---------------------------------------------------------------------- 1006 !! This routine transforms symmetry matrices expressed in the 1007 !! basis of the crystal axis into rotations in cartesian axis. 1008 ! 1009 IMPLICIT NONE 1010 ! 1011 INTEGER :: isym 1012 REAL(DP) :: sa(3,3), sb(3,3) 1013 ! 1014 DO isym = 1,nsym 1015 sa(:,:) = DBLE( s(:,:,isym) ) 1016 sb = MATMUL( bg, sa ) 1017 sr(:,:,isym) = MATMUL( at, TRANSPOSE(sb) ) 1018 ENDDO 1019 ! 1020 END SUBROUTINE s_axis_to_cart 1021 ! 1022 ! 1023 !----------------------------------------------------------------------- 1024 SUBROUTINE find_sym_ifc( nat, tau, ityp ) 1025 !----------------------------------------------------------------------- 1026 !! This routine finds the point group of the crystal, by eliminating 1027 !! the symmetries of the Bravais lattice which are not allowed 1028 !! by the atomic positions (for use in the FD package). 1029 ! 1030 IMPLICIT NONE 1031 ! 1032 INTEGER, INTENT(IN) :: nat 1033 !! number of atoms 1034 INTEGER, INTENT(IN) :: ityp(nat) 1035 !! the type of each atom 1036 REAL(DP), INTENT(IN) :: tau(3,nat) 1037 !! postion of each atom 1038 ! 1039 ! ... local variables 1040 ! 1041 INTEGER :: i 1042 LOGICAL :: sym(48) 1043 ! if true the corresponding operation is a symmetry operation 1044 ! 1045 IF ( .NOT. ALLOCATED(irt) ) ALLOCATE( irt(48,nat) ) 1046 irt(:,:) = 0 1047 ! 1048 ! ... Here we find the true symmetries of the crystal 1049 CALL sgam_at_ifc( nat, tau, ityp, sym ) 1050 ! 1051 ! ... Here we re-order all rotations in such a way that true sym.ops 1052 ! are the first nsym; rotations that are not sym.ops. follow 1053 nsym = copy_sym( nrot, sym ) 1054 ! 1055 ! ... check if inversion (I) is a symmetry. 1056 ! If so, it should be the (nsym/2+1)-th operation of the group 1057 invsym = ALL( s(:,:,nsym/2+1) == -s(:,:,1) ) 1058 ! 1059 CALL inverse_s() 1060 ! 1061 CALL s_axis_to_cart() 1062 ! 1063 RETURN 1064 ! 1065 END SUBROUTINE find_sym_ifc 1066 ! 1067 ! 1068 !----------------------------------------------------------------------- 1069 SUBROUTINE sgam_at_ifc( nat, tau, ityp, sym ) 1070 !----------------------------------------------------------------------- 1071 !! Given the point group of the Bravais lattice, this routine finds 1072 !! the subgroup which is the point group of the considered crystal. 1073 !! Non symmorphic groups are allowed, provided that fractional 1074 !! translations are allowed (nofrac=.false), that the unit cell is 1075 !! not a supercell. 1076 ! 1077 !! On output, the array sym is set to .TRUE.. for each operation 1078 !! of the original point group that is also a symmetry operation 1079 !! of the crystal symmetry point group. 1080 ! 1081 IMPLICIT NONE 1082 ! 1083 INTEGER, INTENT(IN) :: nat 1084 !! number of atoms in the unit cell 1085 INTEGER, INTENT(IN) :: ityp(nat) 1086 !! species of each atom in the unit cell 1087 REAL(DP), INTENT(IN) :: tau(3,nat) 1088 !! cartesian coordinates of the atoms 1089 LOGICAL, INTENT(OUT) :: sym(48) 1090 !! flag indicating if sym.op. isym in the parent group 1091 !! is a true symmetry operation of the crystal 1092 ! 1093 ! ... local variables 1094 ! 1095 INTEGER :: na, kpol, nb, irot, i, j 1096 ! counters 1097 REAL(DP) , ALLOCATABLE :: xau(:,:), rau(:,:) 1098 ! atomic coordinates in crystal axis 1099 LOGICAL :: fractional_translations 1100 REAL(DP) :: ft_(3) 1101 ! 1102 ALLOCATE( xau(3,nat) ) 1103 ALLOCATE( rau(3,nat) ) 1104 ! 1105 ! ... Compute the coordinates of each atom in the basis of 1106 ! the direct lattice vectors. 1107 ! 1108 DO na = 1, nat 1109 xau(:,na) = bg(1,:)*tau(1,na) + bg(2,:)*tau(2,na) + bg(3,:)*tau(3,na) 1110 ENDDO 1111 ! 1112 ! ... check if the identity has fractional translations 1113 ! (this means that the cell is actually a supercell). 1114 ! When this happens, fractional translations are disabled, 1115 ! because there is no guarantee that the generated sym.ops. 1116 ! form a group. 1117 ! 1118 nb = 1 1119 irot = 1 1120 ! 1121 fractional_translations = .NOT. nofrac 1122 ! 1123 DO na = 2, nat 1124 IF ( fractional_translations ) THEN 1125 IF (ityp(nb) == ityp(na) ) THEN 1126 ft_(:) = xau(:,na) - xau(:,nb) - NINT( xau(:,na) - xau(:,nb) ) 1127 ! 1128 sym(irot) = checksym( irot, nat, ityp, xau, xau, ft_ ) 1129 ! 1130 IF ( sym (irot) .AND. & 1131 (ABS(ft_(1)**2 + ft_(2)**2 + ft_(3)**2) < 1.d-8) ) & 1132 CALL errore( 'sgam_at_ifc', 'overlapping atoms', na ) 1133 ENDIF 1134 ENDIF 1135 ENDDO 1136 ! 1137 nsym_ns = 0 1138 ! 1139 DO irot = 1, nrot 1140 DO na = 1, nat 1141 ! rau = rotated atom coordinates 1142 rau(:,na) = s(1,:,irot) * xau(1,na) + & 1143 s(2,:,irot) * xau(2,na) + & 1144 s(3,:,irot) * xau(3,na) 1145 ENDDO 1146 ! 1147 ! ... first attempt: no fractional translation 1148 ft(:,irot) = 0 1149 ft_(:) = 0.d0 1150 ! 1151 sym(irot) = checksym( irot, nat, ityp, xau, rau, ft_ ) 1152 ! 1153 IF (.NOT.sym(irot) .AND. fractional_translations) THEN 1154 nb = 1 1155 ! 1156 DO na = 1, nat 1157 IF (ityp(nb) == ityp(na) ) THEN 1158 ! 1159 ! second attempt: check all possible fractional translations 1160 ! 1161 ft_(:) = rau(:,na) - xau(:,nb) - NINT( rau(:,na) - xau(:,nb) ) 1162 ! 1163 sym(irot) = checksym( irot, nat, ityp, xau, rau, ft_ ) 1164 ! 1165 IF (sym(irot) ) THEN 1166 nsym_ns = nsym_ns + 1 1167 ft(:,irot) = ft_(:) 1168 GOTO 100 1169 ENDIF 1170 ENDIF 1171 ENDDO 1172 ! 1173 ENDIF 1174 ! 1175 100 CONTINUE 1176 ! 1177 ENDDO 1178 ! 1179 DEALLOCATE( rau ) 1180 DEALLOCATE( xau ) 1181 ! 1182 RETURN 1183 ! 1184 END SUBROUTINE sgam_at_ifc 1185 ! 1186 !----------------------------------------------------------------------- 1187 FUNCTION check_grid_sym ( nr1, nr2, nr3 ) RESULT ( compatible ) 1188 !--------------------------------------------------------------------- 1189 !! Check that symmetry operations and FFT grid are compatible 1190 !! Needed to prevent trouble with real-space symmetrization 1191 ! 1192 IMPLICIT NONE 1193 ! 1194 INTEGER, INTENT(IN) :: nr1, nr2, nr3 1195 LOGICAL :: compatible, bad 1196 INTEGER :: isym,i,j 1197 ! 1198 compatible = .true. 1199 DO isym = 1, nsym 1200 ! 1201 bad = ( MOD( s(2,1,isym)*nr1, nr2) /= 0 .OR. & 1202 MOD( s(3,1,isym)*nr1, nr3) /= 0 .OR. & 1203 MOD( s(1,2,isym)*nr2, nr1) /= 0 .OR. & 1204 MOD( s(3,2,isym)*nr2, nr3) /= 0 .OR. & 1205 MOD( s(1,3,isym)*nr3, nr1) /= 0 .OR. & 1206 MOD( s(2,3,isym)*nr3, nr2) /= 0 ) 1207 IF ( bad ) THEN 1208 WRITE( stdout, '(5x,"warning: symmetry operation # ",i2, & 1209 & " not compatible with FFT grid. ")') isym 1210 WRITE( stdout, '(3i4)') ( (s(i,j,isym), j=1,3), i=1,3 ) 1211 compatible = .false. 1212 ENDIF 1213 ! 1214 ENDDO 1215 ! 1216 END FUNCTION check_grid_sym 1217 ! 1218 !----------------------------------------------------------------------- 1219 SUBROUTINE remove_sym( nr1, nr2, nr3 ) 1220 !--------------------------------------------------------------------- 1221 !! Compute ftau used for symmetrization in real space (phonon, exx) 1222 !! ensure that they are commensurated with the FFT grid. 1223 ! 1224 IMPLICIT NONE 1225 ! 1226 INTEGER, INTENT(IN) :: nr1, nr2, nr3 1227 ! 1228 ! ... local variables 1229 ! 1230 LOGICAL :: sym(48) 1231 INTEGER :: isym, nsym_, i, j 1232 REAL(dp) :: ftaux(3) 1233 ! 1234 nsym_ = nsym 1235 sym(1:nsym_) = .TRUE. 1236 nsym_na = 0 1237 ! 1238 DO isym = 1, nsym_ 1239 ! 1240 ! check that the grid is compatible with the S rotation 1241 ! 1242 IF ( MOD( s(2,1,isym)*nr1, nr2) /= 0 .OR. & 1243 MOD( s(3,1,isym)*nr1, nr3) /= 0 .OR. & 1244 MOD( s(1,2,isym)*nr2, nr1) /= 0 .OR. & 1245 MOD( s(3,2,isym)*nr2, nr3) /= 0 .OR. & 1246 MOD( s(1,3,isym)*nr3, nr1) /= 0 .OR. & 1247 MOD( s(2,3,isym)*nr3, nr2) /= 0 ) THEN 1248 sym(isym) = .FALSE. 1249 WRITE( stdout, '(5x,"warning: symmetry operation # ",i2, & 1250 & " not compatible with FFT grid. ")') isym 1251 WRITE( stdout, '(3i4)') ( (s(i,j,isym), j=1,3), i=1,3 ) 1252 sym(isym) = .FALSE. 1253 IF ( ABS(ft(1,isym)) > eps2 .OR. & 1254 ABS(ft(2,isym)) > eps2 .OR. & 1255 ABS(ft(3,isym)) > eps2 ) nsym_ns = nsym_ns-1 1256 ENDIF 1257 ! 1258 ! convert ft to FFT coordinates, check if compatible with FFT grid 1259 ! for real-space symmetrization 1260 ! 1261 ftaux(1) = ft(1,isym) * nr1 1262 ftaux(2) = ft(2,isym) * nr2 1263 ftaux(3) = ft(3,isym) * nr3 1264 ! check if the fractional translations are commensurate 1265 ! with the FFT grid, discard sym.op. if not 1266 ! (needed because ph.x symmetrizes in real space) 1267 IF ( ABS(ftaux(1) - NINT(ftaux(1)) ) / nr1 > eps2 .OR. & 1268 ABS(ftaux(2) - NINT(ftaux(2)) ) / nr2 > eps2 .OR. & 1269 ABS(ftaux(3) - NINT(ftaux(3)) ) / nr3 > eps2 ) THEN 1270 ! WRITE( stdout, '(5x,"warning: symmetry operation", & 1271 ! & " # ",i2," not allowed. fractional ", & 1272 ! & "translation:"/5x,3f11.7," in crystal", & 1273 ! & " coordinates")') isym, ft_ 1274 sym(isym) = .FALSE. 1275 nsym_na = nsym_na + 1 1276 nsym_ns = nsym_ns - 1 1277 ENDIF 1278 ! 1279 ENDDO 1280 ! 1281 ! ... count symmetries, reorder them exactly as in "find_sym" 1282 ! 1283 nsym = copy_sym( nsym_, sym ) 1284 invsym = ALL( s(:,:,nsym/2+1) == -s(:,:,1) ) 1285 ! 1286 CALL inverse_s() 1287 CALL s_axis_to_cart() 1288 ! 1289 END SUBROUTINE remove_sym 1290 ! 1291 ! 1292 !-------------------------------------------------------------------------- 1293 INTEGER FUNCTION mcm( i, j ) 1294 !------------------------------------------------------------------------ 1295 !! Returns minimum common multiple of two integers 1296 !! if i=0, returns j, and vice versa; if i<0 or j<0, returns -1. 1297 ! 1298 INTEGER, INTENT(IN) :: i,j 1299 INTEGER :: n1,n2,k 1300 ! 1301 IF (i < 0 .OR. j < 0) THEN 1302 mcm = -1 1303 ELSEIF (i == 0 .AND. j == 0) THEN 1304 mcm = 0 1305 ELSE 1306 n1 = MIN(i,j) 1307 n2 = MAX(i,j) 1308 DO k = 1, n1 1309 mcm = k*n2 1310 IF (MOD(mcm,n1) == 0) RETURN 1311 ENDDO 1312 mcm = n2 1313 ENDIF 1314 ! 1315 END FUNCTION mcm 1316 ! 1317 ! 1318END MODULE symm_base 1319