From 4038cab76fe4630e39710819423ce6ec11365302 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 23 Apr 2020 12:17:04 -0400 Subject: [PATCH] Latest working changes to read-CAC --- src/elements.f90 | 2 - src/io.f90 | 97 ++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 77 insertions(+), 22 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 6265b3c..484057f 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -294,8 +294,6 @@ module elements integer :: i - max_ng_node = 0 - do i=1, n select case(trim(adjustl(element_types(i)))) case('fcc') diff --git a/src/io.f90 b/src/io.f90 index 22692e8..964a816 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -776,7 +776,7 @@ module io integer :: i, inod, ibasis, j, k, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, & atom_type_map(10), etype_map(10), etype, lat_type, new_lattice_map(10), & atom_type - real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), new_displace(3) + 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 @@ -950,7 +950,7 @@ module io end if end subroutine read_pycac - subroutine read_cac(file, displace, temp_box_bd) + 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 @@ -958,8 +958,8 @@ module io character(len=100) :: textholder, element_type, esize character(len=2) :: atom_species - integer :: i, j, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10) - real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3), in_ori(3,3) + integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10) + real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), temp_box_bd(6), 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 @@ -968,6 +968,9 @@ module io end if !Open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') + + !Now initialiaze some important variables + max_basis_num = 10 !Read header information read(11, *) textholder @@ -979,9 +982,38 @@ module io !Read box_boundaries read(11,*) textholder - read(11,*) box_bd(1:2), texholder - read(11,*) box_bd(3:4), texholder - read(11,*) box_bd(5:6), texholder + read(11,*) temp_box_bd(1:2), texholder + read(11,*) temp_box_bd(3:4), texholder + read(11,*) temp_box_bd(5:6), texholder + + !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(11,*) textholder @@ -1001,25 +1033,50 @@ module io !Start the reading loop do i = 1, ele_in read(11,*) j, ele, 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) + !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_node(:,1,2) - r_node(:,1,1))/esize + lat_vec(:,2) = (r_node(:,1,4) - r_node(:,1,1))/esize + lat_vec(:,3) = (r_node(:,1,5) - r_node(:,1,1))/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 - - !Now calculate the orientation matrix based on the lattice type. - lat_vec = r_in(:,1,2) - r_in(:,1,1) - lat_vec = lat_vec / norm2(lat_vec) - - !Now figure out if is an existing lattice_type - call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type) - - + 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(in_lattice_type, esize, lat_type, sub_box_num, r_in) + + case('Atom') + read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1) + call add_atom(in_basis_types(ibasis), sub_box_num, r_in(:,1,1)) end select end do - end subroutine read_cac + end subroutine read_lmpcac subroutine set_cac(apos) !This code parses input values