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