将ELFCAR、CHGCAR转化成MS阅读可读的fortran代码

最新推荐文章于 2025-11-07 14:44:46 发布

原创最新推荐文章于 2025-11-07 14:44:46 发布 · 105 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#科技

VASP 专栏收录该内容

28 篇文章

订阅专栏

PROGRAM CHGCAR2GRD

  IMPLICIT none

  CHARACTER(*), PARAMETER :: vasp_form = '(10(1X,E11.5))'
  CHARACTER(*), PARAMETER :: dmol_form = '(1p,e12.5)'
  REAL(8) :: pi
  CHARACTER*80 HEADER
  REAL(8), DIMENSION(3,3) :: A
  INTEGER, DIMENSION(10) :: num, iatom
  REAL(8) :: scale
  REAL(8), DIMENSION(3) :: length, angle
  REAL(8), DIMENSION(:,:,:), POINTER :: work1, work2
  INTEGER :: nspin, nx, ny, nz, sum, number, i, j, k, l

  pi= 4.*atan(1.)
  nspin=1
  do i=1,10
      num(i)=0
  end do

  OPEN(unit=16,file='ELFCAR')
  READ(16,*)
 
  OPEN(unit=15,file='vasp.grd')
 
  WRITE(15,'(A)') 'VASP charge density'
  WRITE(15,'(A)') dmol_form
 
  READ(16,*) scale
  DO i=1,3
     READ(16,*) (A(j,i),j=1,3)
     length(i)=sqrt(A(1,i)**2+A(2,i)**2+A(3,i)**2)*scale
  END DO
 
  angle(1)=acos(dot_product(A(:,2),A(:,3))/(length(2)*length(3)))*180/pi
  angle(2)=acos(dot_product(A(:,3),A(:,1))/(length(3)*length(1)))*180/pi
  angle(3)=acos(dot_product(A(:,1),A(:,2))/(length(1)*length(2)))*180/pi

  WRITE(15,'(6F8.3)') (length(i),i=1,3),(angle(i),i=1,3)
 
  !number=0
  !sum=0
  !print*,"How many kinds of the atoms in you CHGCAR?"
  !print*,"The number of the kinds is:"
  !read*,number
  !read(16, *, end=20) (iatom(i), i=1, number)
  !20 continue
  READ(16,'(A)') HEADER  
  i=0; j=0; k=0; l=0
  READ(HEADER,*,END=12) i,j,k,l 
  12  sum=i+j+k+l
  !do i=1,number
  !sum=sum+iatom(i)
  !end do
  read (16,*)
  do i = 1, sum
      read (16,*)
  end do
  read (16,*)
 
  READ(16,*) nx,ny,nz
  WRITE(15,'(3I5)') nx-1,ny-1,nz-1
  WRITE(15,'(7I5)') 1,0,nx-1,0,ny-1,0,nz-1

  IF(NSPIN==1) THEN
     ALLOCATE(work1(nx,ny,nz))
     READ(16,vasp_form) (((work1(i,j,k),i=1,nx),j=1,ny),k=1,nz)
     WRITE(15,dmol_form) (((work1(i,j,k),i=1,nx),j=1,ny),k=1,nz)
     DEALLOCATE(work1)
  ELSE
     ALLOCATE(work1(nx,ny,nz),work2(nx,ny,nz))
     READ(16,vasp_form) (((work1(i,j,k),i=1,nx),j=1,ny),k=1,nz)
     READ(16,*)
     READ(16,vasp_form) (((work2(i,j,k),i=1,nx),j=1,ny),k=1,nz)
     work1=work1-work2
     WRITE(15,dmol_form) (((work1(i,j,k),i=1,nx),j=1,ny),k=1,nz)
     DEALLOCATE(work1,work2)
  END IF

  CLOSE(15)
  CLOSE(16)

END PROGRAM CHGCAR2GRD