elk_3Dplot2VESTA

! ---------------------------------------------
program elk_3Dplot2VESTA ! (plot3dFile, title)
    ! Fortran 90 version
    ! compile
    !   gfortran elk_3Dplot2VESTA
    ! using
    !   ./a.out
    ! Reference: Eike Schwier, "3D Visualization of electron density in VESTA", 2012-11-16.

    character(9):: plot3dFile = "RHO3D.OUT" ! plot3dFile
    character(12):: title = "generic_name" ! title
    real(8):: aB = 1.0 ! 0.52917720859 ! stay in bohr for VESTA compatibility with GEOMETRY.OUT

    real(8):: a1, a2, a3 ! components of 1. lattice vector
    real(8):: b1, b2, b3 ! components of 2. lattice vector
    real(8):: c1, c2, c3 ! components of 3. lattice vector

    real(8):: dunit
    character(1):: selec_unit

    dunit = 1.0
    write(6,*) "1: e/Bohr^3,  2: e/Angstrom^3"
    read(5,*) selec_unit
    if( selec_unit == "2" ) then
      dunit = 6.748334 !1/(0.52917720859**3)
      write(6,*) "2: e/Angstrom^3 ( = 6.748334 * e/Bohr^3)"
    else
      write(6,*) "1: e/Bohr^3"
    end if

    call elk_read_geo(a1,a2,a3,b1,b2,b3,c1,c2,c3)

    call elk_3Ddata_inout(a1,a2,a3,b1,b2,b3,c1,c2,c3, aB, plot3dFile, title, dunit)
   
end program elk_3Dplot2VESTA


subroutine elk_read_geo(a1,a2,a3,b1,b2,b3,c1,c2,c3)
    ! ################################################
    ! BEGIN - read GEOMETRY.OUT
    ! ################################################
    
    real(8):: a1, a2, a3 ! components of 1. lattice vector
    real(8):: b1, b2, b3 ! components of 2. lattice vector
    real(8):: c1, c2, c3 ! components of 3. lattice vector
    
    character(12):: GEOfile = "GEOMETRY.OUT"
    character(12):: line

    integer:: g
    open(1, file = GEOfile)
    do g=0, 20

        read(1,'(12a)') line
        if(line(1:5) == "atoms") exit
        if(line(1:4) == "avec") then
            read(1,'(3(f14.9, 4x))' ) a1, a2, a3
            read(1,'(3(f14.9, 4x))' ) b1, b2, b3
            read(1,'(3(f14.9, 4x))' ) c1, c2, c3
        end if
    end do       
    
    return
 
end subroutine elk_read_geo


subroutine elk_3Ddata_inout(a1,a2,a3,b1,b2,b3,c1,c2,c3, aB, plot3dFile, title, dunit)
    ! ################################################
    ! BEGIN - read "3D data".OUT
    ! ################################################

    character(9):: plot3dFile
    character(12):: title
    real(8):: aB

    real(8):: a1, a2, a3 ! components of 1. lattice vector
    real(8):: b1, b2, b3 ! components of 2. lattice vector
    real(8):: c1, c2, c3 ! components of 3. lattice vector

    real(8):: dunit

    character(12):: dumpStr
    integer:: numX, numY, numZ
    real(8):: dumpVar, intensity
    integer:: g

    Open (2, file=plot3dFile)
    read(2,'(3(1x,i5),12a)') numX, numY, numZ, dumpStr ! get number of grid points in the x,y and z direction of plot 3D

    Open (3, file="RHO3D.GRD")
    write(6,*) "Output file: RHO3D.GRD"
    write(3,*) title
    write(3,'(3(f7.5,1x), 3(f6.3, 1x))') &
      sqrt(a1**2.0+a2**2.0+a3**2.0)*aB, sqrt(b1**2.0+b2**2.0+b3**2.0)*aB, sqrt(c1**2.0+c2**2.0+c3**2.0)*aB, &
      elk_vec2ang(b1,b2,b3,c1,c2,c3), elk_vec2ang(a1,a2,a3,c1,c2,c3), elk_vec2ang(a1,a2,a3,b1,b2,b3)
    write(3,'(3(i5,1x))') numX, numY, numZ

    do g=1, int(numX*numY*numZ)

        read(2,'(4e18.14)') dumpVar, dumpVar, dumpVar, intensity
        write(3,*) intensity * dunit

    end do

end subroutine elk_3Ddata_inout


function elk_vec2ang(x1,x2,x3,y1,y2,y3)

    real(8):: x1,x2,x3,y1,y2,y3
    real(8):: pi

    pi = acos(-1.0d0)
    elk_vec2ang = (360/2/pi)*acos( ( x1*y1 + x2*y2 + x3*y3 ) / ( sqrt(x1**2.0+x2**2.0+x3**2.0) * sqrt(y1**2.0+y2**2.0+y3**2.0) ) )

end function elk_vec2ang
! ---------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数
[無料でホームページを作成] [通報・削除依頼]