Working read from pycac.out file format

src/Makefile

@@ -1,7 +1,7 @@
 FC=ifort
 FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
 #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
-MODES=mode_create.o mode_merge.o mode_convert.o
+MODES=mode_create.o mode_merge.o mode_convert.o 
 OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o sorts.o
 

src/call_mode.f90

@@ -1,4 +1,4 @@
-subroutine call_mode(arg_pos,mode)
+subroutine call_mode(arg_pos)
     !This code is used to parse the command line argument for the mode information and calls the required 
     !mode module.
 
@@ -10,7 +10,6 @@ subroutine call_mode(arg_pos,mode)
     implicit none    
 
     integer, intent(out) :: arg_pos
-    character(len=100), intent(in) :: mode
 
     select case(mode)
     case('--create')

src/call_option.f90

@@ -39,6 +39,8 @@ subroutine call_option(option, arg_pos)
         call run_delete(arg_pos)
     case('-set_cac')
         arg_pos=arg_pos +3
+    case('-set_types')
+        arg_pos = arg_pos + 3 + atom_types
     case('-redef_box')
         call redef_box(arg_pos)
     case default

src/elements.f90

@@ -13,6 +13,8 @@ module elements
     character(len=100), allocatable :: type_ele(:) !Element type 
     integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
+    !Element result data structures
+    real(kind=8), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:)
 
     integer, save :: ele_num !Number of elements
     integer, save :: node_num !Total number of nodes
@@ -22,6 +24,8 @@ module elements
     integer, allocatable :: sbox_atom(:), tag_atom(:)
     real(kind =dp),allocatable :: r_atom(:,:) !atom position
     integer :: atom_num=0 !Number of atoms
+    !Atom result data structures information
+    real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:)
 
     !Mapping atom type to provided name
     character(len=2), dimension(10) :: type_to_name
@@ -669,17 +673,17 @@ module elements
     esize = size_ele(ie)
     select case(iface)
         case(1)
-            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /)
+            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /)
         case(2)
-            pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
+            pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
         case(3)
-            pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
+            pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
         case(4)
-            pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
+            pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
         case(5)
-            pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
+            pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
         case(6)
-            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /)
+            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /)
     end select
 
     !Now transform it to real space and adjust it to the position of the element in the first node.
@@ -736,4 +740,57 @@ module elements
 
     end subroutine lattice_map
 
+    subroutine alloc_dat_arrays(n,m)
+        !This subroutine used to provide initial allocation for the atom and element data arrays
+        integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
+        integer :: allostat
+
+        !Allocate element arrays
+        if (n > 0) then 
+            allocate(force_node(3,max_basisnum, max_ng_node, n), &
+                     virial_node(6,max_basisnum, max_ng_node, n), &
+                     energy_node(max_basisnum,max_ng_node,n),    &
+                     stat=allostat)
+            if(allostat > 0) then 
+                print *, "Error allocating element data arrays in mode_metric becaus of:", allostat
+                stop
+            end if
+
+        end if
+
+        if (m > 0) then
+            allocate(force_atom(3, m), &
+                     virial_atom(6, m), &
+                     energy_atom(m),    &
+                     stat=allostat)
+            if(allostat > 0) then 
+                print *, "Error allocating atom data arrays in mode_metric becaus of:", allostat
+                stop
+            end if
+        end if
+
+    end subroutine 
+
+    subroutine add_atom_data(ia, eng, force, virial)
+        !Function which sets the atom data arrays 
+        integer, intent(in) :: ia
+        real(kind=dp), intent(in) :: eng, force(3), virial(6)
+
+        energy_atom(ia) = eng
+        force_atom(:,ia) = force(:) 
+        virial_atom(:,ia) = virial(:) 
+        return
+    end subroutine add_atom_data
+
+    subroutine add_element_data(ie, eng, force, virial)
+        !Function which sets the element data arrays
+        integer, intent(in) :: ie
+        real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
+                                     force(3,max_basisnum, max_ng_node), &
+                                     virial(6,max_basisnum,max_ng_node)
+        energy_node(:,:,ie)   = eng
+        force_node(:,:,:,ie)  = force
+        virial_node(:,:,:,ie) = virial
+        return
+    end subroutine add_element_data
 end module elements

src/io.f90

@@ -106,7 +106,7 @@ module io
                 call write_lmpcac(outfiles(i))
             case default
                 print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
-                         " is not accepted for writing. Please select from: xyz and try again"
+                         " is not accepted for writing. Please select from: xyz, lmp, vtk, mb, restart, cac and try again"
                 stop 
 
             end select
@@ -599,9 +599,24 @@ module io
                 infilenum=infilenum+1
                 infiles(infilenum) =  temp_infile
                 exit
+            case('out')
+                if(atom_types == 0) then 
+                    print *, "Please run -set_types command prior to running code requiring reading in pycac_*.out files"
+                    stop 3
+                end if
+                select case(trim(adjustl(mode)))
+                case('--convert','--metric')
+                    infilenum = infilenum+1
+                    infiles(infilenum) = temp_infile
+                    exit
+                case default
+                    print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode))
+                    stop 3
+                end select
+
             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, restart."
