Merge pull request #6 from aselimov/ft-write-lammpscac

Ft write lammpscac

src/Makefile

@@ -1,6 +1,7 @@
+ 
 FC=ifort
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
-FFLAGS=-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)
 
@@ -8,7 +9,7 @@ OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box
 .SUFFIXES: .c .f .f90 .F90 .o
 
 cacmb: $(OBJECTS) 
-	$(FC) $(FFLAGS) $(OBJECTS) -o $@
+	$(FC) $(FFLAGS) $(OBJECTS) parameters.o -o $@
 
 .f90.o: 
 	$(FC) $(FFLAGS) -c $< 
@@ -30,4 +31,4 @@ main.o io.o build_subroutines.o: elements.o
 call_mode.o : $(MODES)
 $(MODES) io.o: atoms.o box.o
 $(MODES) main.o : io.o
-testfuncs.o elements.o mode_create.o: subroutines.o
+testfuncs.o elements.o mode_create.o: subroutines.o

src/elements.f90

@@ -101,10 +101,28 @@ module elements
 
         real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat
 
+        integer :: inod, i
+        real(kind=dp), dimension(3,max_ng_node) :: adjustVar
+
+        adjustVar(:,:) = 0.0_dp
+
         select case(trim(ele_type))
 
         case('fcc')
-            cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, fcc_cell))
+            if(lmpcac) then 
+                do inod = 1, 8
+                    do i = 1,3
+                        if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                            adjustVar(i,inod) = -0.5_dp
+                        else 
+                            adjustVar(i, inod) = 0.5_dp
+                        end if
+                    end do
+                end do
+                adjustVar(:,1:8) = matmul(fcc_mat, adjustVar(:,1:8))
+            end if
+            cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8)
+            cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
         case default
             print *, "Element type ", trim(ele_type), " currently not accepted"
             stop

src/io.f90

@@ -9,6 +9,7 @@ module io
 
     integer :: outfilenum = 0, infilenum = 0
     character(len=100) :: outfiles(10), infiles(10)
+    
 
    public 
    contains
@@ -62,9 +63,14 @@ module io
                 outfilenum=outfilenum+1
                 outfiles(outfilenum) =  temp_outfile
                 exit
+            case('cac')
+                lmpcac = .true. 
+                outfilenum=outfilenum+1
+                outfiles(outfilenum) =  temp_outfile
+                exit
             case default    
-                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", &
-                         "please input a filename with extension from following list: xyz, lmp, vtk."
+                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
+                         "please input a filename with extension from following list: xyz, lmp, vtk, cac."
                 read(*,*) temp_outfile
 
             end select
@@ -72,7 +78,6 @@ module io
 
     end subroutine get_out_file 
 
-
     subroutine write_out
         !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
 
@@ -94,6 +99,8 @@ module io
                 call write_mb(outfiles(i))
             case('restart')
                 call write_pycac(outfiles(i))
+            case('cac')
+                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"
@@ -193,6 +200,71 @@ module io
         end do
     end subroutine write_lmp
 
+    subroutine write_lmpcac(file)
+        !This subroutine writes out a .lmp style dump file
+        character(len=100), intent(in) :: file
+        integer :: write_num, i, inod, ibasis
+        real(kind=dp) :: mass
+
+1       format(i16, '   Eight_Node', 4i16)
+2       format(i16, '   Atom', 4i16)
+3       format(3i16,3f23.15)
+
+        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
+
+        !Comment line
+        write(11, '(a)') '# CAC input file made with cacmb'
+        write(11, '(a)')
+        !Calculate total atom number
+        write_num = atom_num + ele_num
+
+        !Write total number of atoms + elements
+        write(11, '(i16, a)') write_num, ' cac elements'
+        !Write number of atom types 
+        write(11, '(i16, a)') atom_types, ' atom types'
+
+        write(11,'(a)') ' '
+        !Write box bd
+        write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
+        write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
+        write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'
+
+        !Masses
+        write(11, '(a)') 'Masses'
+
+        write(11, '(a)') ' '
+        do i =1, atom_types
+            call atommass(type_to_name(i),mass)
+            write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i)
+        end do
+        write(11, '(a)') ' '
+
+        write(11, '(a)') 'CAC Elements'
+        write(11, '(a)') ' '
+
+        !Write element nodal positions
+        do i = 1, ele_num
+            select case(trim(adjustl(type_ele(i))))
+            case('fcc')
+                !The first entry is the element specifier 
+                write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
+                do ibasis = 1, basisnum(lat_ele(i))
+                    do inod = 1, 8
+                        !Nodal information for every node
+                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
+                    end do
+                end do
+            end select
+        end do
+
+        do i = 1, atom_num
+            !Element specifier dictating that it is an atom
+            write(11,2) ele_num+i, 1, 1, 1, 1
+            !Write the atomic information
+            write(11,3) 1, 1, type_atom(i), r_atom(:,i)
+        end do
+    end subroutine write_lmpcac
+
     subroutine write_vtk(file)
         !This subroutine writes out a vtk style dump file
         integer :: i, j, inod, ibasis
