Merge pull request #4 from aselimov/bg--fixes-to-mode-merge-reading

Fixes to file reading to ensure that mode_merge works correctly

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)
 

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

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)