Dear WIEN2k users,
Unfortunately since WIEN2k_23.1 a severe bug was introduced in
SRC_symmetso/angle.f
It may lead to a wrong symmetry reduction due to spin-orbit interactions.
The bug is active in magnetic (spinpolarized) spin-orbit calculations at
least for B and F cubic and orthorhombic CXY, CXZ and CYZ lattices,
where the magnetization direction was multiplied by the primitive
Bravais matrix and thus e.g the (0 0 1) direction was moved into (1 1 1)
or (1 1 0), resulting into wrong Euler angles theta and phi and thus to
wrong symmetry.
It should NOT be active for P, H and R lattices.
If you have done SO calculations with WIEN2k_23 or 24, please check
grep THETA case.outsymso
and check if theta (the angle between M and z) and phi agree with
expectations according to the selected M-direction in case.inso.
Please replace the attached angle.f subroutine in SRC_symmetso and
recompile.
Best regards
Peter Blaha
PS: Special thanks to Stefaan, who reported the problem of SO calculations.
Am 26.11.2024 um 16:22 schrieb Peter Blaha:
initso_lapw should reduce the symmetry operations to 16, not 12 (if M=
0 0 1). It makes cubic --> tetragonal
Unfortunately, a quick test shows that in *outputsymso theta and phi
is wrong and thus it reduces to a wrong symmetry.
I'll debug it as quick as possible.
Peter
Am 26.11.2024 um 16:03 schrieb Stefaan Cottenier via Wien:
Dear WIEN2k community,
I am struggling with lapwdm, for calculating the orbital magnetic
moment. This feature worked fine many years ago, but I am not able to
get it working right now.
Tests are done with WIEN2k 23.2, not with the current version (not
available on the resources to which I currently have access). But I
guess this does not matter for this test.
bcc Fe is the test case:
blebleble s-o calc. M|| 0.00 0.00 1.00
B 1 229
RELA
5.410352 5.410352 5.410352 90.000000 90.000000 90.000000
ATOM -1: X=0.00000000 Y=0.00000000 Z=0.00000000
MULT= 1 ISPLIT=-2
Fe1 NPT= 781 R0=.000050000 RMT= 2.19000 Z: 26.00000
LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000
0.0000000 1.0000000 0.0000000
0.0000000 0.0000000 1.0000000
0 NUMBER OF SYMMETRY OPERATIONS
Initialization:
init_lapw -prec 1 -rkmax 7.5 -numk 8000
init_so_lapw
In the latter step, all defaults were accepted and symmetso was
allowed to run. It reduces the number of symmetry operations from 48
to 12. There are several ‘warnings’ in case.outsymso, but I guess
these are normal and indicate the 36 symmetry operations that are
eliminated.
The case is run by:
runsp -so -cc 0.00001
and converges in 12 iterations, without problems.
I then prepare the following case.indmc file:
-12. Emin cutoff energy
1 number of atoms for which density matrix is
calculated
1 4 0 1 2 3 index of 1st atom, number of L's, L1
1 3 r-index, (l,s)index
This should give the orbital contribution to the orbital moment for
the s, p, d and f orbitals of the iron atom (not all of them relevant
for this element and this property, I know).
When running lapwdm, there is an error message in stdout :
x lapwdm -c -so -up
Error: check case.outputdmup, symmetry might be wrong
The output file case.scfdmup has only a single line:
Spin-polarized + s-o calculation, M|| 0.000 0.000 1.000
And the output file case.outputdmup terminates with an extra line
(which is presumably an error message) after heaving dealt with the
second symmetry operation:
2 Euler angles: a,b,c: 270.0 90.0 0.0
# of operation, phase, det:
2 4.71238896548353 0.000000000000000E+000
symm. operation 2 so-det= 0.000000000000000E+000
Did I overlook something, has something changed to the procedure (I
can’t find any hint for this in the usersguide) or is this feature
broken?
Thanks,
Stefaan
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at: http://www.mail-archive.com/
wien@zeus.theochem.tuwien.ac.at/index.html
--
-----------------------------------------------------------------------
Peter Blaha, Inst. f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-158801165300
Email: peter.bl...@tuwien.ac.at
WWW: http://www.imc.tuwien.ac.at WIEN2k: http://www.wien2k.at
-------------------------------------------------------------------------
SUBROUTINE ANGLE(XMS,alat,alpha,lattic,theta,phi)
IMPLICIT REAL*8 (A-H,O-Z)
character*4 lattic
logical ortho
dimension alat(3),alpha(3),br2(3,3)
!*******************************************************************
!
DIMENSION XMS(3),dif(3),help(3)
do i=1,3
dif(i)=xms(i)
enddo
pi=4.d0*atan(1.d0)
alpha0=alpha(1)
beta0=alpha(2)
gamma=alpha(3)
gamma1=alpha(3)*pi/180.d0
beta1=alpha(2)*pi/180.d0
alpha1=alpha(1)*pi/180.d0
cosg1=0.
IF(LATTIC(1:1).EQ.'H') GOTO 10
IF(LATTIC(1:1).EQ.'S') GOTO 20
IF(LATTIC(1:1).EQ.'P') GOTO 20
IF(LATTIC(1:1).EQ.'B') GOTO 30
IF(LATTIC(1:1).EQ.'F') GOTO 40
IF(LATTIC(1:1).EQ.'C') GOTO 50
IF(LATTIC(1:1).EQ.'R') GOTO 80
write(6,*)lattic
STOP 'LATTIC WRONG'
!
!.....HEXAGONAL CASE
10 BR2(1,1)=SQRT(3.d0)/2.d0
BR2(1,2)=-.5d0
BR2(1,3)=0.0d0
BR2(2,1)=0.0d0
BR2(2,2)=1.0d0
BR2(2,3)=0.0d0
BR2(3,1)=0.0d0
BR2(3,2)=0.0d0
BR2(3,3)=1.d0
ORTHO=.FALSE.
GOTO 100
!
!.....PRIMITIVE LATTICE CASE
20 continue
!
cosg1=(cos(gamma1)-cos(alpha1)*cos(beta1))/sin(alpha1)/sin(beta1)
gamma0=acos(cosg1)
BR2(1,1)=1.0d0*sin(gamma0)*sin(beta1)
BR2(1,2)=1.0d0*cos(gamma0)*sin(beta1)
BR2(1,3)=1.0d0*cos(beta1)
BR2(2,1)=0.0d0
BR2(2,2)=1.0d0*sin(alpha1)
BR2(2,3)=1.0d0*cos(alpha1)
BR2(3,1)=0.0d0
BR2(3,2)=0.0d0
BR2(3,3)=1.0d0
ORTHO=.TRUE.
if(alpha(3).ne.90.d0) ortho=.false.
if(alpha(2).ne.90.d0) ortho=.false.
if(alpha(1).ne.90.d0) ortho=.false.
! write(*,*) alpha0,beta0,gamma0,ortho,br2
GOTO 100
!
!.....BC CASE (DIRECT LATTICE)
30 CONTINUE
BR2(1,1)=-0.5d0
BR2(1,2)=0.5d0
BR2(1,3)=0.5d0
BR2(2,1)=0.5d0
BR2(2,2)=-0.5d0
BR2(2,3)=0.5d0
BR2(3,1)=0.5d0
BR2(3,2)=0.5d0
BR2(3,3)=-0.5d0
ORTHO=.TRUE.
GOTO 100
!
!.....FC CASE (DIRECT LATTICE)
40 CONTINUE
BR2(1,1)=0.0d0
BR2(1,2)=0.5d0
BR2(1,3)=0.5d0
BR2(2,1)=0.5d0
BR2(2,2)=0.0d0
BR2(2,3)=0.5d0
BR2(3,1)=0.5d0
BR2(3,2)=0.5d0
BR2(3,3)=0.0d0
ORTHO=.TRUE.
GOTO 100
!
!.....CXY CASE (DIRECT LATTICE)
50 CONTINUE
IF(LATTIC(2:3).EQ.'XZ') GOTO 60
IF(LATTIC(2:3).EQ.'YZ') GOTO 70
BR2(1,1)=0.5d0
BR2(1,2)=-0.5d0
BR2(1,3)=0.0d0
BR2(2,1)=0.5d0
BR2(2,2)=0.5d0
BR2(2,3)=0.0d0
BR2(3,1)=0.0d0
BR2(3,2)=0.0d0
BR2(3,3)=1.0d0
ORTHO=.TRUE.
GOTO 100
!
!.....CXZ CASE (DIRECT LATTICE)
60 CONTINUE
!.....CXZ ORTHOROMBIC CASE
if(gamma.eq.90.d0) then
BR2(1,1)=0.5d0
BR2(1,2)=0.0d0
BR2(1,3)=-0.5d0
BR2(2,1)=0.0d0
BR2(2,2)=1.0d0
BR2(2,3)=0.0d0
BR2(3,1)=0.5d0
BR2(3,2)=0.0d0
BR2(3,3)=0.5d0
ORTHO=.TRUE.
GOTO 100
ELSE
!.....CXZ MONOCLINIC CASE
write(*,*) 'gamma not equal 90'
SINAB=SIN(gamma1)
COSAB=COS(gamma1)
!
BR2(1,1)=0.5d0*sinab
BR2(1,2)=0.5d0*cosab
BR2(1,3)=-0.5d0
BR2(2,1)=0.0d0
BR2(2,2)=1.0d0
BR2(2,3)=0.0d0
BR2(3,1)=0.5d0*sinab
BR2(3,2)=0.5d0*cosab
BR2(3,3)=0.5d0
ORTHO=.FALSE.
GOTO 100
ENDIF
!
!.....CYZ CASE (DIRECT LATTICE)
70 CONTINUE
BR2(1,1)=1.0d0
BR2(1,2)=0.0d0
BR2(1,3)=0.0d0
BR2(2,1)=0.0d0
BR2(2,2)=0.5d0
BR2(2,3)=0.5d0
BR2(3,1)=0.0d0
BR2(3,2)=-0.5d0
BR2(3,3)=0.5d0
ORTHO=.TRUE.
GOTO 100
!.....RHOMBOHEDRAL CASE
80 BR2(1,1)=1/2.d0/sqrt(3.d0)
BR2(1,2)=-1/2.d0
BR2(1,3)=1/3.d0
BR2(2,1)=1/2.d0/SQRT(3.d0)
BR2(2,2)=1*0.5d0
BR2(2,3)=1/3.d0
BR2(3,1)=-1/SQRT(3.d0)
BR2(3,2)=0.d0
BR2(3,3)=1/3.d0
ORTHO=.FALSE.
GOTO 100
!
100 CONTINUE
999 format(3f15.5)
!
cosgam=cos(alpha(3)/180.d0*pi)
singam=sin(alpha(3)/180.d0*pi)
IF (.not.ortho) THEN
help(1)=dif(1)
help(2)=dif(2)
help(3)=dif(3)
if(lattic(1:1).eq.'R') then
dif(1)=help(1)*BR2(1,1)+help(2)*BR2(2,1)+help(3)*BR2(3,1)
dif(2)=help(1)*BR2(1,2)+help(2)*BR2(2,2)+help(3)*BR2(3,2)
dif(3)=help(1)*BR2(1,3)+help(2)*BR2(2,3)+help(3)*BR2(3,3)
elseif(lattic(1:3).eq.'CXZ') then
dif(1)=help(1)*singam
dif(2)=(help(1)*cosgam*alat(1)+help(2)*alat(2))/alat(2)
dif(3)=help(3)
else
dif(1)=(help(1)*BR2(1,1)*alat(1)+help(2)*BR2(2,1)*alat(2)+ &
help(3)*BR2(3,1)*alat(3)) ! /alat(1)
dif(2)=(help(1)*BR2(1,2)*alat(1)+help(2)*BR2(2,2)*alat(2)+ &
help(3)*BR2(3,2)*alat(3)) ! /alat(2)
dif(3)=(help(1)*BR2(1,3)*alat(1)+help(2)*BR2(2,3)*alat(2)+ &
help(3)*BR2(3,3)*alat(3)) ! /alat(3)
endif
else
dif(1)=xms(1)*alat(1)
dif(2)=xms(2)*alat(2)
dif(3)=xms(3)*alat(3)
endif
dist=0.d0
DO L=1,3
DIST=DIST+dif(L)*dif(L)
enddo
DIST=SQRT(DIST)
THETA=ACOS(dif(3)/dist)
XX=SQRT(dif(1)**2+dif(2)**2)
IF (XX.LT.1D-5) THEN
PHI=0.D0
ELSE
PHI=ACOS(dif(1)/XX)
IF (ABS(dif(2)).GT.1D-5) PHI=PHI*dif(2)/ABS(dif(2))
END IF
write(6,*) 'THETA, PHI', theta,phi
RETURN
END
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
