added the types option to discard; support multiple discards

src/autocp@proc.f90

@@ -69,7 +69,7 @@ module subroutine autocritic(line)
     use tools, only: uniqc
     use tools_io, only: uout, ferror, faterr, lgetword, equal, isexpression_or_word,&
        string, warning, tictac
-    use types, only: realloc
+    use types, only: realloc, discard_cp_expr
     use param, only: pi
     character*(*), intent(in) :: line
 
@@ -100,7 +100,7 @@ module subroutine autocritic(line)
     real*8  :: iniv(4,3), xdum(3)
     integer :: nt
     logical :: ok, dryrun, dochk
-    character(len=:), allocatable :: word, str, discexpr
+    character(len=:), allocatable :: word, str
     real*8 :: gfnormeps
     integer :: ntetrag
     real*8, allocatable :: tetrag(:,:,:)
@@ -120,7 +120,8 @@ module subroutine autocritic(line)
     real*8 :: nucepsh
     type(grhandle) :: gr
     type(mesh) :: meshseed
-    logical :: typeok(4)
+    logical :: typeok(4), laux(4)
+    type(discard_cp_expr), allocatable :: discard(:)
 
     real*8, parameter :: gradeps_check = 1d-4 ! minimum gradeps requirement for addcp (grids)
 
@@ -160,7 +161,6 @@ module subroutine autocritic(line)
 
     ! parse the input
     lp = 1
-    discexpr = ""
     typeok = .true.
     do while (.true.)
        word = lgetword(line,lp)
@@ -344,11 +344,36 @@ module subroutine autocritic(line)
              endif
           end do
        elseif (equal(word,"discard")) then
-          ok = isexpression_or_word(discexpr,line,lp)
+          laux = .true.
+          lpo = lp
+          word = trim(lgetword(line,lp))
+          if (equal(word,"types")) then
+             laux = .false.
+             word = trim(lgetword(line,lp))
+             do i = 1, len(word)
+                if (word(i:i) == "n") then
+                   laux(1) = .true.
+                elseif (word(i:i) == "b") then
+                   laux(2) = .true.
+                elseif (word(i:i) == "r") then
+                   laux(3) = .true.
+                elseif (word(i:i) == "c") then
+                   laux(4) = .true.
+                else
+                   call ferror('autocritic','TYPES input letters can only be one of n, b, r, c',faterr,line,syntax=.true.)
+                   return
+                end if
+             end do
+          else
+             lp = lpo
+          end if
+
+          ok = isexpression_or_word(str,line,lp)
           if (.not. ok) then
              call ferror("autocritic","wrong DISCARD keyword",faterr,line,syntax=.true.)
              return
           end if
+          call add_discard_expression(str,laux)
        elseif (equal(word,"types")) then
           typeok = .false.
           word = trim(lgetword(line,lp))
@@ -362,7 +387,6 @@ module subroutine autocritic(line)
              elseif (word(i:i) == "c") then
                 typeok(4) = .true.
              else
-                write (*,*) "bleh! ->", word(i:i), "<-"
                 call ferror('autocritic','TYPES input letters can only be one of n, b, r, c',faterr,line,syntax=.true.)
                 return
              end if
@@ -560,8 +584,16 @@ module subroutine autocritic(line)
        string(nucepsh*dunit0(iunit),'e',decimal=3), iunitname0(iunit)
     write (uout,'("  CPs are degenerate if any Hessian element abs value is less than: ",A)') &
        string(CP_hdegen,'e',decimal=3)
-    if (len_trim(discexpr) > 0) &
-       write (uout,'("  Discard CP expression: ",A)') trim(discexpr)
+    if (allocated(discard)) then
+       do i = 1, size(discard,1)
+          if (all(discard(i)%typeok)) then
+             write (uout,'("  Discard CP expression: ",A)') trim(discard(i)%s)
+          else
+             write (uout,'("  Discard CP expression: ",A," applies to: n=",A,", b=",A,", r=",A,", c=",A)') &
+                trim(discard(i)%s), (string(discard(i)%typeok(j)),j=1,4)
+          end if
+       end do
+    end if
     if (.not.all(typeok)) then
        write (uout,'("  CP types to keep: nuclei=",A,", bonds=",A,", rings=",A,", cages=",A)') &
           (string(typeok(j)),j=1,4)
@@ -760,7 +792,7 @@ module subroutine autocritic(line)
 
              if (ok) then
                 !$omp critical (addcp)
