t-nissie/00AkaiKKR.md

## 00AkaiKKR.md

      
    Raw
  

              00AkaiKKR.md
            
          
    AkaiKKRについてのメモ

この文章のオリジナルは https://gist.github.com/t-nissie/4edcf8d62144b4d77f5c に置いてある。
使用するパソコンはIntel Core i5程度のノートパソコンで十分。
入手

cpa2002v009c.tgzのAugust 26, 2015版を
http://kkr.phys.sci.osaka-u.ac.jp/jp/
から入手した。
$ md5sum cpa2002v009c.tgz | grep 798c7f8175032d2f55137bed49d15409
MD5     (cpa2002v009c.tgz)   =   798c7f8175032d2f55137bed49d15409

コンパイル

gfortranでコンパイルできるようにパッチcpa2002v009c.patchをあてる。
tar xf cpa2002v009c.tgz
mv cpa2002v009c/ cpa2002v009c.2015-08-26
cd cpa2002v009c.2015-08-26/
patch -p1 -b < /SOMEWHERE/cpa2002v009c.patch
make
make gpd
make spc
cd util/
gfortran -o fmg fmg.f
cd ..

bcc鉄

bcc鉄のSCFの計算

./specx < in/fe >> data/fe.out; cat data/fe.out

bcc鉄のDOSの計算

sed -e 's/go/dos/' -e 's/    4    /    8    /' in/fe > in/fe.dos
./specx < in/fe.dos > data/fe.dos
cat data/fe.info
     5.2700       -2522.8176206   2.16492   ←SCF計算の正確な値
     5.2700       -2518.5344701   0.17646   ←DOSの計算の不正確な値
./gpd data/fe.dos

鉄の強磁性状態とlocal moment disorder (LMD) 状態とを比較して相転移温度Tcを求める

ポテンシャルの準備

まず、さきほどの鉄の計算を元にスピンが反対向きのサイト用のポテンシャルを用意する。
cd util
cat inputsample
./fmg < inputsample
cd ..

本番実行

diff -u in/fe in/fe_lmd
./specx < in/fe_lmd >> data/fe_lmd.out; cat data/fe_lmd.out
./specx < in/fe_lmd >> data/fe_lmd.out; cat data/fe_lmd.out   # 収束しなかったのでもう一回
cat data/fe_lmd.info                             # スピンモーメントが０になっていることを確認

平均場近似 Tc=(2/3)(E_LMD-E_Ferro)/k_B [K], where k_B=6.3336E-6 (Ry/K) で
(/ (* (- -2522.8055945 -2522.8176206) 0.66666666666666666) 6.3336E-6)
1265.851964099398

コバルトでもできるはず。
ニッケルは縦揺らぎ（スピンの伸び縮み）が大きいのであまり正確でないTcがでるかも。
逆？
fcc Ni_x Fe_1-xのDOSの計算

無限鉄結晶に１つだけ水素不純物がある状態

in/fe_hを使う。濃度を0%にして計算する。
cmd2% run "specx < in/fe_h"
cmd2% run "specx < in/fe_h"
cmd2% run "specx < in/fe_h"   # 10e-6に収束するまで計算する。
             :
                             *** type-Fe       H  (z=  1.0) ***
   core charge in the muffin-tin sphere = 0.0000000
   valence charge in the cell (spin up  ) =     0.24232(s)   0.11719(p)   0.03840(d)
   valence charge in the cell (spin down) =     0.27799(s)   0.14372(p)   0.03678(d)
   total charge=   0.85640   valence charge (up/down)=   0.39791   0.45849
   spin moment=  -0.06057  orbital moment=   0.00000
             :

水素の位置をずらすこともできるらしい。
See H. Akai et al.: progress in theoretical physics Suppl. vol.101 p.11 (1990).
hcpコバルト

コバルトのSCF計算

in/co。hcpって指定するけど、本当は六方格子で、unit cell中に２つ原子をおくことに注意。
atmicxにa,b,cをつけ忘れないこと。
cmd2% run "specx < in/co"
cmd2% run "specx < in/co"
cmd2% run "specx < in/co"
cmd2% run "specx < in/co"

コバルトのDOS計算

in/co    編集
cmd2% run "specx < in/co > data/co.dos"
cmd2% ./gpd data/co.dos
【応用編】half-metallic half-Heusler alloy CrAs

Masafumi Shirai et al.: JJAP 39 L1118 (2000).
【応用編】リチウムイオン二次電池の起電力の計算

Positive electrode (LiCoO2): CoO2 + Li+ + e- → LiCoO2
Negative electrode (graphite): Li → Li+ + e-
Coは高価
2000, 6,000, 400 310 30-60 8000
% run "specx < in/licoo2 >> data/licoo2.out" ; grep err data/licoo2.out
% run "specx < in/licoo2 >> data/licoo2.out" ; grep err data/licoo2.out
% run "specx < in/licoo2 >> data/licoo2.out" ; grep err data/licoo2.out
% run "specx < in/licoo2 >> data/licoo2.out" ; grep err data/licoo2.out
% run "specx < in/licoo2 >> data/licoo2.out" ; grep err data/licoo2.out # 収束するまで繰り返す

