You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1105 lines
41 KiB
1105 lines
41 KiB
module io
|
|
|
|
use elements
|
|
use parameters
|
|
use atoms
|
|
use box
|
|
|
|
implicit none
|
|
|
|
integer :: outfilenum = 0, infilenum = 0
|
|
character(len=100) :: outfiles(100), infiles(100), in_lattice_type=''
|
|
logical :: force_overwrite
|
|
real(kind=dp) :: in_lapa=0.0
|
|
public
|
|
contains
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutines for writing out data files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
subroutine get_out_file(filename)
|
|
|
|
implicit none
|
|
|
|
character(len=100), intent(in) :: filename
|
|
character(len=100) :: temp_outfile
|
|
character(len=1) :: overwrite
|
|
logical :: file_exists
|
|
|
|
!If no filename is provided then this function is called with none and prompts user input
|
|
if (filename=='none') then
|
|
print *, "Please specify a filename or extension to output to:"
|
|
read(*,*) temp_outfile
|
|
else
|
|
temp_outfile = filename
|
|
end if
|
|
|
|
!Infinite loop which only exists if user provides valid filetype
|
|
overwrite = 'r'
|
|
do while(.true.)
|
|
|
|
!Check to see if file exists, if it does then ask user if they would like to overwrite the file
|
|
inquire(file=trim(temp_outfile), exist=file_exists)
|
|
if (file_exists.and.(.not.(force_overwrite))) then
|
|
if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)"
|
|
read(*,*) overwrite
|
|
if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then
|
|
print *, "Please specify a new filename with extension:"
|
|
read(*,*) temp_outfile
|
|
else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then
|
|
continue
|
|
else
|
|
print *, "Please pick either y or n"
|
|
read(*,*) overwrite
|
|
end if
|
|
|
|
end if
|
|
|
|
if (scan(temp_outfile,'.',.true.) == 0) then
|
|
print *, "No extension included on filename, please type a full filename that includes an extension."
|
|
read(*,*) temp_outfile
|
|
cycle
|
|
end if
|
|
select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
|
|
case('xyz', 'lmp', 'vtk', 'mb', 'restart')
|
|
outfilenum=outfilenum+1
|
|
outfiles(outfilenum) = temp_outfile
|
|
exit
|
|
case('cac')
|
|
lmpcac = .true.
|
|
outfilenum=outfilenum+1
|
|
outfiles(outfilenum) = temp_outfile
|
|
exit
|
|
case default
|
|
print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
|
|
"please input a filename with extension from following list: xyz, lmp, vtk, cac."
|
|
read(*,*) temp_outfile
|
|
|
|
end select
|
|
end do
|
|
|
|
end subroutine get_out_file
|
|
|
|
subroutine write_out
|
|
!This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
|
|
|
|
integer :: i
|
|
|
|
!Find max esize which will be needed later
|
|
call set_max_esize
|
|
|
|
do i = 1, outfilenum
|
|
|
|
print *, "Writing data out to ", trim(adjustl(outfiles(i)))
|
|
|
|
!Pull out the extension of the file and call the correct write subroutine
|
|
select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
|
|
case('xyz')
|
|
call write_xyz(outfiles(i))
|
|
case('lmp')
|
|
call write_lmp(outfiles(i))
|
|
case('vtk')
|
|
call write_vtk(outfiles(i))
|
|
case('mb')
|
|
call write_mb(outfiles(i))
|
|
case('restart')
|
|
call write_pycac(outfiles(i))
|
|
case('cac')
|
|
call write_lmpcac(outfiles(i))
|
|
case default
|
|
print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
|
|
" is not accepted for writing. Please select from: xyz and try again"
|
|
stop
|
|
|
|
end select
|
|
end do
|
|
end subroutine write_out
|
|
|
|
|
|
|
|
subroutine write_xyz(file)
|
|
!This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
|
|
character(len=100), intent(in) :: file
|
|
|
|
integer :: i, inod, ibasis
|
|
|
|
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
!Write total number of atoms + elements
|
|
write(11, '(i16)') node_num+atom_num
|
|
|
|
!Write comment line
|
|
write(11, '(a)') "#Node + atom file created using cacmb"
|
|
|
|
!Write nodal positions
|
|
do i = 1, ele_num
|
|
do inod = 1, ng_node(lat_ele(i))
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 0, r_node(:,ibasis,inod,i)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
!Write atom positions
|
|
do i = 1, atom_num
|
|
write(11, '(2i16, 3f23.15)') type_atom(i), 1, r_atom(:,i)
|
|
end do
|
|
|
|
!Finish writing
|
|
close(11)
|
|
end subroutine write_xyz
|
|
|
|
subroutine write_lmp(file)
|
|
!This subroutine writes out a .lmp style dump file
|
|
character(len=100), intent(in) :: file
|
|
integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3), interp_num
|
|
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass
|
|
|
|
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
!Comment line
|
|
write(11, '(a)') '# lmp file made with cacmb'
|
|
write(11, '(a)')
|
|
!Calculate total atom number
|
|
write_num = atom_num
|
|
do i = 1,ele_num
|
|
if(type_ele(i) == 'fcc') write_num = write_num + basisnum(lat_ele(i))*size_ele(i)**3
|
|
end do
|
|
!Write total number of atoms + elements
|
|
write(11, '(i16, a)') write_num, ' atoms'
|
|
!Write number of atom types
|
|
write(11, '(i16, a)') atom_types, ' atom types'
|
|
|
|
write(11,'(a)') ' '
|
|
!Write box bd
|
|
write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
|
|
write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
|
|
write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'
|
|
|
|
!Masses
|
|
write(11, '(a)') 'Masses'
|
|
|
|
write(11, '(a)') ' '
|
|
do i =1, atom_types
|
|
call atommass(type_to_name(i),mass)
|
|
write(11, '(i16, f23.15)') i, mass
|
|
end do
|
|
write(11, '(a)') ' '
|
|
|
|
!Write atom positions
|
|
write(11, '(a)') 'Atoms'
|
|
write(11, '(a)') ' '
|
|
do i = 1, atom_num
|
|
write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i)
|
|
end do
|
|
|
|
!Write refined element atomic positions
|
|
interp_num = 0
|
|
do i = 1, ele_num
|
|
call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
|
|
select case(trim(adjustl(type_ele(i))))
|
|
case('fcc','bcc')
|
|
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
|
|
interp_num = interp_num+1
|
|
call apply_periodic(r_interp(:,iatom))
|
|
write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
|
|
end do
|
|
end select
|
|
end do
|
|
end subroutine write_lmp
|
|
|
|
subroutine write_lmpcac(file)
|
|
!This subroutine writes out a .lmp style dump file
|
|
character(len=100), intent(in) :: file
|
|
integer :: write_num, i, inod, ibasis
|
|
real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3)
|
|
|
|
1 format(i16, ' Eight_Node', 4i16)
|
|
2 format(i16, ' Atom', 4i16)
|
|
3 format(3i16,3f23.15)
|
|
|
|
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
!Comment line
|
|
write(11, '(a)') '# CAC input file made with cacmb'
|
|
write(11, '(a)')
|
|
!Calculate total atom number
|
|
write_num = atom_num + ele_num
|
|
|
|
!Write total number of atoms + elements
|
|
write(11, '(i16, a)') write_num, ' cac elements'
|
|
!Write number of atom types
|
|
write(11, '(i16, a)') atom_types, ' atom types'
|
|
|
|
write(11,'(a)') ' '
|
|
!Write box bd
|
|
write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
|
|
write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
|
|
write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'
|
|
|
|
!Masses
|
|
write(11, '(a)') 'Masses'
|
|
|
|
write(11, '(a)') ' '
|
|
do i =1, atom_types
|
|
call atommass(type_to_name(i),mass)
|
|
write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i)
|
|
end do
|
|
write(11, '(a)') ' '
|
|
|
|
write(11, '(a)') 'CAC Elements'
|
|
write(11, '(a)') ' '
|
|
|
|
!Set up the nodal adjustment variables for all the different element types. This adjusts the node centers
|
|
!from the center of the unit cell (as formulated in this code) to the corners of the unit cells
|
|
do inod = 1, 8
|
|
do i = 1,3
|
|
if(is_equal(cubic_cell(i, inod),0.0_dp)) then
|
|
fcc_adjust(i,inod) = -0.5_dp
|
|
else
|
|
fcc_adjust(i, inod) = 0.5_dp
|
|
end if
|
|
end do
|
|
end do
|
|
fcc_adjust = matmul(fcc_mat, fcc_adjust)
|
|
|
|
!Write element nodal positions
|
|
do i = 1, ele_num
|
|
select case(trim(adjustl(type_ele(i))))
|
|
case('fcc')
|
|
!Now orient the current adjustment vector to the correct orientation
|
|
local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i))
|
|
!The first entry is the element specifier
|
|
write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
do inod = 1, 8
|
|
!Nodal information for every node
|
|
rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod)
|
|
write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout
|
|
end do
|
|
end do
|
|
end select
|
|
end do
|
|
|
|
do i = 1, atom_num
|
|
!Element specifier dictating that it is an atom
|
|
write(11,2) ele_num+i, 1, 1, 1, 1
|
|
!Write the atomic information
|
|
write(11,3) 1, 1, type_atom(i), r_atom(:,i)
|
|
end do
|
|
end subroutine write_lmpcac
|
|
|
|
subroutine write_vtk(file)
|
|
!This subroutine writes out a vtk style dump file
|
|
integer :: i, j, inod, ibasis
|
|
character(len=100), intent(in) :: file
|
|
|
|
1 format('# vtk DataFile Version 4.0.1', / &
|
|
'CAC output -- cg', / &
|
|
'ASCII')
|
|
11 format('# vtk DataFile Version 4.0.1', / &
|
|
'CACmb output -- atoms', / &
|
|
'ASCII')
|
|
2 format('DATASET UNSTRUCTURED_GRID')
|
|
3 format('POINTS', i16, ' float')
|
|
4 format(/'CELLS', 2i16)
|
|
5 format(/'CELL_TYPES', i16)
|
|
12 format(/'CELL_DATA', i16)
|
|
16 format(/'POINT_DATA', i16)
|
|
17 format('SCALARS weight float', / &
|
|
'LOOKUP_TABLE default')
|
|
18 format('SCALARS atom_type float', / &
|
|
'LOOKUP_TABLE default')
|
|
|
|
20 format('SCALARS lattice_type float', /&
|
|
'LOOKUP_TABLE default')
|
|
|
|
21 format('SCALARS esize float', /&
|
|
'LOOKUP_TABLE default')
|
|
|
|
!First we write the vtk file containing the atoms
|
|
open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
write(11, 11)
|
|
write(11, 2)
|
|
write(11, 3) atom_num
|
|
do i = 1, atom_num
|
|
write(11, '(3f23.15)') r_atom(:,i)
|
|
end do
|
|
write(11,4) atom_num, atom_num*2
|
|
do i = 1, atom_num
|
|
write(11, '(2i16)') 1, i-1
|
|
end do
|
|
write(11, 5) atom_num
|
|
do i = 1, atom_num
|
|
write(11, '(i16)') 1
|
|
end do
|
|
write(11, 16) atom_num
|
|
write(11, 18)
|
|
do i = 1, atom_num
|
|
write(11, '(i16)') type_atom(i)
|
|
end do
|
|
close(11)
|
|
|
|
open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
write(11,1)
|
|
write(11,2)
|
|
write(11,3) node_num
|
|
do i = 1, ele_num
|
|
do inod=1, ng_node(lat_ele(i))
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
write(11, '(3f23.15)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
|
|
end do
|
|
end do
|
|
end do
|
|
write(11, 4) ele_num, ele_num + node_num
|
|
do i =1, ele_num
|
|
write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1)
|
|
end do
|
|
write(11,5) ele_num
|
|
do i = 1, ele_num
|
|
if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
|
|
end do
|
|
write(11,12) ele_num
|
|
write(11,20)
|
|
do i = 1, ele_num
|
|
write(11, '(i16)') lat_ele(i)
|
|
end do
|
|
write(11,21)
|
|
do i = 1, ele_num
|
|
write(11, '(i16)') size_ele(i)
|
|
end do
|
|
close(11)
|
|
end subroutine
|
|
|
|
subroutine write_pycac(file)
|
|
!This subroutine writes restart files meant to be used with the McDowell Group CAC code.
|
|
!NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the
|
|
!each element has only 1 atom type at the node.
|
|
character(len=100), intent(in) :: file
|
|
integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_num, &
|
|
etype
|
|
real(kind=dp) :: box_vec(3)
|
|
|
|
1 format('time' / i16, f23.15)
|
|
2 format('number of elements' / i16)
|
|
3 format('number of nodes' / i16)
|
|
4 format('element types' / i16)
|
|
5 format('number of atoms' / i16)
|
|
6 format('number of grains' / i16)
|
|
7 format('boundary ' / 3a1)
|
|
8 format('box bound' / 6f23.15)
|
|
9 format('box length' / 3f23.15)
|
|
10 format('box matrix')
|
|
11 format(3f23.15)
|
|
12 format('coarse-grained domain')
|
|
13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis x y z')
|
|
14 format('atomistic domain' / 'ia grain_atom type_atom x y z')
|
|
15 format('maximum lattice periodicity length' / 3f23.15)
|
|
16 format('Number of lattice types and atom types '/ 2i16)
|
|
17 format('lattice type IDs')
|
|
18 format('lattice types for grains')
|
|
19 format('max nodes per element' / i16)
|
|
20 format('max interpo per element' / i16)
|
|
21 format('atom types to elements')
|
|
|
|
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
write(11,1) timestep, total_time
|
|
write(11,2) ele_num
|
|
|
|
!Below writes the header information for the restart file
|
|
|
|
|
|
!First figure out all of the unique element types
|
|
unique_num = 0
|
|
unique_index(:) = 0
|
|
eleloop:do i = 1, ele_num
|
|
do j =1 , unique_num
|
|
if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. &
|
|
( lat_ele(i) == lat_ele(unique_index(j)) ) ) then
|
|
cycle eleloop
|
|
end if
|
|
end do
|
|
unique_num = unique_num + 1
|
|
unique_index(unique_num) = i
|
|
unique_size(unique_num) = size_ele(i)
|
|
end do eleloop
|
|
|
|
!Calculate the max number of atoms per element
|
|
select case(max_ng_node)
|
|
case(8)
|
|
interp_max = (max_esize)**3
|
|
case default
|
|
interp_max = 0
|
|
end select
|
|
write(11,20) interp_max
|
|
write(11,3) node_num
|
|
write(11,19) max_ng_node
|
|
write(11,4) unique_num
|
|
write(11,5) atom_num
|
|
write(11,6) 1 !Grain_num is ignored
|
|
write(11,16) lattice_types, atom_types
|
|
write(11,21)
|
|
do i = 1, atom_types
|
|
write(11,*) i, type_to_name(i)
|
|
end do
|
|
write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3)
|
|
write(11,18)
|
|
write(11,'(2i16)') 1,1 !This is another throwaway line that is meaningless
|
|
write(11,17)
|
|
!This may have to be updated in the future but currently the only 8 node element is fcc
|
|
do i = 1, lattice_types
|
|
select case(ng_node(i))
|
|
case(8)
|
|
write(11, *) i, 'fcc'
|
|
end select
|
|
end do
|
|
write(11,15) 1.0_dp, 1.0_dp, 1.0_dp !Another throwaway line that isn't needed
|
|
write(11,8) box_bd
|
|
write(11,9) box_bd(2)-box_bd(1), box_bd(4) - box_bd(3), box_bd(6)-box_bd(5)
|
|
write(11,10)
|
|
!Current boxes are limited to being rectangular
|
|
do i = 1,3
|
|
box_vec(:) = 0.0_dp
|
|
box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
|
|
write(11,11) box_vec
|
|
end do
|
|
!We write this as box_mat ori and box_mat current
|
|
do i = 1,3
|
|
box_vec(:) = 0.0_dp
|
|
box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
|
|
write(11,11) box_vec
|
|
end do
|
|
|
|
!write the element information
|
|
if(ele_num > 0) then
|
|
write(11,12)
|
|
|
|
do i = 1, unique_num
|
|
write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i)))
|
|
end do
|
|
ip = 0
|
|
write(11,13)
|
|
do i = 1, ele_num
|
|
!Figure out the ele type
|
|
do j = 1, unique_num
|
|
if ( unique_size(j) == size_ele(i)) then
|
|
etype = j
|
|
exit
|
|
endif
|
|
end do
|
|
write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i))
|
|
do inod = 1, ng_node(lat_ele(i))
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
ip = ip + 1
|
|
write(11, '(2i16, 3f23.15)') ip, ibasis, r_node(:, ibasis, inod, i)
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
|
|
!Now write the atomic information
|
|
if(atom_num /= 0) then
|
|
write(11,14)
|
|
do i = 1, atom_num
|
|
write(11, '(3i16, 3f23.15)') i, 1, type_atom(i), r_atom(:,i)
|
|
end do
|
|
end if
|
|
|
|
close(11)
|
|
end subroutine write_pycac
|
|
|
|
subroutine write_mb(file)
|
|
|
|
!This subroutine writes the cacmb formatted file which provides necessary information for building models
|
|
character(len=100), intent(in) :: file
|
|
|
|
integer :: i, j, k, inod, ibasis
|
|
|
|
!Open the .mb file for writing
|
|
open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
|
|
|
|
!First write the box boundary information
|
|
!Write the global box boundaries
|
|
write(11,*) box_bd(:)
|
|
!Write the number of sub_boxes in the system
|
|
write(11,*) sub_box_num
|
|
!For every subbox write the orientation, sub box boundar
|
|
do i = 1, sub_box_num
|
|
write(11,*) sub_box_ori(:,:,i)
|
|
write(11,*) sub_box_bd(:,i)
|
|
end do
|
|
|
|
!Write the number of atom types in the current model and all of their names
|
|
write(11,*) atom_types, (type_to_name(i)//' ', i=1, atom_types)
|
|
!Write the number of lattice_types, basisnum and number of nodes for each lattice type
|
|
write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
|
|
!Now for every lattice type write the basis atom types
|
|
write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)
|
|
!Now for every lattice type write the lattice parameters
|
|
write(11,*) (lapa(i), i = 1, lattice_types)
|
|
|
|
!Now write the numbers of elements and atoms
|
|
write(11,*) atom_num, ele_num
|
|
|
|
!Write out atoms first
|
|
do i = 1, atom_num
|
|
write(11,*) tag_atom(i), type_atom(i), sbox_atom(i), r_atom(:,i)
|
|
end do
|
|
|
|
!Write out the elements, this is written in two stages, one line for the element and then 1 line for
|
|
!every basis at every node
|
|
do i = 1, ele_num
|
|
write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i)
|
|
do inod = 1, ng_node(lat_ele(i))
|
|
do ibasis =1, basisnum(lat_ele(i))
|
|
write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
close(11)
|
|
end subroutine write_mb
|
|
|
|
|
|
!!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine get_in_file(filename)
|
|
|
|
implicit none
|
|
|
|
character(len=100), intent(in) :: filename
|
|
character(len=100) :: temp_infile
|
|
logical :: file_exists
|
|
|
|
!If no filename is provided then this function is called with none and prompts user input
|
|
if (filename=='none') then
|
|
print *, "Please specify a filename or extension to output to:"
|
|
read(*,*) temp_infile
|
|
else
|
|
temp_infile = filename
|
|
end if
|
|
|
|
!Infinite loop which only exists if user provides valid filetype
|
|
do while(.true.)
|
|
|
|
!Check to see if file exists, if it does then ask user if they would like to overwrite the file
|
|
inquire(file=trim(temp_infile), exist=file_exists)
|
|
if (.not.file_exists) then
|
|
print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists"
|
|
read(*,*) temp_infile
|
|
cycle
|
|
end if
|
|
|
|
select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
|
|
case('restart', 'mb', 'cac')
|
|
infilenum=infilenum+1
|
|
infiles(infilenum) = temp_infile
|
|
exit
|
|
case default
|
|
print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", &
|
|
"please input a filename with extension from following list: mb, restart."
|
|
read(*,*) temp_infile
|
|
|
|
end select
|
|
end do
|
|
|
|
end subroutine get_in_file
|
|
|
|
subroutine read_in(i, displace, temp_box_bd)
|
|
!This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
|
|
|
|
integer, intent(in) :: i
|
|
real(kind=dp), dimension(3), intent(in) :: displace
|
|
real(kind=dp), dimension(6), intent(out) :: temp_box_bd
|
|
|
|
!Pull out the extension of the file and call the correct write subroutine
|
|
select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:))))
|
|
case('mb')
|
|
call read_mb(infiles(i), displace, temp_box_bd)
|
|
case('restart')
|
|
call read_pycac(infiles(i), displace, temp_box_bd)
|
|
case('cac')
|
|
call read_lmpcac(infiles(i), displace, temp_box_bd)
|
|
case default
|
|
print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
|
|
" is not accepted for writing. Please select from: mb and try again"
|
|
stop
|
|
|
|
end select
|
|
|
|
end subroutine read_in
|
|
|
|
subroutine read_mb(file, displace, temp_box_bd)
|
|
!This subroutine reads in an mb file for operation
|
|
|
|
character(len=100), intent(in) :: file
|
|
real(kind=dp), dimension(3), intent(in) :: displace
|
|
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
|
|
|
|
integer :: i, j, k, l, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
|
|
new_type_to_type(10), new_lattice_types, sbox, new_lattice_map(10)
|
|
character(len=100) :: etype
|
|
real(kind=dp) :: r(3), newdisplace(3)
|
|
real(kind=dp), allocatable :: r_innode(:,:,:)
|
|
character(len = 2) :: new_type_to_name(10)
|
|
!First open the file
|
|
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
|
|
|
|
!Read in the box boundary and grow the current active box bd
|
|
read(11, *) temp_box_bd(:)
|
|
|
|
do i = 1, 3
|
|
if (abs(displace(i)) > lim_zero) then
|
|
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
|
|
else
|
|
newdisplace(i)=displace(i)
|
|
end if
|
|
temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
|
|
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
|
|
end do
|
|
|
|
call grow_box(temp_box_bd)
|
|
!Read in the number of sub_boxes and allocate the variables
|
|
read(11, *) n
|
|
|
|
if (sub_box_num == 0) then
|
|
call alloc_sub_box(n)
|
|
else
|
|
call grow_sub_box(n)
|
|
end if
|
|
|
|
!Read in subbox orientations and boundaries
|
|
do i = 1, n
|
|
!Read in orientation with column major ordering
|
|
read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3)
|
|
!Read in subbox boundaries
|
|
read(11,*) sub_box_bd(:,sub_box_num+i)
|
|
|
|
do j = 1, 3
|
|
sub_box_bd(2*j-1,sub_box_num+i) = sub_box_bd(2*j-1, sub_box_num+i) + displace(j)
|
|
sub_box_bd(2*j,sub_box_num+i) = sub_box_bd(2*j, sub_box_num+i) + displace(j)
|
|
end do
|
|
end do
|
|
|
|
!Read in the number of atom types and all their names
|
|
read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
|
|
!Now fit these into the global list of atom types, after this new_type_to_type is the actual global
|
|
!type of the atoms within this file
|
|
do i = 1, new_atom_types
|
|
call add_atom_type(new_type_to_name(i), new_type_to_type(i))
|
|
end do
|
|
!Read the number of lattice types, basisnum, and number of nodes for each lattice type
|
|
read(11,*) new_lattice_types, (basisnum(i), i = lattice_types+1, lattice_types+new_lattice_types), &
|
|
(ng_node(i), i = lattice_types+1, lattice_types+new_lattice_types)
|
|
!Define max_ng_node and max_basis_num
|
|
max_basisnum = maxval(basisnum)
|
|
max_ng_node = maxval(ng_node)
|
|
!Read the basis atom types for every lattice
|
|
read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = lattice_types+1, lattice_types+new_lattice_types)
|
|
!Convert the basis_atom types
|
|
do j = lattice_types+1, lattice_types+new_lattice_types
|
|
do i = 1, basisnum(j)
|
|
basis_type(i,j) = new_type_to_type(basis_type(i,j))
|
|
end do
|
|
end do
|
|
!Read the lattice parameters for every lattice type
|
|
read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types)
|
|
|
|
!Now we loop over all new lattice types and check to see if they are exactly the same as any old lattice types
|
|
!If they are then we map the new type to the old type.
|
|
k = lattice_types + 1
|
|
new_lattice_map(:) = 0
|
|
new_loop:do i = 1, new_lattice_types
|
|
old_loop:do j = 1, lattice_types
|
|
!First check all the lattice level variables
|
|
if ((basisnum(lattice_types + i) == basisnum(j)).and. &
|
|
(ng_node(lattice_types + i) == ng_node(j))) then
|
|
!Now check the basis level variables
|
|
do ibasis =1, basisnum(i)
|
|
if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then
|
|
cycle old_loop
|
|
end if
|
|
end do
|
|
new_lattice_map(i) = j
|
|
cycle new_loop
|
|
end if
|
|
end do old_loop
|
|
new_lattice_map(i) = k
|
|
k = k+1
|
|
end do new_loop
|
|
|
|
!Read number of elements and atoms and allocate arrays
|
|
read(11, *) in_atoms, in_eles
|
|
call grow_ele_arrays(in_eles, in_atoms)
|
|
allocate(r_innode(3,max_basisnum, max_ng_node))
|
|
|
|
print *, "Read in ", in_eles, " elements and ", in_atoms, " atoms from ", trim(adjustl(file))
|
|
print *, "New box dimensions are: ", box_bd
|
|
|
|
!Read the atoms
|
|
do i = 1, in_atoms
|
|
read(11,*) j, type, sbox, r(:)
|
|
r = r+newdisplace
|
|
call add_atom(j, new_type_to_type(type), sbox+sub_box_num, r)
|
|
end do
|
|
|
|
!Read the elements
|
|
do i = 1, in_eles
|
|
read(11, *) j, type, size, sbox, etype
|
|
do inod = 1, ng_node(type)
|
|
do ibasis =1, basisnum(type)
|
|
read(11,*) k, l, r_innode(:, ibasis, inod)
|
|
r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
|
|
end do
|
|
end do
|
|
call add_element(j, etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode)
|
|
end do
|
|
|
|
!Close the file being read
|
|
close(11)
|
|
|
|
!Only increment the lattice types if there are elements, if there are no elements then we
|
|
!just overwrite the arrays
|
|
lattice_types = maxval(new_lattice_map)
|
|
|
|
|
|
sub_box_num = sub_box_num + n
|
|
|
|
call set_max_esize
|
|
end subroutine read_mb
|
|
|
|
subroutine read_pycac(file, displace, temp_box_bd)
|
|
!Add subroutine for reading in restart files from PyCAC. This code currently only
|
|
!works for 8 node elements.
|
|
character(len=100), intent(in) :: file
|
|
real(kind=dp), dimension(3), intent(in) :: displace
|
|
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
|
|
|
|
integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
|
|
atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), &
|
|
atom_type
|
|
real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3)
|
|
character(len=100) :: textholder, in_lattype_map(10)
|
|
character(len=2) :: atomic_element
|
|
!First open the file
|
|
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
|
|
|
|
!Disregard unneeded information
|
|
do i = 1, 3
|
|
read(11,*) textholder
|
|
end do
|
|
|
|
!Read element number
|
|
read(11,*) in_eles
|
|
|
|
!Discard info and read ng_max_node
|
|
do i = 1, 5
|
|
read(11,*) textholder
|
|
end do
|
|
read(11,*) max_ng_node
|
|
|
|
!Read element types (only needed inside this subroutine)
|
|
read(11,*) textholder
|
|
read(11,*) ele_types
|
|
|
|
!Read in atom num
|
|
read(11,*) textholder
|
|
read(11,*) in_atoms
|
|
|
|
!read in lattice_types and atom types
|
|
do i = 1,3
|
|
read(11,*) textholder
|
|
end do
|
|
read(11,*) in_lat_num, in_atom_types
|
|
|
|
|
|
!Read define atom_types by name
|
|
read(11,*) textholder
|
|
do i = 1, in_atom_types
|
|
read(11,*) j, atomic_element
|
|
call add_atom_type(atomic_element, atom_type_map(i))
|
|
end do
|
|
|
|
!Read in the boundary
|
|
read(11,*) textholder
|
|
read(11,*) box_bc
|
|
|
|
!Disregard useless info
|
|
do i = 1, 3
|
|
read(11,*) textholder
|
|
end do
|
|
|
|
!Read in lat_type map
|
|
do i = 1, in_lat_num
|
|
read(11,*) j, in_lattype_map(i)
|
|
ng_node(lattice_types+i) = 8 !Only cubic elements are currently supported in pyCAC
|
|
end do
|
|
|
|
!Disregard useless info
|
|
do i =1 , 3
|
|
read(11,*) textholder
|
|
end do
|
|
|
|
!Read box boundaries and displace them if necessary
|
|
read(11,*) temp_box_bd(:)
|
|
do i = 1, 3
|
|
if (abs(displace(i)) > lim_zero) then
|
|
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
|
|
else
|
|
newdisplace(i)=displace(i)
|
|
end if
|
|
temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
|
|
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
|
|
end do
|
|
|
|
call grow_box(temp_box_bd)
|
|
|
|
!Allocate sub_box
|
|
if (sub_box_num == 0) then
|
|
call alloc_sub_box(1)
|
|
else
|
|
call grow_sub_box(1)
|
|
end if
|
|
|
|
!Because orientations and other needed sub_box information isn't really
|
|
!present within the restart file we just default a lot of it.
|
|
sub_box_ori(:,:,sub_box_num+1) = identity_mat(3)
|
|
sub_box_bd(:, sub_box_num+1) = temp_box_bd
|
|
|
|
!Read in more useless info
|
|
do i = 1, 10
|
|
read(11,*) textholder
|
|
end do
|
|
|
|
!Now read the ele_type to size and lat map
|
|
do i = 1, ele_types
|
|
read(11,*) j, etype_map(i)
|
|
etype_map(i) = etype_map(i) + 1
|
|
end do
|
|
|
|
|
|
!Now set up the lattice types. In this code it assumes that lattice_type 1 maps to
|
|
!atom type 1 because it only allows 1 atom per basis in pyCAC at the moment.
|
|
do i = 1, in_lat_num
|
|
basisnum(lattice_types+i) = 1
|
|
basis_type(1,lattice_types+i) = atom_type_map(i)
|
|
end do
|
|
|
|
!Figure out the lattice type maps in case we have repeated lattice_types
|
|
k = lattice_types + 1
|
|
new_lattice_map(:) = 0
|
|
new_loop:do i = 1, in_lat_num
|
|
old_loop:do j = 1, lattice_types
|
|
!First check all the lattice level variables
|
|
if ((basisnum(lattice_types+i) == basisnum(j)).and. &
|
|
(ng_node(lattice_types+i) == ng_node(j))) then
|
|
!Now check the basis level variables
|
|
do ibasis =1, basisnum(i)
|
|
if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then
|
|
cycle old_loop
|
|
end if
|
|
end do
|
|
new_lattice_map(i) = j
|
|
cycle new_loop
|
|
end if
|
|
end do old_loop
|
|
new_lattice_map(i) = k
|
|
k = k+1
|
|
end do new_loop
|
|
|
|
!Read more useless data
|
|
read(11,*) textholder
|
|
|
|
!set max values and allocate variables
|
|
max_basisnum = maxval(basisnum)
|
|
max_ng_node = maxval(ng_node)
|
|
call grow_ele_arrays(in_eles, in_atoms)
|
|
|
|
!Now start reading the elements
|
|
if(in_eles > 0) then
|
|
read(11,*) textholder
|
|
do i = 1, in_eles
|
|
read(11,*) j, etype, k, lat_type
|
|
do inod = 1, 8
|
|
read(11, *) k, l, r_in(:,1,inod)
|
|
r_in(:,1,inod) = r_in(:,1,inod) + newdisplace
|
|
end do
|
|
call add_element(j, in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in)
|
|
end do
|
|
end if
|
|
|
|
if(in_atoms > 0) then
|
|
|
|
if (in_eles > 0) then
|
|
!Read useless data
|
|
read(11,*) textholder
|
|
read(11,*) textholder
|
|
end if
|
|
|
|
do i = 1, in_atoms
|
|
read(11,*) j, k, atom_type, r_in_atom(:)
|
|
r_in_atom = r_in_atom + newdisplace
|
|
call add_atom(j,atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
|
|
end do
|
|
!Close file
|
|
close(11)
|
|
|
|
lattice_types = maxval(new_lattice_map)
|
|
|
|
sub_box_num = sub_box_num + 1
|
|
|
|
call set_max_esize
|
|
end if
|
|
end subroutine read_pycac
|
|
|
|
subroutine read_lmpcac(file, displace, temp_box_bd)
|
|
!This subroutine is used to read .cac files which are used with the lammpsCAC format
|
|
character(len=100), intent(in) :: file
|
|
real(kind=dp), dimension(3), intent(in) :: displace
|
|
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
|
|
|
|
character(len=100) :: textholder, element_type
|
|
character(len=2) :: atom_species
|
|
integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10), esize, &
|
|
lat_type
|
|
real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), newdisplace(3)
|
|
|
|
!First check to make sure that we have set the needed variables
|
|
if(is_equal(in_lapa,0.0_dp).or.(in_lattice_type=='')) then
|
|
print *, "Please use set_cac to set needed parameters to read in .cac file"
|
|
stop 3
|
|
end if
|
|
!Open the file
|
|
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
|
|
|
|
!Now initialiaze some important variables if they aren't defined
|
|
if (max_basisnum==0) max_basisnum = 10
|
|
if (max_ng_node==0) max_ng_node=8
|
|
|
|
!Read header information
|
|
read(11, *) textholder
|
|
|
|
!Read number of elements
|
|
read(11, *) ele_in, textholder
|
|
read(11, *) type_in, textholder
|
|
|
|
!Read box_boundaries
|
|
read(11,*) temp_box_bd(1:2), textholder
|
|
read(11,*) temp_box_bd(3:4), textholder
|
|
read(11,*) temp_box_bd(5:6), textholder
|
|
|
|
!Shift the box boundaries if needed
|
|
do i = 1, 3
|
|
if (abs(displace(i)) > lim_zero) then
|
|
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
|
|
else
|
|
newdisplace(i)=displace(i)
|
|
end if
|
|
temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
|
|
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
|
|
end do
|
|
|
|
!Grow box boundaries
|
|
call grow_box(temp_box_bd)
|
|
|
|
!Allocate sub_box
|
|
if (sub_box_num == 0) then
|
|
call alloc_sub_box(1)
|
|
else
|
|
call grow_sub_box(1)
|
|
end if
|
|
|
|
!Because orientations and other needed sub_box information isn't really
|
|
!present within the .cac file we just default a lot of it.
|
|
sub_box_ori(:,:,sub_box_num+1) = identity_mat(3)
|
|
sub_box_bd(:, sub_box_num+1) = temp_box_bd
|
|
sub_box_num = sub_box_num + 1
|
|
|
|
|
|
!Read useless information
|
|
read(11,*) textholder
|
|
|
|
!Read atomic masses
|
|
do i = 1, type_in
|
|
read(11,*) j, mass, textholder
|
|
call ATOMMASSSPECIES(mass, atom_species)
|
|
call add_atom_type(atom_species, type_map(i))
|
|
end do
|
|
|
|
!Read useless info
|
|
read(11,*) textholder
|
|
|
|
!Start the reading loop
|
|
do i = 1, ele_in
|
|
read(11,*) j, element_type, in_basis, esize
|
|
!Check to see if we need to grow the max_basis_num
|
|
select case(trim(adjustl(element_type)))
|
|
case('Eight_Node')
|
|
!Read in all the data
|
|
do j = 1, 8*in_basis
|
|
read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,ibasis,inod)
|
|
end do
|
|
|
|
!Now calculate the lattice vectors and shift the nodal points from the corners to the center of the unit cell
|
|
!Please check the nodal numbering figure in the readme in order to understand which nodes are used for the
|
|
!calculation
|
|
lat_vec(:,1) = (r_in(:,1,2) - r_in(:,1,1))/(2*esize)
|
|
lat_vec(:,2) = (r_in(:,1,4) - r_in(:,1,1))/(2*esize)
|
|
lat_vec(:,3) = (r_in(:,1,5) - r_in(:,1,1))/(2*esize)
|
|
|
|
!Now shift all the nodal positions
|
|
select case(trim(adjustl(in_lattice_type)))
|
|
case('fcc','FCC')
|
|
do ibasis = 1, in_basis
|
|
r_in(:,ibasis,1) = r_in(:,ibasis,1) + lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,2) = r_in(:,ibasis,2) - lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,3) = r_in(:,ibasis,3) - lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,4) = r_in(:,ibasis,4) + lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,5) = r_in(:,ibasis,5) + lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,6) = r_in(:,ibasis,6) - lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,7) = r_in(:,ibasis,7) - lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace
|
|
r_in(:,ibasis,8) = r_in(:,ibasis,8) + lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace
|
|
end do
|
|
case default
|
|
print *, in_lattice_type, " is not an accepted lattice type. Please select from: fcc"
|
|
end select
|
|
!Now map it to either an existing or new lattice type
|
|
call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type)
|
|
|
|
!Now add the element
|
|
call add_element(0,in_lattice_type, esize, lat_type, sub_box_num, r_in(:,1:max_basisnum,1:max_ng_node))
|
|
|
|
case('Atom')
|
|
read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1)
|
|
call add_atom(0,in_basis_types(ibasis), sub_box_num, r_in(:,1,1))
|
|
end select
|
|
end do
|
|
|
|
end subroutine read_lmpcac
|
|
|
|
subroutine set_cac(apos)
|
|
!This code parses input values
|
|
integer, intent(in) :: apos
|
|
integer :: arglen, arg_pos
|
|
|
|
character(len=100) :: textholder
|
|
|
|
arg_pos = apos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) then
|
|
print *, "Missing lattice parameter for set_input_lat"
|
|
end if
|
|
read(textholder,*) in_lapa
|
|
print *, in_lapa
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) then
|
|
print *, "Missing lattice type for set_input_lat"
|
|
end if
|
|
read(textholder,*) in_lattice_type
|
|
print *, in_lattice_type
|
|
|
|
end subroutine set_cac
|
|
end module io
|