From b52d7761e09795fc9fd050615d1edcb68a26dc75 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 2 Jun 2020 16:42:12 -0400 Subject: [PATCH] First working version of reading in .cac format. Works for two separate orientations --- src/call_option.f90 | 2 +- src/elements.f90 | 5 +++-- src/io.f90 | 54 +++++++++++++++++++++++++-------------------- src/main.f90 | 2 -- 4 files changed, 34 insertions(+), 29 deletions(-) diff --git a/src/call_option.f90 b/src/call_option.f90 index b8cc3b8..368e542 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -34,7 +34,7 @@ subroutine call_option(option, arg_pos) case('-delete') call run_delete(arg_pos) case('-set_cac') - arg_pos = arg_pos+3 + arg_pos=arg_pos +3 case default print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.' stop 3 diff --git a/src/elements.f90 b/src/elements.f90 index 484057f..6252427 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -646,7 +646,7 @@ module elements if ((basisnum(j) == in_bnum).and.(ng_node(j) == in_ngnodes).and.(is_equal(lapa(j),in_lapa))) then !Now check lattice level variables do ibasis = 1, basisnum(j) - if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle old_loop + if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle lat_loop end do lat_type = j exit lat_loop @@ -657,9 +657,10 @@ module elements if( lat_type == 0) then lattice_types = lattice_types + 1 basisnum(lattice_types) = in_bnum - basis_types(:,lattice_types) = in_btypes + basis_type(:,lattice_types) = in_btypes ng_node(lattice_types) = in_ngnodes lapa(lattice_types) = in_lapa + lat_type = lattice_types end if end subroutine lattice_map diff --git a/src/io.f90 b/src/io.f90 index 964a816..ab90f85 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -591,7 +591,7 @@ module io end if select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('restart', 'mb') + case('restart', 'mb', 'cac') infilenum=infilenum+1 infiles(infilenum) = temp_infile exit @@ -618,6 +618,8 @@ module io call read_mb(infiles(i), displace, temp_box_bd) case('restart') call read_pycac(infiles(i), displace, temp_box_bd) + case('cac') + call read_lmpcac(infiles(i), displace, temp_box_bd) case default print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), & " is not accepted for writing. Please select from: mb and try again" @@ -956,10 +958,11 @@ module io real(kind=dp), dimension(3), intent(in) :: displace real(kind = dp), dimension(6), intent(out) :: temp_box_bd - character(len=100) :: textholder, element_type, esize + character(len=100) :: textholder, element_type character(len=2) :: atom_species - 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) + integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10), esize, & + lat_type + real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), 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 @@ -969,22 +972,21 @@ module io !Open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') - !Now initialiaze some important variables - max_basis_num = 10 + !Now initialiaze some important variables if they aren't defined + if (max_basisnum==0) max_basisnum = 10 + if (max_ng_node==0) max_ng_node=8 !Read header information read(11, *) textholder - read(11, *) textholder !Read number of elements read(11, *) ele_in, textholder read(11, *) type_in, textholder !Read box_boundaries - read(11,*) textholder - read(11,*) temp_box_bd(1:2), texholder - read(11,*) temp_box_bd(3:4), texholder - read(11,*) temp_box_bd(5:6), texholder + read(11,*) temp_box_bd(1:2), textholder + read(11,*) temp_box_bd(3:4), textholder + read(11,*) temp_box_bd(5:6), textholder !Shift the box boundaries if needed do i = 1, 3 @@ -1016,7 +1018,6 @@ module io !Read useless information read(11,*) textholder - read(11,*) textholder !Read atomic masses do i = 1, type_in @@ -1026,13 +1027,11 @@ module io end do !Read useless info - do i = 1, 3 - read(11,*) textholder - end do + read(11,*) textholder !Start the reading loop do i = 1, ele_in - read(11,*) j, ele, element_type, in_basis, esize + read(11,*) j, 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') @@ -1044,10 +1043,10 @@ module io !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 - + lat_vec(:,1) = (r_in(:,1,2) - r_in(:,1,1))/(2*esize) + lat_vec(:,2) = (r_in(:,1,4) - r_in(:,1,1))/(2*esize) + lat_vec(:,3) = (r_in(:,1,5) - r_in(:,1,1))/(2*esize) + !Now shift all the nodal positions select case(trim(adjustl(in_lattice_type))) case('fcc','FCC') @@ -1068,7 +1067,7 @@ module io 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) + call add_element(in_lattice_type, esize, lat_type, sub_box_num, r_in(:,1:max_basisnum,1:max_ng_node)) case('Atom') read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1) @@ -1082,17 +1081,24 @@ module io !This code parses input values integer, intent(in) :: apos integer :: arglen, arg_pos + + character(len=100) :: textholder arg_pos = apos + 1 - call get_command_argument(arg_pos, in_lapa, arglen) + call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) then print *, "Missing lattice parameter for set_input_lat" end if + read(textholder,*) in_lapa + print *, in_lapa arg_pos = arg_pos + 1 - call get_command_argument(arg_pos, in_lattice_type, arglen) + call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) then print *, "Missing lattice type for set_input_lat" end if - end subroutine set_input_lat(arg_pos) + read(textholder,*) in_lattice_type + print *, in_lattice_type + + end subroutine set_cac end module io diff --git a/src/main.f90 b/src/main.f90 index 8d153ce..53fa4fa 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -61,10 +61,8 @@ program main case('-wrap') wrap_flag=.true. - !This gives necessary information in order to correctly read .cac files case('-set_cac') call set_cac(i) - end select end do !Determine if a mode is being used and what it is. The first argument has to be the mode