LiCoO3: -3096.74712993 Ry
% run "specx < in/li >> data/li.out" ; grep err data/li.out
% run "specx < in/li >> data/li.out" ; grep err data/li.out
% run "specx < in/li >> data/li.out" ; grep err data/li.out

Li: -14.8370012
% run "specx < in/vccoo2 >> data/vccoo2.out" ; grep err data/vccoo2.out

VcCoO3: -3081.6365328 Ry
起電力は？
-V = Etot(LiCoO2) - Etot(CoO2) - Etot(Li)
= -3096.74712993 + 14.8370012 + 3081.6365328
= -0.27359 Ry
= -3.72 eV
比較するときはmxlは同じにする。
% diff -u in/licoo2 in/vccoo2
--- in/licoo2	2015-08-25 11:54:34.000000000 +0900
+++ in/vccoo2	2016-03-04 09:30:25.143318378 +0900
@@ -1,5 +1,5 @@
-c----------------------LiCoO2--------------------------------
-      go  data/licoo2
+c----------------------VcCoO2--------------------------------
+      go  data/vccoo2
 c------------------------------------------------------------
 c   brvtyp     a         c/a   b/a   alpha   beta   gamma
      rhb    9.373   ,      ,    ,   32.97 ,     ,    ,
@@ -14,8 +14,8 @@
       3
 c------------------------------------------------------------
 c   type     ncmp    rmt    field   mxl  anclr   conc
-    Li          1      1    0.000   2
-                                           3     100
+    Vc          1      1    0.000   2
+                                           0     100
     Co          1      1    0.000   2
                                           27     100
     O           1      1    0.000   2
@@ -26,7 +26,7 @@
 c------------------------------------------------------------
 c   atmicx                   atmtyp
      0           ,  0           ,  0       , Co
-     0.5a        ,  0.5b        ,  0.5c    , Li
+     0.5a        ,  0.5b        ,  0.5c    , Vc
      0.26a       ,  0.26b       ,  0.26c   , O
      0.74a       ,  0.74b       ,  0.74c   , O
 c------------------------------------------------------------
LiAlO2の計算

Oの2pを積分に入れるか入れないかを注意。
-798.8252139
100%
-5.90
10%
-3.86
GaAs

GaAsのSCF計算

./specx < in/gaas >> data/gaas.out; grep err data/gaas.out

GaAsのDOSの計算

下記のgp5d.fを使った。
sed -e 's/go/dos/' -e 's/    4    /    8    /' in/gaas > in/gaas.dos
./specx < in/gaas.dos >> data/gaas.dos
SOMEWHERE/gp5d data/gaas.dos
cp data/gaas.plt data/gaas.gp

data/gaas.gpを適宜編集してPDOSを描く。
grep -A 201 'DOS of component 1' gaas-nmag.dos | tail -n 201 | ruby -e 'load "interpolations.rb"; load "integrations.rb"; d = DataSet.new(STDIN,0,1); printf "%.3f\n", (-1.1250..0.0).simpson(100){|x| d.linear_interpolation(x)}'

Total 9.0
Ga    s 0.407  p 0.331  d 3.320
As    s 0.680  p 0.927  d 0.020

希薄磁性半導体 (DMS)

(GaMn)AsのSCF計算

(GaMn)AsのDOSの計算

Hyperfine

keyword: Fermi contact interaction
【付録】gp5d.f

gp5d.f は gpd.f replacement。特徴は

独自サブルーチンlftfil()を標準内部関数に置き換え: subroutine lftfil(file) → adjustl(file)
独自サブルーチンchleng()を標準内部関数に置き換え: subroutine chleng(a,ln) → len_trim(file)
インラインデータとしてGNUPLOT 5の新機能の$u << EOFの導入。PDOSのプロットに便利。
不正でないEPSファイルの出力: set terminal postscript landscape → set terminal postscript eps
PostScript標準フォント（欧文基本14フォント）への置き換え: arial → Times-Roman
シバンの書き換え: #!/usr/bin/gnuplot→ #!/usr/bin/env gnuplot
-*-Gnuplot-*- for Emacs


## cpa2002v009c.patch
diff -ur cpa2002v009c.downloaded/in/fe cpa2002v009c/in/fe
--- cpa2002v009c.downloaded/in/fe	2015-08-25 12:37:17.000000000 +0900
+++ cpa2002v009c/in/fe	2016-03-02 13:07:09.000000000 +0900
@@ -1,5 +1,5 @@
 c----------------------Fe------------------------------------
-    spc   data/fe
+    go   data/fe
 c------------------------------------------------------------
 c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     bcc      5.27  ,     ,      ,      ,       ,      ,
@@ -38,13 +38,12 @@
 c  "end" is not necessary unless other data used only
 c   for "spc" exist below.
 #end
-
 c  the following lines are used to get energy dispersion curves
- 300
- 0 0 0
- 0 1 0
- 0.5 0.5 0
- 0.5 0.5 0.5
- 0 0 0
- 0.5 0.5 0
+# 300
+# 0 0 0
+# 0 1 0
+# 0.5 0.5 0
+# 0.5 0.5 0.5
+# 0 0 0
+# 0.5 0.5 0

