Better xTB restart

src/xtb/GNUmakefile

@@ -3,7 +3,8 @@ $(error USE_TBLITE must be set in order to use xTB methods)
 endif
 
 OBJ = xtb_energy.o xtb_input.o xtb_energy_gradient.o \
-      xtb_rdinput.o xtb_guess.o xtb_spinpol.o xtb_print_param.o
+      xtb_rdinput.o xtb_guess.o xtb_spinpol.o xtb_print_param.o \
+      xtb_restart_read.o xtb_restart_write.o
 
 LIBRARY = libnwxtb.a
 

src/xtb/xtb.fh

@@ -1,7 +1,9 @@
       double precision acc, rcharge, kt, broyden
       logical gfn1, gfn2, gff, guess
       integer uhf, verb, nspin
+      character*255 restart_out
       common /cxtb/
      $    acc, rcharge, kt, broyden,
      $    gfn1, gfn2, gff, guess,
-     $    uhf, verb, nspin
+     $    uhf, verb, nspin,
+     $    restart_out

src/xtb/xtb_energy.F

@@ -31,7 +31,7 @@
 
       integer rtdb
       integer geom, i, iat
-      logical converged, status, oprint, ovecs, oparam
+      logical converged, status, oprint, ovecs, oparam, restart
       double precision walltime, cputime
 
       integer istat
@@ -216,17 +216,17 @@
         !
         ! Compute starting charges 
         !
