VASP 自旋极化计算

vasp 计算摘记

(2010-11-17 03:57:01)

转载

标签： vasp

分类： 工作问题

1， vasp soc计算
a）非自洽计算SOC
      首先进行自旋极化计算（当然，如果体系没有磁矩，如Bi2Se3，初始磁矩取为零）。
      然后进行非自洽SOC计算：
               ICHARG = 11
        
                Spin orbit
                LSORBIT = .TRUE.
               SAXIS   =  0 0 1   ：因为之前进行了自旋极化计算，这里不需要重新指定MAGMOM。
               NBANDS  =  160 ：vasp要求NBNADS在自旋极化计算和soc计算时有精确两倍关系。
b) 全自洽计算SOC

2， vasp 考察体系的dipole potential
      在做全自洽计算时，打开选项
       LVTOT = T
     
       此后会生成LOCPOT文件。此文件可以用下面小程序读取：
c!
c!   Program calculates the line plot or plane-averaged plot of the
c!   3D grid files, CHGCAR, LOCPOT, PARCHG
c!  
        dimension coord(5000,3),vec1(3),a1(3),a2(3),a3(3),
     &            veci(3),vec0(3),nntype(10),vec2(3)
        character aaa
        character*20 ifile,ofile
        real, allocatable:: ev(:)

        write(6,*)'Enter the name of potential file:
     &          CHGCAR,LOCPOT,PARCHG'
        read(5,*)ifile
        write(6,*)ifile
        if(ifile.eq.'CHGCAR')then
            ndiv = 5
            write(6,*)'Total Charge density: Datatype=',ndiv
                elseif(ifile.eq.'LOCPOT')then
            ndiv = 5
            write(6,*)'Local potential: Datatype=',ndiv
                else
            ndiv = 10
            write(6,*)'partial charge density: Datatype=',ndiv
        endif

        open(unit=13,file=ifile,status='OLD',iostat=ierror)
        read(13,*)
        read(13,*)
        read(13,*)(a1(k),k=1,3)
        read(13,*)(a2(k),k=1,3)
        read(13,*)(a3(k),k=1,3)
        write(6,*)'Enter number of atom types'
        read(5,*)ntype
        read(13,*)(nntype(k),k=1,ntype)

C***************************************************************
C***************************************************************
C***************************************************************
C************   PLEASE input Z0 and Lattice-3  *****************
C***************************************************************
C***************************************************************
C***************************************************************
        natom = 0
        do i = 1,ntype
           natom = natom + nntype(i)
        enddo

        read(13,*)
            vec1(:) = 0.0
        do i = 1,natom
           read(13,*)(coord(i,k),k=1,3)
               vec1(:) = vec1(:) + coord(i,:)/natom
        enddo

                write(6,*)'Check the Origin '
                vec2(:) = vec1(1)*a1(:) + vec1(2)*a2(:) + vec1(3)*a3(:)
                write(6,*)vec2
                write(6,*)
     &           'Check the Origin: Is it correct for your purpose ???'
                write(6,*)
                write(6,*)
                write(6,*)

        read(13,*)
        read(13,*)n1,n2,n3
        write(6,120)n1,n2,n3

                allocate(ev(n1*n2*n3))
            ev(:) = 0.0

c       do i = 1,n1
c          write(6,*)i,a1*(i-1)/n1
c       enddo
c       do j = 1,n2
c          write(6,*)j,a2*(j-1)/n2
c       enddo
 
 118    format(a3,3f10.5)
 119    format('##',a3,3f10.5)
 120    format('##',3I5)
        write(6,119)'a1',(a1(k),k=1,3)
        write(6,119)'a2',(a2(k),k=1,3)
        write(6,119)'a3',(a3(k),k=1,3)
        Nline = n1*n2*n3/ndiv
        if(Nline*ndiv .eq. n1*n2*n3)then
          do iline = 1,Nline
             indg = (iline-1)*ndiv
             read(13,*)(ev(indg +k),k=1,ndiv)
          enddo
        else
          nrem = n1*n2*n3 - Nline*ndiv
          do iline = 1,Nline
             indg = (iline-1)*ndiv
             read(13,*)(ev(indg +k),k=1,ndiv)
          enddo
             indg = (Nline-0)*ndiv
             read(13,*)(ev(indg +k),k=1,nrem)
        endif

C==================================================
         write(*,*)"move the origin to zero? (yes=1/no=0)"
         read(*,*) ishift
          if(ishift .eqv. 0) then
            write(*,*)"OK, No shift !"
            vec2(:)=0
           else
             write(*,*) "Well, shift !"
            endif
C==================================================


C ======Line Plot=========================================
C       n1p = n1/2 + 1
C       n2p = n2/2 + 1
C       xplot = (n1p -1)*a1(1)/float(n1)
C       yplot = (n2p -1)*a2(2)/float(n2)
        x0 =vec2(1)
                y0 =vec2(2)
                z0 =vec2(3)
                n1p = nint(n1*x0/a1(1)) + 1
                n2p = nint(n2*y0/a2(2)) + 1
                n3p = nint(n3*z0/a3(3)) + 1
                xplot = (n1p -1)*a1(1)/float(n1)
                yplot = (n2p -1)*a2(2)/float(n2)
                zplot = (n3p -1)*a3(3)/float(n3)