diff -ur cpa2002v009c.downloaded/makefile cpa2002v009c/makefile
--- cpa2002v009c.downloaded/makefile	2015-08-25 11:58:18.000000000 +0900
+++ cpa2002v009c/makefile	2016-03-03 17:11:42.000000000 +0900
@@ -13,18 +13,18 @@
 ##########################
 # In the case of ifort
 ##########################
-fort = ifort
-#flag = -mp1 -i-dynamic -mcmodel=medium
-flag = -O2 -mp1 -mcmodel=medium
-omp = -openmp
-nomp = -openmp-stubs
+# fort = ifort
+# #flag = -mp1 -i-dynamic -mcmodel=medium
+# flag = -O2 -mp1 -mcmodel=medium
+# omp = -openmp
+# nomp = -openmp-stubs
 #
 ##########################
 # In the case of gfortran
 ##########################
-#fort = gfortran
-#flag = -O2
-#omp = -fopenmp
+fort = gfortran
+flag = -O2 -Wall
+omp = -fopenmp
 #omp =
 #
 program = specx
diff -ur cpa2002v009c.downloaded/source/spmain.f cpa2002v009c/source/spmain.f
--- cpa2002v009c.downloaded/source/spmain.f	2015-08-26 17:53:59.000000000 +0900
+++ cpa2002v009c/source/spmain.f	2015-08-25 11:57:25.000000000 +0900
@@ -181,10 +181,7 @@
    12 write(*,1310)(atmicp(l,j),l=1,3),atmtyp(j)
  1310 format('   position=',3f13.8,'  type=',a8)
       if(ids .eq. 5) then
-      call xtoken(token,*18)
-      go to 19
-   18 call errtrp(1,'spmain','fspin data not found')
-   19 read(token,'(g15.0)')fspin
+      read(*,*)fspin
       write(*,'(1x)')
       write(*,'(a,f13.6)')'   fixed spin moment=',fspin
       endif

## gp5d.f
      program gp5d
c-----------------------------------------------------------------------
c     construct DOS data from cpa98 output.
c     coded by H.Akai, 22 Feb. 1998, Osaka
c     modifed by H. Akai, 27 Nov. 2010
c     modifed by T. Nishimatsu, 7 Mar. 2016
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      parameter(msex=201,icmpx=10,maxlin=10000,mmxl=4)
      real*8 dat(mmxl*icmpx+3,msex,2)
      character file*256,a*120,fmt*80,fmtx*80,magtyp*4
      logical exist
      data big/1d30/
c--- get input file
      file=' '
      call getarg(1,file)
      if(file .eq. ' ') then
      write(*,'(a)')'Usage: gp5d file'
      stop
      endif
      file = adjustl(file)
c
c     If a directory of the same name exists, inquire cannot
c     ditect existence of a file correctly. Therefore one should
c     use 'open' statement together with error return.
      inquire(file=file,exist=exist)
      if(.not. exist) then
      write(*,'(a/a)')'File not found','Usage: gp5d file'
      stop
      endif
      open(12,file=file)
c---  create output file
      ip=index(file,'.')
      if(ip .ge. 2) file=file(1:ip-1)
      file = adjustl(file)
      if(len_trim(file) .le. 0) then
      write(*,'(a/a)')'Illegal file name','Usage: gp5d file'
      stop
      endif
c
c--- read-in files (unit=12) and find DOS data.
c
      mse=0
      mxl=0
      icmp=0
      magtyp=' '
      do 50 line=1,maxlin
c     --- get the number of meshes, l_max, and number of components
      read(12,'(a)',end=10)a
      if(a(5:28) .eq. 'meshr  mse    ng     mxl') then
      read (12,*)meshr,mse,ng,mxl
      endif
      if(a(4:8) .eq. 'ntyp='.and. a(13:17) .eq. 'natm='
     &   .and. a(22:27) .eq. 'ncmpx=') then
      read(a(28:29),'(i2)')icmp
      endif
c     --- get magtyp
      if(a(68:74) .eq. 'magtyp=') then
      magtyp=a(75:78)
      endif
c     --- the last data before the DOS data begin are adopted
      if(a(2:19) .eq. 'DOS of component 1') go to 100
   50 continue
      write(*,'(a)')'Seems not to contain DOS data.'
      stop
c 100 write(*,'(a,3i3)')' mse,mxl,icmp=',mse,mxl,icmp
  100 continue
c     --- check the sizes ---
      if(mse .gt. msex) then
      write(*,'(a,i3)')' ***error...mse>',msex
      stop
      endif
      if(mxl .gt. mmxl) then
      write(*,'(a,i2)')' ***error...mxl>',mmxl
      stop
      endif
      if(icmp .gt. icmpx) then
      write(*,'(a,i2)')' ***error...icmp>',icmpx
      stop
      endif
      if(magtyp .ne. 'mag' .and. magtyp .ne. 'nmag') then
      write(*,'(a,a4)')' ***error...illegal magtyp=',magtyp
      stop
      endif
      if(magtyp .eq. 'mag') then
      ismx=2
      else
      ismx=1
      endif
