TABLE OF CONTENTS
ABINIT/symlatt [ Functions ]
NAME
symlatt
FUNCTION
From the unit cell vectors (rprimd) and the corresponding metric tensor, find the Bravais lattice and its symmetry operations (ptsymrel). 1) Find the shortest possible primitive vectors for the lattice 2) Determines the holohedral group of the lattice, and the axes to be used for the conventional cell (this is a delicate part, in which the centering of the reduced cell must be taken into account) The idea is to determine the basis vectors of the conventional cell from the reduced cell basis vectors. 3) Generate the symmetry operations of the holohedral group
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
msym=default maximal number of symmetries rprimd(3,3)=dimensional primitive translations for real space (bohr) tolsym=tolerance for the symmetries
OUTPUT
bravais(11): bravais(1)=iholohedry
bravais(2)=center
bravais(3:11)=coordinates of rprim in the axes
of the conventional bravais lattice (*2 if center/=0)
nptsym=number of point symmetries of the Bravais lattice
ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
of the Bravais lattice in real space in terms
of primitive translations.
NOTES
WARNING: bravais(1) might be given a negative value in another routine, if the cell is non-primitive. The holohedral groups are numbered as follows (see international tables for crystallography (1983), p. 13) iholohedry=1 triclinic 1bar iholohedry=2 monoclinic 2/m iholohedry=3 orthorhombic mmm iholohedry=4 tetragonal 4/mmm iholohedry=5 trigonal 3bar m iholohedry=6 hexagonal 6/mmm iholohedry=7 cubic m3bar m Centering center=0 no centering center=-1 body-centered center=-3 face-centered center=1 A-face centered center=2 B-face centered center=3 C-face centered
PARENTS
ingeo,inqpt,m_ab7_symmetry,m_effective_potential_file,m_tdep_sym m_use_ga,symanal,symbrav,thmeig
CHILDREN
holocell,matr3inv,smallprim,symrelrot,wrtout
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 76 subroutine symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym) 77 78 use defs_basis 79 use m_errors 80 use m_profiling_abi 81 82 !This section has been created automatically by the script Abilint (TD). 83 !Do not modify the following lines by hand. 84 #undef ABI_FUNC 85 #define ABI_FUNC 'symlatt' 86 use interfaces_14_hidewrite 87 use interfaces_32_util 88 use interfaces_41_geometry, except_this_one => symlatt 89 !End of the abilint section 90 91 implicit none 92 93 !Arguments ------------------------------------ 94 !scalars 95 integer,intent(in) :: msym 96 integer,intent(out) :: nptsym 97 real(dp),intent(in) :: tolsym 98 !arrays 99 integer,intent(out) :: bravais(11),ptsymrel(3,3,msym) 100 real(dp),intent(in) :: rprimd(3,3) 101 102 !Local variables------------------------------- 103 !scalars 104 integer,parameter :: mgen=4 105 integer :: center,fact,found,foundc,ia,ib,icase,igen,iholohedry,ii,index,isym 106 integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign 107 real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca 108 real(dp) :: scalarprod,scb,trace,val 109 character(len=500) :: message 110 !arrays 111 integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/) 112 integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3) 113 integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3) 114 integer :: order(mgen) 115 real(dp) :: axes(3,3),axesinvt(3,3),cell_base(3,3),coord(3,3),metmin(3,3) 116 real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3) 117 118 !************************************************************************** 119 120 identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 121 nvecta(1)=2 ; nvectb(1)=3 122 nvecta(2)=1 ; nvectb(2)=3 123 nvecta(3)=1 ; nvectb(3)=2 124 125 !-------------------------------------------------------------------------- 126 !Reduce the input vectors to a set of minimal vectors 127 call smallprim(metmin,minim,rprimd) 128 129 !DEBUG 130 !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1) 131 !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2) 132 !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3) 133 !ENDDEBUG 134 135 !-------------------------------------------------------------------------- 136 !Examine the angles and vector lengths 137 ang90(:)=0 138 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1 139 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1 140 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1 141 equal(:)=0 142 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1 143 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1 144 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1 145 146 !DEBUG 147 !write(std_out,*)' ang90=',ang90(:) 148 !write(std_out,*)' equal=',equal(:) 149 !ENDDEBUG 150 151 !----------------------------------------------------------------------- 152 !Identification of the centering 153 154 foundc=0 155 !Default values 156 fact=1 ; center=0 157 cell_base(:,:)=minim(:,:) 158 159 !Examine each holohedral group 160 !This search is ordered : should not be happy with tetragonal, 161 !while there is FCC ... 162 do index=1,6 163 164 ! If the holohedry is already found, exit 165 if(foundc==1)exit 166 167 ! Initialize the target holohedry 168 iholohedry=list_holo(index) 169 170 ! DEBUG 171 ! write(std_out,*)' symlatt : trial holohedry',iholohedry 172 ! ENDDEBUG 173 174 orthogonal=0 175 if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1 176 177 ! Now, will examine different working hypothesis. 178 ! The set of these hypothesis is thought to cover all possible cases ... 179 180 ! Working hypothesis : the basis is orthogonal 181 if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then 182 fact=1 ; center=0 183 cell_base(:,:)=minim(:,:) 184 ! Checks that the basis vectors are OK for the target holohedry 185 call holocell(cell_base,0,foundc,iholohedry,tolsym) 186 end if 187 188 ! Select one trial direction 189 do itrial=1,3 190 191 ! If the holohedry is already found, exit 192 if(foundc==1)exit 193 194 ia=nvecta(itrial) ; ib=nvectb(itrial) 195 196 ! This is in case of hexagonal holohedry 197 if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then 198 reduceda=metmin(ib,ia)/metmin(ia,ia) 199 fact=1 ; center=0 200 if(abs(reduceda+0.5d0)<tolsym)then 201 cell_base(:,1)=minim(:,ia) 202 cell_base(:,2)=minim(:,ib) 203 cell_base(:,3)=minim(:,itrial) 204 ! Checks that the basis vectors are OK for the target holohedry 205 call holocell(cell_base,0,foundc,iholohedry,tolsym) 206 else if(abs(reduceda-0.5d0)<tolsym)then 207 cell_base(:,1)=minim(:,ia) 208 cell_base(:,2)=minim(:,ib)-minim(:,ia) 209 cell_base(:,3)=minim(:,itrial) 210 ! Checks that the basis vectors are OK for the target holohedry 211 call holocell(cell_base,0,foundc,iholohedry,tolsym) 212 end if 213 end if 214 215 ! Working hypothesis : the conventional cell is orthogonal, 216 ! and the two other vectors are axes of the conventional cell 217 if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then 218 219 ! Compute the reduced coordinate of trial vector in the basis 220 ! of the two other vectors 221 reduceda=metmin(itrial,ia)/metmin(ia,ia) 222 reducedb=metmin(itrial,ib)/metmin(ib,ib) 223 cell_base(:,ia)=minim(:,ia) 224 cell_base(:,ib)=minim(:,ib) 225 if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. & 226 & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then 227 if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib 228 if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia 229 fact=2 230 cell_base(:,itrial)= & 231 & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0 232 call holocell(cell_base,0,foundc,iholohedry,tolsym) 233 else if( abs(abs(reduceda)-0.5d0)<tolsym .and.& 234 & abs(abs(reducedb)-0.5d0)<tolsym ) then 235 fact=2 ; center=-1 236 cell_base(:,itrial)= & 237 & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0 238 call holocell(cell_base,0,foundc,iholohedry,tolsym) 239 end if 240 end if 241 242 ! Working hypothesis : the conventional cell is orthogonal, and 243 ! the trial vector is one of the future axes, and the face perpendicular to it is centered 244 if(foundc==0 .and. iholohedry==3 .and. & 245 & ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then 246 fact=2 ; center=itrial 247 cell_base(:,ia)=minim(:,ia)+minim(:,ib) 248 cell_base(:,ib)=minim(:,ia)-minim(:,ib) 249 cell_base(:,itrial)=minim(:,itrial) 250 ! Checks that the basis vectors are OK for the target holohedry 251 call holocell(cell_base,0,foundc,iholohedry,tolsym) 252 end if 253 254 ! DEBUG 255 ! write(std_out,*)' after test_b, foundc=',foundc 256 ! ENDDEBUG 257 258 ! Working hypothesis : the conventional cell is orthogonal, and 259 ! the trial vector is one of the future axes 260 if(foundc==0 .and. orthogonal==1)then 261 ! Compute the projection of the two other vectors on the trial vector 262 reduceda=metmin(itrial,ia)/metmin(itrial,itrial) 263 reducedb=metmin(itrial,ib)/metmin(itrial,itrial) 264 ! If both projections are half-integer, one might have found an axis 265 if( abs(abs(reduceda)-0.5d0)<tolsym .and.& 266 & abs(abs(reducedb)-0.5d0)<tolsym ) then 267 vecta(:)=minim(:,ia)-reduceda*minim(:,itrial) 268 vectb(:)=minim(:,ib)-reducedb*minim(:,itrial) 269 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 270 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 271 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 272 ! Note the order of selection : body-centered is prefered 273 ! over face centered, which is correct for the tetragonal case 274 if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then 275 ! The lattice is body centered 276 fact=2 ; center=-1 277 cell_base(:,ia)=vecta(:)+vectb(:) 278 cell_base(:,ib)=vecta(:)-vectb(:) 279 cell_base(:,itrial)=minim(:,itrial) 280 call holocell(cell_base,0,foundc,iholohedry,tolsym) 281 else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then 282 ! The lattice is face centered 283 fact=2 ; center=-3 284 cell_base(:,ia)=2.0d0*vecta(:) 285 cell_base(:,ib)=2.0d0*vectb(:) 286 cell_base(:,itrial)=minim(:,itrial) 287 call holocell(cell_base,0,foundc,iholohedry,tolsym) 288 end if 289 end if 290 end if 291 292 ! DEBUG 293 ! write(std_out,*)' after test_c, foundc=',foundc 294 ! ENDDEBUG 295 296 ! Working hypothesis : the conventional cell is orthogonal, 297 ! and body centered with no basis vector being an axis, 298 ! in which case the basis vectors must be equal (even for orthorhombic) 299 if(foundc==0 .and. orthogonal==1 .and. & 300 & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then 301 ! Compute the combination of the two other vectors 302 vecta(:)=minim(:,ia)+minim(:,ib) 303 vectb(:)=minim(:,ia)-minim(:,ib) 304 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 305 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 306 ! Project the trial vector on the first of the two vectors 307 reduceda=( minim(1,itrial)*vecta(1)+ & 308 & minim(2,itrial)*vecta(2)+ & 309 & minim(3,itrial)*vecta(3) )/norm2a 310 reducedb=( minim(1,itrial)*vectb(1)+ & 311 & minim(2,itrial)*vectb(2)+ & 312 & minim(3,itrial)*vectb(3) )/norm2b 313 if( abs(abs(reduceda)-0.5d0)<tolsym )then 314 ! The first vector is an axis 315 fact=2 ; center=-1 316 cell_base(:,ia)=vecta(:) 317 vecta(:)=minim(:,itrial)-reduceda*vecta(:) 318 vectb(:)=0.5d0*vectb(:) 319 cell_base(:,ib)=vecta(:)+vectb(:) 320 cell_base(:,itrial)=vecta(:)-vectb(:) 321 call holocell(cell_base,0,foundc,iholohedry,tolsym) 322 else if( abs(abs(reducedb)-0.5d0)<tolsym )then 323 ! The second vector is an axis 324 fact=2 ; center=-1 325 cell_base(:,ib)=vectb(:) 326 vectb(:)=minim(:,itrial)-reducedb*vectb(:) 327 vecta(:)=0.5d0*vecta(:) 328 cell_base(:,ia)=vectb(:)+vecta(:) 329 cell_base(:,itrial)=vectb(:)-vecta(:) 330 call holocell(cell_base,0,foundc,iholohedry,tolsym) 331 end if 332 end if 333 334 ! Working hypothesis : the conventional cell is orthogonal, 335 ! and face centered, in the case where two minimal vectors are equal 336 if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then 337 ! Compute the combination of these two vectors 338 vecta(:)=minim(:,ia)+minim(:,ib) 339 vectb(:)=minim(:,ia)-minim(:,ib) 340 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 341 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 342 ! Project the trial vector on the two vectors 343 reduceda=( minim(1,itrial)*vecta(1)+ & 344 & minim(2,itrial)*vecta(2)+ & 345 & minim(3,itrial)*vecta(3) )/norm2a 346 reducedb=( minim(1,itrial)*vectb(1)+ & 347 & minim(2,itrial)*vectb(2)+ & 348 & minim(3,itrial)*vectb(3) )/norm2b 349 if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. & 350 & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then 351 fact=2 ; center=-3 352 cell_base(:,itrial)= & 353 & (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0 354 cell_base(:,ia)=vecta(:) 355 cell_base(:,ib)=vectb(:) 356 call holocell(cell_base,0,foundc,iholohedry,tolsym) 357 end if 358 end if 359 360 ! Working hypothesis : the conventional cell is orthogonal, 361 ! face centered, but no two vectors are on the same "square" 362 if(foundc==0 .and. orthogonal==1)then 363 ! Compute the combination of these two vectors 364 vecta(:)=minim(:,ia)+minim(:,ib) 365 vectb(:)=minim(:,ia)-minim(:,ib) 366 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 367 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 368 ! The trial vector length must be equal to one of these lengths 369 if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then 370 fact=2 ; center=-3 371 cell_base(:,ia)=vecta(:)+minim(:,itrial) 372 cell_base(:,ib)=vecta(:)-minim(:,itrial) 373 ! Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib) 374 norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2 375 norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2 376 reduceda=( cell_base(1,ia)*vectb(1)+ & 377 & cell_base(2,ia)*vectb(2)+ & 378 & cell_base(3,ia)*vectb(3) )/norm2a 379 reducedb=( cell_base(1,ib)*vectb(1)+ & 380 & cell_base(2,ib)*vectb(2)+ & 381 & cell_base(3,ib)*vectb(3) )/norm2b 382 if( abs(abs(reduceda)-0.5d0)<tolsym .and. & 383 & abs(abs(reducedb)-0.5d0)<tolsym )then 384 cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib) 385 call holocell(cell_base,0,foundc,iholohedry,tolsym) 386 end if 387 else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then 388 fact=2 ; center=-3 389 cell_base(:,ia)=vectb(:)+minim(:,itrial) 390 cell_base(:,ib)=vectb(:)-minim(:,itrial) 391 ! Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib) 392 norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2 393 norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2 394 reduceda=( cell_base(1,ia)*vecta(1)+ & 395 & cell_base(2,ia)*vecta(2)+ & 396 & cell_base(3,ia)*vecta(3) )/norm2a 397 reducedb=( cell_base(1,ib)*vecta(1)+ & 398 & cell_base(2,ib)*vecta(2)+ & 399 & cell_base(3,ib)*vecta(3) )/norm2b 400 if( abs(abs(reduceda)-0.5d0)<tolsym .and. & 401 & abs(abs(reducedb)-0.5d0)<tolsym )then 402 cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib) 403 call holocell(cell_base,0,foundc,iholohedry,tolsym) 404 end if 405 end if 406 end if 407 408 ! Working hypothesis : the cell is rhombohedral, and 409 ! the three minimal vectors have same length and same absolute scalar product 410 if(foundc==0 .and. iholohedry==5 .and. & 411 & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then 412 if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and. & 413 & abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2) )then 414 fact=1 ; center=0 415 cell_base(:,:)=minim(:,:) 416 ! One might have to change the sign of one of the vectors 417 sign12=1 ; sign13=1 ; sign23=1 418 if(metmin(1,2)<0.0d0)sign12=-1 419 if(metmin(1,3)<0.0d0)sign13=-1 420 if(metmin(2,3)<0.0d0)sign23=-1 421 sumsign=sign12+sign13+sign23 422 if(sumsign==-1)then 423 if(sign12==1)cell_base(:,3)=-cell_base(:,3) 424 if(sign13==1)cell_base(:,2)=-cell_base(:,2) 425 if(sign23==1)cell_base(:,1)=-cell_base(:,1) 426 else if(sumsign==1)then 427 if(sign12==-1)cell_base(:,3)=-cell_base(:,3) 428 if(sign13==-1)cell_base(:,2)=-cell_base(:,2) 429 if(sign23==-1)cell_base(:,1)=-cell_base(:,1) 430 end if 431 call holocell(cell_base,0,foundc,iholohedry,tolsym) 432 end if 433 end if 434 435 ! DEBUG 436 ! write(std_out,*)' after test_3a, foundc=',foundc 437 ! write(std_out,*)' after test_3a, itrial=',itrial 438 ! write(std_out,*)' after test_3a, equal(:)=',equal(:) 439 ! ENDDEBUG 440 441 ! Working hypothesis : the cell is rhombohedral, one vector 442 ! is parallel to the trigonal axis 443 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then 444 vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib) 445 norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2 446 reduceda=( minim(1,itrial)*vecta(1)+ & 447 & minim(2,itrial)*vecta(2)+ & 448 & minim(3,itrial)*vecta(3) )/norm2trial 449 reducedb=( minim(1,itrial)*vectb(1)+ & 450 & minim(2,itrial)*vectb(2)+ & 451 & minim(3,itrial)*vectb(3) )/norm2trial 452 ! DEBUG 453 ! write(std_out,*)' reduceda,reducedb=',reduceda,reducedb 454 ! ENDDEBUG 455 if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and. & 456 & abs(abs(reducedb)-1.0d0/3.0d0)<tolsym ) then 457 ! Possibly change of sign to make positive the scalar product with 458 ! the vector parallel to the trigonal axis 459 if(reduceda<zero)vecta(:)=-vecta(:) 460 if(reducedb<zero)vectb(:)=-vectb(:) 461 ! Projection on the orthogonal plane 462 vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial) 463 vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial) 464 ! These two vectors should have an angle of 120 degrees 465 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 466 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 467 ! DEBUG 468 ! write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod 469 ! ENDDEBUG 470 if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then 471 fact=1 ; center=0 472 if(scalarprod>0.0d0)vectb(:)=-vectb(:) 473 ! Now vecta and vectb have an angle of 120 degrees 474 cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:) 475 cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:) 476 cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:) 477 ! DEBUG 478 ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1) 479 ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2) 480 ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3) 481 ! ENDDEBUG 482 call holocell(cell_base,0,foundc,iholohedry,tolsym) 483 end if 484 end if 485 end if 486 487 ! Working hypothesis : the cell is rhombohedral, one vector 488 ! is in the plane perpendicular to the trigonal axis 489 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then 490 vecta(:)=minim(:,ia)+minim(:,ib) 491 vectb(:)=minim(:,ia)-minim(:,ib) 492 norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2 493 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 494 norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2 495 reduceda=( cell_base(1,itrial)*vecta(1)+ & 496 & cell_base(2,itrial)*vecta(2)+ & 497 & cell_base(3,itrial)*vecta(3) )/norm2trial 498 reducedb=( cell_base(1,itrial)*vectb(1)+ & 499 & cell_base(2,itrial)*vectb(2)+ & 500 & cell_base(3,itrial)*vectb(3) )/norm2trial 501 if(abs(norm2trial-norm2a)<tolsym*norm2a .and. & 502 & abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial )then 503 fact=1 ; center=0 504 cell_base(:,1)=minim(:,ia) 505 cell_base(:,2)=-minim(:,ib) 506 cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial) 507 call holocell(cell_base,0,foundc,iholohedry,tolsym) 508 else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. & 509 & abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial )then 510 fact=1 ; center=0 511 cell_base(:,1)=minim(:,ia) 512 cell_base(:,2)=minim(:,ib) 513 cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial) 514 call holocell(cell_base,0,foundc,iholohedry,tolsym) 515 end if 516 end if 517 518 ! Working hypothesis : the cell is rhombohedral, two vectors 519 ! are in the plane perpendicular to the trigonal axis 520 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then 521 vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib) 522 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 523 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 524 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 525 if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then 526 ! This is in order to have 120 angle between vecta and vectb 527 if(scalarprod>0.0d0)vectb(:)=-vectb(:) 528 reduceda=( cell_base(1,itrial)*vecta(1)+ & 529 & cell_base(2,itrial)*vecta(2)+ & 530 & cell_base(3,itrial)*vecta(3) )/norm2a 531 reducedb=( cell_base(1,itrial)*vectb(1)+ & 532 & cell_base(2,itrial)*vectb(2)+ & 533 & cell_base(3,itrial)*vectb(3) )/norm2b 534 fact=1 ; center=0 535 cell_base(:,1)=minim(:,itrial) 536 if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then 537 cell_base(:,2)=minim(:,itrial)-vecta(:) 538 cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:) 539 call holocell(cell_base,0,foundc,iholohedry,tolsym) 540 else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then 541 cell_base(:,2)=minim(:,itrial)-vecta(:) 542 cell_base(:,3)=minim(:,itrial)+vectb(:) 543 call holocell(cell_base,0,foundc,iholohedry,tolsym) 544 else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then 545 cell_base(:,2)=minim(:,itrial)+vectb(:) 546 cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:) 547 call holocell(cell_base,0,foundc,iholohedry,tolsym) 548 else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then 549 cell_base(:,2)=minim(:,itrial)+vecta(:) 550 cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:) 551 call holocell(cell_base,0,foundc,iholohedry,tolsym) 552 else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then 553 cell_base(:,2)=minim(:,itrial)+vecta(:) 554 cell_base(:,3)=minim(:,itrial)-vectb(:) 555 call holocell(cell_base,0,foundc,iholohedry,tolsym) 556 else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then 557 cell_base(:,2)=minim(:,itrial)-vectb(:) 558 cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:) 559 call holocell(cell_base,0,foundc,iholohedry,tolsym) 560 end if 561 end if 562 end if 563 564 ! Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees 565 if(foundc==0 .and. iholohedry==2 .and. & 566 & ang90(ia)==1 .and. ang90(ib)==1 ) then 567 fact=1 ; center=0 568 cell_base(:,1)=minim(:,ia) 569 cell_base(:,2)=minim(:,itrial) 570 cell_base(:,3)=minim(:,ib) 571 ! Checks that the basis vectors are OK for the target holohedry 572 call holocell(cell_base,0,foundc,iholohedry,tolsym) 573 end if 574 575 ! Monoclinic holohedry, one-face-centered cell 576 ! Working hypothesis, two vectors have equal length. 577 do icase=1,5 578 if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then 579 vecta(:)=cell_base(:,ia)+cell_base(:,ib) 580 vectb(:)=cell_base(:,ia)-cell_base(:,ib) 581 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 582 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 583 ! The minim(:,trial) vector belongs to the 584 ! plane parallel to the cell_base(:,ia),cell_base(:,ib) plane 585 ! In that plane, must try minim(:,itrial), 586 ! as well as the 4 different combinations of 587 ! minim(:,itrial) with the vectors in the plane 588 if(icase==1)vectc(:)=minim(:,itrial) 589 if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia) 590 if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib) 591 if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia) 592 if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib) 593 norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2 594 sca=vectc(1)*vecta(1)+& 595 & vectc(2)*vecta(2)+& 596 & vectc(3)*vecta(3) 597 scb=vectc(1)*vectb(1)+& 598 & vectc(2)*vectb(2)+& 599 & vectc(3)*vectb(3) 600 ! DEBUG 601 ! write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb 602 ! ENDDEBUG 603 if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then 604 fact=2 ; center=3 605 ! The itrial direction is centered 606 cell_base(:,3)=vectc(:) 607 if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then 608 cell_base(:,2)=vecta(:) 609 cell_base(:,1)=vectb(:) 610 call holocell(cell_base,0,foundc,iholohedry,tolsym) 611 else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then 612 cell_base(:,2)=vectb(:) 613 cell_base(:,1)=vecta(:) 614 call holocell(cell_base,0,foundc,iholohedry,tolsym) 615 end if 616 end if 617 end if 618 end do ! icase=1,5 619 620 ! Monoclinic holohedry, one-face-centered cell, but non equivalent. 621 ! This case, one pair of vectors is orthogonal 622 if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then 623 vecta(:)=minim(:,ia) 624 vectb(:)=minim(:,ib) 625 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 626 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 627 ! Project the trial vector on the two vectors 628 reduceda=( minim(1,itrial)*vecta(1)+ & 629 & minim(2,itrial)*vecta(2)+ & 630 & minim(3,itrial)*vecta(3) )/norm2a 631 reducedb=( minim(1,itrial)*vectb(1)+ & 632 & minim(2,itrial)*vectb(2)+ & 633 & minim(3,itrial)*vectb(3) )/norm2b 634 if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then 635 fact=2 ; center=3 636 if(abs(abs(reduceda)-0.5d0)<tolsym)then 637 cell_base(:,2)=vecta(:) 638 cell_base(:,3)=vectb(:) 639 cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:)) 640 call holocell(cell_base,0,foundc,iholohedry,tolsym) 641 else if(abs(abs(reducedb)-0.5d0)<tolsym)then 642 cell_base(:,2)=vectb(:) 643 cell_base(:,3)=vecta(:) 644 cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:)) 645 call holocell(cell_base,0,foundc,iholohedry,tolsym) 646 end if 647 end if 648 end if 649 650 ! Monoclinic holohedry, one-face-centered cell, but non equivalent. 651 ! This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh 652 if(foundc==0 .and. iholohedry==2)then 653 ! Try to find a vector that belongs to the mediator plane, or the binary vector. 654 ! There must be one such vector, if centered monoclinic and no pair of vectors of equal length, 655 ! either among the three vectors, or among one of their differences or sums. 656 ! And there must be, among the two other vectors, one vector whose projection 657 ! on this vector is half the length of this vector. 658 vecta(:)=minim(:,ia) 659 vectb(:)=minim(:,ib) 660 ! Try the different possibilities for the vector on which the projection will be half ... 661 do ii=1,5 662 if(ii==1)vectc(:)=minim(:,itrial) 663 if(ii==2)vectc(:)=minim(:,itrial)+vecta(:) 664 if(ii==3)vectc(:)=minim(:,itrial)-vecta(:) 665 if(ii==4)vectc(:)=minim(:,itrial)+vectb(:) 666 if(ii==5)vectc(:)=minim(:,itrial)-vectb(:) 667 norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2 668 ! Project the two vectors on the trial vector 669 reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial 670 reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial 671 found=0 672 if(abs(abs(reduceda)-0.5d0)<tolsym)then 673 vin1(:)=vectc(:) 674 vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:)) 675 vext(:)=vectb(:) 676 found=1 677 else if(abs(abs(reducedb)-0.5d0)<tolsym)then 678 vin1(:)=vectc(:) 679 vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:)) 680 vext(:)=vecta(:) 681 found=1 682 end if 683 if(found==1)exit 684 end do 685 ! Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector. 686 ! One of them must be the binary vector if any. 687 ! On the other hand, vext is out-of-plane. Might belong to the mediator plane or not. 688 ! If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2. 689 ! The binary axis must be stored in cell_base(:,2) for conventional C-cell 690 if(found==1)then 691 found=0 692 693 ! Test vin1 being the binary axis 694 norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2 695 reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial 696 if(abs(reduceda)<tolsym)then ! vin1 is the binary axis and vext is in the mediator plane 697 found=1 698 cell_base(:,1)=vin2(:) 699 cell_base(:,2)=vin1(:) 700 cell_base(:,3)=vext(:) 701 else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection 702 found=1 703 cell_base(:,1)=vin2(:) 704 cell_base(:,2)=vin1(:) 705 cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half 706 else 707 ! Test vin2 being the binary axis 708 norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2 709 reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial 710 if(abs(reduceda)<tolsym)then ! vin2 is the binary axis and vext is in the mediator plane 711 found=1 712 cell_base(:,1)=vin1(:) 713 cell_base(:,2)=vin2(:) 714 cell_base(:,3)=vext(:) 715 else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection 716 found=1 717 cell_base(:,1)=vin1(:) 718 cell_base(:,2)=vin2(:) 719 cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half 720 end if 721 end if 722 723 if(found==1)then 724 fact=2 ; center=3 725 call holocell(cell_base,0,foundc,iholohedry,tolsym) 726 end if 727 end if 728 end if 729 730 end do ! Do-loop on three different directions 731 end do ! Do-loop on different target holohedries 732 733 if(foundc==0)then 734 iholohedry=1 ; fact=1 ; center=0 735 cell_base(:,:)=minim(:,:) 736 end if 737 738 !DEBUG 739 !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc 740 !write(std_out,*)' center=',center 741 !write(std_out,*)' iholohedry=',iholohedry 742 !ENDDEBUG 743 744 !-------------------------------------------------------------------------- 745 !Final check on the Bravais lattice, using the basis vectors 746 747 !Recompute the metric tensor 748 if(foundc==1)then 749 do ii=1,3 750 metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+& 751 & cell_base(2,:)*cell_base(2,ii)+& 752 & cell_base(3,:)*cell_base(3,ii) 753 end do 754 end if 755 756 !Examine the angles and vector lengths 757 ang90(:)=0 758 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1 759 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1 760 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1 761 equal(:)=0 762 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1 763 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1 764 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1 765 766 !DEBUG 767 !write(std_out,*)' symlatt : recompute the metric tensor ' 768 !write(std_out,*)' ang90=',ang90 769 !write(std_out,*)' equal=',equal 770 !ENDDEBUG 771 772 !The axes will be aligned with the previously determined 773 !basis vectors, EXCEPT for the tetragonal cell, see later 774 axes(:,:)=cell_base(:,:) 775 776 found=0 777 !Check orthogonal conventional cells 778 if(ang90(1)+ang90(2)+ang90(3)==3)then 779 780 ! Cubic system 781 if(equal(1)+equal(2)+equal(3)==3)then 782 ! However, one-face centered is not admitted 783 if(center==0 .or. center==-1 .or. center==-3)then 784 iholohedry=7 ; found=1 785 if(center==0)then 786 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cP (primitive cubic)' 787 else if(center==-1)then 788 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cI (body-centered cubic)' 789 else if(center==-3)then 790 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cF (face-centered cubic)' 791 end if 792 end if 793 end if 794 795 ! Tetragonal system 796 if(found==0 .and. & 797 & (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then 798 ! However, one-face centered or face-centered is not admitted 799 if(center==0 .or. center==-1)then 800 iholohedry=4 ; found=1 801 if(equal(1)==1)then 802 axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3) 803 else if(equal(2)==1)then 804 axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3) 805 else if(equal(3)==1)then 806 axes(:,:)=cell_base(:,:) 807 end if 808 if(center==0)then 809 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tP (primitive tetragonal)' 810 else if(center==-1)then 811 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tI (body-centered tetragonal)' 812 end if 813 end if 814 end if 815 816 ! Orthorhombic system 817 if(found==0)then 818 iholohedry=3 ; found=1 819 axes(:,:)=cell_base(:,:) 820 if(center==0)then 821 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oP (primitive orthorhombic)' 822 else if(center==-1)then 823 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oI (body-centered orthorhombic)' 824 else if(center==1 .or. center==2 .or. center==3)then 825 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oC (one-face-centered orthorhombic)' 826 else if(center==-3)then 827 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oF (face-centered orthorhombic)' 828 end if 829 end if 830 831 else 832 833 ! Hexagonal system 834 if(found==0 .and. ang90(1)==1 .and. ang90(2)==1 .and. equal(3)==1 .and. (2*metmin(2,1)+metmin(1,1))<tolsym*metmin(1,1))then 835 iholohedry=6 ; found=1 836 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hP (primitive hexagonal)' 837 end if 838 839 ! Rhombohedral system 840 if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and. & 841 & abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2) .and. & 842 & abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then 843 iholohedry=5 ; found=1 844 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hR (rhombohedral)' 845 end if 846 847 ! Monoclinic system 848 if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then 849 iholohedry=2 ; found=1 850 if(center==0)then 851 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mP (primitive monoclinic)' 852 else if(center==3)then 853 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mC (one-face-centered monoclinic)' 854 end if 855 end if 856 857 ! Triclinic system 858 if(found==0)then 859 iholohedry=1 ; found=1 860 write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is aP (primitive triclinic)' 861 end if 862 863 end if 864 865 call wrtout(std_out,message,'COLL') 866 867 !-------------------------------------------------------------------------- 868 !Make sure that axes form a right-handed coordinate system 869 !(Note : this should be done in the body of the routine, 870 !by making changes that leave the sign of the mixed product of the three 871 !vectors invariant) 872 determinant=axes(1,1)*axes(2,2)*axes(3,3) & 873 & +axes(1,2)*axes(2,3)*axes(3,1) & 874 & +axes(1,3)*axes(3,2)*axes(2,1) & 875 & -axes(1,1)*axes(3,2)*axes(2,3) & 876 & -axes(1,3)*axes(2,2)*axes(3,1) & 877 & -axes(1,2)*axes(2,1)*axes(3,3) 878 if(determinant<0.0d0)then 879 axes(:,:)=-axes(:,:) 880 end if 881 882 !-------------------------------------------------------------------------- 883 !Prefer symmetry axes on the same side as the primitive axes, 884 !when the changes are allowed 885 do 886 do ia=1,3 887 scprods(ia,:)=axes(1,ia)*rprimd(1,:)+& 888 & axes(2,ia)*rprimd(2,:)+& 889 & axes(3,ia)*rprimd(3,:) 890 norm2trial=sum(axes(:,ia)**2) 891 scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial) 892 end do 893 do ia=1,3 894 norm2trial=sum(rprimd(:,ia)**2) 895 scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial) 896 end do 897 898 ! One should now try all the generators of the 899 ! proper rotations of each Bravais lattice, coupled with change of 900 ! signs of each vector. This is not done systematically in what follows ... 901 ! Here, the third axis is left unchanged 902 if(iholohedry/=5)then 903 if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then 904 axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2) 905 cycle 906 end if 907 end if 908 ! The first (or second) axis is left unchanged 909 if(iholohedry/=5 .and. iholohedry/=6)then 910 if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then 911 axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3) 912 cycle 913 end if 914 if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then 915 axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3) 916 cycle 917 end if 918 end if 919 ! Permutation of the three axis 920 if(iholohedry==5 .or. iholohedry==7)then 921 trace=scprods(1,1)+scprods(2,2)+scprods(3,3) 922 if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then 923 vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3) 924 axes(:,3)=axes(:,2); axes(:,2)=vecta(:) 925 cycle 926 end if 927 if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then 928 vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2) 929 axes(:,2)=axes(:,3); axes(:,3)=vecta(:) 930 cycle 931 end if 932 ! This case is observed when the three new vectors 933 ! are pointing opposite to the three original vectors 934 ! One takes their opposite, then switch to of them, then process 935 ! them again in the loop 936 if(sum(scprods(:,:))<tolsym)then 937 axes(:,1)=-axes(:,1) 938 vecta(:)=-axes(:,2) 939 axes(:,2)=-axes(:,3) 940 axes(:,3)=vecta(:) 941 cycle 942 end if 943 end if 944 exit 945 ! Other cases might be coded ... 946 end do 947 948 !-------------------------------------------------------------------------- 949 950 !DEBUG 951 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',& 952 !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 953 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',& 954 !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3) 955 !ENDDEBUG 956 957 !Compute the coordinates of rprimd in the system defined by axes(:,:) 958 call matr3inv(axes,axesinvt) 959 do ii=1,3 960 coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ & 961 & rprimd(2,ii)*axesinvt(2,:)+ & 962 & rprimd(3,ii)*axesinvt(3,:) 963 end do 964 965 !Check that the coordinates are integers, or half-integer in 966 !the case there is a centering, and generate integer coordinates 967 do ii=1,3 968 do jj=1,3 969 val=coord(ii,jj)*fact 970 if(abs(val-nint(val))>fact*tolsym)then 971 write(message,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')& 972 & 'One of the coordinates of rprimd in axes is non-integer,',ch10,& 973 & 'or non-half-integer (if centering).',ch10,& 974 & 'coord=',coord(:,1),ch10,& 975 & ' ',coord(:,2),ch10,& 976 & ' ',coord(:,3),ch10,& 977 & 'fact=',fact 978 MSG_BUG(message) 979 end if 980 icoord(ii,jj)=nint(val) 981 end do 982 end do 983 984 !Store the bravais lattice characteristics 985 bravais(1)=iholohedry 986 bravais(2)=center 987 bravais(3:5)=icoord(1:3,1) 988 bravais(6:8)=icoord(1:3,2) 989 bravais(9:11)=icoord(1:3,3) 990 991 !-------------------------------------------------------------------------- 992 !Initialize the set of symmetries 993 !Bravais lattices are always invariant under identity and inversion 994 995 !Identity and inversion 996 ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:) 997 nptsym=2 998 999 !Keep this for IFCv70 compiler 1000 if(nptsym/=2)then 1001 write(message,'(a,a,a,a)')ch10,& 1002 & ' symlatt : BUG -',ch10,& 1003 & ' Crazy error, compiler bug ' 1004 call wrtout(std_out,message,'COLL') 1005 end if 1006 1007 !-------------------------------------------------------------------------- 1008 !Initialize some generators 1009 !gen6 is defined in a coordinated system with gamma=120 degrees 1010 gen6(:,:)=0 ; gen6(3,3)=1 ; gen6(1,1)=1 ; gen6(1,2)=-1 ; gen6(2,1)=1 1011 gen3(:,:)=0 ; gen3(1,2)=1 ; gen3(2,3)=1 ; gen3(3,1)=1 1012 gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1 1013 gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1 1014 gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1 1015 1016 !-------------------------------------------------------------------------- 1017 1018 !Define the generators for each holohedry (inversion is already included) 1019 if(iholohedry==6)then 1020 ngen=2 1021 gen(:,:,1)=gen2xy(:,:) ; order(1)=2 1022 gen(:,:,2)=gen6(:,:) ; order(2)=6 1023 else if(iholohedry==5)then 1024 ngen=2 1025 gen(:,:,1)=gen2xy(:,:) ; order(1)=2 1026 gen(:,:,2)=gen3(:,:) ; order(2)=3 1027 else 1028 gen(:,:,1)=gen2y(:,:) ; order(1)=2 1029 gen(:,:,2)=gen2z(:,:) ; order(2)=2 1030 gen(:,:,3)=gen2xy(:,:) ; order(3)=2 1031 gen(:,:,4)=gen3(:,:) ; order(4)=3 1032 if(iholohedry<=4)ngen=iholohedry-1 1033 if(iholohedry==7)ngen=4 1034 end if 1035 1036 !Build the point symmetry operations from generators, in the reduced system 1037 !of coordinates defined by axes(:,:) 1038 if(ngen/=0)then 1039 do igen=1,ngen 1040 do isym=1+nptsym,order(igen)*nptsym 1041 jsym=isym-nptsym 1042 do ii=1,3 1043 ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ & 1044 & gen(:,2,igen)*ptsymrel(2,ii,jsym)+ & 1045 & gen(:,3,igen)*ptsymrel(3,ii,jsym) 1046 end do 1047 end do 1048 nptsym=order(igen)*nptsym 1049 1050 end do 1051 end if 1052 1053 !-------------------------------------------------------------------------- 1054 1055 !Transform symmetry matrices in the system defined by rprimd 1056 call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym) 1057 1058 end subroutine symlatt