-        if (guess) then
-          call xtb_guess(mol,wfn%qat,wfn%qsh,param,mol%nat,
+        status = rtdb_get(rtdb,'xtb:restart',mt_log,1,restart)
+        if (guess.and..not.status) then
+          call xtb_guess(mol,wfn%qat,wfn%qsh,wfn%dpat,wfn%qpat,
+     $                   param,mol%nat,
      $                   calc%bas%nsh,calc%bas%nsh_at,atnum4)
-          !
-          ! Copy charges to second spin-channel, if available
-          !
           if (nspin.eq.2) then
             wfn%qat(:,2) = wfn%qat(:,1)
             wfn%qsh(:,2) = wfn%qsh(:,1)
           endif
-
+        elseif(status) then
+          call xtb_restart_read(mol, calc, wfn)
         endif
 
 
@@ -275,12 +275,14 @@
             write(luout,*)
           endif
         endif
+        call xtb_restart_write(mol, calc, wfn)
+        status = rtdb_put(rtdb,'xtb:restart',mt_log,1,.true.)
+        xtb_energy = .not.ctx%failed()
       endif
 
       !
       ! Check if something failed during the calculation
       !
-      xtb_energy = .not.ctx%failed()
       call ga_brdcst(101,xtb_energy,ma_sizeof(mt_log,1,mt_byte),0)
 
       ! 

src/xtb/xtb_energy_gradient.F

@@ -41,7 +41,7 @@
 
       integer iuhf8
 
-      logical oskel
+      logical oskel, restart
       integer nactive
 
       character(len=16) aname,elem
@@ -239,12 +239,13 @@ c
         !
         ! Compute starting charges 
         !
-        status = rtdb_get(rtdb,'xtb:qat',mt_dbl,nat,wfn%qat)
-        status = status .and. rtdb_get(rtdb,'xtb:qsh',mt_dbl,
-     $                       int(calc%bas%nsh),wfn%qsh)
+        status = rtdb_get(rtdb,'xtb:restart',mt_log,1,restart)
         if (guess.and..not.status) then
-          call xtb_guess(mol,wfn%qat,wfn%qsh,param,mol%nat,
+          call xtb_guess(mol,wfn%qat,wfn%qsh,wfn%dpat,wfn%qpat,
+     $                   param,mol%nat,
      $                   calc%bas%nsh,calc%bas%nsh_at,atnum4)
+        elseif(status) then
+          call xtb_restart_read(mol, calc, wfn)
         endif
 
         !
@@ -336,12 +337,9 @@ c
           call util_flush(luout)
         endif
 
-        status = rtdb_put(rtdb,'xtb:qat',mt_dbl,nat,wfn%qat)
-        status = status .and. rtdb_put(rtdb,'xtb:qsh',mt_dbl,
-     $                       int(calc%bas%nsh),wfn%qsh)
-        if (.not.status)
-     $    call errquit('xtb_energy_gradients: could not put charges',0,
-     $                  RTDB_ERR)   
+        !call xtb_restart_write(wfn%density, int(calc%bas%nao))
+        call xtb_restart_write(mol, calc, wfn)
+        status = rtdb_put(rtdb,'xtb:restart',mt_log,1,.true.)
 
         dipole(1:3) = sum(wfn%dpat(:,:,1),dim=2)
         status = rtdb_put(rtdb,'task:dipole',mt_dbl,3,dipole)

src/xtb/xtb_guess.F

@@ -1,4 +1,5 @@
-      subroutine xtb_guess(mol,qat,qshl,param,nat,nshl,nshl_at,atnum)
+      subroutine xtb_guess(mol,qat,qshl,dpat,qpat,param,nat,nshl,
+     $                     nshl_at,atnum)
       use mctc_io, only: structure_type
       use tblite_param, only : param_record
       use multicharge, only : mchrg_model_type, new_eeq2019_model,
@@ -10,11 +11,11 @@
       integer*4 :: nat, nshl
       integer*4 :: nshl_at(nat)
       integer*4 :: atnum(nat)
-      double precision :: qat(nat), qshl(nshl)
+      double precision :: qat(nat), qshl(nshl), dpat(3,nat), qpat(6,nat)
 
       type(mchrg_model_type) :: model
       integer :: i, j, ksh
-      double precision :: tot
+      double precision :: tot, vec(3), tr
       double precision, allocatable :: rcov(:), cn(:), trans(:,:)
 
 
@@ -39,6 +40,20 @@
         enddo
       enddo
 
+
+      do i=1,nat
+        dpat(:,i) = qat(i)*mol%xyz(:,i)
+        qpat(1,i) = qat(i)*mol%xyz(1,i)**2
+        qpat(2,i) = qat(i)*mol%xyz(1,i)*mol%xyz(2,i)
+        qpat(3,i) = qat(i)*mol%xyz(2,i)**2
+        qpat(4,i) = qat(i)*mol%xyz(3,i)*mol%xyz(1,i)
+        qpat(5,i) = qat(i)*mol%xyz(3,i)*mol%xyz(2,i)
+        qpat(6,i) = qat(i)*mol%xyz(3,i)**2
+      enddo
+
+      dpat = 0d0
+      qpat = 0d0
+
       deallocate(rcov,cn,trans)
 
       end subroutine

src/xtb/xtb_rdinput.F

@@ -21,6 +21,17 @@
       if (.not.rtdb_get(rtdb,'xtb:broyden',mt_dbl,1,broyden)) 
      $   broyden=0.4d0
 
+      ! Figure restart file name
+      if(.not.rtdb_cget(rtdb,'xtb:restart_out',1,restart_out))
+     $  restart_out = ' '
+
+      if (restart_out.eq.' ') then
+        call util_file_name('xtbrestart',.false.,.false.,restart_out)
+      endif  
+
+      !Resolve restart file names to full paths
+      call util_file_name_resolve(restart_out,.false.)
+
       if (method.eq.'gfn2') then
         gfn2 = .true.
         gfn1 = .false.
@@ -51,5 +62,7 @@
      $  call errquit('xtb_rdinput: cannot put in rtdb',7,RTDB_ERR) 
       if (.not.rtdb_put(rtdb,'xtb:broyden',mt_dbl,1,broyden)) 
      $  call errquit('xtb_rdinput: cannot put in rtdb',8,RTDB_ERR) 
+      if(.not.rtdb_cput(rtdb,'xtb:restart_out',1,restart_out))
+     $  call errquit('xtb_rdinput: cannot put in rtdb',9,RTDB_ERR) 
 
       end subroutine

src/xtb/xtb_restart_read.F

@@ -0,0 +1,24 @@
+      subroutine xtb_restart_read(mol, calc, wfn)
+      use mctc_io, only : structure_type
+      use tblite_wavefunction_type, only : wavefunction_type
+      use tblite_xtb_calculator, only: xtb_calculator
+      implicit none
+#include "xtb.fh"
+      type(structure_type) :: mol
+      type(wavefunction_type) :: wfn
+      type(xtb_calculator) :: calc
+      integer :: ispin,unitno
+
+      parameter(unitno=66)
+
+      open(unit=unitno,status='old',form='unformatted',file=restart_out)
+      do ispin=1,nspin
+        call sread(unitno, wfn%qsh(1,ispin), int(calc%bas%nsh))
+        call sread(unitno, wfn%qat(1,ispin), int(mol%nat))
+        call sread(unitno, wfn%dpat(1,1,ispin), int(3*mol%nat))
+        call sread(unitno, wfn%qpat(1,1,ispin), int(6*mol%nat))
+      enddo
+      close(unit=unitno)
+      return
+
+      end

src/xtb/xtb_restart_write.F

@@ -0,0 +1,27 @@
+      subroutine xtb_restart_write(mol, calc, wfn)
+      use mctc_io, only : structure_type
+      use tblite_xtb_calculator, only: xtb_calculator
+      use tblite_wavefunction_type, only : wavefunction_type
+      implicit none
+#include "xtb.fh"
+
+      type(structure_type) :: mol
+      type(xtb_calculator) :: calc
+      type(wavefunction_type) :: wfn
+
+      integer i, unitno, ispin
+      parameter(unitno=66)
+
+      open(unit=unitno,status='unknown',form='unformatted',
+     $     file=restart_out)
+      do ispin=1,nspin
+c        do i=1,nao
+c          call swrite(unitno, dmat(1,i,ispin), int(calc%bas%nao))
+c        enddo
+        call swrite(unitno, wfn%qsh(1,ispin), int(calc%bas%nsh))
+        call swrite(unitno, wfn%qat(1,ispin), int(mol%nat))
+        call swrite(unitno, wfn%dpat(1,1,ispin), int(3*mol%nat))
+        call swrite(unitno, wfn%qpat(1,1,ispin), int(6*mol%nat))
+      enddo
+      close(unit=unitno)
+      end