We recently fixed one problem in SRC_symmetso connected with the
reduction of symmetry and splitting of equivalent atoms into
no-equivalent ones.
Replace clmchange.f with the attached file and recompile.
Regards
On 05/18/2015 01:36 PM, Vojtech Chlan wrote:
Dear WIEN2k community,
I am facing a problem with disturbance of convergence when the symmetry
of the structure is lowered (e.g., by initso_lapw).
It is well known in spin-polarized calculations that introducing the
spin-orbit interaction may reduce symmetry - in dependence on whether
the chosen direction of magnetization is compatible with present
elements of symetry or not. In spin-polarized cases where the symmetry
is not reduced, e.g. uniaxial structure with magnetization parallel to
the axis, the process is usually quite smooth: After well converged
runsp_lapw and initso call, the subsequent "runsp_lapw -so" calculation
starts almost converged and usually converges within a few iterations
(at least when no really heavy elements are present so that the
spin-orbit coupling is rather weak).
However in cases where the symmetry is reduced during initso (symmetso),
the continuation by runsp_lapw -so is not so smooth. According to my
experience there is often a substantial jump in the charge distance, the
magnetic moment may start to collapse etc. In fact, the jump in
convergence is present regardless the -so switch (or other parameters
that are usually changed during initso, e.g., Emax).
To illustrate this behaviour I did a few tests on barium hexagonal
ferrite (SG #194, P63/mmc):
Firstly, it was fully converged with quite standard parameters - using
PBE-GGA+U (U=4.5eV, J=0), with RKM=6.0, 100 k-points; wien version 14.2
was used.
Now, since there is the hexagonal axis (direction 001), starting initso
and setting the magnetization in 001 naturally does not reduce symmetry.
As expected, everything is nice and smooth when one starts runsp_lapw
-so .... there is only :DIS = 0.015 in the first iteration and the
calculation converges within a few iterations.
But when I set magnetization in direction 100, which kills half of the
symmetry operations (and 11 non-equivalent atoms become 15), the first
iteration starts with a jump in :DIS = 2.31. The -so switch is
irrelevant, the jump is there even for runsp_lapw without the -so switch.
My (naive) understanding of the origin of this jump is that it arises
from the change of basis for those atom sorts that became nonequivalent.
The inclusion of magnetization reduces also the local point group
symmetries of atoms (also often accompanied by change in rotation
matrices), which sumsequently changes their lists of LM expansions in
case.in2 file. The increase of LMs in expansion then manifests after
first iteration also in CLM files and one gets a jump in convergence.
When such changes concern only a few atoms in the structure or the
change in basis is small, it seems the calculation can often be
converged despite the jump. However when more and more atoms have their
expansions changed the jump becomes higher, eventually, for large
structures and in cases where the magnetization strips the structure of
almost all symmetries the jump becomes irrecoverable: the calculations
(with or without -so) crash typically in second iteration (when the new
LMs are first mixed) or even in first (in SELECT), in some cases the
runs survive without a crash but the potentials go all crazy and for
example magnetic moments collapse (anyway the convergence fails). I have
partially learned to live with that and partially learned to circumvent
this by reducing the symmetry not all at once, but in steps (and
re-converge in between) and by using other tricks to avoid crash and
maintain convergence. However, currently I try to switch on the
spin-orbit in a system too large where this simply does not help.
I can be all wrong, blaming the change in LMs, but I could not find any
other cause for the jump. And I believe I cannot touch the LM expansions
as they are given by the point group symmetry of the site.
I made a few naive attempts to fool some routines (mixer, clmextrapol)
into translating the clm files into ones with a new set of LMs, but
without proper knowledge of what I was doing, I naturally only ended up
with nonsenses and segmentation faults.
Is there any possibility to "smoothly renormalize" the density for a
changed set of LMs?
Well, perhaps I am blind to some completely different and obvious
solution, so any help would be appreciated.
Best regards
Vojtech
_______________________________________________
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
--
P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300 FAX: +43-1-58801-165982
Email: bl...@theochem.tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/staff/tc_group_e.php
--------------------------------------------------------------------------
subroutine clmchange(ifile,ifw,nat,ieq,jrj,iord,iz,tau,ksym)
implicit real*8(a-h,o-z)
INCLUDE 'param.inc'
complex*16 imag,ro2
integer iz(3,3,48),ksym(48),ksym2(48),kv(3),kvtest(3,48)
dimension TAU(3,noper)
! this program serves to read clm/potential files in FLAPW WIEN program
! then to construct the new files in the enlarged part of BZ
! when s-o coupling and magnetization are present
! jri(nat) is number of radial points for atom nat
CHARACTER*19 nkktext
character a(111)
dimension jrj(*),ieq(*)
character, allocatable :: al(:,:)
real*8,allocatable ::rho(:,:),ro1(:,:)
integer,allocatable :: lmmax(:),kv1(:,:)
allocate(lmmax(nat))
!
imag=(0.d0,1.d0)
tpi=2.d0*acos(-1.d0)
! read(5,*)nat,(jri(j),j=1,nat)
! read(5,*)cclm1
333 format(' number of atoms',i3,'number of points',2i7)
444 format(' file 1 vsn*',f5.1)
read(ifile,1980,end=999,err=999)(a(i),i=1,79)
a(78)='s'
a(79)='o'
write(6,*) 'Generating clm file',ifile
write(ifw,1980)(a(i),i=1,79)
read(ifile,1980)(a(i),i=1,79)
write(ifw,1980)(a(i),i=1,79)
read(ifile,1980)a(1)
write(ifw,1980)a(1)
jatso=0
do 79 jatom=1,nat
do 80 jeq=1,ieq(jatom)
if(jeq.gt.1) write(6,*) 'duplicating atom',jatom
jatso=jatso+1
if(jeq.eq.1) read(ifile,2000)(a(i),i=1,15),j
write(ifw,1990)jatso
1990 FORMAT(3X,'ATOM NUMBER=',I3,5X,10A4)
if(jeq.eq.1)then
read(ifile,2000)(a(i),i=1,15),lmmax(jatom)
allocate (al(lmmax(jatom),111),rho(lmmax(jatom),nrad))
endif
write(ifw ,2111)lmmax(jatom)
2111 FORMAT(3X,'NUMBER OF LM',I3)
if(jeq.eq.1)then
do 35 lm1=1,lmmax(jatom)
read(ifile,2000)a(1)
read(ifile,2000)a(1)
read(ifile,1980)(al(lm1,i),i=1,25)
read(ifile,1980)a(1)
read(ifile,2020)(rho(lm1,j),j=1,jrj(jatom))
35 continue
endif
do 36 lm1=1,lmmax(jatom)
write(ifw,2000)a(1)
write(ifw,2000)a(1)
write(ifw ,1980)(al(lm1,i),i=1,25)
write(ifw, 1980)a(1)
write(ifw,2020)(rho(lm1,j),j=1,jrj(jatom))
36 continue
do 88 i=1,6
if(jeq.eq.1)read(ifile,1980)a(1)
write(ifw,1980)a(1)
88 continue
80 continue
deallocate (al,rho)
79 continue
if(ifile.le.26) then
deallocate(lmmax)
return
endif
read(ifile,1980)(a(i),i=1,79)
write(ifw,1980)(a(i),i=1,79)
read(ifile,1980)a(1)
write(ifw,1980)a(1)
!
! read(ifile,2030)(a(i),i=1,13),nwav1
read(ifile,'(a19)') nkktext
read(nkktext,'(9x,i10)',err=6767) nwav1
goto 6768
6767 read(nkktext,'(13x,i6)') nwav1
6768 continue
allocate (ro1(2,nwav1*iord),kv1(3,nwav1*iord))
do 140 j=1,nwav1
140 read(ifile,2070)(kv1(jx,j),jx=1,3),(ro1(k,j),k=1,2)
345 format(i5,2e15.7)
2070 format(3x,3i5,2es19.12)
2030 format(13a1,i6)
1980 format(80a1)
2000 format(15a1,i3)
2020 format(3x,4es19.12)
! **** clm file read *****
nw=nwav1
write(6,*)' nw k j knew rhonew kold'
! print *, nwav1,iord,(ksym(i),i=1,iord)
do 20 j=1,nwav1
jsplit=1
! do 20 j=1,10
do 22 k=1,iord
ksym2(k)=ksym(k)
do 22 i1=1,3
kvtest(i1,k)=0
do 22 i2=1,3
kvtest(i1,k)=kvtest(i1,k)+iz(i1,i2,k)*kv1(i2,j)
22 continue
do 23 k1=1,iord
if(ksym(k1).eq.0)then
do k2=1,iord
if(ksym2(k2).eq.1) then
if(kvtest(1,k1).eq.kvtest(1,k2).and. &
kvtest(2,k1).eq.kvtest(2,k2).and. &
kvtest(3,k1).eq.kvtest(3,k2)) goto 23
endif
enddo
! vector kv is new, generate its little-star and block equal in ksym2
do k=1,iord
if(ksym(k).eq.1) then
do i1=1,3
kv(i1)=0
do i2=1,3
kv(i1)=kv(i1)+iz(i1,i2,k)*kvtest(i2,k1)
enddo
enddo
do k2=1,iord
if(kv(1).eq.kvtest(1,k2).and. &
kv(2).eq.kvtest(2,k2).and. &
kv(3).eq.kvtest(3,k2)) ksym2(k2)=1
enddo
endif
enddo
nw=nw+1
jsplit=jsplit+1
TK = 0.0D+0
DO l = 1, 3
TK = TK + TAU(l,k1)*kv1(l,J)*TPI
enddo
do 25 l=1,3
kv1(l,nw)=kvtest(l,k1)
25 continue
ro2=cmplx(ro1(1,j),ro1(2,j))* EXP(DCMPLX(TK)*(-IMAG))
ro1(1,nw)=dble(ro2)
ro1(2,nw)=aimag(ro2)
if(ifw.eq.55) write(6,500)nw,j,k,(kv1(l,nw),l=1,3), &
(ro1(l,nw),l=1,2),(kv1(l,j),l=1,3)
500 format(3i6,3x,3i4,2f10.6,3i4)
endif
23 continue
ro1(1:2,nw-jsplit+2:nw)=ro1(1:2,nw-jsplit+2:nw)/jsplit
ro1(1:2,j)=ro1(1:2,j)/jsplit
21 continue
20 continue
write(ifw,2031) nw
2031 format(' ',i10,' NUMBER OF PW')
do j=1,nw
write(ifw,2070)(kv1(jx,j),jx=1,3),(ro1(k,j),k=1,2)
enddo
deallocate(lmmax)
deallocate (ro1,kv1)
return
!
999 write(6,*) ' Warning: clmfile',ifile,' skipped'
deallocate(lmmax)
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
