Latest working version of cacmb

src/Makefile

@@ -2,8 +2,8 @@
 .DEFAULT_GOAL := all
 
 FC=mpif90
-FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal 
+#FFLAGS=-Wall -mcmodel=large -O0 -g -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal 
-FFLAGS=-O2 -mcmodel=large
+FFLAGS=-O3 -g
 
 OBJDIR=obj
 SRCS := $(wildcard *.f90)

src/elements.f90

@@ -33,7 +33,8 @@ module elements
     integer :: atom_types = 0
 
     !Variables for creating elements based on primitive cells
-    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3)
+    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), &
+                     cube20(3,20)   
     integer :: cubic_faces(4,6)
 
     !Below are lattice type arrays which provide information on the general form of the elements.
@@ -72,6 +73,29 @@ module elements
                               0.0_dp, 1.0_dp, 1.0_dp /), &
                              shape(fcc_cell))
 
+        !Now initialize the 20 node cube element
+        cube20 = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 1.0_dp, 1.0_dp, &
+                             0.0_dp, 1.0_dp, 1.0_dp, &
+                             0.5_dp, 0.0_dp, 0.0_dp, &
+                             1.0_dp, 0.5_dp, 0.0_dp, &
+                             0.5_dp, 1.0_dp, 0.0_dp, &
+                             0.0_dp, 0.5_dp, 0.0_dp, &
+                             0.0_dp, 0.0_dp, 0.5_dp, &
+                             1.0_dp, 0.0_dp, 0.5_dp, &
+                             1.0_dp, 1.0_dp, 0.5_dp, &
+                             0.0_dp, 1.0_dp, 0.5_dp, &
+                             0.5_dp, 0.0_dp, 1.0_dp, &
+                             1.0_dp, 0.5_dp, 1.0_dp, &
+                             0.5_dp, 1.0_dp, 1.0_dp, &
+                             0.0_dp, 0.5_dp, 1.0_dp /), &
+                             shape(cube20))
+
         !Now we create a list containing the list of vertices needed to describe the 6 cube faces
         cubic_faces(:,1) = (/ 1, 2, 3, 4 /) 
         cubic_faces(:,2) = (/ 1, 2, 6, 5 /)
@@ -384,6 +408,8 @@ module elements
                 ng_node(i) = 8
             case('bcc')
                 ng_node(i) = 8
+            case('20fcc')
+                ng_node(i) = 20
             end select
 
             if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
@@ -462,6 +488,34 @@ module elements
                             end do
                         end do
 
+                    end do
+                end do
+            end do 
+        case('20fcc')
+            !Now loop over all the possible sites
+            do it = 1, esize
+                t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
+                do is =1, esize
+                    s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
+                    do ir = 1, esize
+                        r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
+                        call quad_rhombshape(r,s,t,a_shape)
+
+                        do ibasis = 1, bnum
+                            ia = ia + 1
+                            do inod = 1, 20
+                                type_interp(ia) = basis_type(ibasis,lat_type_temp)
+                                r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)
+
+                                if(present(data_interp)) then 
+                                    !If data is present then interpolate data arrays as well
+                                    data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod)
+                                    data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod)
+                                    data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod)
+                                end if
+                            end do
+                        end do
+
                     end do
                 end do
             end do 
@@ -492,6 +546,42 @@ module elements
         return
     end subroutine rhombshape
 
