Bug-fix (in Japanease) for SiF4 and so

Posted on : May 30, 2003 (Fri) 13:15:46

by H.akai

3月以来、KKR-CPAに手をつける暇が全く無く、Iwasawaさんのご指摘のSiF4でatrmrot
がエラーを出す現象もチェックすることができませんでした。今日、少し時間があっ
たので、チェックしてみましたところ、アルゴリズムに思わぬ落とし穴がありまし
た。結晶の点群をきめるのに、結晶のgeometrical form factorを計算して、回転した
結晶のgeometrical form factorがもとの結晶のそれと一致するかで、その回転対称性
をもつか否か判定しているのですが（これは結晶の　回転＋並進　のうち並進部分を
除くためです）　geometrical form factorの絶対値だけを使って判定するために（つ
まり位相部分を落とす）空間反転したものとしないものが区別ができません。そのよ
うな理由でエラーを出しておりました。atmrot.fを修正いたしました。数行ほどの変
更ですが、少なくとも上記の問題は取り除かれております。atmrot.fを添付されてい
るものと置き換えて使用下さい。次回のバージョンアップでは修正版をアップロード
いたします。
大阪大学大学院理学研究科
赤井久純


At 14:38 14/03/03 +0900, you wrote:
>/0071/ 2003/03/14(Fri) 14:38:42
>Misako Iwasawa
>--
>初めまして。
>SiF4の入力データを作ろうと試みています。
>ナーバルリサーチの結晶構造のデータを参考に、
>
>bccにして、
>
> 0.00000a 0.00000b 0.00000c Si
> 0.33000a 0.33000b 0.33000c F
>-0.33000a 0.00000b 0.00000c F
> 0.00000a -0.33000b 0.00000c F
> 0.00000a 0.00000b -0.33000c F
>
>としましたが、
>
>***err in atmrot...procedure fails
>
>となってしまいました。どこが間違っているのでしょうか？
>ご教示下さい。
>


