From 7a7ea160e23962ff51c8f3ce95bbf7e0689b8856 Mon Sep 17 00:00:00 2001 From: Daniel Mejia-Rodriguez Date: Mon, 11 Apr 2022 17:44:46 -0400 Subject: [PATCH] Better xTB restart --- src/xtb/GNUmakefile | 3 ++- src/xtb/xtb.fh | 4 +++- src/xtb/xtb_energy.F | 18 ++++++++++-------- src/xtb/xtb_energy_gradient.F | 20 +++++++++----------- src/xtb/xtb_guess.F | 21 ++++++++++++++++++--- src/xtb/xtb_rdinput.F | 13 +++++++++++++ src/xtb/xtb_restart_read.F | 24 ++++++++++++++++++++++++ src/xtb/xtb_restart_write.F | 27 +++++++++++++++++++++++++++ 8 files changed, 106 insertions(+), 24 deletions(-) create mode 100644 src/xtb/xtb_restart_read.F create mode 100644 src/xtb/xtb_restart_write.F diff --git a/src/xtb/GNUmakefile b/src/xtb/GNUmakefile index be8510584d..868d041b65 100644 --- a/src/xtb/GNUmakefile +++ b/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 diff --git a/src/xtb/xtb.fh b/src/xtb/xtb.fh index 15635ca2c3..dc6920fda1 100644 --- a/src/xtb/xtb.fh +++ b/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 diff --git a/src/xtb/xtb_energy.F b/src/xtb/xtb_energy.F index cdf54651d1..d2b419b809 100644 --- a/src/xtb/xtb_energy.F +++ b/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) ! diff --git a/src/xtb/xtb_energy_gradient.F b/src/xtb/xtb_energy_gradient.F index 23b0b0f624..00b5e6477c 100644 --- a/src/xtb/xtb_energy_gradient.F +++ b/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) diff --git a/src/xtb/xtb_guess.F b/src/xtb/xtb_guess.F index bb1c71d2b6..0e737c6f02 100644 --- a/src/xtb/xtb_guess.F +++ b/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 diff --git a/src/xtb/xtb_rdinput.F b/src/xtb/xtb_rdinput.F index 4a9ec21570..fc61218160 100644 --- a/src/xtb/xtb_rdinput.F +++ b/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 diff --git a/src/xtb/xtb_restart_read.F b/src/xtb/xtb_restart_read.F new file mode 100644 index 0000000000..0d681e57f3 --- /dev/null +++ b/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 diff --git a/src/xtb/xtb_restart_write.F b/src/xtb/xtb_restart_write.F new file mode 100644 index 0000000000..e3f39237b7 --- /dev/null +++ b/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