aselimov/CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

First attempts at figuring out the element format for lammpscac

Browse Source
This commit is contained in:
Alex Selimov 2019-11-29 22:23:22 -05:00
parent 54aa50b605
commit c1150ae4b1
2 changed files with 49 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/io.f90
View File

@ -8,6 +8,7 @@ module io
integer :: outfilenum = 0, infilenum = 0 integer :: outfilenum = 0, infilenum = 0
character(len=100) :: outfiles(10), infiles(10) character(len=100) :: outfiles(10), infiles(10)
logical lmpcac
public public
contains contains
@ -57,7 +58,12 @@ module io
cycle cycle
end if end if
select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
case('xyz','lmp','vtk','cac') case('xyz','lmp','vtk')
outfilenum=outfilenum+1
outfiles(outfilenum) = temp_outfile
exit
case('cac')
lmpcac = .true.
outfilenum=outfilenum+1 outfilenum=outfilenum+1
outfiles(outfilenum) = temp_outfile outfiles(outfilenum) = temp_outfile
exit exit
@ -98,18 +104,18 @@ module io
cycle cycle
end if end if
if (scan(temp_outfile,'.',.true.) == 0) then if (scan(temp_infile,'.',.true.) == 0) then
print *, "No extension included on filename, please type a full filename that includes an extension." print *, "No extension included on filename, please type a full filename that includes an extension."
read(*,*) temp_infile read(*,*) temp_infile
cycle cycle
end if end if
select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
case('cac') case('cac')
infilenum=infilenum+1 infilenum=infilenum+1
infiles(infilenum) = temp_infile infiles(infilenum) = temp_infile
exit exit
case default case default
print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", & print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), " not currently accepted. ", &
"please input a filename with extension from following list: cac." "please input a filename with extension from following list: cac."
read(*,*) temp_infile read(*,*) temp_infile

46
src/mode_create.f90
View File

@ -11,7 +11,7 @@ module mode_create
character(len=100) :: name, element_type character(len=100) :: name, element_type
real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3) orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), adjustVar(3,8)
integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6) integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6)
logical :: dup_flag, dim_flag logical :: dup_flag, dim_flag
@ -91,14 +91,30 @@ module mode_create
continue continue
else else
call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) if(lmpcac) then
do inod = 1, 8
do i = 1,3
if(is_equal(cubic_cell(i, inod),0.0_dp)) then
adjustVar(i,inod) = -0.5_dp*esize
else
adjustVar(i, inod) = 0.5_dp*esize
end if
end do
end do
else
adjustVar(:,:)=0.0_dp
end if
! call cell_init(lattice_parameter, esize, element_type, orient, cell_mat)
!If the user doesn't pass any build instructions than we just put the cell mat into the element_array !If the user doesn't pass any build instructions than we just put the cell mat into the element_array
call alloc_ele_arrays(1,0) call alloc_ele_arrays(1,0)
!Add the basis atoms to the unit cell !Add the basis atoms to the unit cell
do inod = 1, max_ng_node do inod = 1, max_ng_node
do ibasis = 1, basisnum(1) do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:) r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, &
matmul(fcc_mat, esize*(cubic_cell(:,inod)+adjustVar(:,inod)))) &
+ basis_pos(:,ibasis,1)
end do end do
end do end do
do i = 1,3 do i = 1,3
@ -276,7 +292,7 @@ module mode_create
!Internal variables !Internal variables
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, & integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, &
type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o
real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3) real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3), adjustVar(3,8)
real(kind=dp), allocatable :: resize_lat_array(:,:) real(kind=dp), allocatable :: resize_lat_array(:,:)
logical, allocatable :: lat_points(:,:,:) logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8) logical :: node_in_bd(8)
@ -284,6 +300,21 @@ module mode_create
!Do some value initialization !Do some value initialization
max_esize = esize max_esize = esize
!If we are writing out to lammpscac format then we have to adjust the nodal positions
if(lmpcac) then
do inod = 1, 8
do i = 1,3
if(is_equal(cubic_cell(i, inod),0.0_dp)) then
adjustVar(i,inod) = -0.25_dp
else
adjustVar(i, inod) = 0.25_dp
end if
end do
end do
else
adjustVar(:,:)=0.0_dp
end if
!First find the bounding lattice points (min and max points for the box in each dimension) !First find the bounding lattice points (min and max points for the box in each dimension)
numlatpoints = 1 numlatpoints = 1
do i = 1, 3 do i = 1, 3
@ -301,7 +332,6 @@ module mode_create
continue continue
end select end select
!Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case
!and one for the regular case !and one for the regular case
if (esize==2) then if (esize==2) then
@ -415,8 +445,10 @@ module mode_create
if(all(node_in_bd)) then if(all(node_in_bd)) then
lat_ele_num = lat_ele_num+1 lat_ele_num = lat_ele_num+1
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) do inod = 1, 8
r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
end do
!Now set all the lattice points contained within an element to false !Now set all the lattice points contained within an element to false
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))