-                call sy%addcp(sy%iref,x0,discexpr,cpeps,nuceps,nucepsh,max(gradeps_check,gfnormeps),&
+                call sy%addcp(sy%iref,x0,discard,cpeps,nuceps,nucepsh,max(gradeps_check,gfnormeps),&
                    typeok=typeok)
                 !$omp end critical (addcp)
              end if
@@ -823,6 +855,30 @@ module subroutine autocritic(line)
        write (uout,*)
     end if
 
+  contains
+    subroutine add_discard_expression(str,type)
+      character(len=*), intent(in) :: str
+      logical, intent(in) :: type(4)
+
+      integer :: n
+      type(discard_cp_expr), allocatable :: temp(:)
+
+      if (.not.allocated(discard)) then
+         allocate(discard(1))
+         n = 1
+      else
+         n = size(discard,1)
+         allocate(temp(n+1))
+         temp(1:n) = discard
+         call move_alloc(temp,discard)
+         n = n + 1
+      end if
+
+      discard(n)%s = str
+      discard(n)%typeok = type
+
+    end subroutine add_discard_expression
+
   end subroutine autocritic
 
   !> Report the results of the critical point search.

src/systemmod.f90

@@ -19,7 +19,7 @@
 module systemmod
   use grid1mod, only: grid1
   use hashmod, only: hash
-  use types, only: integrable, pointpropable
+  use types, only: integrable, pointpropable, discard_cp_expr
   use fieldmod, only: field
   use crystalmod, only: crystal
   use types, only: thread_info
@@ -226,11 +226,11 @@ module subroutine grdall(s,xpos,lprop,pmask)
        real*8, intent(out) :: lprop(s%npropi)
        logical, intent(in), optional :: pmask(s%npropi)
      end subroutine grdall
-     module subroutine addcp(s,id,x0,discexpr,cpeps,nuceps,nucepsh,gfnormeps,itype,typeok)
+     module subroutine addcp(s,id,x0,discard,cpeps,nuceps,nucepsh,gfnormeps,itype,typeok)
        class(system), intent(inout) :: s
        integer, intent(in) :: id
        real*8, intent(in) :: x0(3)
-       character*(*), intent(in) :: discexpr
+       type(discard_cp_expr), allocatable, intent(in) :: discard(:)
        real*8, intent(in) :: cpeps
        real*8, intent(in) :: nuceps
        real*8, intent(in) :: nucepsh

src/systemmod@proc.f90

@@ -1375,13 +1375,13 @@ end subroutine grdall
   !> length, discard the critical point if it gives a non-zero value
   !> for this expression. If typeok is present, only accept a CP type
   !> if the corresponding element in typeok is true.
-  module subroutine addcp(s,id,x0,discexpr,cpeps,nuceps,nucepsh,gfnormeps,itype,typeok)
+  module subroutine addcp(s,id,x0,discard,cpeps,nuceps,nucepsh,gfnormeps,itype,typeok)
     use tools_io, only: faterr, ferror
     use types, only: scalar_value
     class(system), intent(inout) :: s
     integer, intent(in) :: id
     real*8, intent(in) :: x0(3) !< Position of the CP, in Cartesian coordinates
-    character*(*), intent(in) :: discexpr !< Discard expression
+    type(discard_cp_expr), allocatable, intent(in) :: discard(:) !< Discard expressions
     real*8, intent(in) :: cpeps !< Discard CPs closer than cpeps from other CPs
     real*8, intent(in) :: nuceps !< Discard CPs closer than nuceps from nuclei
     real*8, intent(in) :: nucepsh !< Discard CPs closer than nucepsh from hydrogen
@@ -1390,32 +1390,30 @@ module subroutine addcp(s,id,x0,discexpr,cpeps,nuceps,nucepsh,gfnormeps,itype,ty
     logical, intent(in), optional :: typeok(4) !< Filter out by CP type
 
     real*8 :: fval
-    logical :: ok, isexpr, istypeok
+    logical :: ok
     character(len=:), allocatable :: errmsg
     type(scalar_value) :: res
-    integer :: idx
+    integer :: idx, i
 
-    ! flags
-    istypeok = present(typeok)
-    if (istypeok) istypeok = .not.all(typeok)
-    isexpr = (len_trim(discexpr) > 0)
-
-    ! if type filter and discard expression are both on, pre-filter
-    if (istypeok .and. isexpr) then
+    ! if discard expression is on
+    if (allocated(discard)) then
+       ! pre-filter by type here
        call s%f(id)%grd(x0,2,res)
        idx = (res%s+5)/2
        if (idx >= 1 .and. idx <= 4) then
           if (.not.typeok(idx)) return
        end if
-    end if
 
-    ! apply the discard expression
-    if (isexpr) then
-       fval = s%eval(discexpr,errmsg,x0)
-       if (len_trim(errmsg) > 0) &
-          call ferror("addcp","Invalid DISCARD expression: " // trim(errmsg),faterr)
-       ok = (abs(fval) < 1d-30)
-       if (.not.ok) return
+       ! apply the discard expressions
+       do i = 1, size(discard,1)
+          if (discard(i)%typeok(idx)) then
+             fval = s%eval(discard(i)%s,errmsg,x0)
+             if (len_trim(errmsg) > 0) &
+                call ferror("addcp","Invalid DISCARD expression: " // trim(errmsg),faterr)
+             ok = (abs(fval) < 1d-30)
+             if (.not.ok) return
+          end if
+       end do
     end if
 
     ! add the CP and write down the info

src/types.f90

@@ -37,6 +37,7 @@ module types
   public :: gpathp
   public :: basindat
   public :: int_result
+  public :: discard_cp_expr
 
   ! overloaded procedures
   interface realloc
@@ -307,6 +308,12 @@ module types
      real*8, allocatable :: fa3(:,:,:,:,:,:) ! three-center DIs
   end type int_result
 
+  !> Expressions for DISCARD in CP search
+  type discard_cp_expr
+     character(len=:), allocatable :: s
+     logical :: typeok(4)
+  end type discard_cp_expr
+
   interface
      module subroutine scalar_value_clear(s)
        class(scalar_value), intent(inout) :: s