-------------------- 以下は atmrot.f の修正版
--------------------------------
subroutine atmrot(irotat,natm,g,atmicp,itype,protat,isymop,wk
& ,gsf,ftype,gpt,ngpt)
c-----------------------------------------------------------------------
c Construct the atomic rotation matrix. It also gives a table
c describing the rotation of atoms in the unit cell. For example,
c if lineary aligned five atoms 1,2,3,4,5 change their position
c such that it looks like 4,5,1,2,3 after the rotation, the vector
c 'irot' has values (4,5,1,2,3). If the operation iop map the
c atomic arrangement to itself, isymop(iop) is set to 1. Also If it
c does not but map to its inversion, isymop(iop) is set to -1.
c Otherwise, isymop(iop) is set to 0.
c Since the procedure performed here is a bit complicated, the
c style of the program seems to be rather tedious, namely mainly
c constructed with use of "go to" statement. This however is
c the most transparent way to realize the present algorithm.
c Coded by H.Akai, 23 Dec. 96, Osaka
c Revised on 7 Aug. 1997, Duisburg
c Revised on 24 May 2003. Osaka
c-----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
parameter(nwkx=20)
real*8 protat(9,24),atmicp(3,natm),g(3,3),wk(3,natm),c(3)
& ,gpt(3,ngpt),ftype(natm),dsplmt(3),gsf(ngpt)
integer itype(natm),irotat(natm,24),isymop(24),iwk(nwkx)
logical eqvlat,str
c
c write(*,'(/1x,a)')'iop g/u irotat'
do 10 iop=1,24
c --- check only the symmetry operations that are compatible with
c those of the Bravais lattice.
if(isymop(iop) .eq. 0) go to 10
c --- first, reset "isymop" that will be later set when the local
c symmetry is compatible with the symmetry operation "iop".
ind=isymop(iop)
isymop(iop)=0
c --- rotate the cartesian cordinates of all the atoms.
call rotatm(atmicp,wk,protat(1,iop),natm)
str=eqvlat(wk,natm,gpt,ftype,gsf,ngpt)
if(.not. str) go to 10
c write(*,'(1x,a,i2)')'iop=',iop
c write(*,'((1x,9f7.3))')((wk(i,j),i=1,3),j=1,natm)
do 20 ip=0,1
c --- if ind=-1, proper rotation not allowed.
if(ip .eq. 0 .and. ind .eq. -1) go to 20
p=(-1d0)**ip
c --- since the inversion is not included in the above rotation
c it must be considered separately.
do 30 j=1,natm
do 30 i=1,3
30 wk(i,j)=p*wk(i,j)
c write(*,*)'iop,ip,str=',iop,ip,str
c write(*,'((1x,9f7.3))')((wk(i,j),i=1,3),j=1,natm)
c --- for rotated ia-th atom, check if any unrotated atoms come
c to the same position.
c --- First, we determine the shift of the coordinates.
ity=0
do 40 ic=1,natm
it=itype(ic)
do 50 i=1,ity
if(it .eq. iwk(i)) go to 40
50 continue
c write(*,'(1x,a,2i2)')'ic, it=',ic,it
ity=ity+1
if(ity .gt. nwkx) call errtrp(1,'atmrot','iwk overflows')
iwk(ity)=it
call clrari(irotat(1,iop),natm)
do 60 id=ic,natm
if(itype(id) .ne. it) go to 60
c --- try (ic,id) pair.
irotat(ic,iop)=id
do 70 i=1,3
70 dsplmt(i)=wk(i,id)-atmicp(i,ic)
c --- now dsplacement compatible to (ic,id) pair is set.
c write(*,'(1x,a,i2,a,3f10.5)')'target=',id,' dsp=',dsplmt
c --- all atoms for which ia do 80 ia=ic+1,natm
if(itype(ia) .ne. it) go to 80
do 90 ib=ic,natm
if(ib .eq. id .or. itype(ib) .ne. it) go to 90
c --- for this displacement, (ia,ib) pair is examined.
do 100 j=1,3
100 c(j)=0d0
do 110 i=1,3
d=wk(i,ib)-atmicp(i,ia)-dsplmt(i)
c --- inner product d*g gives coeffecients c's that satisfy
c d=c(1)*r(1)+c(2)*r(2)+c(3)*r(3), where r's are the primitive
c lattice vectors.
do 110 j=1,3
110 c(j)=c(j)+d*g(i,j)
c write(*,'(1x,a,2i3,a,3f10.5)')'pair=',ia,ib,' shift=',c
do 120 i=1,3
c --- check if all c's are integer. if not, ia and ib-th atom do not
c sit at the same place.
icc=c(i)+10000.5d0
cc=dble(icc-10000)
c --- if any of c's is not integer, skip the remaining check and
c see the next atom of new ib.
if(abs(c(i)-cc) .gt. 1d-3) go to 90
120 continue
irotat(ia,iop)=ib
go to 80
90 continue
go to 60
80 continue
go to 130
60 continue
go to 20
130 continue
c --- Now, one of the plausible "dsplmt" is obtained. We check
c if this choice explains the position of all other atoms.
c If it is not the case, we should try different ones.
c write(*,'(1x,5(5i2,1x))')(irotat(i,iop),i=1,natm)
c write(*,'(1x,a,3f10.5)')'dsplmt',dsplmt
do 140 ia=1,natm
if(irotat(ia,iop) .ne. 0) go to 140
do 150 ib=1,natm
c --- (ia,ib) pair is now tried.
c --- check only when ia-th and ib-th atoms are of the same type.
if(itype(ia) .ne. itype(ib)) go to 150
do 160 j=1,3
160 c(j)=0d0
do 170 i=1,3
d=wk(i,ib)-atmicp(i,ia)-dsplmt(i)
c --- inner product d*g gives coeffecients c's that satisfy
c d=c(1)*r(1)+c(2)*r(2)+c(3)*r(3), where r's are the primitive
c lattice vectors.
do 170 j=1,3
170 c(j)=c(j)+d*g(i,j)
do 180 i=1,3
c --- check if all c's are integer. if not, ia and ib-th atom do not
c sit at the same place.
icc=c(i)+10000.5d0
cc=dble(icc-10000)
c --- if any of c's is not integer, skip the remaining check and
c see the next atom of new ib.
if(abs(c(i)-cc) .gt. 1d-3) go to 150
180 continue
c ia and ib-th atoms are sitting at the same place.
irotat(ia,iop)=ib
c --- since the atom sitting the same place is found, we can skip
c the remaining check for ib atom.
go to 140
150 continue
c call errtrp(2,'atmrot','no partner available')
go to 40
140 continue
c --- since all the atoms are found to have their partner,
c "isymop" is now set. This is done either direct or indirect
c rotation compatible with the symmetry.
isymop(iop)=(-1)**ip
c write(*,'(1x,2i3,3x,5(5i3,1x))')
c & iop,isymop(iop),(irotat(i,iop),i=1,natm)
go to 10
40 continue
20 continue
call errtrp(1,'atmrot','procedure fails')
10 continue
end
-------------------------------------　ここまで -------------------------------
----------------------------------
Hisazumi Akai
Department of Physics, Osaka University
1-1 Machikaneyama,Toyonaka
Osaka 560-0043, Japan

Tel: +81 6 6850 5738
Fax: +81 6 6850 5741
e-mail: akai@phys.sci.osaka-u.ac.jp