TABLE OF CONTENTS
ABINIT/symfind [ Functions ]
NAME
symfind
FUNCTION
Symmetry finder. From the symmetries of the Bravais lattice (ptsymrel), select those that leave invariant the system, and generate the corresponding tnons vectors. The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972)
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
berryopt = 4/14, 6/16, 7/17: electric or displacement field efield=cartesian coordinates of the electric field gprimd(3,3)=dimensional primitive translations for reciprocal space msym=default maximal number of symmetries natom=number of atoms in cell. noncoll=1 if non-collinear magnetism is activated else 0 nptsym=number of point symmetries of the Bravais lattice nucdipmom(3,natom) (optional) array of nuclear dipole moments nzchempot=if non-zero, means that a z-spatially varying chemical potential is added ptsymrel(3,3,1:msym)= nptsym point-symmetry operations of the Bravais lattice in real space in terms of primitive translations. spinat(3,natom)=initial spin of each atom, in unit of hbar/2. tolsym=tolerance for the symmetries typat(natom)=integer identifying type of atom. use_inversion=1 if inversion can be included in set of symmetries xred(3,natom)=reduced coordinates of atoms in terms of real space primitive translations
OUTPUT
nsym=actual number of symmetries symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations symrel(3,3,1:msym)= nsym symmetry operations in real space in terms of primitive translations tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would be 0 0 0 each for a symmorphic space group)
PARENTS
ingeo,inqpt,m_ab7_symmetry,m_effective_potential_file,m_tdep_sym m_use_ga,thmeig
CHILDREN
wrtout
SOURCE
59 #if defined HAVE_CONFIG_H 60 #include "config.h" 61 #endif 62 63 #include "abi_common.h" 64 65 subroutine symfind(berryopt,efield,gprimd,jellslab,msym,natom,noncoll,nptsym,nsym,& 66 & nzchempot,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,& 67 & nucdipmom) 68 69 use defs_basis 70 use m_errors 71 use m_profiling_abi 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'symfind' 77 use interfaces_14_hidewrite 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: berryopt,jellslab,msym,natom,noncoll,nptsym,nzchempot,use_inversion 85 integer,intent(out) :: nsym 86 real(dp),intent(in) :: tolsym 87 !arrays 88 integer,intent(in) :: ptsymrel(3,3,msym),typat(natom) 89 integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i 90 real(dp),intent(in) :: efield(3),gprimd(3,3),spinat(3,natom),xred(3,natom) 91 real(dp),optional :: nucdipmom(3,natom) 92 real(dp),intent(inout) :: tnons(3,msym) !vz_i 93 94 !Local variables------------------------------- 95 ! TRUE if antiferro symmetries are used with non-collinear magnetism 96 !scalars 97 integer :: found3,foundcl,iatom,iatom0,iatom1,iatom2,iatom3,iclass,iclass0,ii 98 integer :: isym,jj,kk,natom0,nclass,ntrial,printed,trialafm,trialok 99 real(dp) :: spinatcl2,spinatcl20,det 100 logical,parameter :: afm_noncoll=.true. 101 logical :: test_sameabscollin,test_sameabsnoncoll,test_samespin 102 character(len=500) :: message 103 !arrays 104 integer,allocatable :: class(:,:),natomcl(:),typecl(:) 105 real(dp) :: diff(3),efieldrot(3),sxred0(3),symnucdipmom2(3) 106 real(dp) :: symspinat2(3),symxred2(3),trialnons(3) 107 real(dp),allocatable :: spinatcl(:,:),spinatred(:,:) 108 109 !************************************************************************** 110 111 !DEBUG 112 ! write(std_out,*)' symfind : enter' 113 ! call flush(6) 114 ! write(std_out,*)' symfind : nzchempot= ',nzchempot 115 ! write(std_out,*)' ptsymrel matrices are :' 116 ! do isym=1,nptsym 117 ! write(std_out,'(i4,4x,9i4)' )isym,ptsymrel(:,:,isym) 118 ! end do 119 ! write(std_out,*)' symfind : natom=',natom 120 ! do iatom=1,natom 121 ! write(std_out,*)' atom number',iatom 122 ! write(std_out,*)' typat =',typat(iatom) 123 ! write(std_out,*)' spinat =',spinat(:,iatom) 124 ! write(std_out,*)' xred =',xred(:,iatom) 125 ! end do 126 ! call flush(6) 127 !ENDDEBUG 128 129 !Find the number of classes of atoms (type and spinat must be identical, 130 !spinat might differ by a sign, if aligned with the z direction) 131 !natomcl(iclass) will contain the number of atoms in the class 132 !typecl(iclass) will contain the type of the atoms in the class 133 !spinatcl(1:3,iclass) will contain the spinat of the atoms in the class 134 !class(1:natomclass(iclass),iclass) will contain the index of the 135 !atoms belonging to the class 136 ABI_ALLOCATE(class,(natom+3,natom)) 137 ABI_ALLOCATE(natomcl,(natom)) 138 ABI_ALLOCATE(typecl,(natom)) 139 ABI_ALLOCATE(spinatcl,(3,natom)) 140 141 !Initialise with the first atom 142 nclass=1 143 natomcl(1)=1 144 typecl(1)=typat(1) 145 spinatcl(:,1)=spinat(:,1) 146 class(1,1)=1 147 if(natom>1)then 148 do iatom=2,natom 149 ! DEBUG 150 ! write(std_out,*)' symfind : examine iatom=',iatom 151 ! ENDDEBUG 152 foundcl=0 153 do iclass=1,nclass 154 ! Compare the typat and spinat of atom iatom with existing ones. 155 ! Admit either identical spinat, or z-aligned spinat with same 156 ! absolute magnitude 157 if( typat(iatom)==typecl(iclass)) then 158 test_samespin= & 159 & abs(spinat(1,iatom)-spinatcl(1,iclass))<tolsym .and. & 160 & abs(spinat(2,iatom)-spinatcl(2,iclass))<tolsym .and. & 161 & abs(spinat(3,iatom)-spinatcl(3,iclass))<tolsym 162 test_sameabscollin= & 163 & noncoll==0 .and.& 164 & abs(spinat(1,iatom))<tolsym .and. abs(spinatcl(1,iclass))<tolsym .and.& 165 & abs(spinat(2,iatom))<tolsym .and. abs(spinatcl(2,iclass))<tolsym .and.& 166 & abs(abs(spinat(3,iatom))-abs(spinatcl(3,iclass)))<tolsym 167 test_sameabsnoncoll= & 168 & noncoll==1 .and. afm_noncoll .and. & 169 & abs(spinat(1,iatom)+spinatcl(1,iclass))<tolsym .and. & 170 & abs(spinat(2,iatom)+spinatcl(2,iclass))<tolsym .and. & 171 & abs(spinat(3,iatom)+spinatcl(3,iclass))<tolsym 172 if( test_samespin .or. test_sameabscollin .or. test_sameabsnoncoll) then 173 ! DEBUG 174 ! write(std_out,*)' symfind : find it belongs to class iclass=',iclass 175 ! write(std_out,*)' symfind : spinat(:,iatom)=',spinat(:,iatom) 176 ! write(std_out,*)' symfind : spinatcl(:,iclass)=',spinatcl(:,iclass) 177 ! write(std_out,*)' symfind : test_samespin,test_sameabscollin,test_sameabsnoncoll=',& 178 ! & test_samespin,test_sameabscollin,test_sameabsnoncoll 179 ! ENDDEBUG 180 natomcl(iclass)=natomcl(iclass)+1 181 class(natomcl(iclass),iclass)=iatom 182 foundcl=1 183 exit 184 end if 185 end if 186 end do 187 ! If no class with these characteristics exist, create one 188 if(foundcl==0)then 189 nclass=nclass+1 190 natomcl(nclass)=1 191 typecl(nclass)=typat(iatom) 192 spinatcl(:,nclass)=spinat(:,iatom) 193 class(1,nclass)=iatom 194 end if 195 end do 196 end if 197 198 !DEBUG 199 !write(std_out,*)' symfind : found ',nclass,' nclass of atoms' 200 !do iclass=1,nclass 201 !write(std_out,*)' class number',iclass 202 !write(std_out,*)' natomcl =',natomcl(iclass) 203 !write(std_out,*)' typecl =',typecl(iclass) 204 !write(std_out,*)' spinatcl=',spinatcl(:,iclass) 205 !write(std_out,*)' class =',(class(iatom,iclass),iatom=1,natomcl(iclass)) 206 !end do 207 !ENDDEBUG 208 209 !Select the class with the least number of atoms, and non-zero spinat if any 210 !It is important to select a magnetic class of atom, if any, otherwise 211 !the determination of the initial (inclusive) set of symmetries takes only 212 !non-magnetic symmetries, and not both magnetic and non-magnetic ones, see later. 213 iclass0=1 214 natom0=natomcl(1) 215 spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2 216 if(nclass>1)then 217 do iclass=2,nclass 218 spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2 219 if( (natomcl(iclass)<natom0 .and. (spinatcl20<tolsym .or. spinatcl2>tolsym)) & 220 & .or. (spinatcl20<tolsym .and. spinatcl2>tolsym) )then 221 iclass0=iclass 222 natom0=natomcl(iclass) 223 spinatcl20=spinatcl2 224 end if 225 end do 226 end if 227 228 printed=0 229 230 !DEBUG 231 !write(std_out,*)' symfind : has selected iclass0=',iclass0 232 !write(std_out,*)' # iatom xred spinat ' 233 !do iatom0=1,natomcl(iclass0) 234 !iatom=class(iatom0,iclass0) 235 !write(std_out,'(2i4,6f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom) 236 !end do 237 !ENDDEBUG 238 239 !If non-collinear spinat have to be used, transfer them in reduced coordinates 240 if (noncoll==1) then 241 ABI_ALLOCATE(spinatred,(3,natom)) 242 do iatom=1,natom 243 do ii=1,3 244 spinatred(ii,iatom)=gprimd(1,ii)*spinat(1,iatom) & 245 & +gprimd(2,ii)*spinat(2,iatom) & 246 & +gprimd(3,ii)*spinat(3,iatom) 247 end do 248 end do 249 end if 250 251 !Big loop over each symmetry operation of the Bravais lattice 252 nsym=0 253 do isym=1,nptsym 254 255 ! ji: Check whether symmetry operation leaves efield invariant 256 if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 257 & berryopt==14 .or. berryopt==16 .or. berryopt==17) then 258 efieldrot(:) = ptsymrel(:,1,isym)*efield(1) + & 259 & ptsymrel(:,2,isym)*efield(2) + & 260 & ptsymrel(:,3,isym)*efield(3) 261 diff(:)=efield(:)-efieldrot(:) 262 if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle 263 end if 264 265 if (use_inversion==0) then 266 det=ptsymrel(1,1,isym)*ptsymrel(2,2,isym)*ptsymrel(3,3,isym)+& 267 & ptsymrel(2,1,isym)*ptsymrel(3,2,isym)*ptsymrel(1,3,isym)+& 268 & ptsymrel(1,2,isym)*ptsymrel(2,3,isym)*ptsymrel(3,1,isym) - & 269 & (ptsymrel(3,1,isym)*ptsymrel(2,2,isym)*ptsymrel(1,3,isym)+& 270 & ptsymrel(2,1,isym)*ptsymrel(1,2,isym)*ptsymrel(3,3,isym)+& 271 & ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym)) 272 if(det==-1) cycle 273 end if 274 275 ! jellium slab and spatially varying chemical potential cases: 276 ! (actually, an inversion symmetry/mirror plane perpendicular to z symmetry operation might still be allowed... TO BE DONE !) 277 if (jellslab/=0 .or. nzchempot/=0) then 278 ! check whether symmetry operation produce a rotation only in the xy plane 279 if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. & 280 & ptsymrel(3,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle 281 ! check whether symmetry operation does not change the z 282 if( ptsymrel(3,3,isym)/=1 ) cycle 283 end if 284 285 ! Select a tentative set of associated translations 286 ! First compute the symmetric of the first atom in the smallest class, 287 ! using the point symmetry 288 iatom0=class(1,iclass0) 289 sxred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ & 290 & ptsymrel(:,2,isym)*xred(2,iatom0)+ & 291 & ptsymrel(:,3,isym)*xred(3,iatom0) 292 293 ! From the set of possible images, deduce tentative translations, 294 ! and magnetic factor then test whether it send each atom on a symmetric one 295 ntrial=0 296 do ii=1,natom0 297 iatom1=class(ii,iclass0) 298 299 ! The tentative translation is found 300 trialnons(:)=xred(:,iatom1)-sxred0(:) 301 trialafm=1 302 if(abs(spinat(3,iatom1)-spinat(3,iatom0))>tolsym)trialafm=-1 303 304 if(sum(abs(spinat(:,iatom1)*trialafm-spinat(:,iatom0)))>tolsym)then 305 write(message,'(3a,3i5)')& 306 & 'Problem with matching the spin part within a class.',ch10,& 307 & 'isym,iatom0,iatom1=',isym,iatom0,iatom1 308 MSG_ERROR_CLASS(message, "TolSymError") 309 end if 310 ! jellium slab case: check whether symmetry operation has no translational 311 ! component along z 312 if( jellslab/=0 .and. abs(trialnons(3)) > tolsym ) cycle 313 trialok=1 314 315 ! DEBUG 316 ! write(std_out,*)' isym,trialnons(:),trialafm =',isym,trialnons(:),trialafm 317 ! ENDDEBUG 318 319 ! Loop over all classes, then all atoms in the class, 320 ! to find whether they have a symmetric 321 do iclass=1,nclass 322 do jj=1,natomcl(iclass) 323 324 iatom2=class(jj,iclass) 325 ! Generate the tentative symmetric position of iatom2 326 symxred2(:)=ptsymrel(:,1,isym)*xred(1,iatom2)+ & 327 & ptsymrel(:,2,isym)*xred(2,iatom2)+ & 328 & ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:) 329 ! Generate the tentative symmetric spinat of iatom2 330 if (noncoll==0) then 331 symspinat2(:)=trialafm*spinat(:,iatom2) 332 else 333 symspinat2(:)=trialafm*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ & 334 & ptsymrel(:,2,isym)*spinatred(2,iatom2)+ & 335 & ptsymrel(:,3,isym)*spinatred(3,iatom2)) 336 end if 337 if(present(nucdipmom)) then 338 ! Generate the tentative symmetric nuclear dipole moment of iatom2 339 symnucdipmom2(:)=(ptsymrel(:,1,isym)*nucdipmom(1,iatom2)+ & 340 & ptsymrel(:,2,isym)*nucdipmom(2,iatom2)+ & 341 & ptsymrel(:,3,isym)*nucdipmom(3,iatom2)) 342 end if 343 344 ! DEBUG 345 ! write(std_out,'(a,3f12.4,a,3f12.4)') ' Send atom at xred=',xred(:,iatom2),' to ',symxred2(:) 346 ! ENDDEBUG 347 348 ! Check whether there exists an atom of the same class at the 349 ! same location, with the correct spinat and nuclear dipole moment 350 do kk=1,natomcl(iclass) 351 352 found3=1 353 iatom3=class(kk,iclass) 354 ! Check the location 355 diff(:)=xred(:,iatom3)-symxred2(:) 356 diff(:)=diff(:)-nint(diff(:)) 357 if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )found3=0 358 ! Check the spinat 359 if (noncoll==0) then 360 diff(:)=spinat(:,iatom3)-symspinat2(:) 361 else 362 diff(:)=spinatred(:,iatom3)-symspinat2(:) 363 end if 364 if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )found3=0 365 if(present(nucdipmom)) then 366 ! Check the nuclear dipole moment 367 diff(:) = nucdipmom(:,iatom3) - symnucdipmom2(:) 368 if(any(diff>tolsym))found3=0 369 end if 370 371 if(found3==1)exit 372 373 ! End loop over iatom3 374 end do 375 376 if(found3==0)then 377 trialok=0 378 exit 379 end if 380 381 ! End loop over iatom2 382 end do 383 384 if(trialok==0)exit 385 386 ! End loop over all classes 387 end do 388 389 if(trialok==1)then 390 nsym=nsym+1 391 if(nsym>msym)then 392 write(message,'(3a,i0,4a)')& 393 & 'The number of symmetries (including non-symmorphic translations)',ch10,& 394 & 'is larger than maxnsym=',msym,ch10,& 395 & 'Action: increase maxnsym in the input, or take a cell that is primitive, ',ch10,& 396 & 'or at least smaller than the present one.' 397 MSG_ERROR(message) 398 end if 399 ntrial=ntrial+1 400 symrel(:,:,nsym)=ptsymrel(:,:,isym) 401 symafm(nsym)=trialafm 402 tnons(:,nsym)=trialnons(:)-nint(trialnons(:)-tolsym) 403 end if 404 405 ! End the loop on tentative translations 406 end do 407 408 ! End big loop over each symmetry operation of the Bravais lattice 409 end do 410 411 ABI_DEALLOCATE(class) 412 ABI_DEALLOCATE(natomcl) 413 ABI_DEALLOCATE(spinatcl) 414 ABI_DEALLOCATE(typecl) 415 if (noncoll==1) then 416 ABI_DEALLOCATE(spinatred) 417 end if 418 419 !DEBUG 420 write(message,'(a,I0,a)')' symfind : exit, nsym=',nsym,ch10 421 write(message,'(2a)') trim(message),' symrel matrices, symafm and tnons are :' 422 call wrtout(std_out,message,'COLL') 423 do isym=1,nsym 424 write(message,'(i4,4x,3i4,2x,3i4,2x,3i4,4x,i4,4x,3f8.4)' ) isym,symrel(:,:,isym),& 425 & symafm(isym),tnons(:,isym) 426 call wrtout(std_out,message,'COLL') 427 end do 428 429 !stop 430 !ENDDEBUG 431 432 end subroutine symfind