aselimov/cacmb
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'development' into ft-write-lammpscac

Browse source
This commit is contained in:
Alex Selimov 2020-01-13 17:21:44 -05:00
commit d09ebfa7e0
10 changed files with 750 additions and 120 deletions
Download patch file Download diff file Expand all files Collapse all files

137
src/mode_create.f90
View file

@ -6,24 +6,27 @@ module mode_create
use io
use subroutines
use elements
use box
implicit none
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), &
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)
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), adjustVar(3,8)
integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
basis_pos(3,10)
logical :: dup_flag, dim_flag
real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
public
contains
subroutine create()
subroutine create(arg_pos)
! Main subroutine which controls execution
character(len=100) :: textholder
integer, intent(out) :: arg_pos
integer :: i, ibasis, inod
real(kind=dp), allocatable :: r_node_temp(:,:,:)
@ -43,29 +46,22 @@ module mode_create
lat_atom_num = 0
!First we parse the command
call parse_command()
call parse_command(arg_pos)
! Before building do a check on the file
if (outfilenum == 0) then
textholder = 'none'
call get_out_file(textholder)
end if
!Now we setup the unit element and call other init subroutines
call def_ng_node(1, element_type)
allocate(r_node_temp(3,max_basisnum,max_ng_node))
!Get the inverse orientation matrix we will need later
call matrix_inverse(orient,3,orient_inv)
if(dup_flag) then
!We initialize the cell with a lattice_parameter of 1 because we will add the lattice parameter later
call cell_init(1.0_dp, esize, element_type, orient, cell_mat)
!Define box vertices
do i = 1, 8
box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:)
box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter)
end do
call matrix_inverse(orient,3,orient_inv)
!Now get the rotated box vertex positions in lattice space. Should be integer units
box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
!Find the new maxlen
@ -74,21 +70,25 @@ module mode_create
box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
end do
!and then call the build function with the correct transformation matrix
select case(trim(adjustl(element_type)))
case('fcc')
call build_with_rhomb(box_lat_vert, fcc_mat)
case default
print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
"element type"
stop 3
end select
!Now that it is multiply by the lattice parameter
box_bd = box_bd*lattice_parameter
else if(dim_flag) then
continue
!As a note everything is defined so that the lattice parameter is multiplied in at the end
!so we have to divide all the real Angstroms units by the lattice parameter
!Define box_vertices
do i = 1, 8
box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter
end do
!Now get the rotated box vertex positions in lattice space. Should be integer units
box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
!Now get the box_bd in lattice units
do i = 1, 3
box_bd(2*i) = (box_len(i)+origin(i))/lattice_parameter
box_bd(2*i-1) = origin(i)/lattice_parameter
end do
else
if(lmpcac) then
@ -105,16 +105,14 @@ module mode_create
adjustVar(:,:)=0.0_dp
end if
! call cell_init(lattice_parameter, esize, element_type, orient, cell_mat)
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
call alloc_ele_arrays(1,0)
!Add the basis atoms to the unit cell
do inod = 1, max_ng_node
do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, &
matmul(fcc_mat, (esize+1)*cubic_cell(:,inod)+adjustVar(:,inod))) &
+ basis_pos(:,ibasis,1)
r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:)
end do
end do
do i = 1,3
@ -126,12 +124,25 @@ module mode_create
!If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
if(dup_flag.or.dim_flag) then
!Call the build function with the correct transformation matrix
select case(trim(adjustl(element_type)))
case('fcc')
call build_with_rhomb(box_lat_vert, fcc_mat)
case default
print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
"element type"
stop 3
end select
!Now that it is built multiply by the lattice parameter
box_bd = box_bd*lattice_parameter
!Allocate variables
call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
if(lat_atom_num > 0) then
do i = 1, lat_atom_num
do ibasis = 1, basisnum(1)
call add_atom(basis_type(ibasis,1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis,1))
call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
end do
end do
deallocate(r_atom_lat)
@ -141,7 +152,7 @@ module mode_create
do i = 1, lat_ele_num
do inod= 1, ng_node(1)
do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis,1)
r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
end do
end do
call add_element(element_type, esize, 1, r_node_temp)
@ -149,12 +160,20 @@ module mode_create
end if
end if
!The last thing we do is setup the sub_box_boundaries
call alloc_sub_box(1)
sub_box_num = 1
sub_box_ori(:,:,1) = orient
sub_box_bd(:,1) = box_bd
sub_box_array_bd(1,:,1) = 1
sub_box_array_bd(2,1,1) = atom_num
sub_box_array_bd(2,2,1) = ele_num
end subroutine create
!This subroutine parses the command and pulls out information needed for mode_create
subroutine parse_command()
subroutine parse_command(arg_pos)
integer :: arg_pos, ori_pos, i, j, arglen, stat
integer, intent(out) :: arg_pos
integer :: ori_pos, i, j, arglen, stat
character(len=100) :: textholder
character(len=8) :: orient_string
@ -214,33 +233,31 @@ module mode_create
!If the duplicate command is passed then we extract the information on the new bounds.
case('duplicate')
if(dim_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
dup_flag = .true.
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) duplicate(i)
arg_pos = arg_pos + 1
end do
case('dim')
if(dup_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
dim_flag = .true.
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) box_len(i)
arg_pos = arg_pos + 1
end do
case('origin')
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) origin(i)
arg_pos = arg_pos + 1
end do
!If a filetype is passed then we add name.ext to the outfiles list
case('xyz')
textholder = trim(adjustl(name)) //'.xyz'
call get_out_file(textholder)
case default
!Check to see if it is an option command, if so then mode_create must be finished
if(textholder(1:1) == '-') then
exit
!Check to see if a filename was passed
elseif(scan(textholder,'.',.true.) > 0) then
call get_out_file(textholder)
end if
!If it isn't an option then you have to exit
exit
end select
end do
@ -274,13 +291,17 @@ module mode_create
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Set lattice_num to 1
lattice_types = 1
!If we haven't defined a basis then define the basis and add the default basis atom type and position
if (basisnum(1) == 0) then
max_basisnum = 1
basisnum(1) = 1
call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
basis_pos(:,1,1) = 0.0_dp
basis_pos(:,1) = 0.0_dp
end if
!
end subroutine
subroutine build_with_rhomb(box_in_lat, transform_matrix)
@ -290,9 +311,9 @@ module mode_create
integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
!Internal variables
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
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)
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), &
vlat(3), temp_lat(3,8), m, n, o
real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8)
real(kind=dp), allocatable :: resize_lat_array(:,:)
logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8)
@ -466,9 +487,9 @@ module mode_create
!Now figure out how many lattice points could not be contained in elements
allocate(r_atom_lat(3,count(lat_points)))
lat_atom_num = 0
do ix = 1, bd_in_array(3)
do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do iz = 1, bd_in_array(1)
do ix = 1, bd_in_array(1)
!If this point is a lattice point then save the lattice point as an atom
if (lat_points(ix,iy,iz)) then
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
@ -514,4 +535,4 @@ module mode_create
end subroutine error_message
end module mode_create
end module mode_create