+                         "please input a filename with extension from following list: mb, restart, cac, or out."
                 read(*,*) temp_infile
 
             end select
@@ -624,9 +639,11 @@ module io
             call read_pycac(infiles(i), displace, temp_box_bd)
         case('cac')
             call read_lmpcac(infiles(i), displace, temp_box_bd)
+        case('out')
+            call read_pycac_out(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"
+                     " is not accepted for reading. Please select from: mb,restart,cac,out and try again"
             stop 
 
         end select
@@ -781,7 +798,7 @@ module io
 
         integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
                    atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), &
-                   atom_type, stat
+                   atom_type, stat 
         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
@@ -960,6 +977,102 @@ module io
         end if
     end subroutine read_pycac
 
+    subroutine read_pycac_out(file, displace, temp_box_bd)
+        !This subroutine reads in the pyCAC dump file
+    
+
+        !Arguments
+        character(len=100), intent(in) :: file
+        real(kind=dp), dimension(3), intent(in) :: displace
+        real(kind=dp), dimension(6), intent(out) :: temp_box_bd
+
+        !Internal Variables
+        integer :: i, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
+            id, type_node, ilat, esize, tag, type
+        real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), &
+                         ee(1,8), fe(3,1,8), ve(3,1,8), re(3,1,8)
+        character(len=100) :: textholder, fcc
+
+
+        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') 
+
+        !Now initialize some important variables if they aren't defined
+        if (max_basisnum==0) max_basisnum = 1
+        if (max_ng_node==0) max_ng_node=8
+        fcc="fcc"
+
+        !Skip header comment lines
+        read(11, *) textholder
+        read(11, *) textholder
+
+        !Read atom number and element number and grow element arrays by needed amount
+        read(11,*) textholder, in_atoms, textholder, in_eles
+        call grow_ele_arrays(in_eles, in_atoms)
+        call alloc_dat_arrays(in_eles, in_atoms)
+
+        !Read boundary information
+        read(11,*) textholder, box_bc(1:1), box_bc(2:2), box_bc(3:3), 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 boundaries
+        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
+
+        if(in_atoms > 0 ) then
+            !Read atom header
+            read(11,*) textholder
+            do ia = 1, in_atoms
+                read(11,*) tag, type, ra(:), ea, fa(:), va(:)         
+                call add_atom(tag, type, sub_box_num, ra)
+                call add_atom_data(atom_num, ea, fa, va)
+            end do
+
+        end if
+
+        if(in_eles > 0) then 
+            !Add the lattice_types based on the atom types
+            inbtypes=0
+            do i = 1, maxval(type_atom)
+                inbtypes(1) = i 
+                call lattice_map(1, inbtypes,  8 , 1.0_dp, ilat) !Please check documentation on pycac.out formats
+            end do 
+            !Read element and node headers
+            read(11,*) textholder
+            read(11,*) textholder
+
+            !read element information, currently only 8 node elements with 1 basis
+            do ie =1, in_eles
+                    read(11,*)  tag, lat_type, textholder, textholder, esize
+                do inod =1, 8
+                    read(11,*)  textholder, textholder, textholder, re(:,1,inod), ee(1,inod), fe(:,1,inod), ve(:,1,inod)
+                end do
+                call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re)
+                call add_element_data(ele_num, ee, fe, ve)
+            end do
+        end if
+        call set_max_esize
+        return
+  end subroutine
+
     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
@@ -980,7 +1093,7 @@ module io
         !Open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
         
-        !Now initialiaze some important variables if they aren't defined
+        !Now initialize some important variables if they aren't defined
         if (max_basisnum==0) max_basisnum = 10
         if (max_ng_node==0) max_ng_node=8
 
@@ -1109,4 +1222,25 @@ module io
         print *, in_lattice_type
 
     end subroutine set_cac
+
+    subroutine set_types(apos)
+        !This code 
+        integer, intent(in) :: apos
+        integer :: i, j,arglen, arg_pos, ntypes
+
+        character(len=100) :: textholder
+
+        arg_pos = apos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) stop "Missing numtypes in io"
+        read(textholder,*) ntypes
+
+        do i=1,ntypes
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, textholder, arglen)
+            call add_atom_type(textholder, j)
+        end do
+
+        return
+    end subroutine set_types
 end module io

src/main.f90

@@ -63,6 +63,8 @@ program main
         
         case('-set_cac')
             call set_cac(i)
+        case('-set_types')
+            call set_types(i)
         end select 
     end do 
     !Determine if a mode is being used and what it is. The first argument has to be the mode
@@ -71,7 +73,8 @@ program main
 
     argument = trim(adjustl(argument))
     if (argument(1:2) == '--') then 
-        call call_mode(end_mode_arg, argument)
+        mode = argument
+        call call_mode(end_mode_arg)
     end if
 
     !Check to make sure a mode was called

src/parameters.f90

@@ -12,5 +12,7 @@ module parameters
     
     !Numeric constants
     real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
-
+    
+    !Mode type that is being called 
+    character(len=100) :: mode
 end module parameters