+    subroutine quad_rhombshape(r,s,t, shape_fun)
+        !Shape function for 20 node quadratic rhombohedral elements
+        real(kind=8), intent(in) :: r, s, t
+        real(kind=8), intent(out) :: shape_fun(:)
+
+        !Corner nodes
+        shape_fun(1) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp-t)*(-r-s-t-2)/8.0_dp
+        shape_fun(2) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp-t)*(r-s-t-2)/8.0_dp
+        shape_fun(3) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp-t)*(r+s-t-2)/8.0_dp
+        shape_fun(4) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp-t)*(-r+s-t-2)/8.0_dp
+        shape_fun(5) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp+t)*(-r-s+t-2)/8.0_dp
+        shape_fun(6) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp+t)*(r-s+t-2)/8.0_dp
+        shape_fun(7) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp+t)*(r+s+t-2)/8.0_dp
+        shape_fun(8) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp+t)*(-r+s+t-2)/8.0_dp
+
+        !Side nodes, first node r is zero
+        shape_fun(9) = (1-r*r)*(1-s)*(1-t)/4.0_dp
+        shape_fun(11) = (1-r*r)*(1+s)*(1-t)/4.0_dp
+        shape_fun(17) = (1-r*r)*(1-s)*(1+t)/4.0_dp
+        shape_fun(19) = (1-r*r)*(1+s)*(1+t)/4.0_dp
+
+        !node s is zero
+        shape_fun(10) = (1+r)*(1-s*s)*(1-t)/4.0_dp
+        shape_fun(12) = (1-r)*(1-s*s)*(1-t)/4.0_dp
+        shape_fun(18) = (1+r)*(1-s*s)*(1+t)/4.0_dp
+        shape_fun(20) = (1-r)*(1-s*s)*(1+t)/4.0_dp
+
+        !node t is zero
+        shape_fun(13) = (1-r)*(1-s)*(1-t*t)/4.0_dp
+        shape_fun(14) = (1+r)*(1-s)*(1-t*t)/4.0_dp
+        shape_fun(15) = (1+r)*(1+s)*(1-t*t)/4.0_dp
+        shape_fun(16) = (1-r)*(1+s)*(1-t*t)/4.0_dp
+
+
+    end subroutine quad_rhombshape
+
     subroutine delete_atoms(num, in_index)
         !This subroutine deletes atoms from the atom arrays 
         integer, intent(in) :: num
@@ -624,6 +714,9 @@ do i = 1, atom_num
         else if(trim(adjustl(pos_string)) == '-inf') then
             pos=box_bd(2*i-1)
 
+        else if(trim(adjustl(pos_string)) == 'mid') then 
+            pos=(box_bd(2*i)+box_bd(2*i-1))/2
+
         else if (trim(adjustl(pos_string)) == 'rand') then 
             call random_number(rand)
             pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
@@ -647,7 +740,6 @@ do i = 1, atom_num
                 end do 
             end if
 
-            read(cone, *) face
             if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" 
             !Now get the position
             call offset_pos(ele, face, rand_ele_pos)
@@ -670,6 +762,15 @@ do i = 1, atom_num
             end if
             pos = box_bd(2*i) - pos
 
+        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'mid')>0)) then 
+            !Now extract the number we are reducing from midinity
+            if(index(pos_string,'mid') < index(pos_string,'-')) then 
+                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
+            end if
+            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp - pos
+
         else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
             !Now extract the number we are reducing from infinity
            if(index(pos_string,'inf') < index(pos_string,'+')) then 
@@ -679,6 +780,15 @@ do i = 1, atom_num
             end if
             pos = box_bd(2*i-1) + pos
 
+        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'mid')>0)) then 
+            !Now extract the number we are reducing from midinity
+           if(index(pos_string,'mid') < index(pos_string,'+')) then 
+                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
+            else
+                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
+            end if
+            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp + pos
+
         else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
             !Now extract the number we are reducing from infinity
            if(index(pos_string,'inf') < index(pos_string,'*')) then 
@@ -809,6 +919,7 @@ do i = 1, atom_num
         integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
         integer :: allostat
 
+        print *, max_ng_node
         !Allocate element arrays
         if (n > 0) then 
             if (allocated(force_node)) then 

src/io.f90

@@ -177,6 +177,7 @@ module io
         do i = 1,ele_num
             if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
             if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3
+            if(type_ele(i) == '20fcc') write_num = write_num + size_ele(i)**3
         end do
 
         !Write total number of atoms + elements
@@ -212,7 +213,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','bcc')
+            case('fcc','bcc', '20fcc')
                 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                     interp_num = interp_num+1
                     call apply_periodic(r_interp(:,iatom))
@@ -255,6 +256,7 @@ module io
             do i = 1,ele_num
                 if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
                 if(type_ele(i) == 'bcc') write_num = write_num + size_ele(i)**3
+                if(type_ele(i) == '20fcc') write_num = write_num + size_ele(i)**3
             end do
         end if
         !Write total number of atoms
@@ -311,7 +313,7 @@ module io
                     call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
                 end if
                 select case(trim(adjustl(type_ele(i))))
-                case('fcc','bcc')
+                case('fcc','bcc', '20fcc')
                     if(write_dat) then 
                         do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                             interp_num = interp_num+1
@@ -562,7 +564,12 @@ module io
             ip = 0
             write(11,13)
             do i = 1, ele_num
+                select case(type_ele(i))
+                case('fcc','bcc')
                     write(11, '(4i16)') i, basisnum(lat_ele(i)), 2, (size_ele(i)-1)
