From 5949f04103656972d64f706f0a3c1ca723d772ef Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 16 Oct 2020 19:48:06 -0400 Subject: [PATCH 1/4] Working read from pycac.out file format --- src/Makefile | 2 +- src/call_mode.f90 | 3 +- src/call_option.f90 | 2 + src/elements.f90 | 69 +++++++++++++++++++-- src/io.f90 | 144 ++++++++++++++++++++++++++++++++++++++++++-- src/main.f90 | 5 +- src/parameters.f90 | 4 +- 7 files changed, 213 insertions(+), 16 deletions(-) diff --git a/src/Makefile b/src/Makefile index feb933f..c58e37c 100644 --- a/src/Makefile +++ b/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 diff --git a/src/call_mode.f90 b/src/call_mode.f90 index fb187c1..7eae059 100644 --- a/src/call_mode.f90 +++ b/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') diff --git a/src/call_option.f90 b/src/call_option.f90 index d55e1ee..aeb17b5 100644 --- a/src/call_option.f90 +++ b/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 diff --git a/src/elements.f90 b/src/elements.f90 index 6df859b..a12d263 100644 --- a/src/elements.f90 +++ b/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 diff --git a/src/io.f90 b/src/io.f90 index 7bbf3c2..e1ec449 100644 --- a/src/io.f90 +++ b/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 diff --git a/src/main.f90 b/src/main.f90 index 53fa4fa..b305138 100644 --- a/src/main.f90 +++ b/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 diff --git a/src/parameters.f90 b/src/parameters.f90 index f261552..825245d 100644 --- a/src/parameters.f90 +++ b/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 From 95e2ad0b4da6272945df58a91046c0b7373cf88c Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 19 Oct 2020 15:14:12 -0400 Subject: [PATCH 2/4] Current working changes to control-box code --- README.md | 7 ++++- src/box.f90 | 6 +++++ src/elements.f90 | 7 +++++ src/io.f90 | 65 ++++++++++++++++++++------------------------- src/opt_delete.f90 | 1 + src/subroutines.f90 | 61 +----------------------------------------- 6 files changed, 50 insertions(+), 97 deletions(-) diff --git a/README.md b/README.md index c83516c..036981d 100644 --- a/README.md +++ b/README.md @@ -58,7 +58,7 @@ Default duplicate is `1 1 1`. This is used to replicate the element along each d **Dimensions** ``` -dimensions dimx dimy dimz +dim dimx dimy dimz ``` There is no default dimensions as duplicate is the default option. This command assigns a box with user-assigned dimensions and fills it with the desired element. By default atoms fill in the jagged edges at the boundaries if the dimensions command is included. @@ -129,6 +129,11 @@ wrap This will wrap atomic positions back inside the box. Effectively as if periodic boundary conditions are applied so that atoms which exit from one side of the simulation cell enter back in through the other. +###Mode Metric +``` +--metric cmetric nfiles {filepaths} +``` + ## Options Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files. diff --git a/src/box.f90 b/src/box.f90 index faf268c..83d81a8 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -97,4 +97,10 @@ module box end do return end subroutine in_sub_box + + subroutine reset_box + !This subroutine just resets the box boundary and position + box_bc = "ppp" + box_bd(:) = 0 + end subroutine reset_box end module box diff --git a/src/elements.f90 b/src/elements.f90 index a12d263..ab1a60e 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -793,4 +793,11 @@ module elements virial_node(:,:,:,ie) = virial return end subroutine add_element_data + + subroutine reset_data + !This function resets all of the data arrays for the elements and atoms + atom_num = 0 + ele_num = 0 + node_num = 0 + end subroutine reset_data end module elements diff --git a/src/io.f90 b/src/io.f90 index e1ec449..556ce1f 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -583,49 +583,41 @@ module io temp_infile = filename end if - !Infinite loop which only exists if user provides valid filetype - do while(.true.) + !Check to see if file exists, if it does then ask user if they would like to overwrite the file + inquire(file=trim(temp_infile), exist=file_exists) + if (.not.file_exists) then + print *, "The file ", trim(adjustl(filename)), " does not exist. Please input an existing file to read." + stop 3 + end if - !Check to see if file exists, if it does then ask user if they would like to overwrite the file - inquire(file=trim(temp_infile), exist=file_exists) - if (.not.file_exists) then - print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists" - read(*,*) temp_infile - cycle + select case(temp_infile(scan(temp_infile,'.',.true.)+1:)) + case('restart', 'mb', 'cac') + infilenum=infilenum+1 + infiles(infilenum) = temp_infile + 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(temp_infile(scan(temp_infile,'.',.true.)+1:)) - case('restart', 'mb', 'cac') - 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, cac, or out." - read(*,*) temp_infile - + select case(trim(adjustl(mode))) + case('--convert','--metric') + infilenum = infilenum+1 + infiles(infilenum) = temp_infile + case default + print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode)) + stop 3 end select - end do + + 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, cac, or out." + stop 3 + end select end subroutine get_in_file subroutine read_in(i, displace, temp_box_bd) - !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine + !This subroutine reads in file i integer, intent(in) :: i real(kind=dp), dimension(3), intent(in) :: displace @@ -1067,6 +1059,7 @@ module io end do call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re) call add_element_data(ele_num, ee, fe, ve) + node_num = node_num + 8 end do end if call set_max_esize diff --git a/src/opt_delete.f90 b/src/opt_delete.f90 index 09cefb1..c3bcc9e 100644 --- a/src/opt_delete.f90 +++ b/src/opt_delete.f90 @@ -3,6 +3,7 @@ module opt_delete use parameters use subroutines use elements + use neighbors implicit none diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 1015968..b272c16 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -198,66 +198,6 @@ module subroutines end do end subroutine - subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) - !This subroutine builds a cell list based on rc_off - - !----------------------------------------Input/output variables------------------------------------------- - - integer, intent(in) :: numinlist !The number of points within r_list - - real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of - !the cell list. - - real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells - - integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. - - integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell - - integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. - - integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list - - !----------------------------------------Begin Subroutine ------------------------------------------- - - integer :: i, j, cell_lim, c(3) - real(kind=dp) :: box_len(3) - integer, allocatable :: resize_cell_list(:,:,:,:) - - !First calculate the number of cells that we need in each dimension - do i = 1,3 - 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 - - !Initialize/allocate variables - cell_lim = 10 - allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) - - !Now place points within cell - do i = 1, numinlist - !c is the position of the cell that the point belongs to - do j = 1, 3 - c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 - end do - - !Place the index in the correct position, growing if necessary - num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 - if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then - allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) - resize_cell_list(1:cell_lim, :, :, :) = cell_list - resize_cell_list(cell_lim+1:, :, :, :) = 0 - call move_alloc(resize_cell_list, cell_list) - end if - - cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i - which_cell(:,i) = c - end do - - return - end subroutine build_cell_list - - subroutine check_right_ortho(ori, isortho, isrighthanded) !This subroutine checks whether provided orientations in the form: ! | x1 x2 x3 | @@ -301,4 +241,5 @@ module subroutines return end subroutine check_right_ortho + end module subroutines From 6e085176975f908e1e5421cceb76892b8edc4fed Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 20 Oct 2020 01:03:31 -0400 Subject: [PATCH 3/4] Working continuum metric calculation code --- src/Makefile | 5 +- src/call_mode.f90 | 3 + src/elements.f90 | 10 +- src/functions.f90 | 120 ++++++++++++++++++++++ src/io.f90 | 4 +- src/main.f90 | 12 ++- src/mode_metric.f90 | 245 ++++++++++++++++++++++++++++++++++++++++++++ src/neighbors.f90 | 142 +++++++++++++++++++++++++ 8 files changed, 531 insertions(+), 10 deletions(-) create mode 100644 src/mode_metric.f90 create mode 100644 src/neighbors.f90 diff --git a/src/Makefile b/src/Makefile index c58e37c..6106104 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,9 +1,9 @@ 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 mode_metric.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 +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o neighbors.o $(MODES) $(OPTIONS) call_option.o sorts.o .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -42,3 +42,4 @@ elements.o : sorts.o $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o +opt_delete.o mode_metric.o : neighbors.o diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 7eae059..8d73e25 100644 --- a/src/call_mode.f90 +++ b/src/call_mode.f90 @@ -5,6 +5,7 @@ subroutine call_mode(arg_pos) use mode_create use mode_convert use mode_merge + use mode_metric use parameters implicit none @@ -18,6 +19,8 @@ subroutine call_mode(arg_pos) call convert(arg_pos) case('--merge') call merge(arg_pos) + case('--metric') + call metric(arg_pos) case default print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", & "accepted modes and rerun." diff --git a/src/elements.f90 b/src/elements.f90 index ab1a60e..1254a43 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -747,24 +747,30 @@ module elements !Allocate element arrays if (n > 0) then + if (allocated(force_node)) then + deallocate(force_node, virial_node, energy_node) + end if 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 + print *, "Error allocating element data arrays in mode_metric because of:", allostat stop end if end if if (m > 0) then + if (allocated(force_atom)) then + deallocate(force_atom, virial_atom, energy_atom) + end if 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 + print *, "Error allocating atom data arrays in mode_metric because of:", allostat stop end if end if diff --git a/src/functions.f90 b/src/functions.f90 index 6fed2b7..0a8bc49 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -271,4 +271,124 @@ END FUNCTION StrDnCase norm_dis(1:3) = (rk - rl) norm_dis(4) = norm2(rk-rl) end function + + pure function matinv3(A) result(B) + !! Performs a direct calculation of the inverse of a 3×3 matrix. + real(kind=dp), intent(in) :: A(3,3) !! Matrix + real(kind=dp) :: B(3,3) !! Inverse matrix + real(kind=dp) :: detinv + + if(abs(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) < lim_zero) then + B(:,:) = 0 + return + else + ! Calculate the inverse determinant of the matrix + + detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + ! Calculate the inverse of the matrix + B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2)) + B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1)) + B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1)) + B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2)) + B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1)) + B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1)) + B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2)) + B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1)) + B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1)) + end if + end function + + pure function transpose3(A) result(B) + !!Transposes matrix A + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + integer :: i, j + forall(i =1:3,j=1:3) B(i,j) = A(j,i) + + end function transpose3 + + pure function sqrt3(A) result(B) + !This calculates the square of matrix A fulfilling the equation B*B = A according to + !the algorithm by Franca1989 + + real(kind=dp), intent(in) :: A(3,3) + real(kind=dp) :: B(3,3) + + real(kind=dp) :: Ione, Itwo, Ithree, l, k, phi, Asq(3,3), lambda, Bone, Btwo, Bthree, p + + !Step 1 is calculating the invariants of C + Ione = A(1,1) + A(2,2) + A(3,3) + Asq = matmul(A,A) + Itwo = 0.5_dp *(Ione*Ione - (Asq(1,1) + Asq(2,2) + Asq(3,3))) + Ithree = (A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& + - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& + + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) + + if (Ithree < 0) then + B(:,:)=0.0_dp + return + end if + !Check for an isotropic matrix + k = Ione*Ione - 3*Itwo + if (k < lim_zero) then + lambda = sqrt(Ione/3.0_dp) + B = lambda*identity_mat(3) + else + l = Ione*(Ione*Ione - 9.0_dp/2.0_dp * Itwo) + 27.0_dp/2.0_dp * Ithree + p = l/(k**(1.5_dp)) + + if (p > 1.0 ) then + B(:,:) = 0.0_dp + return + end if + + if ((p< -1).or.(p>1)) then + B(:,:) = 0.0_dp + return + end if + phi = acos(p) + lambda = sqrt(1.0_dp/3.0_dp * (Ione + 2*sqrt(k)*cos(phi/3))) + + !Now calculate invariantes of B + Bthree = sqrt(Ithree) + if((-lambda*lambda + Ione + 2*Ithree/lambda) < 0) then + B(:,:) = 0.0_dp + return + end if + Bone = lambda + sqrt(- lambda*lambda + Ione + 2*Ithree/lambda) + Btwo = (Bone*Bone - Ione)/2.0_dp + + !Now calculate B + if(abs(Bone*Btwo -Bthree) < lim_zero) then + B(:,:) = 0.0_dp + return + end if + B = 1/(Bone*Btwo - Bthree) *(Bone*Bthree*identity_mat(3) + (Bone*Bone - Btwo)*A - matmul(A,A)) + end if + end function sqrt3 + + pure function permutation(i,j,k) result(e) + !Calculates the permutation tensor + integer, intent(in) :: i,j,k + integer :: e + + if ( ((i==1).and.(j==2).and.(k==3)).or. & + ((i==2).and.(j==3).and.(k==1)).or. & + ((i==3).and.(j==1).and.(k==2))) then + e=1 + else if( ((i==2).and.(j==1).and.(k==3)).or. & + ((i==1).and.(j==3).and.(k==2)).or. & + ((i==3).and.(j==2).and.(k==1))) then + e=-1 + else + e=0 + end if + end function permutation + end module functions diff --git a/src/io.f90 b/src/io.f90 index 556ce1f..9c02ba1 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -984,6 +984,7 @@ module io 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 + character(len=1000) :: line open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') @@ -1033,7 +1034,8 @@ module io !Read atom header read(11,*) textholder do ia = 1, in_atoms - read(11,*) tag, type, ra(:), ea, fa(:), va(:) + read(11,'(a)') line(:) + read(line,*) 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 diff --git a/src/main.f90 b/src/main.f90 index b305138..9f7a484 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -111,11 +111,13 @@ program main if(bound_called) call def_new_box !Check to make sure a file was passed to be written out and then write out - ! Before building do a check on the file - if (outfilenum == 0) then - argument = 'none' - call get_out_file(argument) + ! Before building do a check on the file + if (trim(adjustl(mode)) /= "--metric") then + if ((outfilenum == 0)) then + argument = 'none' + call get_out_file(argument) + end if + call write_out end if - call write_out end program main diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 new file mode 100644 index 0000000..279e710 --- /dev/null +++ b/src/mode_metric.f90 @@ -0,0 +1,245 @@ +module mode_metric + !This mode is used to calculate continuum metrics for the j + + use parameters + use io + use elements + use neighbors + + implicit none + + integer :: nfiles + character(len=100) :: metric_type + real(kind=dp), allocatable :: met(:,:) + + !Save reference positions + integer :: np, npreal, nmet + real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:) + + public + contains + subroutine metric(arg_pos) + !This is the main calling subroutine for the metric code + integer, intent(out) :: arg_pos + character(len=100) :: infile, outfile + + integer :: i, ibasis, inod, np_temp, ppos + real(kind=dp), dimension(6) :: temp_box_bd + + !These are the variables containing the cell list information + integer, dimension(3) :: cell_num + integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:) + integer, allocatable :: cell_list(:,:,:,:) + + !Parse the command arguments + call parse_command(arg_pos) + + !Now read the first file + call read_in(1, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + np = atom_num + max_basisnum*max_ng_node*ele_num + print *,np + allocate(r_zero(3,atom_num+max_basisnum*max_ng_node*ele_num), & + r_curr(3,atom_num+max_basisnum*max_ng_node*ele_num)) + r_zero(:,:) = -huge(1.0_dp) + + !Set up the met variable for the user desired metric + select case(trim(adjustl(metric_type))) + case('def_grad') + allocate(met(9, np)) + case('microrotation') + allocate(met(3,np)) + end select + + !Now set the reference positions + call convert_positions(r_zero, npreal) + + !Now calculate the neighbor list for the reference configuration + call calc_neighbor(5.0_dp, r_zero, np) + + !Reset element and box + call reset_data + call reset_box + + !Now loop over new files + do i = 2, nfiles + call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) + call convert_positions(r_curr, np_temp) + if (npreal /= np_temp) then + print *, "Error in mode_metric where number of points in ", i, "th file is ", np_temp, " and number of points in" & + // "reference file is", npreal + end if + call calc_metric + !Now create the output file num and write out to xyz format + ppos = scan(trim(infiles(i)),".", BACK= .true.) + if ( ppos > 0 ) then + outfile = infiles(i)(1:ppos)//'xyz' + else + outfile = infiles(i)//'.xyz' + end if + call write_metric_xyz(outfile) + call reset_data + call reset_box + end do + end subroutine metric + + subroutine parse_command(arg_pos) + !This subroutine parses the arguments for mode metric + integer, intent(out) :: arg_pos + + integer :: i, arglen + character(len=100) :: textholder + logical :: file_exists + + !First read the metric to be used + call get_command_argument(2,metric_type,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + select case(trim(adjustl(metric_type))) + case("microrotation") + continue + case default + print *, "Mode metric does not accept metric ", metric_type, ". Please select from: microrotation" + stop 3 + end select + + !Now read the number of files to read and allocate the variables + call get_command_argument(3, textholder, arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) nfiles + + !Now read the files to be read + do i = 1, nfiles + call get_command_argument(3+i, textholder, arglen) + call get_in_file(textholder) + end do + + arg_pos = 4+nfiles + return + end subroutine parse_command + + subroutine calc_metric + !This subroutine calculates the continuum metric that we require + integer :: i, j, k, nei, ip, jp + real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), & + U(3,3), R(3,3), Rskew(3,3), oldrij(3) + + !Loop over all points + do ip = 1, np + eta(:,:) = 0.0_dp + omega(:,:) = 0.0_dp + def_grad(:,:) = 0.0_dp + do jp = 1, nei_num(ip) + !Calculate the neighbor vec in current configuration + nei = nei_list(jp, ip) + rij = r_curr(:,nei) - r_curr(:,ip) + oldrij = r_zero(:,nei) - r_zero(:,ip) + + !Calculate eta and omega + do i = 1,3 + do j = 1,3 + omega(i,j) = omega(i,j) + rij(i) * oldrij(j) + eta(i,j) = eta(i,j) + oldrij(i) * oldrij(j) + end do + end do + end do + + eta_inv=matinv3(eta) + def_grad=matmul(omega,eta_inv) + + select case(trim(adjustl(metric_type))) + case('def_grad') + k = 1 + do i = 1,3 + do j = 1, 3 + met(k, ip) = def_grad(i,j) + end do + end do + case('microrotation') + met(:,ip) = 0.0_dp + if(.not.all(def_grad == 0)) then + !Now calculate microrotation + ftf = matmul(transpose3(def_grad), def_grad) + U = sqrt3(ftf) + if(.not.all(abs(U) < lim_zero)) then + R = matmul(def_grad, matinv3(U)) + Rskew = 0.5_dp * ( R - transpose3(R)) + do k =1,3 + do j = 1,3 + do i = 1,3 + met(k,ip) = met(k,ip) -0.5*permutation(i,j,k)*Rskew(i,j) + end do + end do + end do + end if + end if + end select + end do + return + end subroutine + + subroutine convert_positions(rout, npoints) + !This subroutine just converts current atom and element arrays to a single point based form + real(kind=dp), dimension(3,atom_num+max_ng_node*max_basisnum*ele_num), intent(inout) :: rout + integer, intent(out) :: npoints + + integer :: i, inod, ibasis + + npoints=0 + print *, atom_num + max_ng_node*max_basisnum*ele_num + print *, rout(:,1) + + if(atom_num > 0) then + do i = 1, atom_num + rout(:,tag_atom(i)) = r_atom(:,i) + npoints = npoints + 1 + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) & + = r_node(:,ibasis,inod,i) + npoints = npoints + 1 + end do + end do + end do + end if + + end subroutine convert_positions + + subroutine write_metric_xyz(outfile) + character(len=100), intent(in) :: outfile + + integer :: i, inod, ibasis + real(kind = dp) :: r(3), eng + open (unit=11, file=trim(adjustl(outfile)), action='write', position='rewind', status = 'replace') + !Write the header + write(11,*) npreal + + select case(metric_type) + case('def_grad') + write(11,*) "type x y z F11 F12 F13 F21 F22 F23 F31 F32 F33" + case('microrotation') + write(11,*) "type x y z micro1 micro2 micro3" + end select + + if(atom_num > 0) then + do i = 1, atom_num + write(11,*) type_atom(i), r_atom(:,i), met(:,tag_atom(i)) + end do + end if + + if (ele_num > 0) then + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11,*) basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i), & + met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) + end do + end do + end do + end if + end subroutine write_metric_xyz + +end module mode_metric diff --git a/src/neighbors.f90 b/src/neighbors.f90 new file mode 100644 index 0000000..0faebca --- /dev/null +++ b/src/neighbors.f90 @@ -0,0 +1,142 @@ +module neighbors + + use parameters + use elements + use subroutines + use functions + + integer, allocatable :: nei_list(:,:), nei_num(:) + real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:) + public + contains + + subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) + !This subroutine builds a cell list based on rc_off + + !----------------------------------------Input/output variables------------------------------------------- + + integer, intent(in) :: numinlist !The number of points within r_list + + real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of + !the cell list. + + real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells + + integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. + + integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell + + integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. + + integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list + + !----------------------------------------Begin Subroutine ------------------------------------------- + + integer :: i, j, cell_lim, c(3) + real(kind=dp) :: box_len(3) + integer, allocatable :: resize_cell_list(:,:,:,:) + + !First calculate the number of cells that we need in each dimension + do i = 1,3 + 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 + + !Initialize/allocate variables + cell_lim = 10 + allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) + + !Now place points within cell + 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 + + !c is the position of the cell that the point belongs to + do j = 1, 3 + c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 + end do + + !Place the index in the correct position, growing if necessary + num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 + if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then + allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) + resize_cell_list(1:cell_lim, :, :, :) = cell_list + resize_cell_list(cell_lim+1:, :, :, :) = 0 + call move_alloc(resize_cell_list, cell_list) + end if + + cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i + which_cell(:,i) = c + end do + + return + end subroutine build_cell_list + + subroutine calc_neighbor(rc_off, r_list, n) + !This function populates the neighbor list in this module + + real(kind=dp), intent(in) :: rc_off + integer, intent(in) :: n + real(kind=dp), dimension(3,n) :: r_list + + integer :: i, c(3), ci, cj, ck, num_nei, nei + !Variables for cell list code + integer, dimension(3) ::cell_num + integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) + integer :: which_cell(3,n) + + !First reallocate the neighbor list codes + if (allocated(nei_list)) then + deallocate(nei_list,nei_num) + end if + + allocate(nei_list(100,n),nei_num(n)) + + !Now first pass the position list and and point num to the cell algorithm + call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) + + !Now loop over every point and find it's 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 + + !c is the position of the cell that the point + c = which_cell(:,i) + + !Loop over all neighboring cells + do ci = -1, 1, 1 + do cj = -1, 1, 1 + do ck = -1, 1, 1 + !First check to make sure that the neighboring cell exists + if(any((c + (/ ck, cj, ci /)) == 0)) cycle + if( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. & + (c(3) + ci > cell_num(3))) cycle + + do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci) + nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) + + !Check to make sure the atom isn't the same index as the atom we are checking + !and that the neighbor hasn't already been deleted + if((nei /= i)) then + + !Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that + !atom and calculate the initial neighbor vector + if (norm2(r_list(:,nei)-r_list(:,i)) < rc_off) then + + nei_num(i) = nei_num(i) + 1 + nei_list(nei_num(i), i) = nei + + end if + end if + end do + end do + end do + end do + + end do pointloop + + return + end subroutine calc_neighbor + +end module neighbors From 2ea388b82ad3c49ab5cb0c7c7e5699e68369fb4c Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 23 Oct 2020 11:37:28 -0400 Subject: [PATCH 4/4] Current working changes with some updates to comments for accuracy --- src/elements.f90 | 5 ++--- src/mode_metric.f90 | 18 ++++++++++-------- src/opt_disl.f90 | 1 - src/opt_group.f90 | 2 +- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 1254a43..226eb60 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -501,6 +501,7 @@ module elements !We go from largest index to smallest index just to make sure that we don't miss !accidentally overwrite values which need to be deleted do i = num, 1, -1 + node_num = node_num - ng_node(lat_ele(sorted_index(i))) if(sorted_index(i) == ele_num) then r_node(:,:,:,sorted_index(i)) = 0.0_dp type_ele(sorted_index(i)) ='' @@ -509,7 +510,6 @@ module elements sbox_ele(sorted_index(i)) = 0 tag_ele(sorted_index(i)) = 0 else - node_num = node_num - ng_node(lat_ele(sorted_index(i))) r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num) type_ele(sorted_index(i)) = type_ele(ele_num) size_ele(sorted_index(i)) = size_ele(ele_num) @@ -538,8 +538,7 @@ module elements max_bd(:) = -huge(1.0_dp) min_bd(:) = huge(1.0_dp) - - do i = 1, atom_num +do i = 1, atom_num do j = 1, 3 if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 index 279e710..10838c6 100644 --- a/src/mode_metric.f90 +++ b/src/mode_metric.f90 @@ -10,11 +10,11 @@ module mode_metric integer :: nfiles character(len=100) :: metric_type - real(kind=dp), allocatable :: met(:,:) + real(kind=dp) :: rc_off !Save reference positions integer :: np, npreal, nmet - real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:) + real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:), met(:,:) public contains @@ -37,7 +37,6 @@ module mode_metric !Now read the first file call read_in(1, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) np = atom_num + max_basisnum*max_ng_node*ele_num - print *,np allocate(r_zero(3,atom_num+max_basisnum*max_ng_node*ele_num), & r_curr(3,atom_num+max_basisnum*max_ng_node*ele_num)) r_zero(:,:) = -huge(1.0_dp) @@ -101,18 +100,23 @@ module mode_metric stop 3 end select + !Now read the cutoff radius + call get_command_argument(3,textholder,arglen) + if (arglen == 0) stop "Incomplete mode metric command, check documentation" + read(textholder, *) rc_off + !Now read the number of files to read and allocate the variables - call get_command_argument(3, textholder, arglen) + call get_command_argument(4, textholder, arglen) if (arglen == 0) stop "Incomplete mode metric command, check documentation" read(textholder, *) nfiles !Now read the files to be read do i = 1, nfiles - call get_command_argument(3+i, textholder, arglen) + call get_command_argument(4+i, textholder, arglen) call get_in_file(textholder) end do - arg_pos = 4+nfiles + arg_pos = 5+nfiles return end subroutine parse_command @@ -184,8 +188,6 @@ module mode_metric integer :: i, inod, ibasis npoints=0 - print *, atom_num + max_ng_node*max_basisnum*ele_num - print *, rout(:,1) if(atom_num > 0) then do i = 1, atom_num diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 19f8ed8..dc0e72c 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -561,7 +561,6 @@ module opt_disl !Now reset the list for the scanning algorithm delete_num = 0 - !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom !that reside on the same plane. If loop_radius is > 0 then we build a circular vacancy cluster if(loop_radius > 0) then do i = 1, atom_num diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 855e04e..33cf685 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -87,7 +87,7 @@ module opt_group continue case default print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", & - "Please select from atoms, nodes, or both." + "Please select from atoms, elements, or both." end select arg_pos = arg_pos + 1