From c1150ae4b13b1418d9419452b455d41d028310fb Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 29 Nov 2019 22:23:22 -0500 Subject: [PATCH] First attempts at figuring out the element format for lammpscac --- src/io.f90 | 14 ++++++++++---- src/mode_create.f90 | 46 ++++++++++++++++++++++++++++++++++++++------- 2 files changed, 49 insertions(+), 11 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 948fd52..af66da3 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -8,6 +8,7 @@ module io integer :: outfilenum = 0, infilenum = 0 character(len=100) :: outfiles(10), infiles(10) + logical lmpcac public contains @@ -57,7 +58,12 @@ module io cycle end if 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 outfiles(outfilenum) = temp_outfile exit @@ -98,18 +104,18 @@ module io cycle 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." read(*,*) temp_infile cycle end if - select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:)) + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) case('cac') infilenum=infilenum+1 infiles(infilenum) = temp_infile exit 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." read(*,*) temp_infile diff --git a/src/mode_create.f90 b/src/mode_create.f90 index e822798..53880c9 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -11,7 +11,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) + 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) logical :: dup_flag, dim_flag @@ -91,14 +91,30 @@ module mode_create continue 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 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) = 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 do i = 1,3 @@ -276,7 +292,7 @@ module mode_create !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) + 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(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) @@ -284,6 +300,21 @@ 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.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) numlatpoints = 1 do i = 1, 3 @@ -301,7 +332,6 @@ 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 @@ -415,8 +445,10 @@ module mode_create if(all(node_in_bd)) then 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 do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))