Merge read-cac branch into main branch

README.md

@@ -329,7 +329,6 @@ This command will delete all overlapping atoms within a specific cutoff radius `
 This option is primarily used when reading data from non .mb formats. This code simply sets the orientation variable for the specified sub box `sbox`.
 
 ****
-
 ## Position Specification
 
 Specifying positions in cacmb can be done through a variety of ways. Examples of each format is shown below.

src/call_option.f90

@@ -33,6 +33,8 @@ subroutine call_option(option, arg_pos)
         call sbox_ori(arg_pos)
     case('-delete')
         call run_delete(arg_pos)
+    case('-set_cac')
+        arg_pos=arg_pos +3
     case default
         print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
         stop 3 

src/elements.f90

@@ -36,13 +36,13 @@ module elements
     integer :: lattice_types = 0
     integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type
     integer :: max_esize=0 !Maximum number of atoms per side of element
+    real(kind=dp) :: lapa(10)
 
     !These variables contain information on the basis, for simplicities sake we limit 
     !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
     !This can be easily increased with no change to efficiency
     integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type 
     integer :: basis_type(10,10)
-    real(kind=dp) :: lapa(10)
 
     !Additional module level variables we need 
     logical :: wrap_flag
@@ -323,8 +323,6 @@ module elements
 
         integer :: i
 
-        max_ng_node = 0
-
         do i=1, n
             select case(trim(adjustl(element_types(i))))
             case('fcc')
@@ -676,4 +674,39 @@ module elements
     end select
     end subroutine
 
+    subroutine lattice_map(in_bnum, in_btypes, in_ngnodes, in_lapa, lat_type)
+        !This subroutine maps an input lattice type to either a new lattice type or an existing one depending on basis_type and 
+        !number of nodes at the atoms
+
+        integer, intent(in) :: in_ngnodes, in_bnum, in_btypes(10) !Input variables
+        real(kind=dp), intent(in) :: in_lapa
+        integer, intent(out) :: lat_type
+
+        integer j, ibasis
+
+        lat_type = 0
+        lat_loop:do j = 1, lattice_types
+            !Check all the lattice level variables
+            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 lat_loop
+                end do
+                lat_type = j
+                exit lat_loop
+            end if
+        end do lat_loop
+
+        !If it doesn't match an existing lattice type we add it 
+        if( lat_type == 0) then 
+            lattice_types = lattice_types + 1
+            basisnum(lattice_types) = in_bnum
+            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
+
 end module elements

src/io.f90

@@ -8,9 +8,9 @@ module io
     implicit none
 
     integer :: outfilenum = 0, infilenum = 0
-    character(len=100) :: outfiles(100), infiles(100)
+    character(len=100) :: outfiles(100), infiles(100), in_lattice_type=''
     logical :: force_overwrite
-
+    real(kind=dp) :: in_lapa=0.0
    public 
    contains
 
@@ -196,7 +196,7 @@ module io
         do i = 1, ele_num
             call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
             select case(trim(adjustl(type_ele(i))))
-            case('fcc')
+            case('fcc','bcc')
                 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                     interp_num = interp_num+1
                     call apply_periodic(r_interp(:,iatom))
@@ -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"
@@ -776,7 +778,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
-        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
@@ -949,4 +951,154 @@ module io
             call set_max_esize
         end if
     end subroutine read_pycac
+
+    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
+        real(kind = dp), dimension(6), intent(out) :: temp_box_bd
+        
+        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), 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 
+            print *, "Please use set_cac to set needed parameters to read in .cac file"
+            stop 3
+        end if
+        !Open the file
+        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
+        
+        !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 number of elements
+        read(11, *) ele_in, textholder 
+        read(11, *) type_in, textholder
+
+        !Read box_boundaries
+        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
+            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 atomic masses
+        do i = 1, type_in
+            read(11,*) j, mass, textholder
+            call ATOMMASSSPECIES(mass, atom_species) 
+            call add_atom_type(atom_species, type_map(i))
+        end do
+
+        !Read useless info
+        read(11,*) textholder
+
+        !Start the reading loop
+        do i = 1, ele_in
+            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')
+                !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_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')
+                    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
+                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(:,1:max_basisnum,1:max_ng_node))
+                
+            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_lmpcac
+
+    subroutine set_cac(apos)
+        !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, 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, textholder, arglen)
+        if (arglen==0) then 
+            print *, "Missing lattice type for set_input_lat"
+        end if
+        read(textholder,*) in_lattice_type
+        print *, in_lattice_type
+
+    end subroutine set_cac
 end module io

src/main.f90

@@ -60,6 +60,9 @@ program main
         !This lets us know if we need to wrap atomic positions back into the cell
         case('-wrap')
             wrap_flag=.true.
+        
+        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