[Wien] lapwdm not working (anymore)?
Peter Blaha
peter.blaha at tuwien.ac.at
Tue Nov 26 17:06:38 CET 2024
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 at 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 at zeus.theochem.tuwien.ac.at/index.html
>
--
-----------------------------------------------------------------------
Peter Blaha, Inst. f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-158801165300
Email: peter.blaha at tuwien.ac.at
WWW: http://www.imc.tuwien.ac.at WIEN2k: http://www.wien2k.at
-------------------------------------------------------------------------
-------------- next part --------------
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
More information about the Wien
mailing list