started implementation of the edit keyword

src/crystalmod.f90

@@ -203,6 +203,7 @@ module crystalmod
      procedure :: cell_delaunay !< Transform to the Delaunay primitive cell
      procedure :: reorder_atoms !< reorder the atoms in the crystal/molecule
      procedure :: wholemols !< Re-assign atomic types to have an asymmetric unit with whole molecules
+     procedure :: delete_atoms !< Delete a list of atoms
 
      ! symmetry (symmetry)
      procedure :: sitesymm !< Determine the local-symmetry group symbol for a point
@@ -511,6 +512,12 @@ module subroutine wholemols(c,ti)
        class(crystal), intent(inout) :: c
        type(thread_info), intent(in), optional :: ti
      end subroutine wholemols
+     module subroutine delete_atoms(c,nat,iat,ti)
+       class(crystal), intent(inout) :: c
+       integer, intent(in) :: nat
+       integer, intent(in) :: iat(nat)
+       type(thread_info), intent(in), optional :: ti
+     end subroutine delete_atoms
      module function sitesymm(c,x0,eps0,leqv,lrotm)
        class(crystal), intent(in) :: c
        real*8, intent(in) :: x0(3)

src/crystalmod@edit.f90

@@ -697,4 +697,57 @@ module subroutine wholemols(c,ti)
 
   end subroutine wholemols
 
+  !> Delete the atoms with IDs in the array iat(1:nat) from the
+  !> structure.
+  module subroutine delete_atoms(c,nat,iat,ti)
+    use crystalseedmod, only: crystalseed
+    use types, only: realloc
+    class(crystal), intent(inout) :: c
+    integer, intent(in) :: nat
+    integer, intent(in) :: iat(nat)
+    type(thread_info), intent(in), optional :: ti
+
+    type(crystalseed) :: seed
+    logical, allocatable :: useatoms(:), usespcs(:)
+    integer :: i
+
+    ! make seed from this crystal
+    call c%makeseed(seed,copysym=.false.)
+
+    ! flag the atoms and species to use
+    allocate(useatoms(c%ncel),usespcs(c%nspc))
+    useatoms = .true.
+    do i = 1, nat
+       useatoms(iat(i)) = .false.
+    end do
+    usespcs = .false.
+    do i = 1, c%ncel
+       if (useatoms(i)) usespcs(c%atcel(i)%is) = .true.
+    end do
+
+    ! re-do the atom and species info
+    seed%nat = 0
+    do i = 1, c%ncel
+       if (useatoms(i)) then
+          seed%nat = seed%nat + 1
+          seed%x(:,seed%nat) = c%atcel(i)%x
+          seed%is(seed%nat) = c%atcel(i)%is
+       end if
+    end do
+    call realloc(seed%x,3,seed%nat)
+    call realloc(seed%is,seed%nat)
+    seed%nspc = 0
+    do i = 1, c%nspc
+       if (usespcs(i)) then
+          seed%nspc = seed%nspc + 1
+          seed%spc(seed%nspc) = c%spc(i)
+       end if
+    end do
+    call realloc(seed%spc,seed%nspc)
+
+    ! build the new crystal
+    call c%struct_new(seed,crashfail=.true.,ti=ti)
+
+  end subroutine delete_atoms
+
 end submodule edit

src/global@proc.F90

@@ -43,7 +43,7 @@ module subroutine critic_main()
     use struct_drivers, only: struct_crystal_input, struct_molcell,&
        struct_atomlabel, struct_sym, struct_sym, struct_charges, struct_write,&
        struct_powder, struct_rdf, struct_compare, struct_comparevc, struct_environ,&
-       struct_econ,&
+       struct_econ, struct_edit,&
        struct_coord, struct_polyhedra, struct_packing, struct_vdw, struct_identify,&
        struct_makemols_neighcrys, struct_molreorder, struct_molmove,&
        struct_kpoints, struct_bz, struct_newcell, struct_amd
@@ -93,6 +93,12 @@ module subroutine critic_main()
           end if
           call struct_crystal_input(subline,ismoli,.true.,.not.quiet,s0=sy)
 
+          ! newcell
+       elseif (equal(word,'edit')) then
+          call check_structure_defined(ok)
+          if (.not.ok) cycle
+          call struct_edit(sy,.not.quiet)
+
           ! newcell
        elseif (equal(word,'newcell')) then
           call check_structure_defined(ok)

src/struct_drivers.f90

@@ -41,6 +41,7 @@ module struct_drivers
   public :: struct_polyhedra
   public :: struct_packing
   public :: struct_vdw
+  public :: struct_edit
   public :: struct_newcell
   public :: struct_molcell
   public :: struct_identify
