Allow printing GFN-FF topology lists (#545)

Signed-off-by: Thomas Rose <39367840+Thomas3R@users.noreply.github.com>

src/gfnff/topology.f90

@@ -22,7 +22,7 @@ module xtb_gfnff_topology
    implicit none
    private
 
-   public :: TGFFTopology
+   public :: TGFFTopology, TPrintTopo
 
 
    !> Topology information for a given system
@@ -96,6 +96,23 @@ module xtb_gfnff_topology
 
    end type TGFFTopology
 
+   ! logicals for GFN-FF topology list printout
+   type :: TPrintTopo
+     logical :: nb      = .false.
+     logical :: bpair   = .false.
+     logical :: alist   = .false.
+     logical :: blist   = .false.
+     logical :: tlist   = .false.
+     logical :: vtors   = .false.
+     logical :: vbond   = .false.
+     logical :: vangl   = .false.
+     logical :: warning = .false.
+
+   contains
+
+   procedure :: any
+
+   end type TPrintTopo
 
 contains
 
@@ -124,5 +141,12 @@ subroutine zero(self)
 
 end subroutine zero
 
+function any(self) result(tf)
+  class(TPrintTopo), intent(in) :: self
+  logical :: tf
+
+  tf = self%nb.or.self%bpair.or.self%alist.or.self%blist.or. &
+     & self%tlist.or.self%vtors.or.self%vbond.or.self%vangl
+end function any
 
 end module xtb_gfnff_topology

src/main/json.f90

@@ -238,3 +238,105 @@ subroutine write_json_reduced_masses(ijson,freqres)
    write(ijson,'(3x,f15.8,",")') (freqres%rmass(i),i=1,freqres%n3true-1)
    write(ijson,'(3x,f15.8,"],")') freqres%rmass(freqres%n3true)
 end subroutine write_json_reduced_masses
+
+subroutine write_json_gfnff_lists(n, topo, printTopo)
+  use xtb_gfnff_topology, only : TGFFTopology
+  use xtb_gfnff_topology, only : TPrintTopo
+  include 'xtb_version.fh'
+  !> gfnff topology lists
+  type(TGFFTopology), intent(in) :: topo
+  !> topology printout booleans
+  type(TPrintTopo), intent(in) :: printTopo
+  character(len=:),allocatable :: cmdline
+  integer :: iunit, j, n
+
+  call open_file(iunit,'gfnff_lists.json','w')
+  ! header
+  write(iunit,'("{")')
+  ! lists printout
+  if(printTopo%nb)then ! nb(20,n)
+    write(iunit,'(3x,''"nb":'',"[")')
+    do j=1, n-1
+      write(iunit,'(3x,"[",*(i7,:,","))',advance='no') topo%nb(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(i7,:,","),"]",/)',advance='no') topo%nb(:,n)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%bpair)then ! bpair(n*(n+1)/2) packed symmetric matrix
+    write(iunit,'(3x,''"bpair":'',"[")')
+    write(iunit,'(3x,*(i7,:,","))',advance='no') topo%bpair
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%alist)then ! alist(3,nangl)
+    write(iunit,'(3x,''"alist":'',"[")')
+    do j=1, topo%nangl-1
+      write(iunit,'(3x,"[",*(i8,:,","))',advance='no') topo%alist(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(i8,:,","),"]",/)',advance='no') topo%alist(:,topo%nangl)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%blist)then ! blist(2,nbond)
+    write(iunit,'(3x,''"blist":'',"[")')
+    do j=1, topo%nbond-1
+      write(iunit,'(3x,"[",*(i8,:,","))',advance='no') topo%blist(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(i8,:,","),"]",/)',advance='no') topo%blist(:,topo%nbond)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%tlist)then ! tlist(5,ntors)
+    write(iunit,'(3x,''"tlist":'',"[")')
+    do j=1, topo%ntors-1
+      write(iunit,'(3x,"[",*(i8,:,","))',advance='no') topo%tlist(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(i8,:,","),"]",/)',advance='no') topo%tlist(:,topo%ntors)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%vtors)then ! vtors(2,ntors)
+    write(iunit,'(3x,''"vtors":'',"[")')
+    do j=1, topo%ntors-1
+      write(iunit,'(3x,"[",*(f25.15,:,","))',advance='no') topo%vtors(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(f25.15,:,","),"]",/)',advance='no') topo%vtors(:,topo%ntors)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%vbond)then ! vbond(3,nbond)
+    write(iunit,'(3x,''"vbond":'',"[")')
+    do j=1, topo%nbond-1
+      write(iunit,'(3x,"[",*(f25.15,:,","))',advance='no') topo%vbond(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(f25.15,:,","),"]",/)',advance='no') topo%vbond(:,topo%nbond)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  if(printTopo%vangl)then ! vangl(2,nangl)
+    write(iunit,'(3x,''"vangl":'',"[")')
+    do j=1, topo%nangl-1
+      write(iunit,'(3x,"[",*(f25.15,:,","))',advance='no') topo%vangl(:,j)
+      write(iunit,'("],")')
+    enddo
+    write(iunit,'(3x,"[",*(f25.15,:,","),"]",/)',advance='no') topo%vangl(:,topo%nangl)
+    write(iunit,'("]")')
+    write(iunit,'(3x,"],")')
+  endif
+  ! footer
+  call get_command(length=l)
+  allocate( character(len=l) :: cmdline )
+  call get_command(cmdline)
+  write(iunit,'(3x,''"program call":'',1x,''"'',a,''",'')') cmdline
+  write(iunit,'(3x,''"method": "GFN-FF"'',",")')
+  write(iunit,'(3x,a)') '"xtb version": "'//version//'"'
+  write(iunit,'("}")')
+  call close_file(iunit)
+
+end subroutine write_json_gfnff_lists