c     write(*,*)'ismx=',ismx
c     --- Now the format statement can be constructed
      write(fmtx,'(a,i1,a)')'(1x,f7.4,3x,',mxl,'f10.4)'
      write(fmt,'(a,i2,a)')'((f8.4,f7.2,f8.4,',mxl*icmp,'f7.2))'
c     write(*,*)'format=',fmtx
c     --- read in DOS data
      do 20 is=1,ismx
      do 30 line=1,maxlin
      if(a(2:19) .eq. 'DOS of component 1') then
c     write(*,*)'DOS           1 met'
      read(12,fmtx)((dat(l,k,is),l=3,mxl+3),k=1,mse)
c     write(*,fmtx)((dat(l,k),l=3,mxl+3,is),k=1,mse)
      do 40 i=2,icmp
c     write(*,*)'DOS',i,' met'
      i0=mxl*(i-1)+4
      read(12,'(//)')
   40 read(12,fmtx)
     &    (dmy,(dat(l,k,is),l=i0,i0+mxl-1),k=1,mse)
      endif
      if(a(2:10) .eq. 'total DOS') then
c     write(*,*)'total DOS met'
      read(12,'(1x,f12.7,3x,f13.5)')(dat(1,k,is),dat(2,k,is),k=1,mse-1)
      dat(1,mse,is)=dat(1,mse-1,is)
      dat(2,mse,is)=dat(2,mse-1,is)
      go to 20
      endif
   30 read(12,'(a)',end=10)a
   20 read(12,'(a)',end=10)a
      if(ismx .eq. 1) then
      do 60 i=1,mxl+32
      do 60 k=1,mse
   60 dat(i,k,2)=dat(i,k,1)
      endif
c     --- fix axis design
      xmi=big
      xmx=-big
      ymi=big
      ymx=-big
      do 70 is=1,2
      do 70 k=1,mse
      xmi=min(xmi,dat(1,k,is))
      xmx=max(xmx,dat(1,k,is))
      ymi=min(ymi,dat(2,k,is))
   70 ymx=max(ymx,dat(2,k,is))
      ymi=0d0
c     call ezaxis(xmi,xmx,nex,ntickx,1d0)
      call ezaxis(ymi,ymx,ney,nticky,0d0)
c     xmi=xmi*10d0**nex
c     xmx=xmx*10d0**nex
c     xincr=(xmx-xmi)/dble(ntickx-1)
      ymi=ymi*10d0**ney
      ymx=ymx*10d0**ney
      yincr=(ymx-ymi)/dble(nticky-1)