@@ -122,6 +123,10 @@ module subroutine struct_vdw(s,line)
        type(system), intent(in) :: s
        character*(*), intent(in) :: line
      end subroutine struct_vdw
+     module subroutine struct_edit(s,verbose)
+       type(system), intent(inout) :: s
+       logical, intent(in) :: verbose
+     end subroutine struct_edit
      module subroutine struct_newcell(s,line,verbose)
        type(system), intent(inout) :: s
        character*(*), intent(in) :: line

src/struct_drivers@proc.f90

@@ -1331,7 +1331,7 @@ module subroutine struct_compare(s,line)
     use tools_math, only: crosscorr_triangle, rmsd_walker, umeyama_graph_matching,&
        ullmann_graph_matching, emd, synthetic_powder
     use tools_io, only: getword, equal, faterr, ferror, uout, string, ioj_center,&
-       ioj_left, string, lower, lgetword, fopen_read, fclose, getline, isreal,&
+       ioj_left, string, lower, lgetword, fopen_read, fclose, isreal,&
        file_read_xy
     use types, only: realloc
     use param, only: isformat_unknown, maxzat, pi
@@ -3103,6 +3103,99 @@ module subroutine struct_vdw(s,line)
 
   end subroutine struct_vdw
 
+  !> Edit a molecular or crystal structure
+  module subroutine struct_edit(s,verbose)
+    use systemmod, only: system
+    use tools_io, only: uout, ucopy, uin, getline, lgetword, equal, ferror, faterr,&
+       getword, isinteger, string
+    use types, only: realloc
+    type(system), intent(inout) :: s
+    logical, intent(in) :: verbose
+
+    character(len=:), allocatable :: line, word, errmsg
+    integer :: lp
+    integer :: idx, nat, i
+    integer, allocatable :: iat(:)
+    logical :: ok
+
+    if (verbose) then
+       write (uout,'("* Edit the crystal or molecular structure (EDIT)")')
+       write (uout,*)
+    end if
+
+    ! transform to the primitive?
+    errmsg = "Invalid syntax"
+    do while (getline(uin,line,ucopy=ucopy))
+       lp = 1
+       word = lgetword(line,lp)
+       if (equal(word,"delete")) then
+          word = lgetword(line,lp)
+
+          nat = 0
+          allocate(iat(10))
+          if (equal(word,"atoms") .or. equal(word,"atom")) then
+             ! read the list of atoms
+             do while (.true.)
+                word = getword(line,lp)
+                if (len_trim(word) == 0) exit
+                ok = isinteger(idx,word)
+                if (.not.ok) then
+                   errmsg = "invalid atom number in delete atoms"
+                   goto 999
+                end if
+                if (idx < 1 .or. idx > s%c%ncel) then
+                   errmsg = "atomic ID " // string(idx) // " outside valid range"
+                   goto 999
+                end if
+                nat = nat + 1
+                if (nat > size(iat,1)) call realloc(iat,2*nat)
+                iat(nat) = idx
+             end do
+          elseif (equal(word,"molecules") .or. equal(word,"molecule")) then
+             ! read the list of molecules
+             do while (.true.)
+                word = getword(line,lp)
+                if (len_trim(word) == 0) exit
+                ok = isinteger(idx,word)
+                if (.not.ok) then
+                   errmsg = "invalid molecule number in delete atoms"
+                   goto 999
+                end if
+                if (idx < 1 .or. idx > s%c%nmol) then
+                   errmsg = "molecule ID " // string(idx) // " outside valid range"
+                   goto 999
+                end if
+                do i = 1, s%c%ncel
+                   if (s%c%idatcelmol(i) == idx) then
+                      nat = nat + 1
+                      if (nat > size(iat,1)) call realloc(iat,2*nat)
+                      iat(nat) = i
+                   end if
+                end do
+             end do
+
+             ! transform
+             call s%c%delete_atoms(nat,iat(1:nat))
+             call s%reset_fields()
+             if (verbose) &
+                call s%report(.true.,.true.,.true.,.true.,.true.,.true.,.false.)
+          else
+             goto 999
+          end if
+       elseif (equal(word,"end") .or. equal(word,"endedit")) then
+          exit
+       elseif (len_trim(word) > 0) then
+          goto 999
+       end if
+    end do
+
+    return
+999 continue
+
+    call ferror('struct_edit',errmsg,faterr,line,syntax=.true.)
+
+  end subroutine struct_edit
+
   !> Build a new crystal from the current crystal by cell transformation
   module subroutine struct_newcell(s,line,verbose)
     use systemmod, only: system