+                case('20fcc')
+                    write(11, '(4i16)') i, basisnum(lat_ele(i)), 3, (size_ele(i)-1)
+                end select
                 do inod = 1, ng_node(lat_ele(i))
                     do ibasis = 1, basisnum(lat_ele(i))
                         ip = ip + 1
@@ -624,7 +631,7 @@ module io
         !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
         !every basis at every node
         do i = 1, ele_num
-            write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i)
+            write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), trim(adjustl(type_ele(i)))
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis =1, basisnum(lat_ele(i))
                     write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
@@ -1003,7 +1010,7 @@ module io
         integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
             id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip, atom_type_map(100)
         real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), &
-                         ee(10,8), fe(3,10,8), ve(6,10,8), re(3,10,8), atomic_masses(10)
+                         ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), atomic_masses(10)
         character(len=100) :: textholder, fcc
         character(len=1000) :: line
         character(len=2) :: atomic_element
@@ -1013,7 +1020,7 @@ module io
 
         !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
+        if (max_ng_node==0) max_ng_node=20
         fcc="fcc"
 
         !Skip header comment lines
@@ -1096,6 +1103,11 @@ module io
                     inbtypes(ibasis) = atom_type_map(inbtypes(ibasis))
                 end do
                 call lattice_map(bnum, inbtypes, n, 1.0_dp, lat_type)
+                if (n == 8)then 
+                    fcc="fcc"
+                else if (n == 20) then 
+                    fcc = "20fcc"
+                end if
                 call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re)
                 call add_element_data(ele_num, ee, fe, ve)
             end do
@@ -1135,6 +1147,8 @@ module io
         read(11,*) temp_box_bd(3:4), textholder, textholder
         read(11,*) temp_box_bd(5:6), textholder, textholder
 
+        call grow_ele_arrays(0,atom_in)
+
         print *, "Read in ", atom_in, " atoms from ", trim(adjustl(file))
         !Shift the box boundaries if needed
         do i = 1, 3

src/mode_create.f90

@@ -67,7 +67,7 @@ module mode_create
             end do
             !Now get the rotated box vertex positions in lattice space. Should be integer units and get the new maxlen
             select case(trim(adjustl(element_type)))
-            case('fcc')
+            case('fcc', '20fcc')
                 box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
                 maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
             case('bcc')
@@ -93,7 +93,7 @@ module mode_create
             end do
             !Now get the rotated box vertex positions in lattice space. Should be integer units
             select case(trim(adjustl(element_type)))
-            case('fcc')
+            case('fcc', '20fcc')
                 box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
             case('bcc')
                 box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1
@@ -131,9 +131,11 @@ module mode_create
             !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)
+                call build_with_rhomb(box_lat_vert, fcc_mat, 8)
             case('bcc')
-                call build_with_rhomb(box_lat_vert, bcc_mat)
+                call build_with_rhomb(box_lat_vert, bcc_mat, 8)
+            case('20fcc')
+                call build_with_rhomb(box_lat_vert, fcc_mat, 20)
             case default
                 print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",&
                          "element type"
@@ -280,7 +282,7 @@ module mode_create
 
         !Check special periodicity relations
         select case(trim(adjustl(element_type)))
-        case('fcc')
+        case('fcc', '20fcc')
             do i = 1,3
                 !Check if one of the directions is 110
                 if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.&
@@ -331,19 +333,20 @@ module mode_create
         !
     end subroutine
     
-    subroutine build_with_rhomb(box_in_lat, transform_matrix)
+    subroutine build_with_rhomb(box_in_lat, transform_matrix, nn)
         !This subroutine returns all the lattice points in the box in r_lat
 
         !Inputs
         integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
         real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
+        integer, intent(in) :: nn
         !Internal variables
-        integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), efill_size, &
+        integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,nn), rzero(3), efill_size, &
-                   vlat(3), temp_lat(3,8), m, n, o, j, k, nump_ele, efill_temp_lat(3,8), filzero(3), bd_ele_lat(6),&
+                   vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), &
-                   efill_ele(3,8), ebd(6)
+                   bd_ele_lat(6), efill_ele(3,nn), ebd(6)
-        real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6)
+        real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6)
         logical, allocatable :: lat_points(:,:,:)
-        logical :: node_in_bd(8), add, lat_points_ele(esize,esize,esize), intersect_bd(3)
+        logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3)
 
         !Do some value initialization
         max_esize = esize
@@ -393,7 +396,12 @@ module mode_create
             !If we are working with elements we have to use more complex code
 
             !Initialize finite element 