src/prog/main.F90

@@ -84,6 +84,7 @@ module xtb_prog_main
    use xtb_disp_dftd3param
    use xtb_disp_dftd4
    use xtb_gfnff_param, only : gff_print
+   use xtb_gfnff_topology, only : TPrintTopo
    use xtb_gfnff_convert, only : struc_convert
    use xtb_scan
    use xtb_kopt
@@ -193,6 +194,7 @@ subroutine xtbMain(env, argParser)
    integer :: TID, OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
    integer :: nproc
 
+   type(TPrintTopo) :: printTopo ! gfnff topology printout list
 
    xenv%home = env%xtbhome
    xenv%path = env%xtbpath
@@ -201,12 +203,13 @@ subroutine xtbMain(env, argParser)
    ! ------------------------------------------------------------------------
    !> read the command line arguments
    call parseArguments(env, argParser, xcontrol, fnv, acc, lgrad, &
-      & restart, gsolvstate, strict, copycontrol, coffee)
+      & restart, gsolvstate, strict, copycontrol, coffee, printTopo)
 
    nFiles = argParser%countFiles()
    select case(nFiles)
    case(0)
       if (.not.coffee) then
+         if(printTopo%warning) call env%error("Eventually the input file was given to wrtopo as an argument.",source)
          call env%error("No input file given, so there is nothing to do", source)
       else
          fname = 'coffee'
@@ -854,7 +857,12 @@ subroutine xtbMain(env, argParser)
          call close_file(ich)
       end select
    endif
-
+   if(printTopo%any()) then
+     select type(calc)
+       type is(TGFFCalculator)
+         call write_json_gfnff_lists(mol%n,calc%topo,printTopo)
+     end select
+   endif
    if ((runtyp.eq.p_run_opt).or.(runtyp.eq.p_run_ohess).or. &
       (runtyp.eq.p_run_omd).or.(runtyp.eq.p_run_screen).or. &
       (runtyp.eq.p_run_metaopt).or.(runtyp.eq.p_run_bhess)) then
@@ -1084,7 +1092,7 @@ end subroutine xtbMain
 
 !> Parse command line arguments and forward them to settings
 subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
-      & restart, gsolvstate, strict, copycontrol, coffee)
+      & restart, gsolvstate, strict, copycontrol, coffee, printTopo)
    use xtb_mctc_global, only : persistentEnv
 
    !> Name of error producer
@@ -1117,6 +1125,9 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
    !> Debugging with a lot of caffeine
    logical, intent(out) :: coffee
 
+   !> topology printout list
+   type(TPrintTopo), intent(out) :: printTopo
+
    !> Print the gradient to file
    logical, intent(out) :: lgrad
 
@@ -1125,7 +1136,7 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
 
 !$ integer :: omp_get_num_threads, nproc
    integer :: nFlags
-   integer :: idum
+   integer :: idum, ndum
    real(wp) :: ddum
    character(len=:), allocatable :: flag, sec
 
@@ -1551,11 +1562,74 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
             call env%error("No input file for RMSD bias provided", source)
          end if
 
+      case('--wrtopo')
+         call args%nextArg(sec)
+         if (allocated(sec)) then
+           call setWRtopo(sec,printTopo)
+           if(printTopo%warning) call env%error("A wrtopo argument has been misspelled.",source)
+         else
+           call env%error("The wrtopo keyword is missing an argument.",source)
+         endif
       end select
       call args%nextFlag(flag)
    end do
 
 end subroutine parseArguments
 
+! set booleans for requested topology list printout
+subroutine setWRtopo(sec,printTopo)
+   ! command line argument
+   character(len=*), intent(in) :: sec
+   ! type holds booleans of to be printed topology lists
+   type(TPrintTopo), intent(inout) :: printTopo
+   ! seperator for lists is ","
+   character, parameter :: sep = ","
+   ! current and old position of seperator
+   integer :: curr_pos, old_pos
+   integer :: lenSec, i
+
+   curr_pos = 0
+   old_pos = 0
+   lenSec = len(sec)
+   do i=1, lenSec
+     curr_pos = scan(sec(curr_pos+1:lenSec),sep)+old_pos
+     if(curr_pos.ne.old_pos) then
+       call selectList(sec(old_pos+1:curr_pos-1),printTopo)
+     else
+       call selectList(sec(old_pos+1:lenSec),printTopo)
+       exit
+     endif
+     old_pos=curr_pos
+   enddo
+
+end subroutine setWRtopo
+
+subroutine selectList(secSplit, printTopo)
+   ! part of command line argument
+   character(len=*), intent(in) :: secSplit
+   ! holds booleans of to be printed topology lists
+   type(TPrintTopo), intent(inout) :: printTopo
+
+   select case(secSplit)
+   case("nb")
+     printTopo%nb = .true.
+   case("bpair")
+     printTopo%bpair = .true.
+   case("alist")
+     printTopo%alist = .true.
+   case("blist")
+     printTopo%blist = .true.
+   case("tlist")
+     printTopo%tlist = .true.
+   case("vtors")
+     printTopo%vtors = .true.
+   case("vbond")
+     printTopo%vbond = .true.
+   case("vangl")
+     printTopo%vangl = .true.
+   case default
+     printTopo%warning = .true.
+   end select
+end subroutine selectList
 
 end module xtb_prog_main