From b0941e44820615d26889e6eca87ac6b1b5fef55c Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 10 Dec 2019 19:33:46 -0500 Subject: [PATCH] Fixes to file reading to ensure that mode_merge works correctly --- src/Makefile | 4 +-- src/io.f90 | 28 ++++++++++++++----- src/mode_create.f90 | 66 ++++++++++++++++++++++++++++++--------------- 3 files changed, 68 insertions(+), 30 deletions(-) diff --git a/src/Makefile b/src/Makefile index 22c84eb..c84186e 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,6 +1,6 @@ FC=ifort -FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -#FFLAGS=-c -mcmodel=large -Ofast +#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin +FFLAGS=-mcmodel=large -Ofast MODES=mode_create.o mode_merge.o mode_convert.o OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) diff --git a/src/io.f90 b/src/io.f90 index 20fac52..e6c4645 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -394,10 +394,12 @@ module io real(kind=dp), dimension(3), intent(in) :: displace real(kind = dp), dimension(6), intent(out) :: temp_box_bd - integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles + integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, & + new_type_to_type(10), new_lattice_types character(len=100) :: etype real(kind=dp) :: r(3), newdisplace(3) real(kind=dp), allocatable :: r_innode(:,:,:) + character(len = 2) :: new_type_to_name(10) !First open the file open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') @@ -433,15 +435,26 @@ module io sub_box_num = sub_box_num + n !Read in the number of atom types and all their names - read(11, *) atom_types, (type_to_name(i), i = 1, atom_types) + 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 + !type of the atoms within this file + do i = 1, new_atom_types + call add_atom_type(new_type_to_name(i), new_type_to_type(i)) + end do !Read the number of lattice types, basisnum, and number of nodes for each lattice type - read(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) + read(11,*) new_lattice_types, (basisnum(i), i = lattice_types+1, lattice_types+new_lattice_types), & + (ng_node(i), i = lattice_types+1, lattice_types+new_lattice_types) !Define max_ng_node and max_basis_num max_basisnum = maxval(basisnum) max_ng_node = maxval(ng_node) !Read the basis atom types for every lattice - read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types) - + read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = lattice_types+1, lattice_types+new_lattice_types) + !Convert the basis_atom types + do j = lattice_types+1, lattice_types+new_lattice_types + do i = 1, basisnum(j) + basis_type(i,j) = new_type_to_type(basis_type(i,j)) + end do + end do !Read number of elements and atoms and allocate arrays read(11, *) in_atoms, in_eles call grow_ele_arrays(in_eles, in_atoms) @@ -450,7 +463,7 @@ module io !Read the atoms do i = 1, in_atoms read(11,*) j, type, r(:) - call add_atom(type, r+newdisplace) + call add_atom(new_type_to_type(type), r+newdisplace) end do !Read the elements @@ -462,11 +475,12 @@ module io r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do - + type = type + lattice_types call add_element(etype, size, type, r_innode) end do !Close the file being read close(11) + lattice_types = lattice_types + new_lattice_types end subroutine read_mb end module io \ No newline at end of file diff --git a/src/mode_create.f90 b/src/mode_create.f90 index b82cc82..08e201a 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -12,8 +12,8 @@ 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) - 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) + 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 @@ -58,16 +58,15 @@ module mode_create allocate(r_node_temp(3,max_basisnum,max_ng_node)) - if(dup_flag) then + !Get the inverse orientation matrix we will need later + call matrix_inverse(orient,3,orient_inv) - !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) + if(dup_flag) then - + !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 @@ -76,21 +75,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 call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) @@ -112,6 +115,19 @@ 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 @@ -208,13 +224,21 @@ 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)