From 4868a30e5740d01661808951637161eeda13686a Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 20 Mar 2020 10:48:36 -0400 Subject: [PATCH] Added read_pycac subroutine --- src/io.f90 | 200 +++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 185 insertions(+), 15 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 4aab723..90efe0f 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -389,7 +389,7 @@ module io 10 format('box matrix') 11 format(3f23.15) 12 format('coarse-grained domain') -13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis type x y z') +13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis x y z') 14 format('atomistic domain' / 'ia grain_atom type_atom x y z') 15 format('maximum lattice periodicity length' / 3f23.15) 16 format('Number of lattice types and atom types '/ 2i16) @@ -521,11 +521,10 @@ module io write(11,*) box_bd(:) !Write the number of sub_boxes in the system write(11,*) sub_box_num - !For every subbox write the orientation, sub box boundary, and sub_box_array_bds + !For every subbox write the orientation, sub box boundar do i = 1, sub_box_num write(11,*) sub_box_ori(:,:,i) write(11,*) sub_box_bd(:,i) - write(11,*) ((sub_box_array_bd(j,k,i), j = 1, 2), k = 1, 2) end do !Write the number of atom types in the current model and all of their names @@ -590,13 +589,13 @@ module io end if select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('xyz', 'lmp', 'vtk', 'mb') + case('restart', 'mb') infilenum=infilenum+1 infiles(infilenum) = temp_infile exit case default print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", & - "please input a filename with extension from following list: mb." + "please input a filename with extension from following list: mb, restart." read(*,*) temp_infile end select @@ -615,6 +614,8 @@ module io select case(trim(adjustl(infiles(i)(scan(infiles(i),'.',.true.)+1:)))) case('mb') call read_mb(infiles(i), displace, temp_box_bd) + case('restart') + call read_pycac(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" @@ -643,7 +644,6 @@ module io !Read in the box boundary and grow the current active box bd read(11, *) temp_box_bd(:) - print *, "displace", displace do i = 1, 3 if (abs(displace(i)) > lim_zero) then newdisplace(i) = displace(i) - temp_box_bd(2*i-1) @@ -654,7 +654,6 @@ module io temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) end do - print *, "newdisplace", newdisplace call grow_box(temp_box_bd) !Read in the number of sub_boxes and allocate the variables read(11, *) n @@ -676,14 +675,8 @@ module io sub_box_bd(2*j-1,sub_box_num+i) = sub_box_bd(2*j-1, sub_box_num+i) + displace(j) sub_box_bd(2*j,sub_box_num+i) = sub_box_bd(2*j, sub_box_num+i) + displace(j) end do - !Read in sub_box_array_bd - read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2) end do - !Add the existing element boundaries - 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 - !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 @@ -715,8 +708,8 @@ module io new_loop:do i = 1, new_lattice_types old_loop:do j = 1, lattice_types !First check all the lattice level variables - if ((basisnum(i) == basisnum(j)).and. & - (ng_node(i) == ng_node(j))) then + if ((basisnum(lattice_types + i) == basisnum(j)).and. & + (ng_node(lattice_types + i) == ng_node(j))) then !Now check the basis level variables do ibasis =1, basisnum(i) if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then @@ -770,4 +763,181 @@ module io call set_max_esize end subroutine read_mb + + subroutine read_pycac(file, displace, temp_box_bd) + !Add subroutine for reading in restart files from PyCAC. This code currently only + !works for 8 node elements. + character(len=100), intent(in) :: file + real(kind=dp), dimension(3), intent(in) :: displace + real(kind = dp), dimension(6), intent(out) :: temp_box_bd + + 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) + character(len=100) :: textholder, in_lattype_map(10) + character(len=2) :: atomic_element + !First open the file + open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') + + !Disregard unneeded information + do i = 1, 3 + read(11,*) textholder + end do + + !Read element number + read(11,*) in_eles + + !Discard info and read ng_max_node + do i = 1, 5 + read(11,*) textholder + end do + read(11,*) max_ng_node + + !Read element types (only needed inside this subroutine) + read(11,*) textholder + read(11,*) ele_types + + !Read in atom num + read(11,*) textholder + read(11,*) in_atoms + + !read in lattice_types and atom types + do i = 1,3 + read(11,*) textholder + end do + read(11,*) in_lat_num, in_atom_types + + + !Read define atom_types by name + read(11,*) textholder + do i = 1, in_atom_types + read(11,*) j, atomic_element + call add_atom_type(atomic_element, atom_type_map(i)) + end do + + !Read in the boundary + read(11,*) textholder + read(11,*) box_bc + + !Disregard useless info + do i = 1, 3 + read(11,*) textholder + end do + + !Read in lat_type map + do i = 1, in_lat_num + read(11,*) j, in_lattype_map(i) + ng_node(lattice_types+i) = 8 !Only cubic elements are currently supported in pyCAC + end do + + !Disregard useless info + do i =1 , 3 + read(11,*) textholder + end do + + !Read box boundaries and displace them if necessary + read(11,*) temp_box_bd(:) + 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 + + 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 restart 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 + + !Read in more useless info + do i = 1, 10 + read(11,*) textholder + end do + + !Now read the ele_type to size and lat map + do i = 1, ele_types + read(11,*) j, etype_map(i) + etype_map(i) = etype_map(i) + 1 + end do + + + !Now set up the lattice types. In this code it assumes that lattice_type 1 maps to + !atom type 1 because it only allows 1 atom per basis in pyCAC at the moment. + do i = 1, in_lat_num + basisnum(lattice_types+i) = 1 + basis_type(1,lattice_types+i) = atom_type_map(i) + end do + + !Figure out the lattice type maps in case we have repeated lattice_types + k = lattice_types + 1 + new_lattice_map(:) = 0 + new_loop:do i = 1, in_lat_num + old_loop:do j = 1, lattice_types + !First check all the lattice level variables + if ((basisnum(lattice_types+i) == basisnum(j)).and. & + (ng_node(lattice_types+i) == ng_node(j))) then + !Now check the basis level variables + do ibasis =1, basisnum(i) + if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then + cycle old_loop + end if + end do + new_lattice_map(i) = j + cycle new_loop + end if + end do old_loop + new_lattice_map(i) = k + k = k+1 + end do new_loop + + !Read more useless data + read(11,*) textholder + read(11,*) textholder + + !set max values and allocate variables + max_basisnum = maxval(basisnum) + max_ng_node = maxval(ng_node) + call grow_ele_arrays(in_eles, in_atoms) + + !Now start reading the elements + do i = 1, in_eles + read(11,*) j, etype, k, lat_type + do inod = 1, 8 + read(11, *) j, k, r_in(:,1,inod) + r_in(:,1,inod) = r_in(:,1,inod) + newdisplace + end do + call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in) + end do + + !Read useless data + read(11,*) textholder + read(11,*) textholder + + do i = 1, in_atoms + read(11,*) j, k, atom_type, r_in_atom(:) + r_in_atom = r_in_atom + newdisplace + call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom) + end do + !Close file + close(11) + + lattice_types = maxval(new_lattice_map) + + sub_box_num = sub_box_num + 1 + + call set_max_esize + end subroutine read_pycac end module io