Originally posted by first_yg at 2011-06-03 08:58:38:
能把你的程序發(fā)給我一個嗎，我最不會寫程序

     

      implicit real*8 (a-h,o-z)
      dimension dump(20),num(2),x(100),y(100),z(100),chg(100,100,800)
      dimension chgt(1000),a(3,3)

      write(*,*) 'Choose input file (LOCPOT=1, CHGCAR=2):'
      read(*,*) input
      if (input.eq.1) then
      open(15,file='LOCPOT')
      write(6,*) ' Enter fermi energy: '
      read(*,*) ef
      elseif (input.eq.2) then
      open(15,file='CHGCAR')
      ef=0.0
      elseif ((input.ne.1).and.(input.ne.2)) then
      write(*,*) ' INPUT ERROR, input must equal to 1 or 2 '
      stop
      endif


      write(*,*) 'Spin polarized calculation? (no=1, yes=2):'
      read (*,*) ispin
      if (ispin.eq.1) then
       if (input.eq.1) open(16,file='workfn.dat')
       if (input.eq.2) open(16,file='chgave.dat')
      elseif (ispin.eq.2) then
       if (input.eq.1) then
       open(16,file='workfn-up.dat')
       open(26,file='workfn-dn.dat')
       elseif (input.eq.2) then
       open(16,file='chgave-up.dat')
       open(26,file='chgave-dn.dat')
       endif
      elseif ((ispin.ne.1).and.(ispin.ne.2)) then
      write(*,*) ' INPUT ERROR, ispin must equal to 1 or 2 '
      stop
      endif


      open(19,file='atom.dat')


      read(15,1) dump
      read(15,*) scale
      do i=1,3
      read(15,*) (a(i,j),j=1,3)
      enddo
1    format(20a4)
      aa=sqrt(a(3,1)**2+a(3,2)**2+a(3,3)**2)
      read(15,2) num(1),num(2)
c     write(6,2) num(1),num(2)
2    format(2i4)
      ity=1
      if (num(2) .gt. 0) ity=2
      natm=0
      do i=1,ity
      natm=natm+num(i)
      enddo
      read(15,1) dump
c     write(6,1) dump
      zero=0.0
      do i=1,natm
      read(15,*) x(i),y(i),z(i)
      enddo
      do i=1,natm
      x(i)=z(i)
      do j=i,natm
      if (z(j) .lt. x(i))  then
      x(i)=z(j)
      z(j)=z(i)
      z(i)=x(i)
      endif
      enddo
c     if (a(3,3) .eq. 0) a(3,3)=a(3,2)
      write(19,30) z(i)*aa*scale,zero
30   format(2f10.5)
      enddo
      nat=natm/2+1
      do i=1,nat
c     write(6,4) i,z(i)
   4  format(i5,f10.5)
      enddo
      read(15,1) dump

      do 600 is = 1,ispin
      read(15,*) nx,ny,nz
c     write(6,*) nx,ny,nz
      nn=nx*ny
      read(15,*) (((chg(j,k,iz),j=1,nx),k=1,ny),iz=1,nz)
c     write(6,5) (((chg(j,k,iz),j=1,nx),k=1,ny),iz=1,nz)
   5  format(5(e18.11,1x))
      do i=1,nz
      chgt(i)=0.0
      do j=1,nx
      do k=1,ny
      chgt(i)=chgt(i)+chg(j,k,i)
      enddo
      enddo
      x1=float(i-1)/float(nz)
      chgt(i)=chgt(i)/float(nn)-ef
      write(6+10*is,10) x1*aa*scale,chgt(i)               
10   format(2f12.5)
       enddo

      if (input.eq.1)then
      emax=-9999.0
      do i=1,nz
      if (chgt(i) .gt. emax) emax=chgt(i)
      enddo
      workfn=emax
c     write(6,201) ef
      write(6,301) workfn
      endif
c201  format(' fermi energy =',f10.5)
301  format(' workfunction =',f10.5)

      if ((is.eq.1).and.(ispin.eq.2))then
       if (input.eq.1) then
       read(15,*) (ttt,i=1,natm)

       elseif (input.eq.2) then
       do nn=1,natm
       read(15,'(24x,2i4)') nn1,nn2
       read(15,*) (ttt,i=1,nn2)
       enddo
       read(15,*) (ttt,i=1,natm)
       endif
      endif
600  continue
      stop  
      end

Originally posted by 胡同學(xué) at 2011-05-29 09:53:33:
計算吸附得出的結(jié)果中，功函數(shù)這一項得出work function:left -6.303  right  -6.086eV   這是表示什么意思？  功函數(shù)結(jié)果不是正數(shù)嗎？  怎么得出負(fù)數(shù)結(jié)果，為什么還有左右之分？  請教哪位大牛幫我解釋一下，先表 ...

Originally posted by cavediger at 2011-06-01 18:46:59:
不知道你這是用什么得出來的結(jié)果，但有兩個功函數(shù)是可以理解的，因為你算的是吸附體系，如果只是表面一側(cè)有吸附物的不對稱體系，則相對于上下表面的功函數(shù)會不同，大概你用了什么后處理程序，輸出為左右了.....

Originally posted by 胡同學(xué) at 2011-06-01 20:59:11:
我用VASP算出的，以MedeA平臺，您能說得具體點嗎？我是一邊吸附的，為什么會得出負(fù)數(shù)的結(jié)果，功函數(shù)應(yīng)該是正數(shù)�。�

Originally posted by cavediger at 2011-06-01 22:43:06:
沒有用過MEDIAa，一般是通過LOCALPOT，找出正空能級再減掉費米能級就是功函數(shù)了......

Originally posted by first_yg at 2011-06-02 22:04:14:
請問如何找出真空能級呢？
我在INCAR 中加入IDIPOL =3，可以得到LOCALPOT文件，但是不會處理結(jié)果，麻煩詳細(xì)說一下好嗎？多謝！

Originally posted by first_yg at 2011-06-02 22:28:25:
說錯了，是用LVTOT=.TRUE.
得到LOCPOT文件，然后呢

Originally posted by cavediger at 2011-06-03 08:21:22:
1.自己寫個簡單小程序就可以了，打開LOCPOT文件，里面是SLAB內(nèi)的靜電勢數(shù)據(jù)，對每一個垂直于Z方向面上的靜電勢取平均，然后對Z作圖，真空區(qū)的靜電勢是一個平臺，平臺的高低就是真空能級了。
2.你可以在網(wǎng)上搜 ...

Originally posted by vasp0000 at 2011-06-03 10:37:24:
left和right在outcar中就會有的，結(jié)構(gòu)是對稱的，沒有偶極矩的一般是左右相等。一些情況（例如加電場等）便會左右不等。對于不等的情況，可以通過調(diào)整結(jié)構(gòu)，使上下表面對稱、或者試著添加原子層數(shù)或者真空等。