@@ -249,7 +321,7 @@ module io
         do i = 1, ele_num
             do inod=1, ng_node(lat_ele(i))  
                 do ibasis = 1, basisnum(lat_ele(i))
-                    write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
+                    write(11, '(3f23.15)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
                 end do 
             end do 
         end do
@@ -601,4 +673,4 @@ module io
         !just overwrite the arrays
         if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
     end subroutine read_mb
-end module io
+end module io

src/main.f90

@@ -27,6 +27,17 @@ program main
     ! Command line parsing
     arg_num = command_argument_count()
 
+    !First check if we are writing out to lammpscac format by looping over all arguments
+    do i = 1, arg_num
+        call get_command_argument(i, argument)
+        select case(argument(scan(argument,'.',.true.)+1:))
+        case('cac')
+            lmpcac = .true.
+        case default
+            continue
+        end select
+    end do 
+
     !Determine if a mode is being used and what it is. The first argument has to be the mode
     !if a mode is being used
     call get_command_argument(1, argument)

src/mode_create.f90

@@ -12,7 +12,7 @@ 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), duplicate(3)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), adjustVar(3,8)
     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
@@ -299,7 +299,7 @@ module mode_create
         !Internal variables
         integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), &
                    vlat(3), temp_lat(3,8), m, n, o
-        real(kind=dp) :: v(3), temp_nodes(3,1,8)
+        real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8)
         real(kind=dp), allocatable :: resize_lat_array(:,:)
         logical, allocatable :: lat_points(:,:,:)
         logical :: node_in_bd(8)
@@ -307,6 +307,23 @@ module mode_create
         !Do some value initialization
         max_esize = esize
 
+        !If we are writing out to lammpscac format then we have to adjust the nodal positions
+
+        if(lmpcac) then 
+             do inod = 1, 8
+                do i = 1,3
+                    if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                        adjustVar(i,inod) = -0.5_dp
+                    else 
+                        adjustVar(i, inod) = 0.5_dp
+                    end if
+                end do
+            end do
+
+            adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8)))
+        else
+            adjustVar(:,:)=0.0_dp
+        end if
         !First find the bounding lattice points (min and max points for the box in each dimension)
         numlatpoints = 1
         do i = 1, 3
@@ -324,7 +341,6 @@ module mode_create
             continue
         end select
 
-
         !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case 
         !and one for the regular case
         if (esize==2) then  
@@ -437,9 +453,11 @@ module mode_create
                         end do
 
                         if(all(node_in_bd)) then 
-                            lat_ele_num = lat_ele_num+1
-                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
-
+                            lat_ele_num = lat_ele_num+1 
+                            do inod = 1, 8
+                                r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
+                            end do
+                            
                             !Now set all the lattice points contained within an element to false
                             do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
                                 do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
@@ -505,4 +523,4 @@ module mode_create
     end subroutine error_message
 
 
-end module mode_create
+end module mode_create

src/parameters.f90

@@ -5,4 +5,6 @@ module parameters
     integer, parameter :: dp= selected_real_kind(15,307)
     real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
                                 lim_large = huge(1.0_dp)
+    logical, save :: lmpcac
+    
 end module parameters