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

Merge pull request #7 from aselimov/redo-lammpscac

Browse Source 
Changes to how the adjustment to nodal positions is performed for lam…
master
aselimov 5 years ago committed by GitHub
parent 55fbe679e5 849da1d24a
commit fe6101167f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 69 additions and 54 deletions
Show Stats Download Patch File Download Diff File

19
src/elements.f90
Unescape View File

@ -9,7 +9,7 @@ module elements
!Data structures used to represent the CAC elements. Each index represents an element
character(len=100), allocatable :: type_ele(:) !Element type
integer, allocatable :: size_ele(:), lat_ele(:) !Element siz
integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:) !Element size
real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
integer, save :: ele_num !Number of elements
@ -39,6 +39,7 @@ module elements
!This can be easily increased with no change to efficiency
integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type
integer :: basis_type(10,10)
real(kind=dp) :: lapa(10)
public
contains
@ -140,7 +141,7 @@ module elements
!Allocate element arrays
if(n > 0) then
allocate(type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
allocate(type_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
stat=allostat)
if(allostat > 0) then
print *, "Error allocating element arrays in elements.f90 because of: ", allostat
@ -179,12 +180,17 @@ module elements
allocate(temp_int(n+ele_num+buffer_size))
temp_int(1:ele_size) = lat_ele
temp_int(ele_size+1:) = 0
call move_alloc(temp_int(1:ele_size), lat_ele)
call move_alloc(temp_int, lat_ele)
allocate(temp_int(n+ele_num+buffer_size))
temp_int(1:ele_size) = size_ele
temp_int(ele_size+1:) = 0
call move_alloc(temp_int(1:ele_size), size_ele)
call move_alloc(temp_int, size_ele)
allocate(temp_int(n+ele_num+buffer_size))
temp_int(1:ele_size) = lat_ele
temp_int(ele_size+1:) = 0
call move_alloc(temp_int, sbox_ele)
allocate(char_temp(n+ele_num+buffer_size))
char_temp(1:ele_size) = type_ele
@ -210,9 +216,9 @@ module elements
end if
end subroutine
subroutine add_element(type, size, lat, r)
subroutine add_element(type, size, lat, sbox, r)
!Subroutine which adds an element to the element arrays
integer, intent(in) :: size, lat
integer, intent(in) :: size, lat, sbox
character(len=100), intent(in) :: type
real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
@ -220,6 +226,7 @@ module elements
type_ele(ele_num) = type
size_ele(ele_num) = size
lat_ele(ele_num) = lat
sbox_ele(ele_num) = sbox
r_node(:,:,:,ele_num) = r(:,:,:)
node_num = node_num + ng_node(lat)

38
src/io.f90
Unescape View File

@ -204,7 +204,7 @@ module io
!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
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)
@ -242,16 +242,32 @@ module io
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
write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
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
@ -486,6 +502,8 @@ module io
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
@ -498,7 +516,7 @@ module io
!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, *) i, lat_ele(i), size_ele(i), type_ele(i)
write(11, *) 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)
@ -582,7 +600,7 @@ module io
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
new_type_to_type(10), new_lattice_types
new_type_to_type(10), new_lattice_types, sbox
character(len=100) :: etype
real(kind=dp) :: r(3), newdisplace(3)
real(kind=dp), allocatable :: r_innode(:,:,:)
@ -623,8 +641,6 @@ module io
sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num
sub_box_num = sub_box_num + n
!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
@ -646,6 +662,8 @@ module io
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)
!Read number of elements and atoms and allocate arrays
read(11, *) in_atoms, in_eles
call grow_ele_arrays(in_eles, in_atoms)
@ -659,7 +677,7 @@ module io
!Read the elements
do i = 1, in_eles
read(11, *) n, type, size, etype
read(11, *) n, type, size, sbox, etype
do inod = 1, ng_node(type)
do ibasis =1, basisnum(type)
read(11,*) j, k, r_innode(:, ibasis, inod)
@ -667,7 +685,7 @@ module io
end do
end do
type = type + lattice_types
call add_element(etype, size, type, r_innode)
call add_element(etype, size, type, sbox+n, r_innode)
end do
!Close the file being read
@ -676,5 +694,9 @@ module io
!Only increment the lattice types if there are elements, if there are no elements then we
!just overwrite the arrays
if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
sub_box_num = sub_box_num + n
end subroutine read_mb
end module io

11
src/main.f90
Unescape View File

@ -28,17 +28,6 @@ program main
! Command line parsing
arg_num = command_argument_count()
!First check if we are writing out to lammpscac format by looping over all arguments
do i = 1, arg_num
call get_command_argument(i, argument)
select case(argument(scan(argument,'.',.true.)+1:))
case('cac')
lmpcac = .true.
case default
continue
end select
end do
!Determine if a mode is being used and what it is. The first argument has to be the mode
!if a mode is being used
call get_command_argument(1, argument)

49
src/mode_create.f90
Unescape View File

@ -12,7 +12,7 @@ module mode_create
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), duplicate(3), adjustVar(3,8)
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
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
@ -105,7 +105,7 @@ module mode_create
box_bd(2*i) = maxval(r_node_temp(i,:,:))
box_bd(2*i-1) = origin(i)
end do
call add_element(element_type, esize, 1, r_node_temp)
call add_element(element_type, esize, 1, 1, r_node_temp)
end if
!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
@ -141,7 +141,7 @@ module mode_create
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)
call add_element(element_type, esize, 1, 1, r_node_temp)
end do
end if
end if
@ -199,7 +199,20 @@ module mode_create
call get_command_argument(arg_pos, orient_string, arglen)
if (arglen==0) STOP "Missing orientation in orient command of mode create"
arg_pos = arg_pos+1
call parse_ori_vec(orient_string, orient(i,:))
ori_pos=2
do j = 1,3
if (orient_string(ori_pos:ori_pos) == '-') then
ori_pos = ori_pos + 1
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
if (stat>0) STOP "Error reading orient value"
orient(i,j) = -orient(i,j)
ori_pos = ori_pos + 1
else
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
if(stat>0) STOP "Error reading orient value"
ori_pos=ori_pos + 1
end if
end do
end do
@ -234,6 +247,7 @@ module mode_create
exit
end select
end do
!Calculate the lattice periodicity length in lattice units
do i = 1, 3
lattice_space(i) = norm2(orient(i,:))
@ -264,8 +278,9 @@ module mode_create
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Set lattice_num to 1
!Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable
lattice_types = 1
lapa(1) = lattice_parameter
!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
@ -286,7 +301,7 @@ module mode_create
!Internal variables
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) :: v(3), temp_nodes(3,1,8)
real(kind=dp), allocatable :: resize_lat_array(:,:)
logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8)
@ -294,23 +309,6 @@ module mode_create
!Do some value initialization
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.5_dp
else
adjustVar(i, inod) = 0.5_dp
end if
end do
end do
adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8)))
else
adjustVar(:,:)=0.0_dp
end if
!First find the bounding lattice points (min and max points for the box in each dimension)
numlatpoints = 1
do i = 1, 3
@ -328,6 +326,7 @@ module mode_create
continue
end select
!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
if (esize==2) then
@ -441,9 +440,7 @@ module mode_create
if(all(node_in_bd)) then
lat_ele_num = lat_ele_num+1
do inod = 1, 8
r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
end do
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
!Now set all the lattice points contained within an element to false
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))

Write Preview
Loading…
Cancel
Save