+            select case(element_type)
+            case('fcc','bcc')
                 ele(:,:) = (esize-1) * cubic_cell(:,:)
+            case('20fcc')
+                ele(:,:) = (esize-1)*cube20
+            end select
             
             !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. 
             !A value of true means that the coordinate is a lattice point which is within the box_bd
@@ -423,7 +431,7 @@ module mode_create
                 do iy = 1, bd_in_array(2)
                     do ix = 1, bd_in_array(1)
                         node_in_bd(:) = .false.
-                        do inod = 1, 8
+                        do inod = 1, nn
                             vlat = ele(:,inod) + (/ ix, iy, iz /)
 
                             !Check to see if the node resides at a position containing a lattice point within the box
@@ -445,16 +453,22 @@ module mode_create
             !Now build the finite element region
             lat_ele_num = 0
             lat_atom_num = 0
-            allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/esize))
+            allocate(r_lat(3,nn,numlatpoints/esize), elat(numlatpoints/esize))
 
+            select case(element_type)
+            case('fcc','bcc')
                 ele(:,:) = (esize-1) * cubic_cell(:,:)
+            case('20fcc')
+                ele(:,:) = (esize-1)*cube20
+            end select
+ 
             !Redefined the second 3 indices as the number of elements that fit in the bounds
             do i = 1, 3
                 bd_in_array(3+i) = bd_in_array(i)/esize
             end do
 
             !Now start the element at rzero
-            do inod=1, 8
+            do inod=1, nn
                 ele(:,inod) = ele(:,inod) + rzero
             end do
             do iz = -bd_in_array(6), bd_in_array(6)
@@ -463,7 +477,7 @@ module mode_create
                         node_in_bd(:) = .false.
                         temp_nodes(:,:,:) = 0.0_dp
                         temp_lat(:,:) = 0
-                        do inod = 1, 8
+                        do inod = 1, nn
                             vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
                             !Transform point back to real space for easier checking
                             ! v = matmul(orient, matmul(transform_matrix,v))
@@ -537,7 +551,7 @@ module mode_create
                                     ze: do k = 1, (esize-efill_size)
                                         ye: do j = 1, (esize-efill_size)
                                             xe: do i = 1, (esize-efill_size)
-                                                do inod = 1,8
+                                                do inod = 1,nn
                                                     vlat = efill_ele(:,inod) + (/ i, j, k /)
                                                     if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe
                                                      do o = 1,3

src/mode_metric.f90

@@ -193,7 +193,7 @@ module mode_metric
         integer :: i, inod, ibasis
 
         npoints=0
-
+        rout = -huge(1.0_dp)
         if(atom_num > 0) then 
             do i = 1, atom_num
                 rout(:,tag_atom(i))  = r_atom(:,i)

src/neighbors.f90

@@ -41,7 +41,6 @@ module neighbors
             box_len(i) = box_bd(2*i) - box_bd(2*i-1)
             cell_num(i) = int(box_len(i)/(rc_off/2))+1
         end do 
-        print *, box_len, cell_num
 
         !Initialize/allocate variables 
         cell_lim = 10
@@ -49,13 +48,15 @@ module neighbors
 
         !Now place points within cell
         num_in_cell = 0
+        
         do i = 1, numinlist
             !Check to see if the current point is a filler point and if so just skip it
-            if(r_list(1,i) < -huge(1.0_dp)+1) cycle
+            if(r_list(1,i) < box_bd(1)) cycle
 
             !c is the position of the cell that the point belongs to 
+            if(i == 11651203) print *, r_list(:,i), (r_list(1,i) < box_bd(1))
             do j = 1, 3
-                c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1
+                c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 1
             end do
 
             !Place the index in the correct position, growing if necessary
@@ -102,7 +103,7 @@ module neighbors
         pointloop: do i = 1, n
             
             !First check to see if the point is a filler point, if so then skip it 
-            if(r_list(1,i) < -Huge(-1.0_dp)+1) cycle
+            if(r_list(1,i) < box_bd(1)) cycle
 
             !c is the position of the cell that the point
             c = which_cell(:,i)

src/opt_group.f90

@@ -644,7 +644,7 @@ module opt_group
         if(group_ele_num > 0) then 
             orig_atom_num = atom_num
             !Estimate number of atoms we are adding, this doesn't have to be exact
-            add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
+            add_atom_num = group_ele_num*max_basisnum*(max_esize)**3
             call grow_ele_arrays(0,add_atom_num)
             do i = 1, group_ele_num
                 ie = element_index(i)