c     --- construct a gnuplot program
      open(13,file=trim(file)//'.plt',status='unknown')
      write(13,'(a)')'#!/usr/bin/env gnuplot'
      write(13,'(a)')'#-*-Gnuplot-*-'
      write(13,'(a)')'##'
      write(13,'(a)')'$u << EOF'
      write(13,fmt)((dat(l,k,1),l=1,mxl*icmp+3),k=1,mse)
      write(13,'(a)')'EOF'
      write(13,'(a)')'$d << EOF'
      write(13,fmt)((dat(l,k,2),l=1,mxl*icmp+3),k=1,mse)
      write(13,'(a)')'EOF'
      write(13,'(2a)')'set terminal postscript eps noenhanced color',
     &               ' "Helvetica" 20'
      write(13,'(a)')'set output "'//trim(file)//'.eps"'
c     write(13,'(a)')'set term x11'
      write(13,'(a)')'#set yzeroaxis lt -1 lw 0.1'
c     write(13,'(a)')'#set yzeroaxis lt -1 lw 0.1'
      write(13,'(a)')'set border 15 lw 0.1'
      write(13,'(a)')'#set mxtics 2'
      write(13,'(a)')'#set mytics 2'
c     write(13,'(a,e15.7,a,e15.7)')'set xtics',xmi,',',xincr
      write(13,'(a,e15.7,a,e15.7)')'set ytics',ymi,',',yincr
      write(13,'(a)')'set size 1,1'
      write(13,'(a)')'set origin 0,0'
      write(13,'(a)')'set multiplot'
      write(13,'(a)')'set size 1.0,0.35'
c     write(13,'(a)')'set border 14 lw 0.1'
c     write(13,'(a)')'set lmargin at screen 0.1'
c     write(13,'(a)')'set rmargin at screen 0.9'
      write(13,'(a)')'set origin 0.0,0.5'
      write(13,'(a)')'set lmargin 8'
      write(13,'(a)')'set rmargin 3'
      write(13,'(a)')'set tmargin 0'
      write(13,'(a)')'set bmargin 0'
      write(13,'(3a)')'set title "',trim(file),'" font "Times-Roman,20"'
      write(13,'(a)')'set format x ""'
c     write(13,'(a,e15.7,a,e15.7,a)')'set xrange [',xmi,':',xmx,']'
      write(13,'(a,e15.7,a,e15.7,a)')'set xrange [*:*]'
      write(13,'(a,e15.7,a,e15.7,a)')'set yrange [',ymi,':',ymx,']'
      write(13,'(a)')'plot "$u" using 1:2 title "" with lines lt 1 lw 3'
      write(13,'(a)')'set border 11 lw 0.1'
      write(13,'(a)')'set origin 0.0,0.15'
      write(13,'(a)')'unset title'
      write(13,'(2a)')
     & 'set xlabel "Energy relative to Fermi energy (Ry)"',
     & ' font "Times-Roman,20"'
      write(13,'(2a)')
     & 'set label "DOS (1/Ry/unit cell)" at screen 0.04, screen 0.5',
     & ' center rotate font "Times-Roman,20"'
      write(13,'(a,e15.7,a,e15.7,a)')'set yrange [',ymx,':',ymi,']'
      write(13,'(a)')'set format x'
      write(13,'(a)')'plot "$d" using 1:2 title "" with lines lt 1 lw 3'
      write(13,'(a)')'unset multiplot'
      write(13,'(a)')'set output   # Close the EPS file.'
      write(13,'(a)')'!epstopdf '//trim(file)//'.eps'
c     write(13,'(a)')'pause -1'
c     write(13,'(a)')'#    EOF'
      close(13)
      call system('gnuplot '//trim(file)//'.plt')
      stop
   10 write(*,'(a/a)')'Illegal file specified','Usage: gp5d file'
      end program gp5d
      subroutine ezaxis(xmin,xmax,nex,ntick,adjust)
c-----------------------------------------------------------------------
c     I fix a scale of the axis.
c     In output, xmin is maximum integer for which  xmin*10**nex
c     does not exeed original input xmin.
c     Similary, xmax is maximum integer for which  xmin*10**nex
c     is not smaller than original xmax.
c     ntick gives a number of ticks including both end points.
c     Coded by H.Akai, April 28, 1991
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      dimension nw(12),nh(12)
      data nw/ 4, 5, 6, 7, 8,10,12,14,16,20,25,30/
     &    ,nh/ 4, 5, 6, 7, 4, 5, 6, 7, 8, 4, 5, 6/
     &    ,small/1d-3/
c---- I fix scale of the axis.
c---- in some cases containing x=0 axis will give better presntation.
      if(xmax .lt. xmin) then
      write(6,1000)
 1000 format('   ***err in ezaxis...xmin and xmax not consistent')
      stop
      endif
      if(xmax*xmin .gt. 0d0 .and.
     &    (xmax-xmin)/max(abs(xmax),abs(xmin)) .gt. adjust) then
      xmin=min(xmin,0d0)
      xmax=max(xmax,0d0)
      endif
      wx=xmax-xmin
c---- wx indicates x-data spreading.
      nex=int(log10(wx/2.9d0)+100d0+small)-100
      sx=10d0**nex
      xmin0=xmin/sx
      xmax0=xmax/sx
c---- now xmin and xmax is normalized by 10**nex so as to
c---- satisfy 2.9<xmax-xmin<29.
      xmin=dble(int(xmin0+small))
      xmax=dble(int(xmax0-small))
c---- bug easily comes in the following statement. never
c---- use xmin instead of xmin0 in the logical-if, otherwise
c---- the procedure will fail when -1<xmin<0.
c---- xmin is maximum integer for which xmin*10**nex does not exeed
c---- original xmin.
      if(xmin0 .lt. 0d0) xmin=xmin-1d0
c---- xmax is maximum integer for which xmax*10**nex is not smaller
c---- than original xmax.
      if(xmax0 .gt. 0d0) xmax=xmax+1d0
      wx=xmax-xmin
c---- I then choose appropriate design for tick construction.
      iwx=int(wx+small)
      do 30 i=1,12
      ii=i
      if(nw(i) .ge. iwx) go to 40
   30 continue
c---- i adjust positioning such that the center of the data spreading
c---- is not much distant from that of the axis.
c---- however, if one of xmin or xmax is located near zero, xmin=0 or
c---- xmax=0 is the arrangement of the highest priority.
   40 ws=dble(nw(ii))
      wm=ws-wx
      if(xmin .gt. -small .and. xmin .lt. wm+small) then
        wm=xmin
c---- xmin=0 is realized in this case.
      else if(xmax .lt. small .and. xmax .gt. -wm-small) then
        wm=xmax+wm
c---- xmax=0 is realized in this case.
      else
        if(xmax-xmax0 .gt. xmin0-xmin) wm=wm+1d0
        wm=dble(int(wm/2d0+small))
c---- otherwise, ticks are arranged so as to spread more or less
c---- on centered.
      endif
      xmin=xmin-wm
      xmax=xmin+ws
      ntick=nh(ii)+1
c---- the first and the last ticks of axis are fixed.
c---- now I can start plotting.
      return
      end
cLocal variables:
c  compile-command: "gfortran -g -o gp5d gp5d.f"
cEnd:
	diff -ur cpa2002v009c.downloaded/in/fe cpa2002v009c/in/fe
	--- cpa2002v009c.downloaded/in/fe 2015-08-25 12:37:17.000000000 +0900
	+++ cpa2002v009c/in/fe 2016-03-02 13:07:09.000000000 +0900
	@@ -1,5 +1,5 @@
	c----------------------Fe------------------------------------
	- spc data/fe
	+ go data/fe
	c------------------------------------------------------------
	c brvtyp a c/a b/a alpha beta gamma
	bcc 5.27 , , , , , ,
	@@ -38,13 +38,12 @@
	c "end" is not necessary unless other data used only
	c for "spc" exist below.
	#end
	-
	c the following lines are used to get energy dispersion curves
	- 300
	- 0 0 0
	- 0 1 0
	- 0.5 0.5 0
	- 0.5 0.5 0.5
	- 0 0 0
	- 0.5 0.5 0
	+# 300
	+# 0 0 0
	+# 0 1 0
	+# 0.5 0.5 0
	+# 0.5 0.5 0.5
	+# 0 0 0
	+# 0.5 0.5 0

	diff -ur cpa2002v009c.downloaded/makefile cpa2002v009c/makefile
	--- cpa2002v009c.downloaded/makefile 2015-08-25 11:58:18.000000000 +0900
	+++ cpa2002v009c/makefile 2016-03-03 17:11:42.000000000 +0900
	@@ -13,18 +13,18 @@
	##########################
	# In the case of ifort
	##########################
	-fort = ifort
	-#flag = -mp1 -i-dynamic -mcmodel=medium
	-flag = -O2 -mp1 -mcmodel=medium
	-omp = -openmp
	-nomp = -openmp-stubs
	+# fort = ifort
	+# #flag = -mp1 -i-dynamic -mcmodel=medium
	+# flag = -O2 -mp1 -mcmodel=medium
	+# omp = -openmp
	+# nomp = -openmp-stubs
	#
	##########################
	# In the case of gfortran
	##########################
	-#fort = gfortran
	-#flag = -O2
	-#omp = -fopenmp
	+fort = gfortran
	+flag = -O2 -Wall
	+omp = -fopenmp
	#omp =
	#
	program = specx
	diff -ur cpa2002v009c.downloaded/source/spmain.f cpa2002v009c/source/spmain.f
	--- cpa2002v009c.downloaded/source/spmain.f 2015-08-26 17:53:59.000000000 +0900
	+++ cpa2002v009c/source/spmain.f 2015-08-25 11:57:25.000000000 +0900
	@@ -181,10 +181,7 @@
	12 write(*,1310)(atmicp(l,j),l=1,3),atmtyp(j)
	1310 format(' position=',3f13.8,' type=',a8)
	if(ids .eq. 5) then
	- call xtoken(token,*18)
	- go to 19
	- 18 call errtrp(1,'spmain','fspin data not found')
	- 19 read(token,'(g15.0)')fspin
	+ read(,)fspin
	write(*,'(1x)')
	write(*,'(a,f13.6)')' fixed spin moment=',fspin
	endif
	program gp5d
	c-----------------------------------------------------------------------
	c construct DOS data from cpa98 output.
	c coded by H.Akai, 22 Feb. 1998, Osaka
	c modifed by H. Akai, 27 Nov. 2010
	c modifed by T. Nishimatsu, 7 Mar. 2016
	c-----------------------------------------------------------------------
	implicit real*8 (a-h,o-z)
	parameter(msex=201,icmpx=10,maxlin=10000,mmxl=4)
	real8 dat(mmxlicmpx+3,msex,2)
	character file256,a120,fmt80,fmtx80,magtyp*4
	logical exist
	data big/1d30/
	c--- get input file
	file=' '
	call getarg(1,file)
	if(file .eq. ' ') then
	write(*,'(a)')'Usage: gp5d file'
	stop
	endif
	file = adjustl(file)
	c
	c If a directory of the same name exists, inquire cannot
	c ditect existence of a file correctly. Therefore one should
	c use 'open' statement together with error return.
	inquire(file=file,exist=exist)
	if(.not. exist) then
	write(*,'(a/a)')'File not found','Usage: gp5d file'
	stop
	endif
	open(12,file=file)
	c--- create output file
	ip=index(file,'.')
	if(ip .ge. 2) file=file(1:ip-1)
	file = adjustl(file)
	if(len_trim(file) .le. 0) then
	write(*,'(a/a)')'Illegal file name','Usage: gp5d file'
	stop
	endif
	c
	c--- read-in files (unit=12) and find DOS data.
	c
	mse=0
	mxl=0
	icmp=0
	magtyp=' '
	do 50 line=1,maxlin
	c --- get the number of meshes, l_max, and number of components
	read(12,'(a)',end=10)a
	if(a(5:28) .eq. 'meshr mse ng mxl') then
	read (12,*)meshr,mse,ng,mxl
	endif
	if(a(4:8) .eq. 'ntyp='.and. a(13:17) .eq. 'natm='
	& .and. a(22:27) .eq. 'ncmpx=') then
	read(a(28:29),'(i2)')icmp
	endif
	c --- get magtyp
	if(a(68:74) .eq. 'magtyp=') then
	magtyp=a(75:78)
	endif
	c --- the last data before the DOS data begin are adopted
	if(a(2:19) .eq. 'DOS of component 1') go to 100
	50 continue
	write(*,'(a)')'Seems not to contain DOS data.'
	stop
	c 100 write(*,'(a,3i3)')' mse,mxl,icmp=',mse,mxl,icmp
	100 continue
	c --- check the sizes ---
	if(mse .gt. msex) then
	write(,'(a,i3)')' **error...mse>',msex
	stop
	endif
	if(mxl .gt. mmxl) then
	write(,'(a,i2)')' **error...mxl>',mmxl
	stop
	endif
	if(icmp .gt. icmpx) then
	write(,'(a,i2)')' **error...icmp>',icmpx
	stop
	endif
	if(magtyp .ne. 'mag' .and. magtyp .ne. 'nmag') then
	write(,'(a,a4)')' **error...illegal magtyp=',magtyp
	stop
	endif
	if(magtyp .eq. 'mag') then
	ismx=2
	else
	ismx=1
	endif
	c write(,)'ismx=',ismx
	c --- Now the format statement can be constructed
	write(fmtx,'(a,i1,a)')'(1x,f7.4,3x,',mxl,'f10.4)'
	write(fmt,'(a,i2,a)')'((f8.4,f7.2,f8.4,',mxl*icmp,'f7.2))'
	c write(,)'format=',fmtx
	c --- read in DOS data
	do 20 is=1,ismx
	do 30 line=1,maxlin
	if(a(2:19) .eq. 'DOS of component 1') then
	c write(,)'DOS 1 met'
	read(12,fmtx)((dat(l,k,is),l=3,mxl+3),k=1,mse)
	c write(*,fmtx)((dat(l,k),l=3,mxl+3,is),k=1,mse)
	do 40 i=2,icmp
	c write(,)'DOS',i,' met'
	i0=mxl*(i-1)+4
	read(12,'(//)')
	40 read(12,fmtx)
	& (dmy,(dat(l,k,is),l=i0,i0+mxl-1),k=1,mse)
	endif
	if(a(2:10) .eq. 'total DOS') then
	c write(,)'total DOS met'
	read(12,'(1x,f12.7,3x,f13.5)')(dat(1,k,is),dat(2,k,is),k=1,mse-1)
	dat(1,mse,is)=dat(1,mse-1,is)
	dat(2,mse,is)=dat(2,mse-1,is)
	go to 20
	endif
	30 read(12,'(a)',end=10)a
	20 read(12,'(a)',end=10)a
	if(ismx .eq. 1) then
	do 60 i=1,mxl+32
	do 60 k=1,mse
	60 dat(i,k,2)=dat(i,k,1)
	endif
	c --- fix axis design
	xmi=big
	xmx=-big
	ymi=big
	ymx=-big
	do 70 is=1,2
	do 70 k=1,mse
	xmi=min(xmi,dat(1,k,is))
	xmx=max(xmx,dat(1,k,is))
	ymi=min(ymi,dat(2,k,is))
	70 ymx=max(ymx,dat(2,k,is))
	ymi=0d0
	c call ezaxis(xmi,xmx,nex,ntickx,1d0)
	call ezaxis(ymi,ymx,ney,nticky,0d0)
	c xmi=xmi10d0*nex
	c xmx=xmx10d0*nex
	c xincr=(xmx-xmi)/dble(ntickx-1)
	ymi=ymi10d0*ney
	ymx=ymx10d0*ney
	yincr=(ymx-ymi)/dble(nticky-1)
	c --- construct a gnuplot program
	open(13,file=trim(file)//'.plt',status='unknown')
	write(13,'(a)')'#!/usr/bin/env gnuplot'
	write(13,'(a)')'#--Gnuplot--'
	write(13,'(a)')'##'
	write(13,'(a)')'$u << EOF'
	write(13,fmt)((dat(l,k,1),l=1,mxl*icmp+3),k=1,mse)
	write(13,'(a)')'EOF'
	write(13,'(a)')'$d << EOF'
	write(13,fmt)((dat(l,k,2),l=1,mxl*icmp+3),k=1,mse)
	write(13,'(a)')'EOF'
	write(13,'(2a)')'set terminal postscript eps noenhanced color',
	& ' "Helvetica" 20'
	write(13,'(a)')'set output "'//trim(file)//'.eps"'
	c write(13,'(a)')'set term x11'
	write(13,'(a)')'#set yzeroaxis lt -1 lw 0.1'
	c write(13,'(a)')'#set yzeroaxis lt -1 lw 0.1'
	write(13,'(a)')'set border 15 lw 0.1'
	write(13,'(a)')'#set mxtics 2'
	write(13,'(a)')'#set mytics 2'
	c write(13,'(a,e15.7,a,e15.7)')'set xtics',xmi,',',xincr
	write(13,'(a,e15.7,a,e15.7)')'set ytics',ymi,',',yincr
	write(13,'(a)')'set size 1,1'
	write(13,'(a)')'set origin 0,0'
	write(13,'(a)')'set multiplot'
	write(13,'(a)')'set size 1.0,0.35'
	c write(13,'(a)')'set border 14 lw 0.1'
	c write(13,'(a)')'set lmargin at screen 0.1'
	c write(13,'(a)')'set rmargin at screen 0.9'
	write(13,'(a)')'set origin 0.0,0.5'
	write(13,'(a)')'set lmargin 8'
	write(13,'(a)')'set rmargin 3'
	write(13,'(a)')'set tmargin 0'
	write(13,'(a)')'set bmargin 0'
	write(13,'(3a)')'set title "',trim(file),'" font "Times-Roman,20"'
	write(13,'(a)')'set format x ""'
	c write(13,'(a,e15.7,a,e15.7,a)')'set xrange [',xmi,':',xmx,']'
	write(13,'(a,e15.7,a,e15.7,a)')'set xrange [:]'
	write(13,'(a,e15.7,a,e15.7,a)')'set yrange [',ymi,':',ymx,']'
	write(13,'(a)')'plot "$u" using 1:2 title "" with lines lt 1 lw 3'
	write(13,'(a)')'set border 11 lw 0.1'
	write(13,'(a)')'set origin 0.0,0.15'
	write(13,'(a)')'unset title'
	write(13,'(2a)')
	& 'set xlabel "Energy relative to Fermi energy (Ry)"',
	& ' font "Times-Roman,20"'
	write(13,'(2a)')
	& 'set label "DOS (1/Ry/unit cell)" at screen 0.04, screen 0.5',
	& ' center rotate font "Times-Roman,20"'
	write(13,'(a,e15.7,a,e15.7,a)')'set yrange [',ymx,':',ymi,']'
	write(13,'(a)')'set format x'
	write(13,'(a)')'plot "$d" using 1:2 title "" with lines lt 1 lw 3'
	write(13,'(a)')'unset multiplot'
	write(13,'(a)')'set output # Close the EPS file.'
	write(13,'(a)')'!epstopdf '//trim(file)//'.eps'
	c write(13,'(a)')'pause -1'
	c write(13,'(a)')'# EOF'
	close(13)
	call system('gnuplot '//trim(file)//'.plt')
	stop
	10 write(*,'(a/a)')'Illegal file specified','Usage: gp5d file'
	end program gp5d
	subroutine ezaxis(xmin,xmax,nex,ntick,adjust)
	c-----------------------------------------------------------------------
	c I fix a scale of the axis.
	c In output, xmin is maximum integer for which xmin10*nex
	c does not exeed original input xmin.
	c Similary, xmax is maximum integer for which xmin10*nex
	c is not smaller than original xmax.
	c ntick gives a number of ticks including both end points.
	c Coded by H.Akai, April 28, 1991
	c-----------------------------------------------------------------------
	implicit real*8 (a-h,o-z)
	dimension nw(12),nh(12)
	data nw/ 4, 5, 6, 7, 8,10,12,14,16,20,25,30/
	& ,nh/ 4, 5, 6, 7, 4, 5, 6, 7, 8, 4, 5, 6/
	& ,small/1d-3/
	c---- I fix scale of the axis.
	c---- in some cases containing x=0 axis will give better presntation.
	if(xmax .lt. xmin) then
	write(6,1000)
	1000 format(' ***err in ezaxis...xmin and xmax not consistent')
	stop
	endif
	if(xmax*xmin .gt. 0d0 .and.
	& (xmax-xmin)/max(abs(xmax),abs(xmin)) .gt. adjust) then
	xmin=min(xmin,0d0)
	xmax=max(xmax,0d0)
	endif
	wx=xmax-xmin
	c---- wx indicates x-data spreading.
	nex=int(log10(wx/2.9d0)+100d0+small)-100
	sx=10d0**nex
	xmin0=xmin/sx
	xmax0=xmax/sx
	c---- now xmin and xmax is normalized by 10**nex so as to
	c---- satisfy 2.9<xmax-xmin<29.
	xmin=dble(int(xmin0+small))
	xmax=dble(int(xmax0-small))
	c---- bug easily comes in the following statement. never
	c---- use xmin instead of xmin0 in the logical-if, otherwise
	c---- the procedure will fail when -1<xmin<0.
	c---- xmin is maximum integer for which xmin10*nex does not exeed
	c---- original xmin.
	if(xmin0 .lt. 0d0) xmin=xmin-1d0
	c---- xmax is maximum integer for which xmax10*nex is not smaller
	c---- than original xmax.
	if(xmax0 .gt. 0d0) xmax=xmax+1d0
	wx=xmax-xmin
	c---- I then choose appropriate design for tick construction.
	iwx=int(wx+small)
	do 30 i=1,12
	ii=i
	if(nw(i) .ge. iwx) go to 40
	30 continue
	c---- i adjust positioning such that the center of the data spreading
	c---- is not much distant from that of the axis.
	c---- however, if one of xmin or xmax is located near zero, xmin=0 or
	c---- xmax=0 is the arrangement of the highest priority.
	40 ws=dble(nw(ii))
	wm=ws-wx
	if(xmin .gt. -small .and. xmin .lt. wm+small) then
	wm=xmin
	c---- xmin=0 is realized in this case.
	else if(xmax .lt. small .and. xmax .gt. -wm-small) then
	wm=xmax+wm
	c---- xmax=0 is realized in this case.
	else
	if(xmax-xmax0 .gt. xmin0-xmin) wm=wm+1d0
	wm=dble(int(wm/2d0+small))
	c---- otherwise, ticks are arranged so as to spread more or less
	c---- on centered.
	endif
	xmin=xmin-wm
	xmax=xmin+ws
	ntick=nh(ii)+1
	c---- the first and the last ticks of axis are fixed.
	c---- now I can start plotting.
	return
	end
	cLocal variables:
	c compile-command: "gfortran -g -o gp5d gp5d.f"
	cEnd: