Talk:EcalJ

From Ecal

Jump to: navigation, search

注意:メモにつかってください。上のほうが若くなるように書いていってください。 将来的に残していく内容はべつのところへ移動しましょう。混乱もあります。

Contents

sxcf_falの構造

     subroutine sxcf_fal3_scz

      allocate(expikt(natom))
      do 1001 ip = 1,nq         !;write (*,*) ip,'  out of ',nq,'  k-points ' ! call cputid  (0)
         if(sum(irkip(:,:,ip))==0) cycle ! next ip
  

         do 1100 kx = iqini,iqend !kx=1 corresponds to q=0 is omitted.
            if(sum(irkip(kx,:,ip))==0) cycle ! next kx
            if(newaniso()) then
 
               do
                  if(allocated(vcoud)) deallocate(vcoud)
                  allocate( zcousq(ngb0,ngb0),vcoud(ngb0) )
                  if(sum(abs(qvv-qxx))<1d-6) goto 1133
               enddo
 1133          continue
               if(allocated(ppovlz)) deallocate(ppovlz)
               allocate(ppovl(ngc,ngc),ppovlz(ngb,ngb))
               deallocate(zcousq,ppovl)
            endif

            do irot = 1,ngrp
               if( kx <= nqibz) then
                  if(kr==0) cycle ! next irot
               else
                  if( wgt0(kx-nqibz,irot)==0d0 ) cycle ! next irot
               endif
            enddo  
            do 1000 irot = 1,ngrp
               if( kx <= nqibz) then
                  if(kr==0) then
                     cycle      ! next irot
                  end if
               else
                  if( wgt0(kx-nqibz,irot)==0d0 ) then
                     cycle      ! next irot
                  end if
               endif

               allocate( rmelt(ngb, nctot+nbmax, ntqxx),        cmelt(ngb, nctot+nbmax, ntqxx))

               if(exchange) then
                  allocate( zmel (ngb, nctot+nbmax, ntqxx)

                  if(.not.newaniso()) deallocate(vcoul)
                  cycle         ! next irot
               endif            !exchange


               allocate(freq_r(nw_i:nw)) 
               allocate( zmel (ngb, nstate, ntqxx) )
!               deallocate(rmelt,cmelt) ! strange position
               allocate( zw (nblochpmx,nblochpmx), zwz0(nstate,ntqxx,ntqxx) )
               allocate(zmel1(ngb),zmel1_(ntqxx,ngb,nstate))
               allocate(zwz(niw,nstate,ntqxx,ntqxx), zwzi(nstate,ntqxx,ntqxx)) !sf 22may02

 
               if(ua_auto) then
                  allocate(uaa(nstate,ntqxx))
               endif
               if(iSigma_en==5) goto 2005
               allocate(zwzs(npm*nx))

               deallocate(zwzs)
               if(iSigma_en==4) then
                  deallocate(zwz,zwz0,zwzi)
                  deallocate(zmel,zmel1,zmel1_,zw,freq_r) !zw4,
                  cycle ! next irot
               endif
 2005          continue
               if (.not.allocated(zwzs)) allocate(zwzs(npm*nx))
! freq_r, zmel, zw, zwz0, zmel1, zwz, zwzi, zwzs , zmel1_
               if(ua_auto) deallocate(uaa)
               deallocate(zwz,zwz0,zwzi)
! freq_r, zmel, zw, zmel1, zwzs , zmel1_
               if (iSigma_en==0) then
                  deallocate(zmel,zmel1,zmel1_,zw,freq_r,zwzs) !zw4,
!    all  deleted 
                  cycle ! next irot        !no pole contribution for SE(e_f)  sf 23may02
               endif
               deallocate(zwzs)

               allocate( zw3(ngb,ngb,nwxi:nwx))
               deallocate(zw)   !,zw4
! freq_r, zmel,zmel1,  zw3 , zmel1_
               do 2001 itp = 1,ntqxx ! loop over states (q-k,n)
                  do 2011 it = itini,itend ! nt0p corresponds to efp

                        deallocate(zwz4)
                     endif      ! inner iSigma_en=1 or 3
 2011             continue
 2001          continue         !itp
               deallocate(zw3, zmel)
! freq_r,,zmel1, zmel1_
 1050          continue
               deallocate(zmel1,zmel1_)
               deallocate(freq_r)
 1000       continue ! end do irot
 1100    continue               ! end of kx-loop
         if (allocated(vcoul))deallocate(vcoul)
 1001 continue ! end do ip
      if (allocated(expikt))deallocate(expikt)

時間を図るツール

!TIME0
...
!TIME1 "comment"

とdirectiveを書く。nestしていても構わない。 ソース最後(時間を表示させたい部分)に

!TIMESHOW

と書く。

gawk -f $(addtime) -vSTART=1 $(main)hsfp0.sc.m.F 

としてdirectiveをfortran sourceに直すと!TIMESHOWの部分で時間を表示する。

addtime=script/addtime.awk

memory使用量

gawk -f $(septhen) | gawk -f $(alloclist)

と変換するとdirectiveなどなくて自動的にallocate,deallocateの部分にコードを加える。 6番のファイルへallocate,deallocate毎にそれぞれの変数のメモリー使用量を書き出す。

septhen=script/then_separate.awk
alloclist=script/add_alloclist.awk

hx0の並列化

 \langle I | i j \rangle \langle i j | J \rangle \delta( \omega - \omega_{ij} ) where I,J: product basis, j,j: state index, ωij = EiEj, を計算している。microtetrahedronを使っているので実装が面倒になっている。

現在外側のloopでだけ並列化がされている。 内側loopのrotationの回数がk点で異なり各k点のloadが平均化されていないのでk点により時間がかなり異なる。 (hsfpの方は外側のk点とrotationを合わせてloadはほぼ同じになるように並列化されている。)


  1. . hsfpのように一重並列のままrotationも合わせてindexをつくり並列化。対称性により異なるが、cugaseでは最大k点数の数倍にしかならないので並列度は稼げない。しかし、内部の並列性は高い。これは最初にやるべき。
  2. . nbnbが長いloopなので、そこを使い二重並列化。こちらは並列度を稼げるが、同じメモリーを2つ目の並列では持つことになる。symmetrizeはx0kfとは別のところへ持って行き、各k点の計算が終わったらsymmetrizeを行うように変更したほうが流れがすっきりする。
  3. iw-loopで並列化する。rcxqをiwでmemory distributeする。nbnb loopの中にiw loopがある。

iw-loopで並列化するとして

以下はdata distributionせずに

  1. "after matmul zmelt"の部分"は別途並列化する。
  2. melpln2も別途並列する。

zzmelはpsicb_v3 psi2b_v3で作られてzzmelを作りmelpln2のzmeltを加えて、処理が終わるとzzmelは消される。psicb_v3 psi2b_v3の内部もzzmelのdistributionで並列化することが可能。一番最後の軸で切れるか?

  1. symmetrizerの部分もiwで並列化可能。

nbnb

\sum_i^{occ} \sum_j^{unocc} 相当。 -> ip(occ), itp(unocc)になる。

zzmel(nbloch, nctot+nt0, ncc+ntp0) なのでnctot+nt0, ncc+ntp0で分割したいところだが、nbnbとit,itpの対応が意味不明。

rcxq(igb1,igb2,iw,jpm)はigb1,igb2がproduct basisのindex


         do it=1,nctot+nt0
         do itp=1,ncc+ntp0  !! See sxcf_fal2.F L933 around. zmelt1
           zmelt(:,it,itp) =  matmul( zmelt(:,it,itp),ppovlz(:,:) )
         enddo
         enddo

はnbnb相当。

symmetrizeはrcxq(:,:,iw,jpm)のproduct basisで無いと並列化できないが、1,2番目の軸だとtransposeを含むので並列化しにくいやりにくい。


nbnb->it,itpへのloopの作り替えが必要。


program h0fp0の並列化案

subroutine x0kfの並列化

  • zmeltを作るところはnbnbで並列化する。これはzmelt(igb1,it,itp)のit,itpの並列化に相当する。データ並列化も行う。

zmelt(:,it,itp) = matmul( zmelt(:,it,itp),ppovlz(:,:) )があるから。

it,itpはn1b,n2bでindexを返すのでloadが同じになるようにindex分割を行うzmeltのallocationもそれに合わせる。

  • 対称化する前まではrcxqを全てが持つ。
  • 対称化はiwで並列化するしかない。iwでデータ並列化も行う。
  • igwtはtetwt5.Fでつくっている。(<-どこで効いたいのだったか。)


(1)

zzmel: it,itpの並列化

zmelt: zmelt(:,it,itp) = matmul( zmelt(:,it,itp),ppovlz(:,:) )まではit,itpで並列化

rcxq: まずはiwで並列化

symmetrize: do iw=1,nwtで並列化

main h0fp0の並列化

  • iwで並列化するしかない。

hsfp0, correlation mode

もっとも時間がかかるのは

              do ix=1,nx
                 call matzwz2(2, zwix(1:ngb,1:ngb,ix), zmel, ntqxx, nstate,ngb,
    o                 zwzix(1:nstate,1:ntqxx,1:ntqxx)) ! zwz = zmel*(W(0)-v)*zmel

ここでmatzwz2の中身は

       do itpp= 1,ntq
         do itp = 1,ntq
           do  it = 1,nstate
             zwz(it,itp,itpp)=
    &    zdotc(ngb,zmel(1,it,itp),1,CC(1,it,itpp),1 )
           enddo
         enddo
       enddo

という計算。

次に時間がかかるのは

                 do itp=1,ntqxx
                    do it=1,nstate
                       zmel(:,it,itp) =  matmul(zmel(:,it,itp),dconjg(ppovlz(:,:)))
                    enddo
                 enddo

ntqxx is number of bands for \langle i|\Sigma|j\rangle .

この上レベルのloopは

           do 1000 irot = 1,ngrp
        do 1100 kx = iqini,iqend !kx=1 corresponds to q=0 is omitted.
     do 1001 ip = 1,nq         !;write (*,*) ip,'  out of ',nq,'  k-points ' ! call cputid  (0)

となる。

次に時間がかかるのは !  ! -- IPW part.

              allocate( rmelt(ngb, nctot+nbmax, ntqxx), ! nstate= nctot+nband
    &              cmelt(ngb, nctot+nbmax, ntqxx))
              if(debug) print *, ' sxcf_fal1: goto drvmelp2 xxx111'
              allocate(drealzzzmel(nbloch,nt,ntqxx),dimagzzzmel(nbloch,nt,ntqxx))
              drealzzzmel=dreal(zzzmel)
              dimagzzzmel=dimag(zzzmel)
              call drvmelp2( q,             ntqxx, ! q in FBZ
    i              q-qbz_kr,  nbmax, ! q-rk
    i              qibz_k,     ! k in IBZ for e-product basis
    i              isp,ginv,
    i              ngc,ngcmx, ngpmx,nband,itq,
    i              symope, shtv, qbas, qbasinv,qibz,qbz,nqbz,nqibz,
    i              drealzzzmel, dimagzzzmel, nbloch, nt,nctot,
    o              rmelt,cmelt)

そして次が

              do 2001 itp = 1,ntqxx ! loop over states (q-k,n)
                 do 2011 it = itini,itend ! nt0p corresponds to efp

のloop 中に

                       do ix0=1,3
                          ix=ixs+ix0-2
                          do igb2=1,ngb !**** most time consuming part for iSigma_en=3 ******
                             zz2=sum(zmel1(1:ngb)*zw3(1:ngb,igb2, iir*ix)  )
                             do itpp=1,ntqxx
                                zwz4(ix0,itpp)=zwz4(ix0,itpp)+ zz2*zmel1_(itpp,igb2,it)
                             enddo !itpp
                          enddo !igb2
                       enddo   !ix

という計算がある。

  • ntq: band数max

zwzix(1:nstate,1:ntqxx,1:ntqxx))を最後の軸で切る?

ntqxxの軸で切るのが良さそう。

hsfp0, exchange mode

時間がかかるのは !$OMP do

                    do itp= 1,ntqxx
                       do it = 1,nctot+nbmax
                          do ivc=1,ngb
                             zmeltt(it,itp,ivc) =  sum( zmel(:,it,itp)* ppovlz(:,ivc) )
                          enddo
                       enddo
                    enddo

!$OMP end do !$OMP do

                    do 992 itpp= 1,ntqxx
                       do 993 itp = 1,ntqxx
                          if(diagonly.and.(itpp/=itp)) cycle
                          do 994 it  = 1,nctot+nbmax
                             w3p(it,itp,itpp) = 0d0
                             do ivc=1,ngb
                                if(ivc==1.and.kx==iqini) then
                                   vc= wklm(1)* fpi*sqrt(fpi) /wk(kx)
                                else
                                   vc= vcoud(ivc)
                                endif
                                w3p(it,itp,itpp) = w3p(it,itp,itpp)
    &                                + vc * zmeltt(it,itp,ivc)*dconjg(zmeltt(it,itpp,ivc))
                             enddo
994                       continue
993                    continue
992                 continue

ntxqqの軸で切るのが良さそう。

sugase2

lsxC,nosmp

   1 end_program___________________   785.22   785.21   785.21   785.20   785.22    1    1
   3 before_2000_loop______________     0.14     0.14     0.63     0.14     0.63    1    1
   5 before_sxcf_fal3_scz__________     0.00     0.00     0.00     0.00        0    1    1
   7 after_sxcf_fal3_scz___________   718.39   771.18   783.86   494.57   783.86    1    1
   9 reduction_nbandmx_end_________    65.95    13.16     0.00     0.00   289.29    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.00     0.00     0.02     0.00     0.02    1    1
  15 end_3020_loop_________________     0.72     0.72     0.70     0.70     0.72    1    1
  17 before_2000_continue__________     0.01     0.00     0.00     0.00     0.01    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________   718.39   771.18   783.85   494.57   783.85    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.01     0.03     0.09     0.01     0.10    1    2
 106 before_1000___________________     2.66     2.73     1.87     1.62     2.73    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.22     0.13     0.16     0.09     0.22    2    3
 112 bfore_psi2b_v2________________     1.93     2.00     2.03     1.37     2.03    2    3
 114 after_drvmelp2________________     1.37     1.99     3.23     1.37     3.39    2    3
 116 end_exchange_section__________     0.03     0.03     0.03     0.02     0.03    2    3
 118 end_of_matm_vcoult____________     0.45     0.47     0.47     0.32     0.47    2    3
 120 end_of_write_ifsexsp__________   711.67   763.75   775.91   488.31   775.91    2    3
 122 endo_of_zsec_wtt_sum__________     0.05     0.05     0.06     0.04     0.06    2    3
 200 after_readgeig_q_q_rk_________     0.16     0.10     0.21     0.08     0.21    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________     0.30     0.94     2.07     0.30     2.66    2    3
 206 end_drvmelpln2________________     0.51     0.53     0.53     0.36     0.53    2    3

lsxC,4smp

   1 end_program___________________   215.64   215.62   215.62   215.62   215.64    1    1
   3 before_2000_loop______________     0.12     0.47     0.47     0.12     0.47    1    1
   5 before_sxcf_fal3_scz__________     0.00     0.00     0.00     0.00        0    1    1
   7 after_sxcf_fal3_scz___________   197.98   213.55   214.42   136.48   214.42    1    1
   9 reduction_nbandmx_end_________    16.78     0.87     0.00     0.00    77.94    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.00     0.04     0.04     0.00     0.04    1    1
  15 end_3020_loop_________________     0.73     0.69     0.69     0.69     0.73    1    1
  17 before_2000_continue__________     0.01     0.00     0.00     0.00     0.01    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________   197.98   213.54   214.42   136.48   214.42    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.01     0.12     0.09     0.01     0.13    1    2
 106 before_1000___________________     2.71     3.31     1.90     1.70     3.31    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.19     0.14     0.15     0.09     0.19    2    3
 112 bfore_psi2b_v2________________     1.91     1.99     2.02     1.37     2.02    2    3
 114 after_drvmelp2________________     2.02     3.12     2.95     2.02     3.12    2    3
 116 end_exchange_section__________     0.03     0.03     0.03     0.02     0.03    2    3
 118 end_of_matm_vcoult____________     0.46     0.47     0.47     0.32     0.47    2    3
 120 end_of_write_ifsexsp__________   190.60   204.31   206.75   130.05   206.75    2    3
 122 endo_of_zsec_wtt_sum__________     0.05     0.06     0.06     0.04     0.06    2    3
 200 after_readgeig_q_q_rk_________     0.22     0.18     0.11     0.11     0.22    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________     0.86     1.99     1.86     0.86     2.14    2    3
 206 end_drvmelpln2________________     0.54     0.55     0.55     0.36     0.55    2    3

end of write ifsexspとは

                 if(newaniso()) then
                   write(*,'(a,5I10)')'kino: ntqxx,nctot+nbmax,ngb=',ntqxx,nctot+nbmax,ngb
                     allocate(zmeltt(nctot+nbmax,ntqxx,ngb))
!$OMP parallel  private(vc)
!$OMP do
                     do itp= 1,ntqxx
                        do it = 1,nctot+nbmax
                           do ivc=1,ngb
                              zmeltt(it,itp,ivc) =  sum( zmel(:,it,itp)* ppovlz(:,ivc) )
                           enddo
                        enddo
                     enddo
!$OMP end do
!$OMP do
                     do 992 itpp= 1,ntqxx
                        do 993 itp = 1,ntqxx
                           if(diagonly.and.(itpp/=itp)) cycle
                           do 994 it  = 1,nctot+nbmax
                              w3p(it,itp,itpp) = 0d0
                              do ivc=1,ngb
                                 if(ivc==1.and.kx==iqini) then
                                    vc= wklm(1)* fpi*sqrt(fpi) /wk(kx)
c     print *,'wklm(1) vc=',wklm(1),vc
                                 else
                                    vc= vcoud(ivc)
                                 endif
c     zmelt1 =  sum( zmel(:,it,itp)  *ppovlz(:,ivc) )
c     zmelt2 =  sum( zmel(:,it,itpp) *ppovlz(:,ivc) )
                                 w3p(it,itp,itpp) = w3p(it,itp,itpp)
    &                                + vc * zmeltt(it,itp,ivc)*dconjg(zmeltt(it,itpp,ivc))
                              enddo
 994                       continue
 993                    continue
 992                 continue
!$OMP end do
!$OMP end parallel
                     deallocate(zmeltt)
                  else
!$OMP parallel do
                     do itpp= 1,ntqxx
                        do itp = 1,ntqxx
                           if(diagonly.and.(itpp/=itp)) cycle
                           do it  = 1,nctot+nbmax
                              w3p(it,itp,itpp) =dcmplx(
     &                             sum ( dreal(z1p(:,it,itp))*rmelt(:,it,itpp)
     &                             +   dimag(z1p(:,it,itp))*cmelt(:,it,itpp) ) ,
     &                             sum ( dimag(z1p(:,it,itp))*rmelt(:,it,itpp)
     &                             -   dreal(z1p(:,it,itp))*cmelt(:,it,itpp) ) )
                           enddo
                        enddo
                     enddo
!$OMP end parallel do

lsx,nosmp

kino@asahi01:/data1/kino/ECALJ> gawk -f time.awk cugase2_222.nompi/lsx
   1 end_program___________________   196.70   196.69   196.69   196.68   196.70    1    1
   3 before_2000_loop______________     0.30     0.30     0.29     0.29     0.30    1    1
   5 before_sxcf_fal3_scz__________     0.00     0.00     0.00     0.00        0    1    1
   7 after_sxcf_fal3_scz___________   183.60   193.75   195.82   132.34   195.82    1    1
   9 reduction_nbandmx_end_________    12.21     2.06     0.00     0.00    63.48    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.02     0.02     0.00     0.00     0.02    1    1
  15 end_3020_loop_________________     0.55     0.56     0.57     0.55     0.57    1    1
  17 before_2000_continue__________     0.01     0.00     0.00     0.00     0.01    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________   183.59   193.75   195.82   132.34   195.82    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.01     0.01     0.01     0.01     0.01    1    2
 106 before_1000___________________     1.84     1.52     1.22     1.00     1.84    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.08     0.08     0.08     0.06     0.08    2    3
 112 bfore_psi2b_v2________________     1.41     1.46     1.49     1.01     1.49    2    3
 114 after_drvmelp2________________    44.91    47.06    47.49    32.47    47.49    2    3
 116 end_exchange_section__________     0.01     0.01     0.01     0.00     0.01    2    3
 118 end_of_matm_vcoult____________     0.12     0.13     0.13     0.09     0.13    2    3
 120 end_of_write_ifsexsp__________   135.17   143.43   145.35    96.49   145.35    2    3
 122 endo_of_zsec_wtt_sum__________     0.03     0.04     0.04     0.02     0.04    2    3
 200 after_readgeig_q_q_rk_________     0.06     0.09     0.06     0.04     0.09    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________    44.58    46.71    47.16    32.24    47.16    2    3
 206 end_drvmelpln2________________     0.17     0.17     0.17     0.12     0.17    2    3

lsx,4smp

kino@asahi01:/data1/kino/ECALJ> gawk -f time.awk cugase2_222.4mpi/lsx
   1 end_program___________________    84.15    84.14    84.14    84.14    84.15    1    1
   3 before_2000_loop______________     0.39     0.38     0.37     0.37     0.39    1    1
   5 before_sxcf_fal3_scz__________     0.00     0.00     0.00     0.00        0    1    1
   7 after_sxcf_fal3_scz___________    80.03    83.17    82.51    56.50    83.17    1    1
   9 reduction_nbandmx_end_________     3.12     0.00     0.67     0.00    26.68    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.05     0.00     0.00     0.00     0.05    1    1
  15 end_3020_loop_________________     0.54     0.59     0.59     0.54     0.59    1    1
  17 before_2000_continue__________     0.01     0.00     0.00     0.00     0.01    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________    80.03    83.17    82.50    56.50    83.17    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.01     0.01     0.01     0.01     0.01    1    2
 106 before_1000___________________     1.87     1.96     1.06     1.06     1.96    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.08     0.09     0.08     0.06     0.09    2    3
 112 bfore_psi2b_v2________________     1.41     1.47     1.49     1.01     1.49    2    3
 114 after_drvmelp2________________    41.82    42.59    42.66    29.39    42.66    2    3
 116 end_exchange_section__________     0.01     0.01     0.01     0.00     0.01    2    3
 118 end_of_matm_vcoult____________     0.12     0.13     0.13     0.09     0.13    2    3
 120 end_of_write_ifsexsp__________    34.67    36.88    37.04    24.67    37.04    2    3
 122 endo_of_zsec_wtt_sum__________     0.03     0.04     0.04     0.03     0.04    2    3
 200 after_readgeig_q_q_rk_________     0.06     0.07     0.06     0.04     0.07    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________    41.50    42.25    42.33    29.17    42.33    2    3
 206 end_drvmelpln2________________     0.17     0.17     0.17     0.12     0.17    2    3

lx0,nosmp

kino@asahi01:/data1/kino/ECALJ> gawk -f time.awk cugase2_222.nompi/lx0
   1 end_program___________________  2505.04  2418.13  2198.23  1905.78  2505.04    1    1
   3 after_readbzdata______________     0.03     0.12     0.04     0.03     0.12    1    1
   5 before_readngmx_______________     0.01     0.01     0.03     0.01     0.05    1    1
   7 before_initreadeigen__________     0.01     0.01     0.01     0.01     0.01    1    1
   9 before_histogram_division_____     0.00     0.00     0.00     0.00        0    1    1
  11 before_realomega_write3111____     0.00     0.00     0.00     0.00        0    1    1
  13 before_open_WVd_______________     0.01     0.00     0.00     0.00     0.01    1    1
  15 after_write_ifwd______________     0.00     0.01     0.06     0.00     0.06    1    1
  17 endo_of_ppbafp_v2_____________     0.03     0.03     0.03     0.03     0.03    1    1
  19 before_1001_loop______________     0.00     0.00     0.00     0.00        0    1    1
  21 before_1003_loop______________     2.04     1.92     1.92     1.92     2.10    1    1
  23 before_x0kf_v2hz______________     0.62     0.62     0.62     0.62     0.64    1    1
  25 x0vk_v4h______________________  1835.67  2094.57  1875.25  1573.97  2114.87    1    1
  27 before_real_omega_____________    42.51    40.84    40.47    40.47    42.55    1    1
  29 real_omega____________________   305.82   269.38   269.20   269.20   305.82    1    1
  31 end_1015_loop_________________   304.94   268.62   268.55   268.55   304.94    1    1
  33 after_write_trace_msg_________     0.00     0.00     0.00     0.00        0    1    1
  35 after_matcinv_________________   276.37   255.82   255.92            276.37       258
  37 after_zw=zw0__________________    22.93     7.31     7.26             22.93       258
  39 end_of_tr_chkwrite____________     5.62     5.49     5.37              5.62       258
  45 imag_omega____________________    11.85    10.47    10.44    10.44    11.85    1    1
  47 end_eeeeeee___________________   247.79     0.00     0.00     0.00   247.79    1    1
  49 before_closing_hbe.d__________    58.48     0.00     0.00     0.00    58.48    1    1
 100 before_if_smbasis_____________     0.00     0.00     0.00     0.00        0    1    1
 102 after_if_smbasis______________     0.00     0.00     0.00     0.00        0    1    1
 104 before_psicb_v3_______________     0.03     0.05     0.08     0.03     0.23    4    6
 106 after_psicb_v3________________     4.23     6.34     5.30     4.23     6.35    4    6
 108 before_melpln2________________     0.32     0.47     0.61     0.32     0.84    4    6
 110 after_melpln2_________________   140.99   203.91   171.89   137.97   210.98    4    6
 112 after_matmul_zmelt____________   580.53   818.90   680.15   574.80   867.59    4    6
 114 before_jpm_ibib_loop__________     0.00     0.00     0.00     0.00        0    4    6
 116 after_rcxq____________________   703.21   991.15   845.48   669.92   991.15    4    6
 118 before_eibzmoden=2____________     1.63     1.60     1.60     1.36     1.63    1    1
 120 after_iqindx2_________________     0.00     0.00     0.00     0.00        0    1    1
 122 before_1011___________________     1.06     0.99     0.98     0.98     1.06    1    1
 124 call_cputid2_qqqqq11__________     0.00     0.00     0.00     0.00        0    1    1
 126 end_matmmsparse_______________    55.92    13.07    25.81     6.98    55.92    2   16
 128 before_iagain=1_______________     0.26     0.06     0.12     0.03     0.26    2   16
 130 end_ig_itimer_icc,_nrotm______     0.00     0.00     0.00     0.00        0    2   16
 132 after_sym_rcxq________________   343.24    54.01   139.22    30.34   343.24    1    1

lx0,4smp

kino@asahi01:/data1/kino/ECALJ> gawk -f time.awk cugase2_222.4mpi/lx0
   1 end_program___________________  1422.32  1325.94  1263.37  1063.06  1422.32    1    1
   3 after_readbzdata______________     0.02     0.03     0.02     0.02     0.06    1    1
   5 before_readngmx_______________     0.03     0.02     0.03     0.02     0.03    1    1
   7 before_initreadeigen__________     0.02     0.02     0.02     0.01     0.02    1    1
   9 before_histogram_division_____     0.00     0.01     0.01     0.00     0.01    1    1
  11 before_realomega_write3111____     0.00     0.00     0.00     0.00        0    1    1
  13 before_open_WVd_______________     0.02     0.00     0.00     0.00     0.02    1    1
  15 after_write_ifwd______________     0.06     0.07     0.07     0.03     0.07    1    1
  17 endo_of_ppbafp_v2_____________     0.03     0.03     0.03     0.03     0.03    1    1
  19 before_1001_loop______________     0.00     0.00     0.00     0.00        0    1    1
  21 before_1003_loop______________     2.05     1.91     1.91     1.91     2.08    1    1
  23 before_x0kf_v2hz______________     0.62     0.62     0.62     0.62     0.64    1    1
  25 x0vk_v4h______________________   889.34  1002.89   940.22   731.37  1002.89    1    1
  27 before_real_omega_____________    42.49    40.45    40.49    40.45    42.53    1    1
  29 real_omega____________________   312.76   269.32   269.36   269.32   312.76    1    1
  31 end_1015_loop_________________   311.87   268.66   268.71   268.66   311.87    1    1
  33 after_write_trace_msg_________     0.00     0.00     0.00     0.00        0    1    1
  35 after_matcinv_________________   275.68   255.89   255.91            275.68       258
  37 after_zw=zw0__________________    30.76     7.27     7.27             30.76       258
  39 end_of_tr_chkwrite____________     5.42     5.50     5.53              5.53       258
  45 imag_omega____________________    12.26    10.44    10.44    10.44    12.26    1    1
  47 end_eeeeeee___________________    64.93     0.00     0.00     0.00    64.93    1    1
  49 before_closing_hbe.d__________    97.55     0.00     0.00     0.00    97.55    1    1
 100 before_if_smbasis_____________     0.00     0.00     0.00     0.00        0    1    1
 102 after_if_smbasis______________     0.00     0.00     0.00     0.00        0    1    1
 104 before_psicb_v3_______________     0.04     0.05     0.08     0.04     0.16    4    6
 106 after_psicb_v3________________     4.24     6.33     5.30     4.24     6.35    4    6
 108 before_melpln2________________     0.36     0.47     0.48     0.36     0.81    4    6
 110 after_melpln2_________________   132.42   188.27   161.14   129.60   196.30    4    6
 112 after_matmul_zmelt____________   173.15   225.03   193.19   155.41   242.69    4    6
 114 before_jpm_ibib_loop__________     0.00     0.00     0.00     0.00        0    4    6
 116 after_rcxq____________________   425.12   547.54   508.78   365.13   547.54    4    6
 118 before_eibzmoden=2____________     2.38     1.60     2.13     1.60     2.38    1    1
 120 after_iqindx2_________________     0.00     0.00     0.00     0.00        0    1    1
 122 before_1011___________________     1.06     0.98     0.98     0.98     1.06    1    1
 124 call_cputid2_qqqqq11__________     0.00     0.00     0.00     0.00        0    1    1
 126 end_matmmsparse_______________    55.32    13.25    25.93     6.95    55.32    2   16
 128 before_iagain=1_______________     0.26     0.06     0.12     0.03     0.26    2   16
 130 end_ig_itimer_icc,_nrotm______     0.00     0.00     0.00     0.00        0    2   16
 132 after_sym_rcxq________________    90.74    15.28    38.06     9.10    90.74    1    1

after matcinvとは

                     if(iq==1) then
                        ix=1
                        zw0(:,1)=0d0
                        zw0(1,:)=0d0
                     else
                        ix=0
                     endif

!     !  Eqs.(37),(38) in PRB81 125102
                     do igb1=ix+1,ngb
                        do igb2=ix+1,ngb
                           epstilde(igb1,igb2)= -vcousq(igb1)*zxq(igb1,igb2,iw)*vcousq(igb2)
                           if(igb1==igb2) epstilde(igb1,igb2)=1+epstilde(igb1,igb2)
                        enddo
                     enddo
                     epstinv(ix+1:ngb,ix+1:ngb)=epstilde(ix+1:ngb,ix+1:ngb)
                     call matcinv(ngb-ix,epstinv(ix+1:ngb,ix+1:ngb))

before closing hbe.d

       if(newaniso2 .and. mpi__rank == 0 ) then ! MIZUHO-IR only on root node

            if(ixc==1011) then !this is only for test.
            ifw0w0i = iopen('W0W0I',0,-1,0)
            read(ifw0w0i) nw_i,nw,niw,nq0i
            print *,'w0w0i: n=',nw_i,nw,niw,nq0i
            read(ifw0w0i) llw(nw_i:nw,1:nq0i)
            read(ifw0w0i) llwI(1:niw,1:nq0i)
c            read(ifw0w0i) w0(nw_i:nw)
c            read(ifw0w0i) w0i(1:niw)
            ifw0w0i = iclose('W0W0I')
            endif

            print *
