conquest/tools/PostProcessing/output_module.f90

546 lines
19 KiB
Fortran

! Module for output routines
module output
use datatypes
implicit none
real(double), dimension(95) :: atrad
character(len=2) :: pte(103)
data pte / "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ", "Ne", &
"Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ", "Ca", &
"Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", &
"Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y ", "Zr", &
"Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "St", &
"Sb", "Te", "I ", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", &
"Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", &
"Lu", "Hf", "Ta", "W ", "Re", "Os", "Ir", "Pt", "Au", "Hg", &
"Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", &
"Pa", "U ", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", &
"Md", "No", "Lr"/
contains
subroutine write_xyz(ci)
use datatypes
use numbers
use local, ONLY: root_file
use global_module, ONLY: ni_in_cell, atom_coord, species_glob
use units, ONLY: BohrToAng, AngToBohr
use pseudo_tm_info, ONLY: pseudo
implicit none
! Passed variables
character(len=50), OPTIONAL :: ci
! Local variables
integer :: i
character(len=50) :: filename
write(*,fmt='(4x,"Writing out coordinates in XYZ format")')
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//".xyz"
else
filename = trim(root_file)//".xyz"
end if
open(unit=17,file=filename)
write(17,fmt='(i7)') ni_in_cell
write(17,fmt='("Comment")')
do i=1,ni_in_cell
write(17,fmt='(a2,3f10.3)') pte(int(pseudo(species_glob(i))%z)), &
atom_coord(1,i)*BohrToAng, atom_coord(2,i)*BohrToAng, atom_coord(3,i)*BohrToAng
end do
close(17)
return
end subroutine write_xyz
! Output CASTEP cell file
subroutine write_cell(ci)
use datatypes
use numbers
use local, ONLY: root_file
use global_module, ONLY: ni_in_cell, atom_coord, species_glob
use units, ONLY: BohrToAng, AngToBohr
use pseudo_tm_info, ONLY: pseudo
use dimens, ONLY: r_super_x, r_super_y, r_super_z
implicit none
! Passed variables
character(len=50), OPTIONAL :: ci
! Local variables
integer :: i
character(len=50) :: filename
write(*,fmt='(4x,"Writing out coordinates in CASTEP .cell format")')
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//".cell"
else
filename = trim(root_file)//".cell"
end if
open(unit=17,file=filename)
write(17,fmt='("%BLOCK LATTICE_CART")')
write(17,fmt='("Bohr")')
write(17,fmt='(3f17.12)') r_super_x, zero, zero
write(17,fmt='(3f17.12)') zero, r_super_y, zero
write(17,fmt='(3f17.12)') zero, zero, r_super_z
write(17,fmt='("%ENDBLOCK LATTICE_CART")')
write(17,fmt='("%BLOCK POSITIONS_FRAC")')
do i=1,ni_in_cell
write(17,fmt='(a2,3f12.8)') pte(int(pseudo(species_glob(i))%z)), &
atom_coord(1,i)/r_super_x, atom_coord(2,i)/r_super_y, atom_coord(3,i)/r_super_z
end do
write(17,fmt='("%ENDBLOCK POSITIONS_FRAC")')
close(17)
return
end subroutine write_cell
subroutine write_xsf(ci)
use datatypes
use numbers, only: zero
use dimens, only: r_super_x, r_super_y, r_super_z
use global_module, only: ni_in_cell, iprint_init, atom_coord, &
species_glob
use species_module, only: species_label
use GenComms, only: inode, ionode, cq_abort
use units, only: BohrToAng, HaToeV
use timer_module
use local, only: flag_write_spin_moments, flag_write_forces, root_file
! Passed variables
character(len=50), OPTIONAL :: ci
! Local variables
integer :: lun, i, step
character(len=2) :: atom_name
character(len=50) :: filename
real(double), dimension(ni_in_cell) :: spin_moment
real(double), dimension(3,ni_in_cell) :: tot_force
real(double) :: tote, eup, edn
step = 1
write(*,fmt='(4x,"Writing out coordinates in XSF format")')
if(flag_write_forces) write(*,fmt='(4x,"Including forces")')
if(flag_write_spin_moments) write(*,fmt='(4x,"Including spin moments")')
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//".xsf"
else
filename = trim(root_file)//".xsf"
end if
open(unit=17,file=filename)
write(17,'(a)') "CRYSTAL"
write(17,'("PRIMVEC ",i8)') step
write(17,fmt='(3f14.8)') r_super_x*BohrToAng, zero, zero
write(17,fmt='(3f14.8)') zero, r_super_y*BohrToAng, zero
write(17,fmt='(3f14.8)') zero, zero, r_super_z*BohrToAng
write(17,'("PRIMCOORD ",i8)') step
write(17,fmt='(2i8)') ni_in_cell, 1
if(flag_write_forces) then
write(*,fmt='(4x,"Force output not yet implemented; just writing positions")')
do i=1,ni_in_cell
atom_name = adjustr(species_label(species_glob(i))(1:2))
write(17,'(a4,3f16.8)') atom_name, atom_coord(:,i)*BohrToAng
!write(17,'(a4,6f16.8)') atom_name, atom_coord(:,i)*BohrToAng,&
! tot_force(:,i)*HaToeV/BohrToAng
end do
else if(flag_write_spin_moments) then
! Read spins
open(unit=18,file="AtomCharge.dat", status='old',iostat=i)
if(i>0) then
write(*,fmt='(4x,"Failed to open file AtomCharge.dat for spins; just writing positions")')
do i=1,ni_in_cell
atom_name = adjustr(species_label(species_glob(i))(1:2))
write(17,'(a4,3f16.8)') atom_name, atom_coord(:,i)*BohrToAng
end do
else
do i=1,ni_in_cell
read(18,*) tote, eup, edn
spin_moment(i) = eup - edn
end do
do i=1,ni_in_cell
atom_name = adjustr(species_label(species_glob(i))(1:2))
write(17,'(a4,6f16.8)') atom_name, atom_coord(:,i)*BohrToAng, &
zero,zero,spin_moment(i)
end do
end if
else
do i=1,ni_in_cell
atom_name = adjustr(species_label(species_glob(i))(1:2))
write(17,'(a4,3f16.8)') atom_name, atom_coord(:,i)*BohrToAng
! species_glob(i),flag_move_atom(1,i),flag_move_atom(2,i), &
end do
end if
end subroutine write_xsf
! Write OpenDX file for charge or current density
subroutine write_dx_density(ci)
use datatypes
use numbers
use local, ONLY: current, root_file, nptsx, nptsy, nptsz, grid_x, grid_y, grid_z, gpv, &
nrptx, nrpty, nrptz, nsampx, nsampy, nsampz
use dimens, ONLY: volume
use units, only : BohrToAng
implicit none
character(len=50), OPTIONAL :: ci
character(len=50) :: filename
integer :: nx, ny, nz, icount, nrx, nry, nrz
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//".dx"
else
filename = trim(root_file)//".dx"
end if
open(unit=17,file=filename)
write(17,fmt='("#Now give details of grid")')
write(17,fmt='("object 1 class gridpositions counts ",3i4)') nptsz*nrptz,nptsy*nrpty,nptsx*nrptx
write(17,fmt='("origin 0.0000 0.0000 0.0000")')
write(17,fmt='("delta",3f10.4)') zero, zero, grid_z*BohrToAng
write(17,fmt='("delta",3f10.4)') zero, grid_y*BohrToAng, zero
write(17,fmt='("delta",3f10.4)') grid_x*BohrToAng, zero, zero
write(17,*)
write(17,fmt='("object 2 class gridconnections counts ",3i4)') nptsz*nrptz,nptsy*nrpty,nptsx*nrptx
write(17,*)
write(17,fmt='("object 3 class array type float rank 0 items ",i10," data follows")') nptsx*nptsy*nptsz*nrptx*nrpty*nrptz
icount = 0
do nrz=1,nrptz
do nz=1,nptsz,nsampz
do nry=1,nrpty
do ny=1,nptsy,nsampy
do nrx=1,nrptx
do nx=1,nptsx,nsampx
icount = icount + 1
if(mod(icount,5)==0) then
write(17,fmt='(f15.8)') current(nx,ny,nz)/gpv
else
write(17,fmt='(f15.8)',advance='no') current(nx,ny,nz)/gpv
endif
end do
end do
end do
end do
end do
end do
write(17,*)
write(17,fmt='(" object ",a18," class field")') '"electron density"'
write(17,fmt='(" component ",a11," 1")') '"positions"'
write(17,fmt='(" component ",a13," 2")') '"connections"'
write(17,fmt='(" component ",a6," 3")') '"data"'
!printf("Total electrons: %20.12lf\n Normalised: %20.12lf\n",n_elec,n_elec/vol);
!printf("Maximum density: %20.12lf\n",a_max);
close(unit=17)
end subroutine write_dx_density
subroutine write_dx_coords(ci)
use datatypes
use numbers
use local, ONLY: root_file, stm_x_min, stm_y_min, stm_z_min
use dimens, ONLY: r_super_x, r_super_y, r_super_z, atomicnum
use species_module, ONLY: n_species
use global_module, ONLY: ni_in_cell, atom_coord, species_glob
use units, only : BohrToAng
implicit none
character(len=50), OPTIONAL :: ci
character(len=50) :: filename
character(len=2) :: ele
integer :: i, j, totnabs
integer, dimension(15,ni_in_cell) :: neigh
integer, dimension(ni_in_cell) :: nabs
real(double) :: size1, size2, dist2, d2, dx, dy, dz
call assign_atomic_radii
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//"XYZ.dx"
else
filename = trim(root_file)//"XYZ.dx"
end if
open(unit=17,file=filename)
write(17,fmt='("object 1 class array type float rank 0 items",i7," data follows")') ni_in_cell
do i=1,ni_in_cell
write(17,fmt='(f8.4)') atrad(atomicnum(species_glob(i)))
end do
write(17,fmt='("attribute ""dep"" string ""positions"" ")')
write(17,fmt='("object 2 class array type float rank 1 shape 3 items",i7," data follows")') ni_in_cell
nabs = 0
neigh = 0
totnabs = 0
do i=1,ni_in_cell
write(17,fmt='(3f12.6)') BohrToAng*(atom_coord(1,i)-stm_x_min),BohrToAng*(atom_coord(2,i)-stm_y_min),&
BohrToAng*(atom_coord(3,i)-stm_z_min)
size1 = atrad(atomicnum(species_glob(i)))
! Find neighbours
do j=i+1,ni_in_cell
size2 = atrad(atomicnum(species_glob(j)))
dist2 = (size1*size1 + size2*size2)
dx = BohrToAng*(atom_coord(1,i) - atom_coord(1,j))
dy = BohrToAng*(atom_coord(2,i) - atom_coord(2,j))
dz = BohrToAng*(atom_coord(3,i) - atom_coord(3,j))
d2 = dx*dx + dy*dy + dz*dz
if(d2<dist2) then
nabs(i) = nabs(i) + 1
nabs(j) = nabs(j) + 1
totnabs = totnabs + 2 ! i->j and j->i
neigh(nabs(i),i) = j
neigh(nabs(j),j) = i
end if
end do
end do
write(17,fmt='("attribute ""dep"" string ""positions"" ")')
write(17,fmt='("object 3 class array type int rank 1 shape 2 items ",i7," data follows")') totnabs
do i=1,ni_in_cell
do j=1,nabs(i)
write(17,fmt='(2i7)') i-1,neigh(j,i)-1 ! Indexing from 0
end do
end do
write(17,fmt='("attribute ""element type"" string ""lines"" ",/,"attribute ""ref"" string ""positions"" ")')
write(17,fmt='("#",/,"object ""molecule"" class field")')
write(17,fmt='("component ""data"" value 1",/,"component ""positions"" value 2")')
write(17,fmt='("component ""connections"" value 3",/,"attribute ""name"" string ""molecule"" ")')
write(17,fmt='("#",/,"end")')
close(unit=17)
end subroutine write_dx_coords
subroutine write_cube(data,ci)
use datatypes
use numbers
use local, ONLY: root_file, stm_x_min, stm_y_min, stm_z_min, nsampx, nsampy, nsampz
use dimens, ONLY: r_super_x, r_super_y, r_super_z, volume, atomicnum
use local, ONLY: current, root_file, nptsx, nptsy, nptsz, grid_x, grid_y, grid_z, gpv, nrptx, nrpty, nrptz
use global_module, ONLY: ni_in_cell, atom_coord, species_glob
use units, ONLY: BohrToAng, AngToBohr
implicit none
real(double), dimension(nptsx,nptsy,nptsz) :: data
character(len=50), OPTIONAL :: ci
character(len=50) :: filename
integer :: i, j, icount, nrx, nry, nrz, nx, ny, nz, isym, nz_total
logical :: SymSamp
character(len=3), dimension(7) :: asymm = (/'X ','Y ','XY ','Z ','XZ ','YZ ','XYZ'/)
real(double) :: shiftx, shifty, shiftz
! Check sampling points are symmetric or not
SymSamp = .true.
isym = 0
if (mod(nptsx-1,nsampx).ne.0) isym = 1
if (mod(nptsy-1,nsampy).ne.0) isym = isym + 2
if (mod(nptsz-1,nsampz).ne.0) isym = isym + 4
if (isym.ne.0) then
SymSamp = .false.
write(*,*) 'Warning: sampling points not symmetric for ',asymm(isym)
end if
! Open file
if(PRESENT(ci)) then
filename = trim(ci)//".cube"
else
filename = trim(root_file)//".cube"
end if
open(unit=17,file=filename)
! Header - improve this later
write(17,fmt='("Conquest charge")')
if(SymSamp) then
write(17,fmt='("Sampling points symmetric")')
else
write(17,fmt='("Sampling points asymmetric in directions ",a2)') asymm(isym)
end if
! Atoms, origin location
write(17,fmt='(i5,3f12.6,i5)') ni_in_cell*nrptx*nrpty*nrptz, zero, zero, zero, 1
! Numbers of grid points, grid increment (i.e. spacing)
write(17,fmt='(i5,3f12.6)') nptsx*nrptx/nsampx,grid_x,zero,zero
write(17,fmt='(i5,3f12.6)') nptsy*nrpty/nsampy,zero,grid_y,zero
write(17,fmt='(i5,3f12.6)') nptsz*nrptz/nsampz,zero,zero,grid_z
! Atomic coordinates
do nrx = 1,nrptx
shiftx = real(nrx-1,double)*r_super_x
do nry = 1,nrpty
shifty = real(nry-1,double)*r_super_y
do nrz = 1,nrptz
shiftz = real(nrz-1,double)*r_super_z
do i=1,ni_in_cell
write(17,fmt='(i5,4f12.6)') atomicnum(species_glob(i)),real(atomicnum(species_glob(i)),double),&
atom_coord(1,i)-stm_x_min + shiftx, &
atom_coord(2,i)-stm_y_min + shifty, &
atom_coord(3,i)-stm_z_min + shiftz
!write(17,fmt='(i5,4f12.6)') species_glob(i),zero,(atom_coord(j,i)*BohrToAng,j=1,3)
end do
end do
end do
end do
! Charge density, z fastest
nz_total = nptsz*nrptz/nsampz
do nrx=1,nrptx
do nx=1,nptsx,nsampx
do nry=1,nrpty
do ny=1,nptsy,nsampy
icount = 0
do nrz=1,nrptz
do nz=1,nptsz,nsampz
icount = icount + 1
if ((icount==nz_total).or.(mod(icount,6)==0)) then
write(17,fmt='(e13.5)') data(nx,ny,nz)!/gpv
else
write(17,fmt='(e13.5)',advance='no') data(nx,ny,nz)!/gpv
endif
end do
end do
end do
end do
end do
end do
close(unit=17)
end subroutine write_cube
! Taken from density_module.f90
subroutine assign_atomic_radii
use datatypes
use numbers, only: zero
use units, only: AngToBohr
implicit none
atrad = zero
atrad(1)=0.35_double*AngToBohr
atrad(2)=0.28_double*AngToBohr
atrad(3)=1.45_double*AngToBohr ! Li
atrad(4)=1.05_double*AngToBohr
atrad(5)=0.85_double*AngToBohr
atrad(6)=0.70_double*AngToBohr
atrad(7)=0.65_double*AngToBohr
atrad(8)=0.60_double*AngToBohr
atrad(9)=0.50_double*AngToBohr ! F
atrad(10)=0.58_double*AngToBohr
atrad(11)=1.80_double*AngToBohr ! Na
atrad(12)=1.50_double*AngToBohr
atrad(13)=1.25_double*AngToBohr
atrad(14)=1.10_double*AngToBohr
atrad(15)=1.00_double*AngToBohr
atrad(16)=1.00_double*AngToBohr
atrad(17)=1.00_double*AngToBohr ! Cl
atrad(18)=1.06_double*AngToBohr
atrad(19)=2.20_double*AngToBohr ! K
atrad(20)=1.80_double*AngToBohr
atrad(21)=1.60_double*AngToBohr
atrad(22)=1.40_double*AngToBohr
atrad(23)=1.35_double*AngToBohr
atrad(24)=1.40_double*AngToBohr
atrad(25)=1.40_double*AngToBohr
atrad(26)=1.40_double*AngToBohr ! Fe
atrad(27)=1.35_double*AngToBohr
atrad(28)=1.35_double*AngToBohr
atrad(29)=1.35_double*AngToBohr
atrad(30)=1.35_double*AngToBohr
atrad(31)=1.30_double*AngToBohr ! Ga
atrad(32)=1.25_double*AngToBohr
atrad(33)=1.15_double*AngToBohr
atrad(34)=1.15_double*AngToBohr
atrad(35)=1.15_double*AngToBohr ! Br
atrad(36)=1.16_double*AngToBohr
atrad(37)=2.35_double*AngToBohr ! Rb
atrad(38)=2.00_double*AngToBohr
atrad(39)=1.80_double*AngToBohr
atrad(40)=1.55_double*AngToBohr ! Zr
atrad(41)=1.45_double*AngToBohr
atrad(42)=1.45_double*AngToBohr
atrad(43)=1.35_double*AngToBohr
atrad(44)=1.30_double*AngToBohr ! Ru
atrad(45)=1.35_double*AngToBohr
atrad(46)=1.40_double*AngToBohr
atrad(47)=1.60_double*AngToBohr ! Ag
atrad(48)=1.55_double*AngToBohr
atrad(49)=1.55_double*AngToBohr ! In
atrad(50)=1.45_double*AngToBohr
atrad(51)=1.45_double*AngToBohr
atrad(52)=1.40_double*AngToBohr
atrad(53)=1.40_double*AngToBohr ! I
atrad(54)=1.40_double*AngToBohr
atrad(55)=2.60_double*AngToBohr ! Cs
atrad(56)=2.15_double*AngToBohr
atrad(57)=1.95_double*AngToBohr ! La
atrad(58)=1.85_double*AngToBohr
atrad(59)=1.85_double*AngToBohr
atrad(60)=1.85_double*AngToBohr
atrad(61)=1.85_double*AngToBohr
atrad(62)=1.85_double*AngToBohr
atrad(63)=1.85_double*AngToBohr ! Eu
atrad(64)=1.80_double*AngToBohr
atrad(65)=1.75_double*AngToBohr
atrad(66)=1.75_double*AngToBohr
atrad(67)=1.75_double*AngToBohr
atrad(68)=1.75_double*AngToBohr
atrad(69)=1.75_double*AngToBohr
atrad(70)=1.75_double*AngToBohr
atrad(71)=1.75_double*AngToBohr ! Lu
atrad(72)=1.55_double*AngToBohr
atrad(73)=1.45_double*AngToBohr
atrad(74)=1.35_double*AngToBohr ! W
atrad(75)=1.35_double*AngToBohr
atrad(76)=1.30_double*AngToBohr
atrad(77)=1.35_double*AngToBohr ! Ir
atrad(78)=1.35_double*AngToBohr
atrad(79)=1.35_double*AngToBohr
atrad(80)=1.50_double*AngToBohr ! Hg
atrad(81)=1.90_double*AngToBohr
atrad(82)=1.80_double*AngToBohr
atrad(83)=1.60_double*AngToBohr
atrad(84)=1.90_double*AngToBohr ! Po
atrad(85)=1.50_double*AngToBohr
atrad(86)=1.50_double*AngToBohr
atrad(87)=2.60_double*AngToBohr
atrad(88)=2.15_double*AngToBohr ! Ra
atrad(89)=1.95_double*AngToBohr
atrad(90)=1.80_double*AngToBohr
atrad(91)=1.80_double*AngToBohr
atrad(92)=1.75_double*AngToBohr
atrad(93)=1.75_double*AngToBohr
atrad(94)=1.75_double*AngToBohr
atrad(95)=1.75_double*AngToBohr
return
end subroutine assign_atomic_radii
subroutine write_banner
use datestamp, ONLY: datestr, commentver
implicit none
character(len=10) :: today, the_time
write(*,fmt='(/"CONQUEST charge density conversion and STM image simulation"/)')
write(*,fmt='("D. R. Bowler (UCL) and A. Nakata (NIMS)")')
call date_and_time(today, the_time)
write(*,fmt='(/4x,"This job was run on ",a4,"/",a2,"/",a2," at ",a2,":",a2,/)') &
today(1:4), today(5:6), today(7:8), the_time(1:2), the_time(3:4)
write(*,&
'(/4x,"Code compiled on: ",a,/10x,"Version comment: ",/6x,a//)') &
datestr, commentver
end subroutine write_banner
end module output