C---------Along the Z axis at the center point------------     
        open(24,file ='vloc-line-z.dat')
        write(24,*)'# Potential line plot along z-direction'
        write(24,*)'# Orthorombinc cell is assumed'
        DO k = 1,n3
           zi = a3(3)*(k-1)/n3

          do i = n1p,n1p
          do j = n2p,n2p
             iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
            write(24,128)xplot,yplot,zi,ev(iadd1)
          enddo; enddo
        ENDDO
 128    format(3f12.7,3x,f10.5)
        close(unit=24)
C---------Along the X axis at the center point------------     
        open(24,file ='vloc-line-x.dat')
        write(24,*)'# Potential line plot along x-direction'
        write(24,*)'# Orthorombinc cell is assumed'
        DO i = 1,n1
           xi = a1(1)*(i-1)/n1

          do j = n2p,n2p
          do k = n3p,n3p
             iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
            write(24,128)xi,yplot,zplot,ev(iadd1)
          enddo; enddo
        ENDDO
        close(unit=24)
C---------Along the Y axis at the center point------------     
        open(24,file ='vloc-line-y.dat')
        write(24,*)'# Potential line plot along y-direction'
        write(24,*)'# Orthorombinc cell is assumed'
        DO j = 1,n2
           yi = a2(2)*(j-1)/n2

          do i = n1p,n1p
          do k = n3p,n3p
             iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
            write(24,128)xplot,yi,zplot,ev(iadd1)
          enddo; enddo
        ENDDO
        close(unit=24)
C====== Average in the XZ-plane  ==============================
        open(24,file ='vloc_XZaveraged.dat')
        DO j = 1,n2
           yi = a2(2)*(j-1)/n2
           avg = 0.0

        do i = 1,n1
        do k = 1,n3
           iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
           veci(:) = 1.*(i-1)/n1 *a1(:) + 1.*(k-1)/n3 *a3(:)

           avg = avg + ev(iadd1)/(n1*n3*1.0)

        enddo; enddo
           write(24,*)yi,avg
        enddo
                close(unit=24)
C=========================================================
C====== Average in the YZ-plane  =========================
        open(24,file ='vloc_YZaveraged.dat')
        DO i = 1,n1
           xi = a1(1)*(i-1)/n1
           avg = 0.0

        do j = 1,n2
        do k = 1,n3
           iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
           veci(:) = 1.*(j-1)/n2 *a1(:) + 1.*(k-1)/n3 *a3(:)

           avg = avg + ev(iadd1)/(n2*n3*1.0)

        enddo; enddo
           write(24,*)xi,avg
        enddo
                close(unit=24)
C=========================================================
C====== Average in the XY-plane  ==============================
        open(24,file ='vloc_XYaveraged.dat')
        DO k = 1,n3
           zi = a3(3)*(k-1)/n3
           avg = 0.0

        do i = 1,n1
        do j = 1,n2
           iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
           veci(:) = 1.*(i-1)/n1 *a1(:) + 1.*(j-1)/n2 *a2(:)
           xi = veci(1)
           yi = veci(2)

           avg = avg + ev(iadd1)/(n1*n2*1.0)

        enddo; enddo
           write(24,*)(zi-z0),avg
        enddo
C=========================================================
        z0 = 0.0
                xL3 = a3(3)
C***************************************************************
C************   PLEASE input Z0 and Lattice-3  *****************
C***************************************************************

        DO k = 1,n3
           zi = a3(3)*(k-1)/n3
           avg = 0.0

        do i = 1,n1
        do j = 1,n2
           iadd1 = (k-1)*n1*n2 + (j-1)*n1 + i
           veci(:) = 1.*(i-1)/n1 *a1(:) + 1.*(j-1)/n2 *a2(:)
           xi = veci(1)
           yi = veci(2)

           avg = avg + ev(iadd1)/(n1*n2*1.0)

        enddo; enddo
c          write(24,*)(zi-z0+xL3),avg
        enddo
C=========================================================

         stop
         end


  3, atomic orbital projection
   LORBIT = 2


ion s py pz px dxy dyz dz2 dxz dx2 tot
1 0.000 0.000 0.000 0.000 0.395 0.039 0.117 0.039 0.397 0.987
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
tot 0.000 0.000 0.000 0.000 0.395 0.039 0.117 0.039 0.397 0.987
1 0.000 0.000 0.000 0.000 0.000 -0.002 -0.002 0.000 -0.001 -0.005
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
tot 0.000 0.000 0.000 0.000 0.000 -0.002 -0.002 0.000 -0.001 -0.005
1 0.000 0.000 0.000 0.000 0.001 0.000 0.001 0.001 0.000 0.003
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
tot 0.000 0.000 0.000 0.000 0.001 0.000 0.001 0.001 0.000 0.003
1 0.000 0.000 0.000 0.000 0.395 0.039 -0.117 0.039 0.397 0.752
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
tot 0.000 0.000 0.000 0.000 0.395 0.039 -0.117 0.039 0.397 0.753

Are the bolded lines indicate the Mx, My and Mz components of the magnetization. In particular, are the values indicated in green, blue and red referring to total-Mx, total-My and total-Mz, respectively. My spin quantization axis is as the default input [0 0 1]