...
               do ircw=1,2
                  if    (ircw==1) then
                     nini=nw_i
                     nend=nw
                     ifrcwx = iopen('WVR.'//charnum5(iq),0,-1,mrecl)
                  elseif(ircw==2) then;  nini=1;      nend=niw;
                     ifrcwx = iopen('WVI.'//charnum5(iq),0,-1,mrecl)
                  endif
                  do iw=nini,nend
c     if(iq<=nqibz) read(ifrcwx, rec=((iq-iqxini)*(nend-nini+1)+ iw-nini+1 ) ) zw !(1:ngb,1:ngb)
                     read(ifrcwx, rec= iw-nini+1 ) zw !(1:ngb,1:ngb)
                     if( iq==1 ) then
                        if(ircw==1) zw(1,1) = w0(iw)
                        if(ircw==2) zw(1,1) = w0i(iw)
                     endif
c     write(ifrcwx,rec=((iq-iqxini)*(nend-nini+1)+ iw-nini+1 ) ) zw !(1:ngb,1:ngb)
                     write(ifrcwx,rec=iw-nini+1) zw !(1:ngb,1:ngb)
                  enddo
                  if    (ircw==1) then
                     ifrcwx = iclose('WVR.'//charnum5(iq))
                  elseif(ircw==2) then
                     ifrcwx = iclose('WVI.'//charnum5(iq))
                  endif
               enddo
            end do
         end if

before melpln2


C --- IPW set
        call readqg('QGpsi',q+rk(:,k)-qq,ginv, qu1x, ngp1, ngvecpB1)
        call readqg('QGpsi',  rk(:,k)-qq,ginv, qu2x, ngp2, ngvecpB2)

c      ngp1 = ngpn(kp)  ! q+k   ntp0 in FBZ
c      ngp2 = ngpn(k)   ! k     np0  in FBZ
!     ngc              ! q          in IBZ
        ngb  = nbloch + ngc ! This is not ngbb for smbasis()=T. oct2005

c                        q       k        q+k
        allocate( zmelt(ngb,nctot+nt0,ncc+ntp0) )
cooo
c        allocate( zmelt1(ngb,nctot+nt0,ncc+ntp0) )
        allocate( z1p(ngb,ngb) )
C ... read plane wave part of eigenfunction.
        call readgeig(q+rk(:,k)-qq, ngpmx,isp2, qu1, geig1)
        call readgeig(  rk(:,k)-qq, ngpmx,isp1, qu2, geig2)
       if(sum(abs(qu1-qu1x))>tolqu) then
          write(6,"('qu1 :',3d23.10)") qu1
          write(6,"('qu1x:',3d23.10)") qu1x
          stop 'x0kf_v4hz:qu1/=qu1x'
        endif
        if(sum(abs(qu2-qu2x))>tolqu) then
          write(6,"('qu2 :',3d23.10)") qu2
          write(6,"('qu2x:',3d23.10)") qu2x
          stop 'x0kf_v4hz:qu2/=qu2x'
        endif

c     qdiff = q   - qbkp(:) + rk(:,k)
        qdiff = q    - qu1    + qu2
        ! q   - (q+k)   + k  is not zero.
        ! qc  -  q1     + q2
        add   = matmul(qbasinv, qdiff)
        nadd  = idint( add + dsign(.5d0,add))  !  print *,' qdif=',qdiff,qbkp(:),rk(:,k)
        if(sum(abs(add-nadd))>1d-10) stop "sexc: abs(add-nadd))>1d-10"
        zmelt = 0d0

after melpln2

        if(ngc/=0) then !Aug2005
          call melpln2(ngp1, ngvecpB1  ! q1=q+k  ; kp ngp1 1:ntp0 q-point
     &             , ngp2, ngvecpB2  ! q2=k    ; k  ngp2 1:nt0  occupied
     &           , ngc,  nadd,
     &       geig1(1:ngp1, itps:itps+ntp0-1), ntp0, ! q1=q+k  ; kp
     &       geig2(1:ngp2, 1:nt0     ),  nt0, ! q2=k    ; k
     i       shtv,  q, q,  symope,qbas,
C... extensiton to nbloch+ngc
     o       zmelt (nbloch+1:nbloch+ngc, nctot+1:nctot+nt0,ncc+1:ncc+ntp0))
        endif
        zmelt(1:nbloch, 1:nctot+nt0, 1:ncc+ntp0) =
     &  zzmel(1:nbloch, 1:nctot+nt0, 1:ncc+ntp0)
!                         k            q+k
        deallocate(zzmel)
        if(debug) write(*,'("4 z1pp definitions begin",3d13.6,$)') sum(abs(zmelt))
        if(debug) call cputid(0)

after_rcxq

!$OMP parallel private(it,itp,iww1,iww2,zmelt2,imagweight)
        do 25 jpm  = 1, npm !
        do 25 ibib = 1, nbnb(k,jpm) !---  ibib loop
          write(*,'(a,5i8)')'kino: ngb,hw(ibib,k,jpm),ihw(ibib,k,jpm)+nhw(ibib,k,jpm)-1=',
     &      ngb,ihw(ibib,k,jpm),ihw(ibib,k,jpm)+nhw(ibib,k,jpm)-1
          if(n1b(ibib,k,jpm) <= nband) then
            it = nctot + n1b(ibib,k,jpm) !valence
            if(it > nctot + nkmax ) cycle
          else
            it = n1b(ibib,k,jpm) - nband !core
          endif
          if( n2b(ibib,k,jpm) <= nband) then
            itp = ncc + n2b(ibib,k,jpm) - itps + 1 !val
            if(itp > ncc + nkqmax-itps+1 ) cycle
          else
            itp =  n2b(ibib,k,jpm) - itps + 1 - nband !core
          endif

          if(jpm==1) then
            if(n2b(ibib,k,jpm)>nbmx)  then  !nbmx
              if(iww1) then
                print *,' nband_chi0 nbmx=',nbmx
                iww1=.false.
              endif
              cycle
            endif
            if( n1b(ibib,k,jpm) <= nbcut .and. nbcut2<n2b(ibib,k,jpm) ) then
              if(iww2) then
                write(6,"(' nband_chi0 nbcut nbcut2 n2b n1b=',4i6)") nbcut,n2b(ibib,k,jpm),n1b(ibib,k,jpm)
                iww2=.false.
              endif
              cycle
            endif

          else !jpm==2 ------------------------------
            if( n1b(ibib,k,jpm) > nbmx) then  !nbmx
              if(iww1) then
                print *,' nband_chi0 nbmx=',nbmx
                iww1=.false.
              endif
              cycle
            endif
            if( n2b(ibib,k,jpm) <= nbcut .and. nbcut2<n1b(ibib,k,jpm) ) then
              if(iww2) then
                write(6,"(' nband_chi0 nbcut nbcut2 n2b n1b=',4i6)") nbcut,n2b(ibib,k,jpm),n1b(ibib,k,jpm)
                iww2=.false.
              endif
              cycle
            endif
          endif
          if (ihw(ibib,k,jpm)+nhw(ibib,k,jpm)-1 >nwt) stop "x0kf_v4hz: iw>nwt"
!$OMP do private(zmelt2)
            do igb2=1, ngb !....................................
              zmelt2 = zmelt(igb2,it,itp)
            do igb1=1,igb2
              z1p(igb1,igb2) = dconjg(zmelt(igb1,it,itp)) * zmelt2
            enddo
            enddo
!$OMP end do
c$$$        endif

!$OMP do private(imagweight)
          do iw = ihw(ibib,k,jpm),ihw(ibib,k,jpm)+nhw(ibib,k,jpm)-1   !iiww=iw+ihw(ibib,k)-1
            imagweight = whw(jhw(ibib,k,jpm)+iw-ihw(ibib,k,jpm))
            if(eibzmoden==2) imagweight = nwgt(k)*imagweight
            do igb2=1,ngb  !this part dominates cpu time most time consuming...........
            do igb1=1,igb2
              rcxq(igb1,igb2,iw,jpm) =  !here we  sum over ibib (or n, n') and k.
     &        rcxq(igb1,igb2,iw,jpm) + z1p(igb1,igb2)*imagweight !sum over spin in hx0fp0
            enddo !igb1
            enddo !igb2
          enddo ! iw
!$OMP end do

 25     continue
!$OMP end parallel

after_sym_rcxq

!$OMP parallel private(rcxq000,icc,itt,icount,rcxqwww,rcxq00,rcxq0)
        allocate(rcxq0(ngb,ngb),rcxq00(ngb,ngb),rcxq000(ngb,ngb),rcxqwww(ngb,ngb))
!$OMP do
        do iw=1,nwt
        do jpm=1,npm
          rcxq000 = 0d0
          icc=0
          do itimer=1,-1,-2
            if(itimer==1 ) itt=1
            if(itimer==-1) itt=2
            icount=0
            if(itimer==1) then
              rcxqwww = rcxq(:,:,iw,jpm)
            else
              rcxqwww = transpose(rcxq(:,:,iw,jpm))
            endif
            rcxq00 = 0d0
            do ig=1,ngrp
              if(eibzsym(ig,itimer)==1) then
               icount=icount+1
               icc=icc+1
               rcxq0 =0d0

               do irotm1 = 1,nrotm(icc)
                 if(abs(zrr(irotm1,icc))<1d-8) cycle
                 rcxq0(:,i2(irotm1,icc)) =rcxq0(:,i2(irotm1,icc)) + rcxqwww(:,i1(irotm1,icc)) * zrr(irotm1,
icc)
               enddo
               do irotm2 = 1,nrotm(icc)
                 if(abs(zrrc(irotm2,icc))<1d-8) cycle
                 rcxq00(i2(irotm2,icc),:)= rcxq00(i2(irotm2,icc),:) + zrrc(irotm2,icc) * rcxq0(i1(irotm2,ic
c),:)
               enddo
              endif
            enddo

           if(itimer==1) then
               rcxq000=rcxq00
            else
               rcxq000=rcxq000+rcxq00
            endif
          enddo
          rcxq(:,:,iw,jpm) = rcxq000/neibz
        enddo
        enddo
!$OMP end  do
        deallocate(rcxq00,rcxq000,rcxq0,rcxqwww)
!$OMP end parallel

lsc,8mpi x nosmpi

kino@asahi01:/data1/kino/ECALJ> gawk -f time.awk cugase2_222.nompi/lsc.0run
   1 end_program___________________  4609.12  4609.09  4609.10  4609.09  4609.12    1    1
   3 before_2000_loop______________     0.09     0.09     0.09     0.09     0.09    1    1
   5 before_sxcf_fal3_scz__________     0.02     0.00     0.00     0.00     0.02    1    1
   7 after_sxcf_fal3_scz___________  4495.34  4608.18  4544.10  3182.77  4608.18    1    1
   9 reduction_nbandmx_end_________   112.82     0.00    64.08     0.00  1425.41    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.02     0.02     0.00     0.00     0.02    1    1
  15 end_3020_loop_________________     0.80     0.80     0.82     0.80     0.82    1    1
  17 before_2000_continue__________     0.03     0.00     0.00     0.00     0.03    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________  4495.34  4608.18  4544.10  3182.77  4608.18    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.01     0.07     0.01     0.01     0.07    1    2
 106 before_1000___________________     1.98     1.90     0.91     0.90     1.98    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.09     0.15     0.09     0.06     0.15    2    3
 112 bfore_psi2b_v2________________     8.13     8.28     8.35     5.60     8.35    2    3
 114 after_drvmelp2________________   494.17   507.35   500.37   345.08   507.35    2    3
 116 end_exchange_section__________     0.07     0.07     0.07     0.04     0.07    2    3
 124 end_reading_3111______________     0.00     0.01     0.03     0.00     0.03    2    3
 126 end_check_freq_r______________     0.00     0.00     0.00     0.00        0    2    3
 128 matmul_zmel_ppovlz____________   980.84  1006.15   983.98   682.17  1006.15    2    3
 130 before_read_frcw_zw___________     6.61     6.84     6.74     4.62     6.84    2    3
 132 read_ifrwc_zw_before_matzwz2__     0.81     0.10     0.64     0.10     0.81    2    3
 134 after_matzwz2_________________   234.17   241.41   237.43   165.68   241.41    2    3
 136 before_ccc333_k-cycle_________     0.12     0.13     0.13     0.09     0.13    2    3
 138 end_do_xxx33333_k-cycle_______  2349.29  2417.42  2386.47  1664.04  2417.42    2    3
 140 before_1385___________________     0.00     0.00     0.00     0.00        0    2    3
 142 end_1385_loop_________________    25.87    26.85    27.51    18.53    27.51    2    3
 144 end_yyy333_k-cycle____________     0.07     0.07     0.07     0.05     0.07    2    3
 146 end_of_deallocate_zmel_zmel1__     0.15     0.15     0.15     0.11     0.15    2    3
 148 end_of_3001___________________     0.07     0.07     0.07     0.05     0.07    2    3
 150 end_of_4001___________________     0.06     0.06     0.06     0.04     0.06    2    3
 152 before_2001___________________    35.07    24.01    26.41    13.68    35.07    2    3
 154 end_2001_loop_________________   357.35   366.68   364.24   263.54   366.68    2    3
 156 end_of_1100_loop______________     0.41     0.40     0.26     0.26     0.41    2    3
 200 after_readgeig_q_q_rk_________     0.18     0.14     0.06     0.04     0.18    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________   491.18   504.40   497.51   343.13   504.40    2    3
 206 end_drvmelpln2________________     1.80     1.80     1.78     1.23     1.80    2    3

lsc,4smp

kino@asahi01:/data1/kino/ECALJ/cugase2_222.4mpi.run2> gawk -f ../time.awk lsc
   1 end_program___________________  2762.22  2762.21  2762.21  2762.20  2762.22    1    1
   3 before_2000_loop______________     0.83     0.82     0.82     0.82     0.83    1    1
   5 before_sxcf_fal3_scz__________     0.02     0.00     0.00     0.00     0.02    1    1
   7 after_sxcf_fal3_scz___________  2535.85  2760.63  2675.50  1935.31  2760.63    1    1
   9 reduction_nbandmx_end_________   224.76     0.00    85.13     0.00   825.32    1    1
  11 end_legas_exchange____________     0.00     0.00     0.00     0.00        0    1    1
  13 end_of_allocate_zsect_________     0.03     0.04     0.03     0.03     0.04    1    1
  15 end_3020_loop_________________     0.72     0.72     0.72     0.72     0.72    1    1
  17 before_2000_continue__________     0.01     0.00     0.00     0.00     0.01    1    1
  19 end_of_exspectrum_____________     0.00     0.00     0.00     0.00        0    1    1
 100 end_of_sxcf_fal2_scz__________  2535.85  2760.63  2675.50  1935.31  2760.63    1    1
 102 before_1001___________________     0.00     0.00     0.00     0.00        0    1    1
 104 before_1100___________________     0.09     0.12     0.09     0.09     0.12    1    2
 106 before_1000___________________     1.59     1.70     1.08     1.08     1.70    2    3
 108 before_ppbafp_v2______________     0.00     0.00     0.00     0.00        0   16   24
 110 before_expikt_________________     0.13     0.09     0.10     0.06     0.13    2    3
 112 bfore_psi2b_v2________________     7.75     8.09     8.20     5.53     8.20    2    3
 114 after_drvmelp2________________   416.52   441.73   441.22   305.86   441.73    2    3
 116 end_exchange_section__________     0.06     0.06     0.06     0.04     0.06    2    3
 124 end_reading_3111______________     0.10     0.01     0.01     0.01     0.10    2    3
 126 end_check_freq_r______________     0.00     0.00     0.00     0.00        0    2    3
 128 matmul_zmel_ppovlz____________   693.07   779.49   696.65   537.76   779.49    2    3
 130 before_read_frcw_zw___________     2.94     3.06     3.05     1.96     3.06    2    3
 132 read_ifrwc_zw_before_matzwz2__     0.55     0.46     0.29     0.21     0.55    2    3
 134 after_matzwz2_________________   221.03   234.72   234.93   163.22   234.93    2    3
 136 before_ccc333_k-cycle_________     0.12     0.13     0.14     0.09     0.14    2    3
 138 end_do_xxx33333_k-cycle_______   793.33   857.67   857.03   591.11   857.67    2    3
 140 before_1385___________________     0.00     0.00     0.00     0.00        0    2    3
 142 end_1385_loop_________________    23.11    25.03    25.74    17.57    25.74    2    3
 144 end_yyy333_k-cycle____________     0.04     0.05     0.05     0.03     0.05    2    3
 146 end_of_deallocate_zmel_zmel1__     0.15     0.16     0.16     0.11     0.16    2    3
 148 end_of_3001___________________     0.06     0.07     0.07     0.05     0.07    2    3
 150 end_of_4001___________________     0.05     0.06     0.06     0.04     0.06    2    3
 152 before_2001___________________    30.41    34.49    22.55    17.76    34.49    2    3
 154 end_2001_loop_________________   344.38   373.06   383.63   274.26   383.63    2    3
 156 end_of_1100_loop______________     0.37     0.38     0.25     0.25     0.38    2    3
 200 after_readgeig_q_q_rk_________     0.21     0.16     0.13     0.08     0.21    2    3
 202 before_drv_melpln2____________     0.00     0.00     0.00     0.00        0    2    3
 204 after_drv_melpln2_____________   413.74   438.89   438.40   303.92   438.89    2    3
 206 end_drvmelpln2________________     1.69     1.77     1.76     1.21     1.77    2    3

after_drvmelp

              allocate( rmelt(ngb, nctot+nbmax, ntqxx), ! nstate= nctot+nband
     &              cmelt(ngb, nctot+nbmax, ntqxx))
               if(debug) print *, ' sxcf_fal1: goto drvmelp2 xxx111'
               allocate(drealzzzmel(nbloch,nt,ntqxx),dimagzzzmel(nbloch,nt,ntqxx))
               drealzzzmel=dreal(zzzmel)
               dimagzzzmel=dimag(zzzmel)
               call drvmelp2( q,             ntqxx, ! q in FBZ
     i              q-qbz_kr,  nbmax, ! q-rk
     i              qibz_k,     ! k in IBZ for e-product basis
     i              isp,ginv,
     i              ngc,ngcmx, ngpmx,nband,itq,
     i              symope, shtv, qbas, qbasinv,qibz,qbz,nqbz,nqibz,
c     i       dreal(zzzmel), dimag(zzzmel), nbloch, nt,nctot,
     i              drealzzzmel, dimagzzzmel, nbloch, nt,nctot,
     o              rmelt,cmelt)
               if(debug) print *, ' sxcf_fal1: end of drvmelp2'
               deallocate(drealzzzmel,dimagzzzmel)
               if(verbose()>50) call timeshow("5 after drvmelp")
               if(nbcut/=0.and.(.not.exchange)) then
                  do it= nctot+1,nctot+min(nbcut,nbmax)
                     rmelt(:, it,:) =0d0
                     cmelt(:, it,:) =0d0
                  enddo
               endif

matmul zmel ppovlz

               allocate( zmel (ngb, nstate, ntqxx) )
               zmel = dcmplx (rmelt,-cmelt)
               if(newaniso()) then
!$OMP parallel
!$OMP do
                  do itp=1,ntqxx
                     do it=1,nstate
                        zmel(:,it,itp) =  matmul(zmel(:,it,itp),dconjg(ppovlz(:,:)))
                     enddo
                  enddo
!$OMP end do nowait
!$OMP end parallel
               endif
               deallocate(rmelt,cmelt)

end do xxx33333 k-cycle

              allocate( zwix(nblochpmx,nblochpmx,nx),stat=ierr )
              if (ierr.ne.0) then
                 write(6,'(i3,a,i3,a,i5)')
     &    mpi__rank,'failed to allocate zwix',ierr,'size=',nblochpmx**2*nx; call flush(6)
                 call mpi_abort()
                 stop
              endif
               do ix=1,nx
                  nrec=ix
c     nrec=(kx-iqini)*niw+ix
c     if(bzcase()==2) nrec= (kx-1)*niw+ix

c                  read(ifrcwi,rec=nrec) zw ! direct access read Wc(0) = W(0) - v
                  read(ifrcwi,rec=nrec) zwix(:,:,ix)  ! direct access read Wc(0) = W(0) - v
c     read(ifrcwi,rec=((kx-2)*niw+ix)) zw  ! direct access read Wc(0) = W(0) - v
               enddo

!$OMP parallel  private(zwzix)
               write(6,*)mpi__rank,'zwzix allocate'; call flush(6)
               allocate(zwzix(1:nstate,1:ntqxx,1:ntqxx),stat=ierr)
               if (ierr.ne.0) then
                  write(6,'(2i5,a,i5,a,i7)')mpi__rank,omp_get_thread_num(),
     &  'failed to allocate zwzix',ierr,'size=',nstate*ntqxx**2; call flush(6)
                  call mpi_abort()
                  stop
               endif
!$OMP do
               do ix=1,nx
                   write(*,'(i3,a,20i7)') mpi__rank,'before matzwz2',
     &       ngb,ngb,size(zmel,dim=1),size(zmel,dim=2),size(zmel,dim=3),nstate,ntqxx,ntqxx
c                  call matzwz2(2, zw(1:ngb,1:ngb), zmel, ntqxx, nstate,ngb,
                  call matzwz2(2, zwix(1:ngb,1:ngb,ix), zmel, ntqxx, nstate,ngb,
     o                 zwzix(1:nstate,1:ntqxx,1:ntqxx)) ! zwz = zmel*(W(0)-v)*zmel
c     o                 zwz(ix,1:nstate,1:ntqxx,1:ntqxx)) ! zwz = zmel*(W(0)-v)*zmel
                  write(6,*)mpi__rank,'matzwz2 called',ix;call flush(6)
                  zwz(ix,1:nstate,1:ntqxx,1:ntqxx)= zwzix(1:nstate,1:ntqxx,1:ntqxx)
c     call matzwzs(zw(1:ngb,1:ngb), zmel, ntq, nstate,ngb,
c     o    zwz(ix,1:nstate,1:ntq,1:ntq)) ! zwz = zmel*(W(0)-v)*zmel
ccc   print *,'sum check2---',sum(abs(zwz(ix,:,:,:)))
               enddo ! ix
!$OMP end  do
               deallocate(zwzix)
!$OMP end parallel
               deallocate(zwix)
!TIME1 "end do xxx33333 k-cycle"

==end 2001 loop==
<pre>
              do 2001 itp = 1,ntqxx ! loop over states (q-k,n)
                  omg = omega(itp)
                  if (omg >= ef) then
                     itini= nt0m+1
                     itend= nt_max
                     iii=  1
                  else
                     itini= 1
                     itend= nt0p
                     iii= -1
                  endif

                  do 2011 it = itini,itend ! nt0p corresponds to efp
                     esmrx = esmr
                     if(it<=nctot) esmrx = 0d0
                     wfac = wfacx2(omg,ef, ekc(it),esmrx)
                     if(GaussSmear()) then
                        if(wfac<wfaccut) cycle ! next it
c     we = .5d0* abs( weavx2(omg,ef, ekc(it),esmr)- omg )
                        we = .5d0* abs( omg-weavx2(omg,ef, ekc(it),esmr) )
                        if(it<=nctot) then !faleev
                           if(wfac>wfaccut) stop "sxcf: it<=nctot.and.wfac/=0"
                        endif
...
                    if (iSigma_en == 1.or.iSigma_en==5) then
                        zwz3=(0d0,0d0)
!$OMP parallel do private(ix,zz2)
                        do ix0=1,3
                           ix=ixs+ix0-2
                           do igb2=2,ngb !**** most time consuming part for iSigma_en=1 ******
                              zz2=sum(zmel1(1:igb2-1)*zw3(1:igb2-1,igb2,iir*ix)  ) +
     &                             .5d0* zmel1(igb2)*zw3(igb2,igb2,iir*ix)
                              zwz3(ix0)=zwz3(ix0)+zz2*zmel(igb2,it,itp)
                           enddo !igb2
                           zwz3(ix0)=2d0*dreal(zwz3(ix0))+
     &                          zmel1(1)*zw3(1,1, iir*ix)*zmel(1,it,itp)
                        enddo   !ix
!$OMP end parallel do
                        if(npm==1) then
                           zsec(itp,itp,ip) = zsec(itp,itp,ip)
     .                          + wfac*alagr3z2(we,freq_r(ixs-1),zwz3,itp,itp)
                        else
                           zsec(itp,itp,ip) = zsec(itp,itp,ip)
     .                          + wfac*alagr3z(we,freq_r(ixs-1),zwz3)
                        endif
c     this contribution to zsec_nn is real (hermitean)
                     elseif(iSigma_en == 3) then
                        allocate(zwz4(3,ntqxx))
                        zwz4=(0d0,0d0)
CYY!$OMP parallel do private(ix,zz2)
                        do ix0=1,3
                           ix=ixs+ix0-2
                           do igb2=1,ngb !**** most time consuming part for iSigma_en=3 ******
                              zz2=sum(zmel1(1:ngb)*zw3(1:ngb,igb2, iir*ix)  )
                             do itpp=1,ntqxx
                                 zwz4(ix0,itpp)=zwz4(ix0,itpp)+ zz2*zmel1_(itpp,igb2,it)
                              enddo !itpp
                           enddo !igb2
                        enddo   !ix
CYY!$OMP end parallel do
                        do 2021 itpp=1,ntqxx
                           if(diagonly.and.(itpp/=itp)) cycle

                           if(npm==1) then
                              zsec(itp,itpp,ip) = zsec(itp,itpp,ip)
     .                             + wfac*alagr3z2(we,freq_r(ixs-1),zwz4(1,itpp),itp,itpp)
                           else
                              zsec(itp,itpp,ip) = zsec(itp,itpp,ip)
     .                             + wfac*alagr3z(we,freq_r(ixs-1),zwz4(1,itpp))
                           endif
 2021                   continue !itpp
                        deallocate(zwz4)
                     endif      ! inner iSigma_en=1 or 3
c     this contribution to zsec_nn' is not hermitean because W(e_n)
c     and must be made hermitean when zsec will be written on disc
 2011             continue
 2001          continue         !itp

openmp書き換え2

lx0, 8mpi

  1 end_program___________________  77.47 105.99  83.74  77.47 134.89    1    1
   3 after_readbzdata______________   0.73   0.73   0.73   0.73   1.47    1    1
   5 before_readngmx_______________   0.02   0.02   0.02   0.01   0.02    1    1
   7 before_initreadeigen__________   0.01   0.02   0.02   0.01   0.02    1    1
   9 before_histogram_division_____   0.03   0.03   0.03   0.02   0.03    1    1
  11 before_realomega_write3111____   0.00   0.00   0.00   0.00      0    1    1
  13 before_open_WVd_______________   0.46   0.00   0.00   0.00   0.46    1    1
  15 after_write_ifwd______________   0.00   0.00   0.00   0.00      0    1    1
  17 endo_of_ppbafp_v2_____________   0.01   0.01   0.01   0.01   0.01    1    1
  19 before_1001_loop______________   0.01   0.01   0.01   0.01   0.01    1    1
  21 before_1003_loop______________   0.08   0.09   0.09   0.08   0.12    1    1
  23 before_x0kf_v2hz______________   0.63   0.56   0.55   0.55   0.64    1    1
  25 x0vk_v4h______________________  66.06  99.67  76.79  66.06 128.35    1    1
  27 before_real_omega_____________   1.44   1.10   1.15   1.02   1.44    1    1
  29 real_omega____________________   5.67   3.52   4.10   3.05   5.67    1    1
  31 end_1015_loop_________________   5.18   3.48   4.04   3.02   5.18    1    1
  33 after_write_trace_msg_________   0.00   0.00   0.00   0.00      0    1    1
  35 after_matcinv_________________   3.81   2.91   3.29   2.54   3.81  163  163
  37 after_zw=zw0__________________   0.88   0.20   0.33   0.10   0.88  163  163
  39 end_of_tr_chkwrite____________   0.49   0.38   0.42   0.32   0.55  163  163
  45 imag_omega____________________   0.27   0.22   0.23   0.18   0.28    1    1
  47 end_eeeeeee___________________   0.00   0.00   0.00   0.00      0    1    1
  49 before_closing_hbe.d__________   2.03   0.00   0.00   0.00   2.03    1    1
 100 before_if_smbasis_____________   0.00   0.00   0.00   0.00      0    1    1
 102 after_if_smbasis______________   0.00   0.00   0.00   0.00      0    1    1
 104    0.00   0.00   0.00   0.00      0    1    1
 106 before_psicb_v3_______________   0.01   0.02   0.01   0.01   0.04    8   40
 108 after_psicb_v3________________   0.45   1.04   0.68   0.45   1.96    8   40
 110 before_melpln2________________   0.11   0.19   0.13   0.11   0.27    8   40
 112 after_melpln2_________________   9.85  20.93  14.63   9.85  32.50    8   40
 114 after_matmul_zmelt____________   7.83  30.31  28.28   7.83  37.99    8   40
 116 before_jpm_ibib_loop__________   0.00   0.00   0.00   0.00      0    8   40
 118 after_rcxq____________________  16.49  40.79  25.15  16.49  54.68    8   40
 120 before_eibzmoden=2____________   0.07   0.16   0.23   0.07   0.27    1    1
 122 after_iqindx2_________________   0.00   0.00   0.00   0.00      0    1    1
 124 before_1011___________________   0.02   0.02   0.02   0.02   0.02    1    1
 126 call_cputid2_qqqqq11__________   0.01   0.00   0.01   0.00   0.01    1    1
 128 end_matmmsparse_______________   2.83   0.34   0.70   0.09   2.83    2   48
 130 before_iagain=1_______________   0.06   0.01   0.01   0.00   0.06    2   48
 132 end_ig_itimer_icc,_nrotm______   0.00   0.00   0.00   0.00      0    2   48
 134 after_sym_rcxq________________  28.17   5.70   6.77   0.50  28.17    1    1

lsc , 8mpi

  1 end_program___________________ 146.01 146.00 146.00 146.00 146.01    1    1
   3 before_2000_loop______________   1.42   1.42   1.42   1.42   1.42    1    1
   5 before_sxcf_fal3_scz__________   0.00   0.00   0.00   0.00      0    1    1
   7 after_sxcf_fal3_scz___________ 136.98 127.34 127.94 110.92 144.52    1    1
   9 reduction_nbandmx_end_________   7.53  17.18  16.58   0.00  33.60    1    1
  11 end_legas_exchange____________   0.00   0.00   0.00   0.00      0    1    1
  13 end_of_allocate_zsect_________   0.01   0.01   0.01   0.01   0.01    1    1
  15 end_3020_loop_________________   0.05   0.05   0.05   0.05   0.05    1    1
  17 before_2000_continue__________   0.01   0.00   0.00   0.00   0.01    1    1
  19 end_of_exspectrum_____________   0.00   0.00   0.00   0.00      0    1    1
 100 end_of_sxcf_fal2_scz__________ 136.98 127.34 127.94 110.92 144.52    1    1
 102 before_1001___________________   0.00   0.00   0.00   0.00      0    1    1
 104 before_1100___________________   0.01   0.01   0.00   0.00   0.02    1    2
 106 before_1000___________________   0.24   0.16   0.14   0.07   0.24    4   14
 108 before_ppbafp_v2______________   0.00   0.00   0.00   0.00      0   96  336
 110 before_expikt_________________   0.08   0.08   0.08   0.08   0.09   16   16
 112 bfore_psi2b_v2________________   0.76   0.70   0.77   0.70   0.83   16   16
 114 after_drvmelp2________________  40.23  41.17  43.50  38.54  47.60   16   16
 116 end_exchange_section__________   0.00   0.00   0.00   0.00      0   16   16
 124 end_reading_3111______________   0.01   0.01   0.01   0.01   0.01   16   16
 126 end_check_freq_r______________   0.00   0.00   0.00   0.00      0   16   16
 128 matmul_zmel_ppovlz____________  33.80  31.03  23.32  20.22  33.80   16   16
 130 before_read_frcw_zw___________   0.30   0.30   0.52   0.26   0.52   16   16
 132 read_ifrwc_zw_before_matzwz2__   0.02   0.02   0.02   0.01   0.02   16   16
 134 after_matzwz2_________________   4.80   4.23   4.61   4.00   4.85   16   16
 136 before_ccc333_k-cycle_________   0.01   0.01   0.01   0.00   0.01   16   16
 138 matzwz2_read_zw_______________   0.20   0.18   0.20   0.14   0.25   16   16
 140 end_do_xxx33333_k-cycle_______  48.61  42.54  46.68  40.03  50.05   16   16
 142 before_1385___________________   0.00   0.00   0.00   0.00      0   16   16
 144 end_1385_loop_________________   2.25   1.93   2.15   1.76   2.38   16   16
 146 end_yyy333_k-cycle____________   0.01   0.00   0.01   0.00   0.01   16   16
 148 end_of_deallocate_zmel_zmel1__   0.02   0.03   0.02   0.02   0.04   16   16
 150 end_of_3001___________________   0.02   0.02   0.02   0.02   0.02   16   16
 152 end_of_4001___________________   0.02   0.02   0.02   0.01   0.02   16   16
 154 before_2001___________________   1.89   2.01   2.23   1.46   2.34   16   16
 156 end_2001_loop_________________   3.67   2.87   3.60   2.65   4.02   16   16
 158 end_of_1100_loop______________   0.03   0.02   0.02   0.01   0.03   16   16
 160 last_proc_in_sxcf_fal2_scz____   0.00   0.00   0.00   0.00      0    1    2
 200 after_readgeig_q_q_rk_________   0.09   0.09   0.10   0.09   0.10   16   16
 202 before_drv_melpln2____________   0.00   0.00   0.00   0.00      0   16   16
 204 after_drv_melpln2_____________  39.89  40.82  43.14  38.19  47.20   16   16
 206 end_drvmelpln2________________   0.17   0.19   0.18   0.17   0.22   16   16

matmul_zmel_ppovlz____________ だけopenmpを書く。 lx0 8mpi 4 smp

  1 end_program___________________  65.47  78.94  48.22  48.22 107.86    1    1
   3 after_readbzdata______________   0.41   0.41   0.52   0.41   0.52    1    1
   5 before_readngmx_______________   0.03   0.03   0.01   0.01   0.03    1    1
   7 before_initreadeigen__________   0.02   0.02   0.01   0.01   0.02    1    1
   9 before_histogram_division_____   0.05   0.05   0.03   0.03   0.05    1    1
  11 before_realomega_write3111____   0.00   0.00   0.00   0.00      0    1    1
  13 before_open_WVd_______________   0.39   0.00   0.00   0.00   0.39    1    1
  15 after_write_ifwd______________   0.00   0.00   0.00   0.00      0    1    1
  17 endo_of_ppbafp_v2_____________   0.01   0.01   0.01   0.01   0.01    1    1
  19 before_1001_loop______________   0.01   0.01   0.01   0.01   0.02    1    1
  21 before_1003_loop______________   0.08   0.10   0.09   0.08   0.11    1    1
  23 before_x0kf_v2hz______________   0.63   0.55   0.55   0.55   0.63    1    1
  25 x0vk_v4h______________________  55.31  72.76  41.74  41.74 102.24    1    1
  27 before_real_omega_____________   1.24   1.13   1.18   1.02   1.24    1    1
  29 real_omega____________________   5.15   3.65   3.85   3.05   5.15    1    1
  31 end_1015_loop_________________   4.67   3.62   3.80   3.01   4.67    1    1
  33 after_write_trace_msg_________   0.00   0.00   0.00   0.00      0    1    1
  35 after_matcinv_________________   3.36   3.00   3.10   2.58   3.36  163  163
  37 after_zw=zw0__________________   0.89   0.28   0.29   0.10   0.89  163  163
  39 end_of_tr_chkwrite____________   0.41   0.34   0.41   0.31   0.44  163  163
  45 imag_omega____________________   0.26   0.22   0.22   0.19   0.26    1    1
  47 end_eeeeeee___________________   0.00   0.00   0.00   0.00      0    1    1
  49 before_closing_hbe.d__________   1.87   0.00   0.00   0.00   1.87    1    1
 100 before_if_smbasis_____________   0.00   0.00   0.00   0.00      0    1    1
 102 after_if_smbasis______________   0.00   0.00   0.00   0.00      0    1    1
 104    0.00   0.00   0.00   0.00      0    1    1
 106 before_psicb_v3_______________   0.01   0.02   0.01   0.01   0.04    8   40
 108 after_psicb_v3________________   0.42   1.04   0.66   0.42   1.99    8   40
 110 before_melpln2________________   0.09   0.18   0.11   0.09   0.28    8   40
 112 after_melpln2_________________   8.30  20.84  12.11   8.30  32.20    8   40
 114 after_matmul_zmelt____________   2.59   6.62   3.44   2.59  10.72    8   40
 116 before_jpm_ibib_loop__________   0.00   0.00   0.00   0.00      0    8   40
 118 after_rcxq____________________  16.29  37.69  19.92  16.29  57.14    8   40
 120 before_eibzmoden=2____________   0.08   0.16   0.16   0.08   0.17    1    1
 122 after_iqindx2_________________   0.00   0.00   0.00   0.00      0    1    1
 124 before_1011___________________   0.02   0.02   0.02   0.02   0.02    1    1
 126 call_cputid2_qqqqq11__________   0.00   0.00   0.00   0.00      0    1    1
 128 end_matmmsparse_______________   2.85   0.35   0.67   0.09   2.85    2   48
 130 before_iagain=1_______________   0.05   0.01   0.01   0.00   0.05    2   48
 132 end_ig_itimer_icc,_nrotm______   0.00   0.00   0.00   0.00      0    2   48
 134 after_sym_rcxq________________  24.43   5.68   4.46   0.51  24.43    1    1

こちらは変わらない。

lsc, 8mpi 4 smp

   1 end_program___________________ 374.66 374.63 374.62 374.62 374.66    1    1
   3 before_2000_loop______________   1.56   1.56   1.56   1.55   1.56    1    1
   5 before_sxcf_fal3_scz__________   0.00   0.00   0.00   0.00      0    1    1
   7 after_sxcf_fal3_scz___________ 355.70 325.51 351.73 306.34 370.32    1    1
   9 reduction_nbandmx_end_________  17.14  47.33  21.12   2.52  66.51    1    1
  11 end_legas_exchange____________   0.00   0.00   0.00   0.00      0    1    1
  13 end_of_allocate_zsect_________   0.01   0.01   0.02   0.01   0.03    1    1
  15 end_3020_loop_________________   0.20   0.20   0.19   0.18   0.20    1    1
  17 before_2000_continue__________   0.03   0.00   0.00   0.00   0.03    1    1
  19 end_of_exspectrum_____________   0.00   0.00   0.00   0.00      0    1    1
 100 end_of_sxcf_fal2_scz__________ 355.70 325.51 351.73 306.32 370.32    1    1
 102 before_1001___________________   0.00   0.00   0.00   0.00      0    1    1
 104 before_1100___________________   0.00   0.02   0.01   0.00   0.02    1    2
 106 before_1000___________________   0.81   0.54   0.48   0.22   0.81    4   14
 108 before_ppbafp_v2______________   0.00   0.00   0.00   0.00      0   96  336
 110 before_expikt_________________   0.10   0.13   0.15   0.10   0.26   16   16
 112 bfore_psi2b_v2________________   2.56   2.42   2.50   2.38   2.76   16   16
 114 after_drvmelp2________________ 122.74 115.34 126.63 109.66 127.17   16   16
 116 end_exchange_section__________   0.00   0.00   0.00   0.00      0   16   16
 124 end_reading_3111______________   0.21   0.25   0.21   0.19   0.27   16   16
 126 end_check_freq_r______________   0.00   0.00   0.00   0.00      0   16   16
 128 matmul_zmel_ppovlz____________  19.96  17.97  16.86  16.69  20.39   16   16
 130 before_read_frcw_zw___________   0.70   0.70   1.13   0.59   1.13   16   16
 132 read_ifrwc_zw_before_matzwz2__   0.01   0.01   0.09   0.01   0.09   16   16
 134 after_matzwz2_________________  16.80  15.15  16.27  14.18  17.39   16   16
 136 before_ccc333_k-cycle_________   0.02   0.00   0.08   0.00   0.08   16   16
 138 matzwz2_read_zw_______________   0.52   0.47   0.54   0.40   0.63   16   16
 140 end_do_xxx33333_k-cycle_______ 168.27 151.61 163.67 142.72 175.66   16   16
 142 before_1385___________________   0.00   0.00   0.00   0.00      0   16   16
 144 end_1385_loop_________________   8.79   7.52   8.45   6.90   9.20   16   16
 146 end_yyy333_k-cycle____________   0.00   0.00   0.04   0.00   0.04   16   16
 148 end_of_deallocate_zmel_zmel1__   0.11   0.07   0.05   0.04   0.13   16   16
 150 end_of_3001___________________   0.05   0.09   0.05   0.05   0.11   16   16
 152 end_of_4001___________________   0.08   0.03   0.09   0.03   0.09   16   16
 154 before_2001___________________   4.07   4.77   4.77   3.73   5.07   16   16
 156 end_2001_loop_________________   9.65   8.22   9.47   7.36  10.56   16   16
 158 end_of_1100_loop______________   0.22   0.07   0.08   0.02   0.22   16   16
 160 last_proc_in_sxcf_fal2_scz____   0.00   0.00   0.00   0.00      0    1    2
 200 after_readgeig_q_q_rk_________   0.24   0.30   0.15   0.15   0.30   16   16
 202 before_drv_melpln2____________   0.00   0.00   0.00   0.00      0   16   16
 204 after_drv_melpln2_____________ 121.77 114.50 125.69 108.85 126.14   16   16
 206 end_drvmelpln2________________   0.58   0.42   0.54   0.40   0.58   16   16

matmul_zmel_ppovlz____________ は高速化されたが、他が3倍以上鈍くなった。全体もかなり鈍くなった。




bug

2f94f2b7c737526bf44b9173ff57d3b93af72186 mpikが動かない。

Exit -1 zhev_tk2: nev /=nevx something wrong.

ctrls.yco5

HEADER    YCo5, hexagonal 
STRUC     ALAT=9.313 
          PLAT=0.8660254 -.5000000 0.0000000 
               0.0000000 1.0000000 0.0000000 
               0.0000000 0.0000000 0.8060000 
SITE      ATOM=Y   POS=0.00000000 0.00000000 0.00000000 
          ATOM=Co  POS=-.28867513 0.50000000 0.00000000 
          ATOM=Co  POS=0.28867513 -.50000000 0.00000000 
          ATOM=Co2 POS=0.43301270 0.25000000 0.40300000 
          ATOM=Co2 POS=0.43301270 -.25000000 0.40300000 
          ATOM=Co2 POS=0.00000000 0.50000000 0.40300000 
SPEC     ATOM=Y   Z=39 
          ATOM=Co  Z=27 
          ATOM=Co2 Z=27 

subs/rdsigm2.F, rv_p_oqpの構造

意味不明・・・。

57         real(8),pointer :: rv_p_oqp(:) =>NULL()
393        allocate(rv_p_oqp(abs(3*mxkp)))
394        if (3*mxkp<0) rv_p_oqp(:)=0.0d0
404       if (-mxkp<0) wgt_rv(:)=0.0d0
425      .  , nsgrps , ipq_iv , rv_p_oqp , wgt_rv , nqps , mxkp , gstar_iv
459       if (procid .eq. master) then
460         do  isp = 1, nsp
481           goto 8888
510           do  iq1 = 1, nqps
514              call dpsadd ( tmp , rv_p_oqp , 3 , 1 , 3 * iq1 - 2 , - 1d0 )
519             if (.not. latvec(1,tolq,plat,tmp)) then
520                call dpscop ( rv_p_oqp , tmp , 3 , 3 * iq1 - 2 , 1 , 1d0 )
522               if (lssym .ge. 4) then
532               endif
533             endif
562             if (alf(2) .ne. 0d0) then
565                call dpsadd ( tmp , rv_p_oqp , 3 , 1 , 3 * iq1 - 2 , - 1d0 )
568               if (abs(tmp(1))+abs(tmp(2))+abs(tmp(3)) .gt. tolq) then
569                  call dpscop ( rv_p_oqp , tmp , 3 , 3 * iq1 - 2 , 1 , 1d0 )
578               endif
589             endif
614           enddo
621  8888     continue
676         enddo
680       endif
706       if (alf(1) .ne. 1 .or. alf(2) .ne. 0) then
711         if (lssym .ge. 4) j1 = 2
712         if (procid .eq. master) then
719      .         , rsstol , i , rv_p_oqp , nbas , 0 , rotm , iwdummy )
725         endif
726       endif
728       if ( .not. ( cmdopt ( '--wsig' , 6 , 0 , outs ) .or.cmdopt (
730          if (associated(rv_p_oqp)) deallocate(rv_p_oqp)
732       endif
844       if (lwsig .eq. 0) then
845         if (cmdopt('--wsig',6,0,outs) .or. cmdopt('-wsig',5,0,outs)) then
846           if (procid .eq. master) then
893             if (lnwmsh) then
900                rv_p_oqp => sbz%rv_p_oqp  ------------ needless.
905                allocate(rv_p_oqp(abs(3*mxkp)))  ---------- delete allocated array first
906                if (3*mxkp<0) rv_p_oqp(:)=0.0d0
923               if (lfbzout) then
927      .            , 0 , ipq_iv , rv_p_oqp , wgt_rv , nqp , mxkp , 0 , 0 )
936      .            , nsgrp , ipq_iv , rv_p_oqp , wgt_rv , nqp , mxkp , gstar_iv
942               endif
944             elseif (lfbzout) then
947               rv_p_oqp => sham%rv_p_oqsig --- delete allocated array first
954             endif
981             if (lqoffo .ne. 0) then
983               do  i = 1, nqp
984                  call dmsadd ( rv_p_oqp , 1 , qoffo , 1 , 1 , 3 , 1 , 1 , 3 *
987               enddo
989             endif
995      .          , rsstol , nqp , rv_p_oqp , nbas , lrot , rotm , delt_rv )
1006           endif
1007         endif
1008       endif
1018       if (lwsig .ne. 0) then
1022         if (lnwmsh) then
1029            rv_p_oqp => sbz%rv_p_oqp  --- delete allocated array first
1035         endif
1037         if (lfbzout) then
1040            allocate(rv_p_oqp(abs(3*mxkp))) --- delete allocated array first
1041            if (3*mxkp<0) rv_p_oqp(:)=0.0d0
1060      .      , 0 , ipq_iv , rv_p_oqp , wgt_rv , nqp , mxkp , 0 , 0 )
1068         endif
1071         if (lqoffo .ne. 0) then
1073           do  i = 1, nqp
1074              call dmsadd ( rv_p_oqp , 1 , qoffo , 1 , 1 , 3 , 1 , 1 , 3 *
1077           enddo
1079         endif
1083           call getqp ( 1 , - ifiz , nqp , nkxyz , lshft , 0 , rv_p_oqp
1093       endif

subs/rdsigm2.F, rv_p_oqsigの構造

55:         real(8),pointer :: rv_p_oqsig(:) =>NULL()
384:       allocate(rv_p_oqsig(abs(3*mxkp)))
385:       if (3*mxkp<0) rv_p_oqsig(:)=0.0d0
403:     .  , 0 , ipq_iv , rv_p_oqsig , wgt_rv , nqsig , mxkp , 0 , 0 )
409:       sham%rv_p_oqsig => rv_p_oqsig

874                if (lnwmsh) then
925                elseif (lfbzout) then
927:              rv_p_oqp => rv_p_oqsig
934                endif

1303:        real(8),pointer :: rv_p_oqsig(:)
1604:       allocate(rv_p_oqsig(abs(3*mxkp)))
1605:       if (3*mxkp<0) rv_p_oqsig(:)=0.0d0
1622:     .  , 0 , ipq_iv , rv_p_oqsig , wgt_rv , nqsig , mxkp , 0 , 0 )
1628:       sham%rv_p_oqsig => rv_p_oqsig

subs/rdsigm2.F, rv_p_ohrsの構造

49:      real(8),pointer :: rv_p_ohrs(:) =>NULL()
374       if (hreal .eq. 1) then

375:         allocate(rv_p_ohrs(ndhrs**2*nttabs*nsp))
376:         rv_p_ohrs(:)=0d0
377       else

378:         allocate(rv_p_ohrs(2*ndhrs**2*nttabs*nsp))
379:         rv_p_ohrs(:)=0d0
380       endif


452       if (procid .eq. master) then
453          do  isp = 1, nsp
503             do  iq1 = 1, nqps
607             enddo

645             call hft2rs ( i , nk1 , nk2 , nk3 , k1 , k2 , k3 , qoff , isp
654:     .           , plat , slat%rv_p_opos ,iv_a_ontabs , iv_a_oiaxs , ndhrs , rv_p_ohrs ) 
668          enddo
672       endif


690:      sham%rv_p_ohrs => rv_p_ohrs 
698       if (alf(1) .ne. 1 .or. alf(2) .ne. 0) then
701       else
704          if (procid .eq. master) then
706             call chksgr ( j1 , ltrans , kcplx , plat , nsp , ndimh , ifis
709:     .           , hreal , sham%iv_p_oindxo , nttabs , iv_a_oiaxs , rv_p_ohrs , ndhrs 
716          endif
717       endif

741       if (mod(lrsig,4) .ge. 2 .and. nsgrp .gt. 1) then
728:         call mpibc1 ( rv_p_ohrs , ndhrs * * 2 * nttabs * nsp , 4 , mlog 
732:         call mpibc1 ( rv_p_ohrs , 2 * ndhrs * * 2 * nttabs * nsp , 4 
755:         nhrss= size(rv_p_ohrs)
757:         hrss = rv_p_ohrs
773          call rsmsym ( i , plat , mxorb , sham%iv_p_oindxo , ndimh , nbas ,
779:     .        , slat%rv_p_osymgr , istb2_iv , nsgrp , ndhrs , hrss , rv_p_ohrs ) 

804:     .              , hreal , sham%iv_p_oindxo , nttabs , iv_a_oiaxs , rv_p_ohrs , ndhrs 
814          if (hreal .eq. 1) then
815:            call mpibc1 ( rv_p_ohrs , ndhrs * * 2 * nttabs * nsp , 4 , mlog 
818          else
819:            call mpibc1 ( rv_p_ohrs , 2 * ndhrs * * 2 * nttabs * nsp , 4 
822          endif
824       endif

827       if (lwsig .eq. 0) then
828          if (cmdopt('--wsig',6,0,outs) .or. cmdopt('-wsig',5,0,outs)) then
829             if (procid .eq. master) then
838                if (lonesp .eq. 1) then
842                   call siged ( 1 , nbas , nsp , ndhrs , plat , slat%rv_p_opos , ndimh
846:     .                 , sham%iv_p_oindxo , hreal , iv_a_ontabs , iv_a_oiaxs , rv_p_ohrs 
857                   call siged ( 0 , nbas , nsp , ndhrs , plat , slat%rv_p_opos , ndimh
861:     .                 , sham%iv_p_oindxo , hreal , iv_a_ontabs , iv_a_oiaxs , rv_p_ohrs 
868                endif

970                call chksgr ( 1 , ltrans , kcplx , plat , nsp , ndimh , ifis2
973:     .              , hreal , sham%iv_p_oindxo , nttabs , iv_a_oiaxs , rv_p_ohrs , ndhrs 
985             endif
986          endif
987       endif


1297:      real(8),pointer :: rv_p_ohrs(:) 

subs/mksym.Fの構造

89       integer,pointer :: iv_p_oipcp(:) =>NULL()
91       integer,pointer :: iv_p_oipc(:) =>NULL()
97       integer,pointer :: iv_p_onrcp(:) =>NULL()
99       integer,pointer :: iv_p_onrc(:) =>NULL()

151       allocate(iv_p_oipc(abs(nsite)))
152       if (nsite<0) iv_p_oipc(:)=0
309       allocate(iv_p_onrc(abs(nspec)))
310       if (nspec<0) iv_p_onrc(:)=0

313       call icopy ( nsite , sarray%iv_p_oips , 1 , iv_p_oipc , 1 )
320      .     , nsgrp , slat%iv_p_oistab , nspec , slabl , nclass , iv_p_oipc ,
322      .     sarray%iv_p_oics , iv_p_onrc )
370       iv_p_oipcp => iv_p_oipc
370       iv_p_oipcp => iv_p_oipc
429       do i_spackv=1,nbasp
430          call spackv_array_copy_i8_i ( 'p' , ssite ( i_spackv ) %class
431      .        , i_copy_size , i_spackv + 1 - 1 , iv_p_oipc )
433       enddo


436       if (mod(mode/10,10) .eq. 2 .or. mod(mode/10,10) .eq. 4) then

442          allocate(iv_p_onrcp(abs(nclspp)))
442          allocate(iv_p_onrcp(abs(nclspp)))
443          if (nclspp<0) iv_p_onrcp(:)=0
443          if (nclspp<0) iv_p_onrcp(:)=0
446          call pvsym2 ( 3 , 2 * nbasp - nbas , nclspp , sarray%iv_p_oics , iv_p_oipcp
448      .        , nspec , slabl , ssite , sarray%rv_p_oclabl , iv_p_onrcp )
455          iv_p_onrc => iv_p_onrcp
455          iv_p_onrc => iv_p_onrcp

463       else

468          nullify(iv_p_onrcp)
468          nullify(iv_p_onrcp)

469       endif

480       sarray%iv_p_onrcp => iv_p_onrcp
482       sarray%iv_p_oipcp => iv_p_oipcp
488       sarray%iv_p_oipc => iv_p_oipc
493       sarray%iv_p_onrc => iv_p_onrc

fp/lmfp.F,rv_p_oposの構造


113:          real(8),pointer :: rv_p_opos(:) =>NULL()

371       if ( iand(1,int(sctrl%lbas)) .ne. 0 ) then
378:         rv_p_opos => slat%rv_p_opos
396:           call shorps ( nbas , plat , ix , pos2_rv , rv_p_opos )
401                call spackv_array_copy_r8_r8 ( 'p' , ssite ( i_spackv ) %pos
402:     .       , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** ssite(i_spackv)%pos=rv_p_opos
412       endif

448:       rv_p_opos => slat%rv_p_opos
470       if (nitrlx .gt. 0) then
519:           call shorps ( nbas , plat , ix , rv_p_opos , pos2_rv )
528       endif

533     4 continue

557 C     --- Setup for iterations in a self-consistency cycle ---

664 C     ---------------- Re-entry point for a new iteration ---------------
665       iter = 1
666     5 continue

692       elseif (mod(irs(1),4) .ge. 1 .and. mod(irs(1),4) .le. 2) then
727:         rv_p_opos => slat%rv_p_opos
732             call spackv_array_copy_r8_r8 ( 'u' , ssite ( i_spackv ) %pos
733:     .     , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** rv_p_opos=ssite(i_spackv)%pos
776       endif

786 C     ... Write positions after file read, and repack

808 C     --- Optionally re-shorten basis vectors ---
809       if (cmdopt('--shorps',8,0,strn)) then
822:         call shorps ( - nbas , plat , ix , pos2_rv , rv_p_opos )  !** This rv_p_opos must be slat%rv_p_opos
827:         call iopos ( .true. , 1 , 'shorps' , nbas , rv_p_opos )
829:         call shorps ( nbas , plat , ix , pos2_rv , rv_p_opos )
834             call spackv_array_copy_r8_r8 ( 'p' , ssite ( i_spackv ) %pos
835:     .     , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** site(i_spackv)%pos= rv_p_opos
844       endif

947       if (maxit .gt. 0) then
1002          if (lsc .gt. 2) then
1003             irs(1) = 3
1004             goto 5
1005 C     kino      ---------------- SCF (iteration) loop end --------------------
1006          endif


1019          if (nitrlx .gt. 0 .and. lsc .le. 2) then
1021 C     Buglet: this eats up a little memory.  Never released.
1022:           allocate(rv_p_opos(abs(3*nbas)))
1023:           if (3*nbas<0) rv_p_opos(:)=0.0d0
1027                call spackv_array_copy_r8_r8 ( 'u' , ssite ( i_spackv ) %pos
1028:     .       , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** rv_p_opos= ssite(i_spackv)%pos
   1044                   goto 98
1058:     .       , natrlx , nvrelx , frc_rv , p_rv , w_rv , 0 , 0d0 , rv_p_opos
1083:                call fixpos ( rv_p_opos , nbas , fptol , ng , plat , slat%rv_p_osymgr
   1093 C     Write updated positions to bsmv file
1102:             call shorps ( nbas , plat , ix , rv_p_opos , pos2_rv )
1115                call spackv_array_copy_r8_r8 ( 'p' , ssite ( i_spackv ) %pos
1116:     .       , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** site(i_spackv)%pos= rv_p_opos
1149:     .         , rv_p_opos )
1170                   call spackv_array_copy_r8_r8 ( 'p' , ssite ( i_spackv ) %pos
1171:     .         , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** site(i_spackv)%pos= rv_p_opos
1218:             call rdistn ( rv_p_opos , wk_rv , nbas , gam ( 1 ) , gam ( 2
   1253             if (lshr) then
   1254                goto 98
   1255             else
   1256                goto 5
   1257             endif
1258          endif  ! if (nitrlx .gt. 0 .and. lsc .le. 2) then  !*** upto 1258, rv_p_opos must be a local array

1259       endif  ! if (maxit.gt.0

1262       if (cmdopt('--wpos=',7,0,strn) .or.
1264:         allocate(rv_p_opos(abs(3*nbas)))
1265:         if (3*nbas<0) rv_p_opos(:)=0.0d0
1269             call spackv_array_copy_r8_r8 ( 'u' , ssite ( i_spackv ) %pos
1270:     .     , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** rv_p_opos= ssite(i_spackv)%pos
1280:         call rdistn ( rv_p_opos , wk_rv , nbas , gam ( 1 ) , gam ( 2
1294:         if (associated(rv_p_opos)) deallocate(rv_p_opos)
1297       endif ! if (cmdopt('--wpos='   !*** in this section rv_p_opos is a local array

  1310       return

   1313 C     --- Setup to start calculation at new shear ---
   1314  98   continue

1328:       call dgemm ( 'N' , 'N' , 3 , nbas , 3 , 1d0 , xv , 3 , rv_p_opos  !** this rv_p_opos must be a local array defined at L1022.
1334:       rv_p_opos => slat%rv_p_opos
1338:     . , 3 , 0d0 , rv_p_opos , 3 )
1344          call spackv_array_copy_r8_r8 ( 'p' , ssite ( i_spackv ) %pos
1345:     .   , i_copy_size , i_spackv + 1 - 1 , rv_p_opos ) !** ssite(i_spackv)%pos=rv_p_opos

   1402       goto 4

scalar relativistic

scalar relativisticなradial 波動関数は、rseqで解かれる。

bndfp.F->augmat.F->potpus->makrwf->rseqで呼び出されている。(他は?)

scalar relativisticな波動関数は二成分の波動関数である(majorityとminority. rseq.Fではg(nr,2)と書かれる)。境界条件はmajority成分のMT端での値とMT端での微分値(とノードの数)で与えられる。minority成分のそれらについて与える必要はない。

MT外では、majority成分のみを扱っている。一方でMT内部では、pi,sigma,tau積分などにおいて、二成分であるとして取り扱う。 potpus.Fを見よ。

Frozen coreがスピン分極していることに注意

現状では、lmfaでMMOMで原子のスピン分極を指定してfrozen coreを作る(空間的にはclosedであり球対称)。これは場合によっては問題になる可能性がある。たとえば、コアの分極とバレンスの分極(揃っているはず)が原子計算において仮定したものと逆向きになって収束してしまうということもありえる。

汎用性を考えれば、スピン分極に関しては平均化したコアを用いたほうが、いたずらな計算誤差を生まないため、 より良いかもしれない。lmfaを二回繰り返すのがありえる。

  1. コア密度に関してはlmfa実行をMMOMなし(分極なし)のself-consistent計算でおこなう。これでfrozen coreとそのfittingを作っておく。
  2. 上で作ったコア電子密度を固定した状態であたえたMMOMに対してlmfa計算を行う。
  • 磁性の計算のエネルギースケールは非常に小さいので注意を要する。

基底の決めかたに関して(小谷による)

「PMTだと基底を引き抜きながらでないとpwemaxを上げていけない」というのはまずは認めざるを得ない。ただしこれはAPWのみの場合(LAPW)でも同様であり、基底を増やしていくと、十分に収束が確認できたと思える以前になんらかの形で線形独立性が悪くなって破綻することがよくある。

  • R(MT半径)を小さくすると、計算は安定するが、Rをあまり小さくしたくない理由もある。
  1. Froze core近似。通常LFOCA=1のFroze core近似を使うが、これが破綻しうる。以下のGaAs、BiVO3など、MT半径依存性が非常に大きいとき、その原因を調べるとこれだった。小さいMT半径では計算が破綻している---たちが悪いのは一見、別な所へ収束して問題ないように見える点(core spilt chargeが大きいとき。GaAsのAsではR=2.0程度でもまずかった)。たとえば、以下のGaAsの例だと、ある程度以上にAsのMTが小さいと、Asの3dをcore扱いできなくなる。
  2. 通常、pwemax=2程度ではphi,phidotで張るほうが有利なのでMTを大きくしたほうが全エネルギーを下げれる。
  3. Rを小さくすると、smooth Hankel(RSMH=R/2にとっている)が原点近傍で発散するような形状のものにちかづく。このためgmax(実空間メッシュのカットオフ)を大きくとらないといけないということになる。
  • 現在、R(MT半径)の決め方は、charge neutralityを尊重してtouchingにとるのを基本にしている。ctrlgen.pyが内部でlmchkを呼び出している。これがよいかどうかは不明。covalent radiusを尊重してもいいようにも思っている。(場合によればスケールする必要もある)。
  • 将来的には、Rを小さめの値で固定してしまうことも考えうる。EHやEH2、さらにはphi,phidotも固定する。こういう方法であれば、エネルギー差の系統性は得やすいだろうし、基底関数を固定することになるので重なり行列を計算しなおす必要もないし、計算も高速化・安定化しうる。非線形な部分が入らないので、エネルギー最小化の意味できれいなformalismになる(Forceとエネルギー差の整合性も完全になる)。フォノン計算のみでなく、結号エネルギー計算にも有利かもしれない。

参考資料

物質案(kino)

  • water dimer File:water_dimer.png

R依存性の長距離のきれいなカーブが難しい。

  • graphite/NFE

物質案(小谷)

  • Fe bcc/fcc

Rの決め方

こんな感じに見える。--Kino 16:26, 8 February 2011 (JST)

まず、core spilloutやfrozen core approximationの良し悪しはRとは無関係であり、 PMTではmuffin-tin半径をtouching/overlappingに取らないのでcore spilloutは無意味なconcept。

以下は図の説明。

  • valence spillout:Rが小さいときに顕著に見えるが、pwemaxを大きく取らないと結果が収束しないことがある。MT半径内+tailでは十分に波動関数を表現できないという当たり前の理由である。
  • orthogonality:Rが大きすぎるとorthogonalityとphi,phidotの計算が破綻する。Rの上限が求まる。
  • inaccuracy inside MT: MT内はphi,phidotという2つの基底で表現されている。scf毎に計算し直すがこのためRが大きくなる程MT内が不正確になる。大きなmagnetic momentが出る場合などに不正確になりやすい。
  • good parameter:pwemaxが大きいところまで破綻せずに計算できるパラメタの方がいい(図のgood parameter)が、pwemaxが小さいところだけ対象にするのならばもっと広い範囲がパラメタになりうる。問題は物質によらず原子だけによるこの領域があるかということ。

多くの自由度があり一度に決めるのは困難。安定することを最優先してパラメタを決める方法を考える。

  1. Rをある程度決めて
  2. それでもorthogonalityで破綻する部分は軌道削除で調整
  3. EHの値による変化は[-2.0:-0.1]にとるならばさほど大きくはないので安定する値にとって構わない。O2などEHの変化がpwemaxの変化と同程度あるが、EHの変化に対しては安定して計算できるので後で決めれば十分。これらの意味でEHの値は最後に決めればいい問題。valence spilloutと関連し、例えば、d電子のmagneticな問題を考えるには1)pwemaxが小さいならば-EHをある程度大きくしてd電子のtailを表現するか、2)pwemaxを大きくとる、と図のvalence spilloutの線を下に下げられる。

メモ

  • -EH=0.1はpwemax=2ではmoleculeのmagnetic solutionには小さすぎる。
  • bulk closed packの物質がorthogonalityが破綻しやすいのでそこからRとEHを決めるのがいいかもしれない。(Feにて検証中)

論文の取りまとめに関するメモ(小谷)

  • どういうデータをとるかは未定だが上述の文献を参考にする。分子に関しては、結合長、振動数、結合エネルギーなどをみたい。pwemaxが小さくてもやれると主張したい。H2Oの緩和計算もしたい。Gaussianと比較できるものも可能かも。
  • 全エネルギー収束チェックにこだわらないほうがいい。他文献をみても他論文をみても全エネルギーの収束を論文で主張するのは難しいと思う。漸近的に収束に向かう振る舞いをしているようにみえる場合もあるが「どんな場合でも」というのが困難。わずかのエネルギー差を判別する状況、不安定性の問題、ebarについては変分原理でないこと(augmentation関数計算のためのreference energy)などがあり困難。またあえて言うなら、わずかのエネルギー差を問題にするため、LMX、LMXA,kmxaも十分に大きくとっておく必要があるかも。(gmax/2=12/2=6 a.u.は36Ryに対応し十分に大きいので問題ないだろうが)。
  • 全エネルギー収束チェックは、論文に含めない方がいいとはおもう。もしやるにしても、理論的に収束が保証される状況でのwell-definedな収束チェックに限定する必要がある。たとえば収束チェックをやるにしてもone-shotでHarris-Folknerエネルギーをみたほうがいいかも。Harris-Folknerエネルギーは http://pmt.sakura.ne.jp/wiki/images/Mtopw9j.pdf の式49であり、入力の密度を同じに取ってone-shotで計算するには、OPTIONS_HF=Tとし、lmf --rs=1,0 feなどとする(--rsについてはlmf --input。rstを上書きさせないためのオプション)。このときにはaugmentation関数も同じになり、バンドエネルギーの差が正確に計算できる。すくなくともaugmentation関数を固定(FRZWF=Tで)チェックした方がいいと思われる。ただ、論文で主張するのは困難な気がする。たとえばダイヤモンドにしても、十分にいい計算になってるはずだがガタつきが誤解を招くかも。
  • むしろ「酸素などの分子や固体のcohesive energyが2-4Ryのカットオフでもけっこういける」という点をセールスポイントにしたい。H2Oもしたい。このためには、原子計算をスーパーセルで計算しないといけない。
  • PRB59,1758,Kresseのp.1769のsmall moleculeの最後あたりには、MTサイズを小さくとる場合のいいわけが書いてある(「我々のspherical Besselによる擬ポテンシャルは多少コア領域が重なってもエラーをおこさない」など。理由不明。:そこでの引用文献Ref.10がunpublishedになっている)。PRB81,125117でのFig.4のFは大きな基底での計算だがone-shot計算でのHarris-Folknerエネルギー。

現状の問題点に関するメモ(小谷)

  • cohesive energyのための原子計算もスーパーセルでやる。pwemaxに関して、原子エネルギーと、分子(固体)でのエネルギーが同様に変わるのならいちばんいい(これは原子近傍でのみ高い振動数のAPWが効くような場合には期待できる。多少絶対値のエネルギーで損しても小さいMT半径をとるならそれなりに期待できるかもしれない。「All-electronに非常に近い物理量を与え得るが、絶対値は多少ずれてしまう。(PAWなどと同様)All-electronとは多少のずれはいたしかたない」のスタンスでいくのがいいと思われる。「MT半径を小さくとってもそんなにひどくはない」というのはセールスポイントであるともいえる。それで、MT半径はcoreやvalenceのspiloutも考慮すべきだが、計算安定性を重視してすこし小さめにとってもよいかもしれない。
  • MTOをいれた強みを除いては、PAWや小口先生のLAPWなど「solar-Williamsの方法」を採用した方法とlmfの基本はおよそ同じと思われる(PWMODE=12に対応)。ただ弱点は、smooth partがover completeなこと;これが以下の計算の異常終了の「MT内のsmooth partの不定性」と関係すると思われる。
  • gmaxを小さくとりたいという要求もある(電荷のための空間メッシュを荒くとって高速化する)。しかしRをある程度小さくとる以上、難しいかもしれない。いまのgmax=12は、おおむねそんなに悪くないように見える。マークさんはgmax=7程度で計算していたがこれは固体でかなりMTが大きい場合。
  • 「Rの決め方」に対するコメント1。酸素などの深い占有軌道をもつ原子ではEHの値はそれなりには影響する。できれば、やはり原子あたり2-3mRyぐらい以下の収束が望ましい。このときにはmtopara.tmp2を見る必要がある。可能性としては、これのEHと0.1程度ずらしたEH2を使うこと;ただmtopara.tmp2が-2.5(下限値)になっている時には注意がいる;smooth Hankelのfittingが実質的に失敗していて、EHを決めきれていない可能性がある。
  • 「Rの決め方」に対するコメント2。不安定性は二種あるようにみえる。以下の計算の異常終了の「MT内のsmooth partの不定性」を参照。
  • 「Rの決め方」に対するコメント3。遷移金属のRは、全エネルギーに結構影響するのである程度は大きくとりたい(木野さんの言うvalence spilit out)。磁性にも関係する場合は微小なエネルギー差を問題にするし、APWでは波数が相当に必要なので、大きく取るのが有利。あまり大きくとりすぎるもまずいので、2.2ぐらい(?)に固定するのがよさそうに思える(期待)。直接に関係してるかどうか不明だが、PRB59,1758,Kresseのp.1770のpseudopotential genarationではFeなどで、core radiusを2.2にしてる。酸素などは結構小さくてもいけるようだ。アルカリ金属、アルカリ土類などが大きすぎる。このあたりのsp系はcoreのspiloutだけに注意してかなり小さくしてもあまり全エネルギーには効かないだろう(これらにはsubs/writebasis.Fに書いてあるようにpセミコアをPZをいれるのがデフォルトにしてある)。それで、およそ最大値を2.2あたりに固定したらpwemax=4程度までの計算ができるようにも思える(期待してるだけ)。spherical bessel関数のノードの位置でおよそはきまってるのかもしれない。小谷のctrlgen.pyでは最大値3.0になるようにしているが、これはLaSrO3でSrが大きすぎて計算できなかったので拘束を入れたため。
  • 「Rの決め方」に対するコメント4。実際の物質で「できるだけ重なりがないようにしておく」というニーズもあるので、遷移金属やf電子系は[covalent radius http://en.wikipedia.org/wiki/Covalent_radius]かそれをすこし小さくするぐらいのサイズでいいのではないか?(遷移金属などでは小さくはしすぎない)。アイオニックなもので重ならないようにするには、sp電子系(とくに1、2族)を小さく取るのでできないか?最適のEHをさがすには、lmfaの与える値がある程度の指標になる(lmfaを改良できます)。
  • 計算の異常終了:以下は、小谷の認識してる分類です。厳密に区別されるわけではない。
  1. 平面波基底の間での線形独立性の喪失:pwemaxを大きくしてくことを考える。interstitial regionで、新たに加わった平面波基底が,既存の平面波基底で展開できるとき、計算が異常終了する場合。この場合NULLな基底(ハミルトニアンにとっても、重なり行列にとっても固有値がほぼゼロになるものとなる)を含むことになり、計算が異常終了する場合。この場合にはHAN_STABILIZE=1e-8(1e-6から1e-10程度)を加えるととりあえず収束させることができる場合がある(ただしそのような計算を応用に使うのは、コントロールできない要素が入るので好ましくない)。ovlmatをみてわかるようにNULLな基底が存在している。これを引き抜くしかない。あるいは収束チェックにのみ注意深くHAN_STABILIZEを用いること。これを固定しておくなら、pwemaxを大きくしていく極限も期待できる。HAN_STABILIZEなしではやたら深い固有値が表示されたりする(これは実際にはNULLベクトルであり、不確定な固有値をもつもの)。これは通常のLAPWでも起こり得る不安定性である。
  2. Hankel関数が、平面波基底で展開できてしまうことによる線形独立性の喪失:interstitial regionでのHankel関数が、与えたpwemaxに対する平面波基底で展開できてしまう場合もありえる。このとき固体においてはpwemaxが2ー4程度でも、MT半径が大きすぎると、STABILIZEを相当に大きくしないと収束しない(あるいは大きくしても収束せず振動する)。このとき「MT内でsmooth part >>true part (charges:をみる)」というような異常がみられる。この異常はMTを小さくとって解消することを考える。木野氏のRの決め方を参照。この異常はLAPWではまずおこらないだろう(LAPWのsmooth partは平面波で完全系であり、おそらく「MT内でsmooth part<true part」が保てる)。
  3. 電荷の正定値性:現状ではFull meshでないために電荷が正定値でなくなり補正を加えざるをえず破綻することもある。ただこれは、上の二つの理由とも関係しており、直接の要因と言えないかもしれない。収束途上でおこるだけなら最終結果に影響しないが、このために不安定性が発生し落ちる場合もあると思われる。
  • 実用的には、pwemax=2-4ぐらいを用いたい。PAWなどと比べて一桁以上(コードの実装は別にして)本質的に高速で安定な手法でないと意味がないだろう。かつ必要とあれば収束チェックもできないといけない。ただ収束チェックが「自動的に」できなくてもそこはしかたないと思う。
  1. 将来的には、原子種によって「MT半径固定、(augmentation関数も固定)、EH、EH2も固定」(どうしても重なりそうなときだけすこしMTをスケールする)ぐらいを試したいと思う。sp電子系だとMT半径は、そんなにエネルギー損することなく(spiloutもなく)小さくとれるだろう。違う構造でのエネルギー比較などもしやすい。バージョン番号やセッティング番号みたいなもので計算を管理することになる。(「収束やトラブルのチェック」と「ユーザー側」を分離する。再現性を確保する。)
  2. 上述の方法だと固定基底の計算法になる。この場合、基底の係数を動かして最小化するだけの問題になる。iterationにおいて重なり積分などは固定になり、安定し高速化もしやすくなることが期待できる。またforceとエネルギー微分が正確に一致するのでフォノン計算なども安定しうる。GW(RPA)での全エネルギー計算にも有利かもしれない(現状でも、LDAからのone-shotで計算できる。しかし、異方性とか上のバンドとかの処理で問題はある)。
  • NaClをpwemax=4程度で調べてみた。電卓のようにしてくりかえし計算してたのでかなりいい加減でデータなしのメモ。
    • NaのR半径、ClのR半径に関しては、どちらの半径にしてもある程度(2.8ぐらい)以上のMT半径だとおちる。HAN_STABILIZE=1e-8を入れるとある程度安定した。
    • Naを2.9ぐらいにすると、HAN_STABILIZEだけでは安定しなくなる。1e-3にしても不安定なまま収束せずに回りつづける。
    • 「Na,ClともR=2.2にして、FRZWF=T(augmentationのphi,phidotを原子計算lmfaで固定)」でやると、self-consistent計算よりエネルギーが下がった。これは,phi,phidotをPDOS重心にとるという処方箋がかならずしも最小化になってないということ。FRZWF=Tは原子ごとにあたえることもできる。
  • CrやMnの等核ダイマーは非常に小さい。R=1.5a.u程度。このようなMT半径では固体計算においては非効率となるがそれでもこのような値に固定するのがよいのか?

Kinoのlmf TODO list

(高),(低)は優先度。

  • (高)R小、railは大きいEHの方がいいか?EHより大がmin energy?
  • (高)\の場合EH2除いてEH依存性
  • (高)phi,phidotの図示。spillout具合を目で見るため。
  • spilloutの値をthretholdにしてRの下限を決定
  • (高)frequency, bulk modulusとの比較
  • (低)sample:interface
  • (低)sample:surface
  • (低) bulk,moleculeの自動選択
  • 小さいoverlap固有値を与える波動関数の図示
  • core spilloutが大きければ計算を止めるoptionの作成。
  • local orbitalを複数channel入れる。

終了したもの

  • (高)scheme: EH2=EH-0.1固定 --- done, Rを小さくすれば収束させられる。
  • (高)EHが大きい時のgmax依存性check --- done,EHが大きくても収束していた。
  • (高)scheme: 2Bから左でEH=-0.1[2:4],EH2=-2.0[all] ---done, この効果はmarginal, Rを変える効果の方が遥かに大きい。
  • phi,phidot固定: NaClで実行。phiを止めることによりenergyが一意的にに上がったり下がったりということはない。
  • dimerをいろいろやってみてR mol parameterの作成

LiFePO4について

  • LiFePO4においては、以下で説明。これをもう少しテスト。--TakaoKotani 11:46, 4 January 2011 (JST)
  • 「LiのMT半径が大きすぎる」。

やはりLiFePO4の場合、LiのMT半径が大きすぎるのかも?covalent radiusの表の値を見ても、その値より大きい。Liイオン電池用ということだから、Liイオン以外がフレームをつくるような構造か?MT半径2.657*.529177=1.405は、表の1.33より大きい。--- Li:Rだけを変えても不安定--Kino 16:51, 4 January 2011 (JST)

  • LiFePO4ではCharges:を見ていくとLiサイトがはげしく荒れる.LiFeのSを抜くだけではダメのようだ。

とにかく落ちるときは、ことごとく、MT内のスムーズ電荷がやたら大きくなっておちる。LiのSを抜くと、 Feで荒れて落ちる。どちらも小さくしないといけないようだ。Liサイト,Feサイトでs-channelを抜いたのちR=2.2にした(RSMHを小さくするのは忘れた)ら28回で収束した。(EH=-0.1,EH2=-2を用いてみた)。ehf=-15388.511392 ehk=-15388.511198。Liのみ小さくした場合にはFeサイトのCharges:が変になっておちた。pwemax=2.nspin=1,バンドギャップ0.33eV.--TakaoKotani 11:46, 4 January 2011 (JST)

  • なぜ半径も小さくする必要があるのか?原因を特定する必要があるが、まだわからない。
「異常なenvelope関数」=「envelope関数の線形結合で、センターのMT内でのみ大きな値をもつようなもの」が存在しうるこ可能性。ハミルトニアンに対しても重なり行列にたいしてもNULLであり、どんな固有関数にでもくっついてしまう。

ということでよいのか?いずれにせよ、もう少し考えてテストする必要あり。

  • Q=で設定されるLiサイトの初期条件はlmfaのインプリメンテーションのために変になってる.その為、最初にみだれて収束が遅くなりうる。

FRZWF=Tでもphi,phidotをつくりなおそうとするのはなぜ?

たしかにこの必要はない。それで、対症療法ですが、lmfp.FのL583のlpnu = 1 を0にしてください。これでいちおういけるとおもいます。 ただ、計算が落ちたりするのは,更新された電子密度がだいぶと変になってる可能性もある。 charges:とか見ればある程度わかりますが。 ーーー 解説: コードを見直した。たぶん以下であってると思う。 あくまで出発点はrst,atmに入っている電子密度です。くわえて、rstにはpnu(対数微分) が入ってます。(どうもatmには入ってないor使われてないようです)。 FRZWF=T だと、スタート時のポテンシャルとpnuが保持されるので、phi,phidotが保持されるということになってます(若干非効率だがphi,phidotを作る(まったく同じ) 計算を繰り返している)。 *atmから出発する時は、これに含まれる電子密度からポテンシャルをつくり、  lmfのどこかで与えられてるデフォルトのpnuを用いて計算する(defpq.F). *rstから出発する時も同様だが,pnuはrstに含まれるものが使われる。 --TakaoKotani 16:07, 7 January 2011 (JST)


止まる例,96d83f8

Ar atom

Q=2 4 2

 e=    -50.00000  q=**********   nr/nre/kc= 489 489 163   de=  0.00D+00
 RSEQ : nit gt 80 and bad nodes for l=2.  Sought 0 but found44.  e=-5.0000D+01
 Exit -1 RSEQ: bad nodes

止まる例

gaussian/PBEの基底状態に対して実行しています。 全原子やってみないとだめそう。

O2 tutorialに従い

      TETRA=0 
      N=-1
      W=0.001

としている。

計算はすべてPBEを用いた。(つもりだが時々修正しきれていなくてLSDAになっている。)

===注意=== 必読。W=0.01がよいかも。

  1. LDAで収束させてから、GGAで収束させるのはまずい。---現状では、lmfaがGGAもしくはLDAでコアを作る。これがずっとrstファイルに引き継がれる。
  2. MT半径をかえたとき。その時にはlmfaから再度実行しないと全エネルギーがずれる。たとえば、Feで、R=1.9006からR=1.910に変えるだけでも200mRyぐらいずれた。理由はlmfaでfrozenコアの形(interpolation)が固定されてしまうために非整合をおこしていることによる。将来、安全側にメッセージだすようするか自動でlmfaを呼びなおすように直す。
  3. MMOMの設定ミス。例えば,Fe2で、MMOM=0 3 0 0などとして、pにモーメントを入れていた。こういうことをすると対応しきれず、最初のiterationでこける。(eigenvalue=に非常に深い固有値が表示される)。
  4. 対称性。場合によってはSYMGRP eにする必要があるかも。ただ、TETRA=0,NW=-1の設定では http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html にあるように有限温度になっているので、十分に高温であれば(W=が大きければ)対称でもいいかもしれない。(現状では高温にするとき十分に大きいBZ_NEVMXを用いないといけないが(NEVMAXでない。lmf --input),少しバグがあり、修正必要.またNEVMXを自動設定するようにしないといけない。)。--TakaoKotani 11:36, 26 November 2010 (JST)
  5. 温度設定。どうもデフォルトの
TETRA=0
N=-1
のときの温度設定、W=0.001(表示されるようにT=157K)では、低すぎて収束させにくいようだ。とくに磁性のあるときなどの安定性に問題あり。FSMOMを使うか、W=0.01にとるのがベター。分子の場合、デフォルトをW=0.01に変えたほうがいいかも。ctrlgen.pyに説明を加えた。


収束メモ--Kn 01:07, 30 November 2010 (JST)

  1. mode12から mixm.*, moms.* を消して mode11へ (MTO -> PMT)
  2. EHのみから EH+EH2へ。
  3. はじめは小さいcutoffでやり、大きいcutoffへ。
  4. smearing energy(W=)を大きく。

OPTやMDをやるとき途中でscriptを切り替えるというわけには行かないので、安定なアルゴリズムをbuilt-inしておくべき。

magneticな場合。

  1. momentの初期値を小さくしておく


Al2

PBEの基底状態はtriplet。原子位置はPBEの平衡距離に取っている。

% const len_au=0.529177
% const disexp=1.246237*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0 nspin=2                                            
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Al POS= 0 0 {dd}*.5
        ATOM=Al POS= 0 0 -{dd}*.5

ctrl.al2はnk=1, MMOM= 0 2 0 0 と修正する。

lmfの計算結果

 RSEQ (warning) eval for l=1 node=0 did not converge: de=1.39D+08 e=-7.7961D+01
 PHIDX  Z=13  l=1  nod=0  bc=1.37e7 Inf  e(bc)=-77.96  e(bc)-e=0.1
 Exit -1 phidx failed to converge

可能。 MMOM= 0 2 0 0がまずいようだ。lmfa al2|grep confすると、初期電子配置がみれるが、pチャンネルの電子数がマイナスになっている。このためiterationの一度目で異常に深い固有値が出る。MMOM= 0 1 0 0でないとまずい。それで見過ごさないように--skip_qvalcheckを入れないときには、lmfaがExit -1するよう改良した。--TakaoKotani 15:06, 26 November 2010 (JST)

Fe2

% const len_au=0.529177
% const disexp=1.005793*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Fe POS= 0 0 {dd}*.5
        ATOM=Fe POS= 0 0 -{dd}*.5

ctrl.fe2はMMOM=0 0 4 0と修正。

計算結果

 evl(nev=12)=0.456 but ef0=0.464 ... restart with larger efmax or nevmx
 Exit -1 bndfp


可能。 NEVMXのオプションはなくして自動化。TETRA=0でNW=-1だとフェルミ分布モードで計算してるため(温度はW=できまる)、バンド数を多めにとらないといけないが自動化した。--TakaoKotani 23:06, 28 November 2010 (JST)

Ge2

% const len_au=0.529177
% const disexp= 1.211273*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Ge POS= 0 0 {dd}*.5
        ATOM=Ge POS= 0 0 -{dd}*.5

ctrl.ge2はMMOM=0 2 0 0と修正。

 Fit local orbitals to sm hankels, species Ge, rmt=2.289
  l  type    Pnu      Eval        K.E.       Rsm       Eh      Q(r>rmt)    Fit K.E.
 Exit -1 : mtchre : failed to match phi'/phi=-2.052 to envelope, l=2

可能。 3dに対応するMTO(local orbital)でこけてる。デフォルトRS3を0.5にもどした。すべての場合にいけるかどうか不明。 またgmaxに対する収束性が悪くなっている可能性もある---RS3が小さいとかなりきついハンケル関数をつかうことになってしまい実空間メッシュのgmaxが多数必要になりうる(まだチェックできてない)。

Li2

% const len_au=0.529177
% const disexp= 1.367254*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Li POS= 0 0 {dd}*.5
        ATOM=Li POS= 0 0 -{dd}*.5

lmfaに失敗する。

 ###### Why is Exit caused? #######################
 # This version can not treat two valence channels 
 # per l (Q for a l-channl is zero if the l is with PZ).
 # This causes a problem typically in Li; then we 
 # can not treat PZ=1.9 and P=2.2 as valence.
 # To avoid this, use Q=0,1 together. This trick supply an 
 # electron to p channel; this trick works fine.

 Exit -1 freeat, l=0:  nonzero charge Q=1 not allowed  with sc PZ=1.9

可能。 ctrlgen.py(正確にはそのなかのlmfa-freeat.Fかどこか)が書き込むデフォルトがまずい。 PZ=11.9がコメントされてない場合(これが標準)にはP=2.6 Q=0,1.000のコメントをはずして生かす必要がある。 デフォルトでPZを用いるように修正した。

Mn2

% const len_au=0.529177
% const disexp=.824049*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Mn POS= 0 0 {dd}*.5
        ATOM=Mn POS= 0 0 -{dd}*.5

nspin=2, MMOM= 0 0 2 0 に修正。 pwemax=0.5,1, 1.5, 2, 3,4,5,6 とやったが0.5,1,5,6の場合に収束せず。

可能。 W=0.01だと、T=1579Kと表示される。収束する(mmom=1.96,XCFUN=1)が、最後のところでおそい(モーメントがうごいてもエネルギーがあまり変化しないような状況ではかなり手間取る)。その後、温度を下げると、mmom=2で収束する。Fixed Moment法であるFSMOMも併用していいかも。)。--TakaoKotani 23:11, 28 November 2010 (JST)

Mg2

% const len_au=0.529177
% const disexp=1.751421*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0   nspin=1                                          
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Mg POS= 0 0 {dd}*.5
        ATOM=Mg POS= 0 0 -{dd}*.5

lmfa fails.

 mtchr2 mode 1  l = 1  r =  3.000000
 l  it  ir      Rsm         Eh        phi'/phi     target
 1   1  -1    1.500000  -10.000000   -1.946E+00  -2.097E+00
 1   2  -4    1.500000   -0.020000   -5.877E-01  -2.097E+00
 1   3  -5    1.500000   -5.010000   -1.712E+00  -2.097E+00
 exit mtchr2 info =-1
 Exit -1 : ftfalo : failed to match log der=-2.097 to envelope, l=1

可能。 第一に、RS3(local orbitalに関するMTOのRSMH)の設定がまずかった。 SPEC_ATOM_RS3のデフォルト値を、def_r8=rmt(j)/2d0からdef_r8=0.5d0にもどした。 (なぜそもそもrmt(j)/2d0にしたかを失念した。メモが不明。原子によってはRS3を変えないといけないか?)。


第二に、また、深いセミコア(-3Ryにある)をデフォルトでMTOとして用いており(多くの場合に不要だと思うが)、この radial関数のエネルギー微分に問題があった(phidx.F-rseq.Fにおいて正しくない計算がおこっていたので、 warningをだしながらまわっていた)。それを、ちょっと対症療法でfixした(数値微分をとる際のエネルギー差を小さくした)。 これをいちいちチェックするのはめんどうなので、RSEQ_ERRORファイルを 吐くようにした(各iterationの初めに消される。終了時に存在してれば、Exit -1を吐いて終了)。 ただし、そもそもlocal orbitalでphidxで作られるphidotは使われていないような気がする;実際答えはMg2のときにはまったく変わらなかった。--->要確認。

Co2

% const len_au=0.529177
% const disexp=1.216687*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0      nspin=2
# mmom= 0 0 4 0
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Co POS= 0 0 {dd}*.5
        ATOM=Co POS= 0 0 -{dd}*.5

lmf fails

 evl(nev=13)=-0.149 but ef0=-0.133 ... restart with larger efmax or nevmx
 Exit -1 bndfp

可能。

Sc2

% const len_au=0.529177
% const disexp= 1.156531*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0      nspin=2
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Sc POS= 0 0 {dd}*.5
        ATOM=Sc POS= 0 0 -{dd}*.5

MMOM=0 0 2 0 へ変更。

lmf fails

 e=     -7.92523  q=   0.00000   nr/nre/kc= 365  35  11   de=  2.60D+09
 RSEQ (warning) eval for l=2 node=0 did not converge: de=2.60D+09 e=-7.9252D+00
 PHIDX  Z=21  l=2  nod=0  bc=17.04 Inf  e(bc)=-7.925  e(bc)-e=0.1
 Exit -1 phidx failed to converge

可能。 最初のeigenvalue=から異常に深い固有値が出る。これはMMOMの設定が適切でないため(Fe2のときは、小谷自身、MMOMでpチャンネルをpolalizeさせてしまって同様の症状でなやんだ。MTO-onlyにするとそういう状況からでも収束させれるので余計に解決が遅れた)。ただし,MMOM=0 0 1 0からだと、30+5回収束にかかった。lmfa sc2|grep confで初期原子配置を確認する。ちゃんと安全側のメッセージを出すように改良した。原子計算もsupercellで行う必要があるが「frozen coreの分極をどう固定するか?」に注意がいる。

Zn2

% const len_au=0.529177
% const disexp= 1.604071*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Zn POS= 0 0 {dd}*.5
        ATOM=Zn POS= 0 0 -{dd}*.5

pwemax=6に変更

 80  0   1   1 413 283   -188.7429532   -188.7429532   -188.7429532   1.8952D+02
 e=   -188.74295  q=**********   nr/nre/kc= 417 413 283   de=  1.90D+02
 RSEQ (warning) eval for l=0 node=1 did not converge: de=1.90D+02 e=-1.8874D+02
 PHIDX  Z=30  l=0  nod=1  bc=-3.524e14 -Inf  e(bc)=-188.7  e(bc)-e=0.1
 Exit -1 phidx failed to converge
SPEC
 ATOM= Zn Z= 30 R= 3.000
  RSMH=  1.500  1.500  1.500  1.500  EH=  -0.6  -0.6  -0.6  -0.6
  RSMH2=  1.500  1.500  1.500  1.500  EH2=  -1.4  -1.4  -1.4  -1.4
     PZ=0,0,13.9 P=0,0,4.2     KMXA={kmxa}  LMX=3 LMXA=4
     MMOM=0 0 0 0

もしくは

SPEC
 ATOM= Zn Z= 30 R= 3.000
  RSMH=  1.500  1.500  1.500  1.500  EH=  -0.8  -0.8  -0.8  -0.8
  RSMH2=  1.500  1.500  1.500  1.500  EH2=  -1.2  -1.2  -1.2  -1.2
     PZ=0,0,13.9 P=0,0,4.2     KMXA={kmxa}  LMX=3 LMXA=4
     MMOM=0 0 0 0

とするとlmfでNaN。

Kohn-Sham energy:
 sumtv= ***************  sumtc=      6863.363997   ekin=****************
 rhoep=             NaN   utot=              NaN   ehks=             NaN
 mixrho: sum smrho  init = 0.118029D+05 0.561669D-25 0.118029D+05       0
 mixrho: sum smrnew new  = 0.101588+200 0.152258+168 0.101588+200       0

 mixing: mode=A  nmix=2  it 9 of 9  beta=.5
 mixrho: dqsum rmsuns=  0.24080+194 NaN          NaN
 mixrealsmooth= T
 wgtsmooth=  1.539600717839002E-003
 mixrho:  sought 2 iter from file mixm; read 3.  RMS DQ=Infe0  last it=1.36e0
 charges:       old           new         screened      rms diff       lin mix
 smooth     188.799813****************************      Infinity**************
 site    1    9.768879****************************      Infinity**************
 site    2    9.768879****************************      Infinity**************

可能。 W=0.001でやってた---magneticでない場合はこれでもいいのかも。 pwemax=2だとOK。pwemax=6はiterationの2回目でおかしくなる。対処法は、 pwemax=2を収束させた後,pwemax=6とする。やってみるとエネルギーゲインがほとんどない。

pwemax=2 ndimh=415基底  c ehf=-7173.3400728
pwemax=6 ndimh=1817基底 c ehf=-7173.340106 

(pwemax=2から5回で収束,kmxaを8にした.経験的にpwemax+2ぐらいだと安全)。

Kr2

% const len_au=0.529177
% const disexp=1.878011*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0     nspin=1
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Kr POS= 0 0 {dd}*.5
        ATOM=Kr POS= 0 0 -{dd}*.5

lmf fails

 e=   -204.35149  q=**********   nr/nre/kc= 431 429 299   de=  9.19D+01
 RSEQ (warning) eval for l=0 node=1 did not converge: de=9.19D+01 e=-2.0435D+02
 PHIDX  Z=36  l=0  nod=1  bc=-8.457e14 -Inf  e(bc)=-204.4  e(bc)-e=0.1
 Exit -1 phidx failed to convergere

Ne2

% const len_au=0.529177
% const disexp=1.275151*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0     nspin=1
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Ne POS= 0 0 {dd}*.5
        ATOM=Ne POS= 0 0 -{dd}*.5

lmf fails

 e=    -42.74043  q=   0.00000   nr/nre/kc= 315  35  11   de=  4.78D+07
 RSEQ (warning) eval for l=1 node=0 did not converge: de=4.78D+07 e=-4.2740D+01
 PHIDX  Z=10  l=1  nod=0  bc=216700 Inf  e(bc)=-42.74  e(bc)-e=0.1
 Exit -1 phidx failed to converge

すこし問題あり。 最初の初期条件(原子を重ねたもの)が悪くて、dualだと、一回目のiterationがまずい。 一度だけEH=-1のみでEH2抜きでやれば、あとはEH=-1,EH2=-2で収束に向かう。 (iterationの一度目で-44Ryぐらいと異常に深い固有値が出るが、これはiteration2回目以後、浅くなっていく)。 他の場合にしても「一回目(or何回か)はEHのみでまわす。そのあとEH2もいれてまわす」 というのが一般的に有効かもしれない、がテストしていない。


pwemax<=1.5はEH->EH+EH2で収束。pwemaxを増やす方法でpwemax>=2が収束するかチェック。--Kino 13:20, 30 November 2010 (JST)

Ni2

% const len_au=0.529177
% const disexp=1.063304*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0      nspin=3
###  mmom= 0 0 2 0
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Ni POS= 0 0 {dd}*.5
        ATOM=Ni POS= 0 0 -{dd}*.5

lmfa fails

 80  0   4   3 485 473     -1.3927650     -1.3927650     -1.3927650
 e=     -1.39276  q=**********   nr/nre/kc= 507 485 473   de=  2.40D+00
 RSEQ : nit gt 80 and bad nodes for l=0.  Sought 3 but found 4.  e=-1.3928D+00
 Exit -1 RSEQ: bad nodes

可能。 nspin=3?

Cr2

% const len_au=0.529177
% const disexp=0.797734*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0      nspin=1
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Cr POS= 0 0 {dd}*.5
        ATOM=Cr POS= 0 0 -{dd}*.5

pwemax=6へ変更。

lmf fails at iteration 1

   26  0.93D-01   27  0.94D-01   28  0.94D-01   29  0.10D+00   30  0.14D+00
   ... skip larger eigenvalues ...
 eigenvalue=
 Exit -1 zhev: zhegv cannot find all evals 111

可能。

Na2

%const len_au=0.529177
%const disexp=1.545153*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0 0 1 0  0 0 1
SITE    ATOM=Na POS= 0 0  -.5*{dd}                          
        ATOM=Na POS= 0 0 .5*{dd}   

lmf fails.

pwemax=2 => fails to converge.

pwemax=3 or 6

 (warning) DOS window (-1,0) reset to (-5.0265,-1.4297)

 evl(nev=10)=-1.928 but ef0=-1.93 ... restart with larger efmax or nevmx
 Exit -1 bndfp


pwemax=4 or 5

 e=  -1102.50403  q=   0.00003   nr/nre/kc= 337 245  81   de=  9.47D+03
 RSEQ (warning) eval for l=0 node=0 did not converge: de=9.47D+03 e=-1.1025D+03
 PHIDX  Z=11  l=0  nod=0  bc=5.874e42 Inf  e(bc)=-1103  e(bc)-e=0.1
 Exit -1 phidx failed to converge

可能。 収束した。たぶんこれもRS3の問題。"local orbitalsの表示の下にrsmが表示されている。

Ti2

% const len_au=0.529177
% const disexp= 0.948721*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Ti POS= 0 0 {dd}*.5
        ATOM=Ti POS= 0 0 -{dd}*.5

MMOM=0 0 2 0.

EH=-0.2, EH2=-1, pwemax=2。 EH=-1.2, EH2-1.4, pwemax=2 などが収束しない。

Ar2

% const len_au=0.529177
% const disexp=1.719724*2/len_au a=10.0/len_au dis=0.0  dd=(disexp+dis)/a
% const nk=1 da=0     nspin=1
STRUC   ALAT={a} DALAT={da} PLAT= 1 0 0  0 1 0   0 0 1
SITE    ATOM=Ar POS= 0 0 {dd}*.5
        ATOM=Ar POS= 0 0 -{dd}*.5

At pwemax=6, lmf fails.

 Make new boundary conditions for phi,phidot..

 site    1   species   1:Ar      
 l  idmod     ql         ebar        pold        ptry        pfree        pnew
 0     0    2.183352  -19.269242    3.940000    3.026431    3.500000    3.500000
 1     0    3.774141   -1.509275    3.920000    3.150348    3.250000    3.250000
 PHIDX  Z=18  l=2  nod=0  bc=1.627e14 2.042e15  e(bc)=-157.7  e(bc)-e=2.814e-9
 PHIDX  Z=18  l=2  nod=0  bc=1.627e14 2.042e15  e(bc)=-157.7  e(bc)-e=1.368e-10
 PHIDX  Z=18  l=2  nod=0  bc=1.627e14 2.042e15  e(bc)=-157.7  e(bc)-e=-1.002e-10
 PHIDX  Z=18  l=2  nod=0  bc=1.627e14 2.042e15  e(bc)=-157.7  e(bc)-e=-1.002e-10
 PHIDX  Z=18  l=2  nod=0  bc=1.627e14 2.042e15  e(bc)=-157.7  e(bc)-e=-1.002e-10
 Exit -1 phidx failed to converge

コードを読むには

  • ecal/lm7Kについては、まずはlm7K/fp/bndfp.Fを軸にみていく必要があります。

これはlm7K/gwd/bndfp.Fと基本的に同じです(fppで違いが出る)。 このルーチンのI/Oをはっきりさせないといけません.これからですが。

  • コードをプリントせよ。geditをつかうとよい。設定で行番号を添付する。端末より紙がよかったりもする。
  • moduleをつかうときは、use m_rdctrl,only: ncutovl,abc1,

などすべてについてonlyをつけたほうがよいかもしれない、明示的になる。


MPI/MPIK対応について

なんとかしないといけない。lmv7.F-lmfp.F-bndfp.F あたりに組み込めばいいと思われるが、検討も必要かと思う。非MPIバージョンとどうcompromizeさせるか?

(MPI,MPIKと2つあるのをどう統一するのかという意味らしい。)

pointer => allocatable,target

pointer => pointer を多用している。 pointerをつなぎ替え続けると memory leakを起こしやすいので poiner => alloctable,target のみにしたほうがいいだろう。


spherical bessel関数の精度

Markさん作のropbes, numerical recipes(NR)のspherical bessel関数を倍精度に変更したもの, PHASEのspherical bessel関数routineの比較。 j_l(x_np)=0となるx_npをbisection法で求めてみた。--Kino 16:04, 2 March 2010 (JST)

exactというのは三角関数を用いた解析解。 difference=x(YYYYYYY)-x(exact), YYYYYYY=その関数名。 eps=bisection法なので最後のiterationでの大きいxと小さいxの差。 E-14はmachine precisionの誤差。

結論

  • NRはもともと単精度で作っているわりに精度がいい。machine precisionの誤差しかない。
  • phaseのはx>1では三角関数を用いた解析解を使用している、つまりexactと呼んでいる関数なのでdifference=0なのは当然。
  • (少なくともこの検査では)ropbesの精度がもっとも悪い。すぐにE-12程度の誤差があり、np=8程度でE-6の大きなずれが生じる。

mshvmt.Fでsum_G j_l(|G|)の類を行っている。j_l(x)が振動する関数なので、足していくと結果は小さな値になる。この部分の精度が気になるならj_lは解析解を用いて、更に|j_l(|G|)|が小さい方から加えるなどという工夫をすべき。

以下詳細。

ropbes
    l   np         x(YYYYYY)                eps(YYYYY)               x(exact)                eps(exact)               difference
    0    1   52   0.3141592653589793D+01   0.8881784197001252D-15   0.3141592653589793D+01   0.8881784197001252D-15   0.0000000000000000D+00
    0    2   51   0.6283185307179458D+01   0.8881784197001252D-15   0.6283185307179586D+01   0.8881784197001252D-15  -0.1278976924368180D-12
    0    3 1001   0.9424777960771150D+01   0.1776356839400250D-14   0.9424777960769379D+01   0.1776356839400250D-14   0.1771027768882050D-11
    0    4 1001   0.1256637061435136D+02   0.1776356839400250D-14   0.1256637061435917D+02   0.1776356839400250D-14  -0.7812417379682302D-11
    0    5 1001   0.1570796326795107D+02   0.1776356839400250D-14   0.1570796326794897D+02   0.1776356839400250D-14   0.2101430141010496D-11
    0    6 1001   0.1884955592125299D+02   0.3552713678800501D-14   0.1884955592153876D+02   0.3552713678800501D-14  -0.2857696301816759D-09
    0    7 1001   0.2199114858412315D+02   0.3552713678800501D-14   0.2199114857512855D+02   0.3552713678800501D-14   0.8994600619871562D-08
    0    8 1001   0.2513274092212302D+02   0.3552713678800501D-14   0.2513274122871834D+02   0.3552713678800501D-14  -0.3065953286807144D-06
    1    1   50   0.4493409457909068D+01   0.8881784197001252D-15   0.4493409457909063D+01   0.8881784197001252D-15   0.4440892098500626D-14
    1    2   52   0.7725251836937529D+01   0.8881784197001252D-15   0.7725251836937707D+01   0.8881784197001252D-15  -0.1776356839400250D-12
    1    3 1001   0.1090412165943218D+02   0.1776356839400250D-14   0.1090412165942890D+02   0.1776356839400250D-14   0.3277378368693462D-11
    1    4 1001   0.1406619391280759D+02   0.1776356839400250D-14   0.1406619391283147D+02   0.1776356839400250D-14  -0.2388667041941517D-10
    1    5 1001   0.1722075527184640D+02   0.3552713678800501D-14   0.1722075527193077D+02   0.3552713678800501D-14  -0.8436984444415430D-10
    1    6 1001   0.2037130295974702D+02   0.3552713678800501D-14   0.2037130295928756D+02   0.3552713678800501D-14   0.4594546965108748D-09
    1    7 1001   0.2351945252136011D+02   0.3552713678800501D-14   0.2351945249868901D+02   0.3552713678800501D-14   0.2267110232878622D-07
    1    8 1001   0.2666605448873015D+02   0.3552713678800501D-14   0.2666605425881267D+02   0.3552713678800501D-14   0.2299174788333858D-06
NR
    0    1   52   0.3141592653589793D+01   0.8881784197001252D-15   0.3141592653589793D+01   0.8881784197001252D-15   0.0000000000000000D+00
    0    2   51   0.6283185307179586D+01   0.8881784197001252D-15   0.6283185307179586D+01   0.8881784197001252D-15   0.0000000000000000D+00
    0    3 1001   0.9424777960769378D+01   0.1776356839400250D-14   0.9424777960769379D+01   0.1776356839400250D-14  -0.1776356839400250D-14
    0    4 1001   0.1256637061435917D+02   0.1776356839400250D-14   0.1256637061435917D+02   0.1776356839400250D-14  -0.1776356839400250D-14
    0    5 1001   0.1570796326794897D+02   0.1776356839400250D-14   0.1570796326794897D+02   0.1776356839400250D-14   0.5329070518200751D-14
    0    6 1001   0.1884955592153875D+02   0.3552713678800501D-14   0.1884955592153876D+02   0.3552713678800501D-14  -0.7105427357601002D-14
    0    7 1001   0.2199114857512854D+02   0.3552713678800501D-14   0.2199114857512855D+02   0.3552713678800501D-14  -0.7105427357601002D-14
    0    8 1001   0.2513274122871833D+02   0.3552713678800501D-14   0.2513274122871834D+02   0.3552713678800501D-14  -0.1065814103640150D-13
    1    1   50   0.4493409457909063D+01   0.8881784197001252D-15   0.4493409457909063D+01   0.8881784197001252D-15   0.0000000000000000D+00
    1    2   52   0.7725251836937708D+01   0.8881784197001252D-15   0.7725251836937707D+01   0.8881784197001252D-15   0.8881784197001252D-15
    1    3 1001   0.1090412165942890D+02   0.1776356839400250D-14   0.1090412165942890D+02   0.1776356839400250D-14  -0.1776356839400250D-14
    1    4 1001   0.1406619391283147D+02   0.1776356839400250D-14   0.1406619391283147D+02   0.1776356839400250D-14  -0.3552713678800501D-14
    1    5 1001   0.1722075527193077D+02   0.3552713678800501D-14   0.1722075527193077D+02   0.3552713678800501D-14   0.7105427357601002D-14
    1    6 1001   0.2037130295928756D+02   0.3552713678800501D-14   0.2037130295928756D+02   0.3552713678800501D-14  -0.3552713678800501D-14
    1    7 1001   0.2351945249868901D+02   0.3552713678800501D-14   0.2351945249868901D+02   0.3552713678800501D-14   0.7105427357601002D-14
    1    8 1001   0.2666605425881268D+02   0.3552713678800501D-14   0.2666605425881267D+02   0.3552713678800501D-14   0.3552713678800501D-14
phase
    0    1   52   0.3141592653589793D+01   0.8881784197001252D-15   0.3141592653589793D+01   0.8881784197001252D-15   0.0000000000000000D+00
    0    2   51   0.6283185307179586D+01   0.8881784197001252D-15   0.6283185307179586D+01   0.8881784197001252D-15   0.0000000000000000D+00
    0    3 1001   0.9424777960769379D+01   0.1776356839400250D-14   0.9424777960769379D+01   0.1776356839400250D-14   0.0000000000000000D+00
    0    4 1001   0.1256637061435917D+02   0.1776356839400250D-14   0.1256637061435917D+02   0.1776356839400250D-14   0.0000000000000000D+00
    0    5 1001   0.1570796326794897D+02   0.1776356839400250D-14   0.1570796326794897D+02   0.1776356839400250D-14   0.0000000000000000D+00
    0    6 1001   0.1884955592153876D+02   0.3552713678800501D-14   0.1884955592153876D+02   0.3552713678800501D-14   0.0000000000000000D+00
    0    7 1001   0.2199114857512855D+02   0.3552713678800501D-14   0.2199114857512855D+02   0.3552713678800501D-14   0.0000000000000000D+00
    0    8 1001   0.2513274122871834D+02   0.3552713678800501D-14   0.2513274122871834D+02   0.3552713678800501D-14   0.0000000000000000D+00
    1    1   50   0.4493409457909063D+01   0.8881784197001252D-15   0.4493409457909063D+01   0.8881784197001252D-15   0.0000000000000000D+00
    1    2   52   0.7725251836937707D+01   0.8881784197001252D-15   0.7725251836937707D+01   0.8881784197001252D-15   0.0000000000000000D+00
    1    3 1001   0.1090412165942890D+02   0.1776356839400250D-14   0.1090412165942890D+02   0.1776356839400250D-14   0.0000000000000000D+00
    1    4 1001   0.1406619391283147D+02   0.1776356839400250D-14   0.1406619391283147D+02   0.1776356839400250D-14   0.0000000000000000D+00
    1    5 1001   0.1722075527193077D+02   0.3552713678800501D-14   0.1722075527193077D+02   0.3552713678800501D-14   0.0000000000000000D+00
    1    6 1001   0.2037130295928756D+02   0.3552713678800501D-14   0.2037130295928756D+02   0.3552713678800501D-14   0.0000000000000000D+00
    1    7 1001   0.2351945249868901D+02   0.3552713678800501D-14   0.2351945249868901D+02   0.3552713678800501D-14   0.0000000000000000D+00
    1    8 1001   0.2666605425881267D+02   0.3552713678800501D-14   0.2666605425881267D+02   0.3552713678800501D-14   0.0000000000000000D+00

GGAのとりこみに関して調べたこと

以下は旧情報。いまはabinitのGGAの関数を用いている。よくできてるともいえないとおもうが、とにかくうごく。 GaAsなどでは,easypbe.Fと非常にきれいに一致した.とにかくWhite-Birdの方法でないとGGAはimplementできない。 低電子密度でこける。またごくわずかにnegativeな密度になるとき(FFTでGの最大値をカットするため)、こけるので、 わずかに1d-8とかの電子密度を全体に加えて正定値の電子密度としている。(negative smrho_w の表示が出る。場合によって1d-3とか加えながらiterationがまわることがある )。

旧情報

lm72のGGAを取り込みたい。たぶん、

  1. mkpot.F内のcall smvxc2で変数slatを渡すように変更。
  2. smvxcm.F、evxcv.F,evxcp.F置き換え。
  3. locpot.F、vxcnsp.F置き換え。
  4. vxc0sp.F,vxcgga.F,vxcgr2.F,vxcloc.F置き換え。
  5. vxcnlm.F,vxcnls.F追加

で済むと思ってます。置き換えてもこちらで開発を進めたところとは干渉しないと思いますので、 マージする必要ないです。

なぜ、それでよいか?とにかくポイントとしては、mkpot.Fの 第0成分と、第1成分のVxcに関係した部分を置き換えるということです。 locpot以下が、第0成分に関係していて、smvxcm、smvxcm2以下が1,2成分に関係している。 それでそれらからコールされるルーチンで置き換えるべきものを置き換えていくということです。

副作用はないと思いますが,もう一度論理的に確認したほうがよいです。

付加的情報。call graphで関連してる部分.完全ではない。

smvxc2はsmvxc.Fに含まれている。 とにかくlocpot.F以下とsmvxc.F以下を置き換えればよいと思われる。 ただ、副作用に注意しないといけないし、改良を加えてる部分との整合性 に注意する必要がある。

  • 以下のグラフは完全でない。機械的に生成すること。

  1. lxcfunがどのxcを選ぶかのスイッチの説明はlm72/rdctrl2.Fのlxcfにある。下のctrlの読み込みを参照。これはlm72でも同じである。
  2. mkpot.F冒頭のio変数同じ(lm7Kではitestというのが最後にあるがテスト用。いまは使ってない)。
  3. smvxc2,smvxcmについては変数slatを渡すように変更されている。smvxc2については"modeの20番台というのが増えている。これが実際使われているのかどうか調べる”。smvxc2のなかでevxcp,evxcvが入ってるがこれがexchange-correlationを生成するコアルーチン。evxcvは本質的ではないがわずかに置き換わっている。

pnuopt

別件ですが、lm72のfreatsでは、pnuoptなどのたぶん最適なpnuを出すアルゴリズムを入れてあるようです。研究します. とにかく、MTが小さいときの破綻の一つはPのとり方にあるようにおもう。LaでMT半径1.5だと、ノード数をMT内で5個とかとれないですから。

lmfの問題点

  • documentを読んでも内部が分からない。
  • programの書き方が特殊すぎてprogramを読もうと変数を追おうにも変数をsubroutineに渡していなくて流れが追えない。

\rightarrowuserになるなら使用可能。

使用して

  • lmfgw-MPIKはすぐに終了してしまう。
RE kotaniの新しいバージョンがおかしいためだとおもう。MPIKブランチチェックできてません.そもそもマークさんのMPI化は大丈夫なんでしょうか?いかせる部分と汚い部分があるかとおもいます。


  • coreのsmooth Hankelのtailを書いてみた。RMTより遠くでつながるようfitできているが,RMTではもとの関数とfitとで飛びがある。RMTで切ってしまって問題ないのだろうか。ここからRMTの下限はこないのだろうか。
RE ありえるとおもいます.気になってます。lmfは従来、packingのよい固体用なのでMT半径を小さくとることがあまり想定されてません。semi-core用にMTOをくわえる(local orbitalのP=14.9など)でいれることはありえますが。
  • scalar relativisticの式は出てきているのは本当にfとgなのだろうか。
RE しらべないとわかりません。
  • spec.orhoc, pot.oqc は同じらしい。
RE しらべないとわかりません。が、ときどき混乱気味です.
全体的なRE 基本的に木野氏のいうとおりだとおもいます。ある時点で一気に書き直す必要があります。とにかく付け焼き刃的な拡張がなされていってるのでぐちゃぐちゃになってきています。まず小谷が知ってることを書き下す必要があるとおもっています。テストシステムを作っておく必要もあります(経験的にはコンソール出力でなく、logファイルのなかの基底の指定、(全)エネルギー関連の値、固有値の比較をすればテストはできます)。これにはとにかく計算をためこんでいってもらえればいいかとおもいます。make checkもきれいとはいいがたい。testgw.pyはどうでしょう?とにかく「出発するctrlファイル。テストのためのコマンド、比較のためのコマンド、比較をOKとする条件(まだあんまり)」を書けば、テストできるようにしたい。lmfの結果の比較にはlogファイルをつかいます。

lmfの原子数の制限

fp/bndfp

     parameter ( nbx=512, n0=10, nppn=12, nab=9)
     integer osig(3,nbx),otau(3,nbx),oppi(3,nbx),oqkkl(3,nbx),
    .  oeqkkl(3,nbx),orhat1(3,nbx)

となっている。ここからcallされている側ではw(orhat1(1,ibas)として使われます。512以上の原子の計算をするにはnbxの値を増やす必要があります。

計算時間

  • O atom 12.0^3 ang cell, 1 processで3min/iteration(itanium2)、対角化の次元=1579。

総計算時間の全体の90%が対角化のためのlapack routineで、全体の50%程はlapackが呼んでいるzgemmが占めています。

  • glysine: 15.0^3 ang cell, 3min/iteration(core i7), 対角化次元=2052
   5061   10.51%   54.63%  augmbl:fp/augmbl.F:274
   4366    9.07%   63.70%  augmbl:fp/augmbl.F:272

で呼んでいるaugqs3, augqp3がそれぞれ総計算量の10%を占めています。

        do  jlm1 = 1, nlma
          do  jlm2 = 1, nlma
            do  k2 = 0, kmax
              do  k1 = 0, kmax
                g(k1,jlm1) = g(k1,jlm1) +
     .                       ppi(k1,k2,jlm1,jlm2,isp)*b(k2,jlm2,i2)
              enddo
            enddo
          enddo
        enddo

こんな感じのloop

c$omp parallel do private(...) schedule(static,1)
      do  i2 = 1, ndimh

とすると4coreで3倍の速度になる。

-g -traceback, per iterationでmpi(mpikでない)による速度upを見ると

1core :8min.
4core 4mpi : 11min.

かえってのろい。


110原子のcrystral,~10x20x30 ang^3、1 processで16min/iteration(itanium2)で対角化の次元=1317。 時間がかかるsubroutineは

times       %       %total     
    176967   15.57%   15.57%  augmbl:fp/augmbl.F:274
    145358   12.79%   28.35%  augmbl:fp/augmbl.F:272
     84462    7.43%   35.78%  hsibl:fp/hsibl.F:312
     72549    6.38%   42.16%  rsibl1:fp/rsibl.F:500
     71198    6.26%   48.42%  pvdf1:fp/dfrce.F:684
     62440    5.49%   53.92%  flocbl:fp/flocbl.F:92
     61241    5.39%   59.30%  rsibl1:fp/rsibl.F:504
     49933    4.39%   63.69%  flocbl:fp/flocbl.F:128
     49559    4.36%   68.05%  __intel_new_memset:??:0
     41228    3.63%   71.68%  rsibl1:fp/rsibl.F:508
     31855    2.80%   74.48%  rlocbl:fp/rlocbl.F:317
     29066    2.56%   77.04%  suphas:fp/suphas.F:73
     23827    2.10%   79.13%  hsibl3:fp/hsibl.F:607
     19352    1.70%   80.84%  N4_M4_Kgas_1:__tmp_ctrsm_right_ker.c:0
     14953    1.32%   82.15%  mkl_blas_zgemm_pst:??:0
     11173    0.98%   83.13%  hsibl3:fp/hsibl.F:601
      9987    0.88%   84.01%  hsibl3:fp/hsibl.F:606

です。

-g -tracebackをつけると per iterationで

1core : 191min. 
4core 4mpi : 89min.

mpi(mpikではない)はあまり効果的でない。上の値とえらく速度が違うが・・・。

core i7, itanium2のmklはSMP並列化をしても速くなりません。 mklをSMP並列にしても高速化されないのでnon-MPI版は1 processで使えば十分です。

スピンのpwemax依存性

現在のmixingでは分子は収束しないものが多い。スピンが絡む遷移金属原子はほとんど収束しない。 1k点(nk=1)で収束しなくてもnk>=2で収束する場合がある。

CrH LSDA magnetic moment=5, re=1.632ang

llmf.crhp30:c mmom=-3.0000001 ehf=-2093.409852 ehk=-2093.4098522
llmf.crhp50.nk2:c mmom=5.0000013 ehf=-2095.9713457 ehk=-2095.9713483
llmf.crhp70:c mmom=5.0000005 ehf=-2097.4174172 ehk=-2097.4174199
llmf.crhp70.nk2:c mmom=5.0000025 ehf=-2097.4249611 ehk=-2097.4251632

where nk2 means nk=2, others nk=1. p?? means pwemax=??.−> pwemax>=5でないととspin momentが収束しない。

MnH LSDA magnetic moment=4, re=1.568ang

llmf.mnhp20:c mmom=4 ehf=-2313.8013886 ehk=-2313.8014416
llmf.mnhp20.nk2:c mmom=4 ehf=-2313.8000842 ehk=-2313.8001614
llmf.mnhp30:c mmom=4 ehf=-2313.8081504 ehk=-2313.808201
llmf.mnhp30.nk2:c mmom=4 ehf=-2313.8074394 ehk=-2313.8074905
llmf.mnhp50:c mmom=4 ehf=-2313.8121663 ehk=-2313.8122168
llmf.mnhp50.nk2:c mmom=4 ehf=-2313.8116622 ehk=-2313.811719
llmf.mnhp70:c mmom=4 ehf=-2313.8137359 ehk=-2313.8137895

これはspin momentに関してpwemaxの収束が速い。

以上はctrlgen.pyで作成。SYMGRP Eに変更。

reply1:モーメントが実数になってないとすると、順位に微妙な縮退がおこって妙なことになってる可能性があります。バンドギャップが出てない。upとdownの準位がへんな具合に縮退してる。酸素原子のparaでも安定して収束しないですし(結局、なんかしら妙なシーソーのような現象が起こる)、なんかバグがありそうなように思います。ctrls.srhもおいてくれたらいいです。パッキングをよくしていったら安定するかも。とにかく原子間の距離をあけてスーパーセルを小さくするとあんていはするようなんですが。。。ってわけで、現状では「固体用(表面もむずかしいだろう)」のコードだと思わないとしかたない。まずは酸素原子を調べてますが、mixingの問題というよりも、どうもhartreeポテンシャルの生成に問題がありそうな気がする(まだ酸素原子のパラ(異方性あり)を調べてるところです)。 酸素原子の件は解決。

CrH, ctrl.crh

%const  au=0.529177                       #<--   These % sections are to define variables.   
%const  a1= 10.0  da=0
% show vars                               
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT  (in  a.u.)  
                                              # DALAT gives "shift" of ALAT  (in  a.u.)  
SITE                                            
  ATOM=H POS=0 0 0
  ATOM=Cr POS=0 0 1.6318

ファイルの構造

例えば、hbe.dの5つめ、BBVECの2行目。BZDATAの1行目にnqbzが書かれています。nqbzを読みたいと思った場合、どれがもっとも正しいはずなのか分からない。(上の例ではBZDATAのnqbzから他のファイルのnqbzを書き出している。)

データの内部構造が汚い、ということに加えてファイルの構造も汚い。データを一元化した方が綺麗になる。 たとえば、*.dは大体パラメタを書いてると思いますが、

  1. parameterを書かせるファイルが多すぎ。
  2. どこにどのデータが書いてあるのか把握できない。
RE たしかにその通りだと思います。データフローが読みにくくなってます。対症療法でやってきたためです。
ちょっと言い訳:データセットには「配列の次元のデータやインデックスをしていするもの」と「(その上での)データ」がふくまれます。で、「(その上での)データ」を渡したいときに、どうしても前者もいっしょにわたしてやる必要があります。また、似たような変数が違う意味で使われたりで、将来的にどうなっていくかが予想できない(「将来どうデータ構造が変更されるかが予想できないこと」を前提にしないといけないですね)。


たとえば@GWDATAというファイルをつくって 内部が

nqbz 27
nqibz 8
<qbz>
1 0 0
0 1 0
...
</qbz>

と書いておく。

nqbzは必ず@GWDATAから読み込む。@GWDATAは全ての*.dとかBZDATAとかを含んでいるので*.dという無駄なファイルを途中でつくらなくてもいい。

nqbzだけ読みたかったら

call keyvalue('GWDATA',nqbz)

BZDATA相当を読みたかったら

call read_bzdata() !---この内部でcall keyvalue('@GWDATA',nqbz); call keyvalue('@GWDATA',qbz,...);...

として読めばいい。どうよ。--Kn 04:42, 9 October 2009 (JST)

qg4gw mode4

SYMLを読むmode3に追加して、mode4としてGWinput <QPNT>のq点を読み込むコードを書きました。--Kino 17:48, 25 September 2009 (JST)

analyzer

ちょっと変です。デザインしなおしたほうがいいかも。



parallel make

  • *.o をつくってから*.aを作っている。本来は*.aがなくてもlinkできるはずだが、*.oの順序に混乱があるらしいため。*.aを作ると*.oのcall/callerの順序がどうでもよくなる。
  • makefileには順序に意味があるruleがあり、どこに書いてもいいわけではない。そのためinclude fileの位置などに苦労した。PLATFORM=を使用せずにinclude make.inc.{...}だけでなんとかしたかったができなかった。
  • gnumake特有の機能を用いていて、gmakeでしか動かない。
  • makefile中の${MAKE}は特殊な意味をもつ。makeを再帰的に呼ぶのなら${MAKE}を書き換えてはいけないとgnumakeのmanualに書いてある。

(Kino)

GWメモ

  • hqpe.sc.m.Fで、最終的にsigmファイルに<chi_i|sigma-vxc^LDA|chi_j>を書き出している。

chi_iはMTO基底(これはQSGWをPMTでできるようにする前の話。march9th).

         ev_se_ev = matmul(evec_invt(1:nhq,1:ntq)
    &           ,matmul(se(1:ntq,1:ntq,ip),evec_inv(1:ntq,1:nhq)))
...
         sigmv(:,:,ip) =  2d0*ev_se_ev(1:nhq,1:nhq)
         write(ifse_out) sigmv(:,:,ip)

ifse_outがsigmのファイルハンドル。ipはirreducible q点のindex, nhqは ハミルトニアンH_0の次元で、基底の数であり、バンド数.(MTOのみのときは、qによらない)。 nhq=ntq(たぶん)。seのうち実際に計算されてる次元数はemax_sigm以下のバンドのみ (この次元数はq点による:その前のループのntqxxにその数が入ってる)で、 それ以外はゼロが入ってる。結局sigmには、そういうカットオフの入った 状態での<chi_i|sigma|chi_j>が格納されることになる。

2d0はたぶん単位の換算(GWコードではところどころ、au=2Ry unitになってる)。 seが自己エネルギー.evec_invはevecの逆。evecは波動関数の基底(MTOのブロッホ和) に関する係数. 9May2009 kotani

  • sigmはlmfにおいてbndfp.Fのrdsigmで読み込まれる.

ifisがファイルハンドルで、w(osigm)が実体.これをhamfb3 で、全てのq-meshの上で作る(空間群での回転操作)。これがosfz これが、hft2rsに渡され、実-mesh上でのsigmaになるこれが、ohrs ham ohrsに格納される。

  • そのあと、このham ohrsは、hamblsで読み込まれる。これは

GWドライバーの波動関数器であるsugw.Fからもよばれる。 hamblsは任意のqに対してのsigmaを生成し、ハミルトニアン全体を 組み立てるルーチン。 (osigはsigmaとは関係ない。原子サイトでの積分。lmfのpaper参照)。

hambls L322   call bloch(i,q,nl,plat,mxorb,w(oiprmb),1,nttabs,w(oiaxs),

w(ohrs)からqに対する<chi_q |sigma_q| chi_q>をつくる。このとき、 にsum_q' <chi_q|chi_q'><chi_q'|sigma_R|chi_q'><chi_q|chi_q'> exp(-iqR) で作ってないのが問題。このあたり書きかけ。May2009 kotani


  • FTMESH(GMAX)をおおきくとれない.実質的には問題ないが、収束チェックの

観点からすれば、もうすこし改善がいるだろう。

========================== mail from Mark. ==================
As for the problems with estat in the SrTiO3 case, it is the same
problem you asked about some months ago.  ftmesh=36 is too high, and
the Bessel function maker can't handle such high-energy PW's: the
estat potential requires a sum of numbers which is ~14 orders of
magnitude smaller than individual terms.
  The SrTiO3 test case in fp/test also had ftmesh=36 (the ctrl file is a
little different; it happens that different machines get nearly the
same results in the test case).  To be safe I shifted the test case to
ftmesh=24.  The total energy changes slightly, but the constant
potential shift is a little different and essentially independent of
computer architecture.


  • Ylm definition.
  True real harmonics YL (denined in GW manual) is given by
    cy(ilm) * yl(ilm,r) * r**l.
    cy ---> sylmc
    yl ---> sylm
  Another way is
    call ropyln(ng,g(1:ng,1),g(1:ng,2),g(1:ng,3),ltop,ng,yl,g2), where
  the true harmonics YL = yl \times r^l.


lmfプログラムのおおまかなながれについて

ctrlファイルの読み込み

いまは、lmv7.fがメインルーチンです。それが呼んでるrdctrl2内のreadctrlで、 m_rdctrl内のそれを呼び出して、m_rdctrl内の変数領域に貯めていく。 で、それをつかって、rdctrl2内で初期値設定ということになります。 で、初期値設定された後、実質的なメインルーチンである lmfpがL641で呼び出されることになります。いまも確認しましたが、最初のlmv7の L461で、

C --- Read recrd ---
     call rdctrl2(recrd,recln0,nrecs,prgnam,vrsion,vn,vn2,F,
    .  slabl,osbz,osctrl,osham,ospot,oslat,osmix,osspec,ossite,osstr,
    .  osarry,osmove,ostb,sstrn)

のslabl以下のものはすべて、ここで設定されます。osなんとかもこのrdctrl2のなかでallocateされます。 sstrnもここで書かれます。

注意点:rdctrl2のpass2=Tは不要(新しいものではすでに消した)

また、lmfp L641で呼び出す前に、L571でpass2で呼び出して、

L571 call rdctrl2(recrd,recln0,nrecs,prgnam,vrsion,vn,vn2,T,<---最後のTがpass2を意味する。

などどしてますが、これはASAとTBバージョンでしか意味がないので(STARTセクションの読み込み)、 lmfp.Fのレベルでコメントアウトしてください。 対応して、m_rdctrl2.Fでpass2のifblockはコメントアウトしてください。


実質のメインlmfpでのループ構造

で、そのあとL571のlmfp に突入します。lmfp.Fでは、それの

L350 C ... Re-entry for shear distortion  
        call defrr(oshr,1)
       4 continue

L471 C ---------------- Re-entry point for a new iteration ---------------
         iter = 1
         5 continue

が大きいループの出発点になっています。この5 continueが通常の電子状態計算および、 それを囲む原子位置緩和計算のときの出発点です。(goto 5が二個みつかりますが 最初の方が通常の電子状態計算のループ、後の方が原子位置緩和のループです。

で、それを取り囲むようにしてgoto 4もでてきます(end文の前 L931). この4 continueとgoto 4(L352)で囲まれたループについてはよく知りません。 shearを計算できるようですが、どの程度うまくやってるのか。。。 (まえのWelcovich?の論文みたいにやってるのか? ぼくは使ったことないです)。


bndfpが電子状態計算のself-consistencyループのコア

で、とにかく電子状態計算のself-consistencyループは 5 continue(L473)とgoto 5(L680)の間でなされています。で、それのコアになってるのは、 L593のbndfpです。(注意:10 continue --- goto 10は小さいループです。 dmatuはLDA+Uに関係してて、そうでない場合はnlibu=0です)。

bndfpにおいては何から何が計算されて戻ってくるのかよくわからないという状況になっています。 bndfpがコアルーチンなわけで(実際、fortranで書くべき部分はせいぜいこいつだけだと思う)、 「input outputをsubroutineの頭でおこなってサイドエフェクトなし」にしたいものです。 でもまだどれだけ数の引数になるかとか検討できてないです。mixingもbndfpの外でやるように したほうがいいようにもおもっています。実際dmatuのmixingだけはbndfpのそとのchkdmu でやってたりします。


将来的には、bndfpの外側はpythonで置き換えるようなことも 可能かと思います。実際、置き換えても計算速度はほとんどかわらないでしょう。 それに簡素化できます。ただバイナリレベルでの結合とかいうのでなくて ファイルレベルでの結合などの緩い結合がのぞましいでしょう。あくまでpythonは glueとしてつかいたいと思っています。

データフローを機械的に読めるツールもつくりたい

とにかく各サブルーチンの変数をデータベース的に解析して( 書かれてるか読まれてるかノータッチか)、そのあとそれを組み合わせてやれば データフローのrecursiveな解析ができるはずです。 ま、だいたいは、analyze.pyのような感じでやろうとおもってます。

方法論などのノートも書いていかないといかん。ある程度はすでにあるんで、 断片的にでもだしていってそれを徐々に整備していく形(たぶんwikiで) でやっていきたいと思っています。


テストシステムが汚い

とにかく、開発で一番重要になってくるのはテストシステムです。 マークさんのテストはmake checkでひととおりやってくれる。これは いい方法です。ですが、そのコード自体はfp/test/test.fpでかなり汚い。 テスト自体がサンプルになってる必要もある。

これをpythonみたいなもので置き換えたい。 pythonみたいなglue languageとfortranの組み合わせ(ほかでもやってるところ はあります)が結局生き残っていくんじゃないかと思っています。 なんでもfortranではやれない。まずは緩い結合(binaryレベルの結合は危険だと思う) がいいでしょう。たとえば、ctrlgen.pyみたいなものです (小谷)Apr.11,2009

上にtestp.fpとそれが使うrdcmdの説明を追加--Kino 18:04, 15 July 2010 (JST)

Personal tools