Implement script driver (#579)

src/extern/CMakeLists.txt

@@ -18,6 +18,7 @@
 set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 
 list(APPEND srcs
+  "${dir}/driver.f90"
   "${dir}/mopac.f90"
   "${dir}/orca.f90"
   "${dir}/turbomole.f90"

src/extern/driver.f90

@@ -0,0 +1,199 @@
+! This file is part of xtb.
+! SPDX-Identifier: LGPL-3.0-or-later
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
+
+module xtb_extern_driver
+  use xtb_mctc_accuracy, only: wp
+  use xtb_mctc_io, only: stdout
+  use xtb_mctc_filetypes, only: fileType
+  use xtb_mctc_symbols, only: toSymbol
+  use xtb_type_calculator, only: TCalculator
+  use xtb_type_data, only: scc_results
+  use xtb_type_environment, only: TEnvironment
+  use xtb_type_molecule, only: TMolecule, len
+  use xtb_type_param, only: scc_parameter
+  use xtb_type_restart, only: TRestart
+  use xtb_io_writer, only: writeMolecule
+  use xtb_mctc_systools
+  use xtb_mctc_strings
+  use xtb_setparam
+  use xtb_readin
+  use xtb_mctc_convert
+  use xtb_fixparam
+  use xtb_scanparam
+  use xtb_sphereparam
+  use xtb_metadynamic
+  use xtb_constrainpot
+  use xtb_extern_turbomole, only: wrtm, rdtm
+  implicit none
+  private
+
+  public :: TDriverCalculator, newDriverCalculator
+
+  type, extends(TCalculator) :: TDriverCalculator
+    type(qm_external) :: ext
+  contains
+    !> Perform single point calculation
+    procedure :: singlepoint
+    !> Write informative printout
+    procedure :: writeInfo
+  end type TDriverCalculator
+
+contains
+
+!> Create a new calculator for driving the external driver program
+  subroutine newDriverCalculator(self, env, ext)
+    !> Instance of the external driver calculator
+    type(TDriverCalculator), intent(out) :: self
+    !> Calculation environment
+    type(TEnvironment), intent(inout) :: env
+    !> Settings for the external driver calculator
+    type(qm_external), intent(in) :: ext
+
+    self%threadsafe = .false.
+    self%ext = ext
+  end subroutine newDriverCalculator
+
+  subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
+        & energy, gradient, sigma, hlgap, results)
+    !> Source of the generated errors
+    character(len=*), parameter :: source = 'extern_driver_singlepoint'
+    !> Calculator instance
+    class(TDriverCalculator), intent(inout) :: self
+    !> Computational environment
+    type(TEnvironment), intent(inout) :: env
+    !> Molecular structure data
+    type(TMolecule), intent(inout) :: mol
+    !> Wavefunction data
+    type(TRestart), intent(inout) :: chk
+    !> Print level for IO
+    integer, intent(in) :: printlevel
+    !> Restart from previous results
+    logical, intent(in) :: restart
+    !> Total energy
+    real(wp), intent(out) :: energy
+    !> Molecular gradient
+    real(wp), intent(out) :: gradient(:, :)
+    !> Strain derivatives
+    real(wp), intent(out) :: sigma(:, :)
+    !> HOMO-LUMO gap
+    real(wp), intent(out) :: hlgap
+    !> Detailed results
+    type(scc_results), intent(out) :: results
+
+    integer :: i, ich
+    integer :: mode_sp_run = 1
+    real(wp) :: efix
+    real(wp) :: dipole(3)
+    logical :: cache, exist
+    logical, parameter :: ccm = .true.
+    logical :: exitRun
+    character(len=*), parameter :: outfmt = &
+                                   '(9x,"::",1x,a,f23.12,1x,a,1x,"::")'
+    real(wp) :: xyz_cached(3, mol%n)
+    integer :: err
+
+    call mol%update
+
+    energy = 0.0_wp
+    gradient(:, :) = 0.0_wp
+    sigma(:, :) = 0.0_wp
+    hlgap = 0.0_wp
+    efix = 0.0_wp
+    dipole(:) = 0.0_wp
+
+    !$omp critical (turbo_lock)
+    inquire (file='gradient', exist=exist)
+    if (exist) then
+      call rdtm(mol%n, .true., energy, gradient, xyz_cached)
+      cache = all(abs(xyz_cached - mol%xyz) < 1.e-10_wp)
+    end if
+    if (.not. cache) then
+      call wrtm(mol%n, mol%at, mol%xyz)
+
+      write (env%unit, '(72("="))')
+      write (env%unit, '(1x,"*",1x,a)') &
+        "letting driver take over the control..."
+      call execute_command_line('exec 2>&1 '//self%ext%executable, exitstat=err)
+      if (err /= 0) then
+        call env%error('driver returned with non-zero exit status, doing the same', source)
+      else
+        write (env%unit, '(1x,"*",1x,a)') &
+          "successful driver run, taking over control again..."
+      end if
+      write (env%unit, '(72("="))')
+
+      call rdtm(mol%n, .true., energy, gradient, xyz_cached)
+    end if
+    !$omp end critical (turbo_lock)
+
+    call env%check(exitRun)
+    if (exitRun) then
+      call env%error("Electronic structure method terminated", source)
+      return
+    end if
+
+    ! ------------------------------------------------------------------------
+    !  various external potentials
+    call constrain_pot(potset, mol%n, mol%at, mol%xyz, gradient, efix)
+    call constrpot(mol%n, mol%at, mol%xyz, gradient, efix)
+    call cavity_egrad(mol%n, mol%at, mol%xyz, efix, gradient)
+    call metadynamic(metaset, mol%n, mol%at, mol%xyz, efix, gradient)
+    call metadynamic(rmsdset, mol%n, mol%at, mol%xyz, efix, gradient)
+
+    ! ------------------------------------------------------------------------
+    !  fixing of certain atoms
+    !  print*,abs(efix/etot)
+    energy = energy + efix
+    results%e_total = energy
+    results%gnorm = norm2(gradient)
+    results%dipole = dipole
+    if (fixset%n .gt. 0) then
+      do i = 1, fixset%n
+        !print*,i,fixset%atoms(i)
+        gradient(1:3, fixset%atoms(i)) = 0
+      end do
+    end if
+
+    if (printlevel .ge. 2) then
+      ! start with summary header
+      if (.not. set%silent) then
+        write (env%unit, '(9x,53(":"))')
+        write (env%unit, '(9x,"::",21x,a,21x,"::")') "SUMMARY"
+      end if
+      write (env%unit, '(9x,53(":"))')
+      write (env%unit, outfmt) "total energy      ", results%e_total, "Eh   "
+      write (env%unit, outfmt) "gradient norm     ", results%gnorm, "Eh/a0"
+      write (env%unit, outfmt) "HOMO-LUMO gap     ", results%hl_gap, "eV   "
+      write (env%unit, '(9x,53(":"))')
+      write (env%unit, '(a)')
+    end if
+  end subroutine singlepoint
+
+  subroutine writeInfo(self, unit, mol)
+
+    !> Calculator instance
+    class(TDriverCalculator), intent(in) :: self
+
+    !> Unit for I/O
+    integer, intent(in) :: unit
+
+    !> Molecular structure data
+    type(TMolecule), intent(in) :: mol
+
+    call generic_header(unit, "Driver", 49, 10)
+  end subroutine writeInfo
+
+end module xtb_extern_driver

src/extern/meson.build

@@ -16,6 +16,7 @@
 # along with xtb.  If not, see <https://www.gnu.org/licenses/>.
 
 srcs += files(
+  'driver.f90',
   'mopac.f90',
   'orca.f90',
   'turbomole.f90',

src/extern/turbomole.f90

@@ -43,7 +43,7 @@ module xtb_extern_turbomole
    implicit none
    private
 
-   public :: wrtm
+   public :: wrtm, rdtm
    public :: TTMCalculator, newTMCalculator
 
    type, extends(TCalculator) :: TTMCalculator

src/main/setup.f90

@@ -23,6 +23,7 @@ module xtb_main_setup
    use xtb_extern_orca, only : TOrcaCalculator, newOrcaCalculator
    use xtb_extern_mopac, only : TMopacCalculator, newMopacCalculator
    use xtb_extern_turbomole, only : TTMCalculator, newTMCalculator
+   use xtb_extern_driver, only : TDriverCalculator, newDriverCalculator
    use xtb_type_calculator, only : TCalculator
    use xtb_type_environment, only : TEnvironment
    use xtb_type_molecule, only : TMolecule
@@ -66,6 +67,7 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input)
    type(TMopacCalculator), allocatable :: mopac
    type(TTMCalculator), allocatable :: turbo
    type(TOniomCalculator), allocatable :: oniom
+   type(TDriverCalculator), allocatable :: driver
 
    logical :: exitRun
 
@@ -108,7 +110,10 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input)
       allocate(turbo)
       call newTMCalculator(turbo, set%extcode, set%extmode)
       call move_alloc(turbo, calc)
-
+   case(p_ext_driver)
+      allocate(driver)
+      call newDriverCalculator(driver, env, set%ext_driver)
+      call move_alloc(driver, calc)
    case(p_ext_oniom)
       allocate(oniom)
       call newOniomCalculator(oniom, env, mol, input)

src/prog/main.F90

@@ -1347,6 +1347,13 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
       case('--orca')
          call set_exttyp('orca')
 
+      case('--driver')
+         call set_exttyp('driver')
+         call args%nextArg(sec)
+         if (allocated(sec)) then
+            set%ext_driver%executable = sec
+         end if
+
       case('--mopac')
          call set_exttyp('mopac')
 

src/set_module.f90

@@ -912,6 +912,8 @@ subroutine set_exttyp(typ)
       set%mode_extrun = p_ext_xtb
    case('orca')
       set%mode_extrun = p_ext_orca
+   case('driver')
+      set%mode_extrun = p_ext_driver
    case('oniom')
       set%mode_extrun = p_ext_oniom
    case('turbomole')

src/setparam.f90

@@ -112,6 +112,7 @@ module xtb_setparam
    integer, parameter :: p_ext_vtb       = -1
    integer, parameter :: p_ext_eht       =  0
    integer, parameter :: p_ext_xtb       =  1
+   integer, parameter :: p_ext_driver    =  3
    integer, parameter :: p_ext_turbomole =  4
    integer, parameter :: p_ext_orca      =  5
    integer, parameter :: p_ext_mopac     = 12
@@ -386,6 +387,7 @@ module xtb_setparam
 
 !! ------------------------------------------------------------------------
 
+   type(qm_external) :: ext_driver
    type(qm_external) :: ext_orca
    type(qm_external) :: ext_turbo
    type(qm_external) :: ext_mopac