diff --git a/src/autocp@proc.f90 b/src/autocp@proc.f90 index 6db5362b..3e91e78e 100644 --- a/src/autocp@proc.f90 +++ b/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. diff --git a/src/systemmod.f90 b/src/systemmod.f90 index 069c1edf..ba4ef14e 100644 --- a/src/systemmod.f90 +++ b/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 diff --git a/src/systemmod@proc.f90 b/src/systemmod@proc.f90 index ac311f72..ff45b54c 100644 --- a/src/systemmod@proc.f90 +++ b/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 diff --git a/src/types.f90 b/src/types.f90 index f0b981a6..d543b58e 100644 --- a/src/types.f90 +++ b/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