diff --git a/man/xcontrol.7.adoc b/man/xcontrol.7.adoc index 96a62d8fe..8d74a9e5b 100644 --- a/man/xcontrol.7.adoc +++ b/man/xcontrol.7.adoc @@ -353,15 +353,13 @@ $opt *maxdispl*='real':: maximum coordinate displacement in `ancopt(3)` +*average conv*='bool':: + average the energy and gradient before checking for convergence to accelerate + numerically noisy potential energy surfaces (default: false). + *s6*='real':: dispersion scaling in ANC generation -*ts*='bool':: - dummy - -*tsroot*='int':: - dummy - *hessian*=lindh-d2|lindh|swart:: model hessian for generation of ANC used in optimization diff --git a/src/optimizer.f90 b/src/optimizer.f90 index b58319fbb..2ec0d93f0 100644 --- a/src/optimizer.f90 +++ b/src/optimizer.f90 @@ -22,9 +22,23 @@ module xtb_optimizer use xtb_type_environment, only : TEnvironment use xtb_bfgs use xtb_david2 + implicit none logical,private,parameter :: profile = .true. + type :: convergence_log + integer :: nlog + real(wp), allocatable :: elog(:) + real(wp), allocatable :: glog(:) + contains + procedure :: set_eg_log + procedure :: get_averaged_energy + procedure :: get_averaged_gradient + end type convergence_log + interface convergence_log + module procedure new_convergence_log + end interface convergence_log + contains subroutine get_optthr(n,olev,ethr,gthr,maxcycle,acc) @@ -100,6 +114,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, & use xtb_type_anc use xtb_type_restart use xtb_type_calculator + use xtb_type_dummycalc use xtb_type_data use xtb_type_timer @@ -158,6 +173,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, & integer, allocatable :: totsym(:) real(wp),allocatable :: pmode(:,:) real(wp),allocatable :: grmsd(:,:) + type(convergence_log), allocatable :: avconv real(wp) :: U(3,3), x_center(3), y_center(3), rmsdval integer :: modef logical :: restart,ex,converged,linear @@ -202,6 +218,18 @@ subroutine ancopt(env,ilog,mol,chk,calc, & maxopt=maxcycle_in endif if(maxopt.lt.maxmicro) maxmicro=maxopt + if (optset%average_conv) then + select type(calc) + class is(TDummyCalculator) + avconv = load_turbomole_log(maxopt) + if (avconv%nlog > 0 .and. pr) then + write(env%unit, '(a, 1x, i0, 1x, a)') & + "Convergence averaging initialized with", avconv%nlog, "entries" + end if + class default + avconv = convergence_log(maxopt) + end select + end if call axis2(mol%n,mol%at,mol%xyz,aaa,bbb,ccc,dumi,dumj) @@ -351,7 +379,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, & ! now everything is prepared for the optimization call relax(env,iter,molopt,anc,restart,maxmicro,maxdispl,ethr,gthr, & & iii,chk,calc,egap,acc,et,maxiter,iupdat,etot,g,sigma,ilog,pr,fail, & - & converged,timer,optset%exact_rf) + & converged,timer,optset%exact_rf,avconv) call env%check(fail) if (fail) then @@ -426,7 +454,7 @@ end subroutine ancopt subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & & ii,chk,calc, & & egap,acc_in,et,maxiter,iupdat,etot,g,sigma,ilog,pr,fail,converged, & - & timer,exact) + & timer,exact,avconv) use xtb_mctc_blas, only : blas_gemv use xtb_type_molecule @@ -462,6 +490,7 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & real(wp),intent(inout) :: egap logical, intent(out) :: converged logical, intent(in) :: exact + type(convergence_log), intent(inout), optional :: avconv type(scc_results) :: res integer :: nvar1,npvar,npvar1 @@ -554,6 +583,18 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & return endif + if (present(avconv)) then + call avconv%set_eg_log(energy, gnorm) + energy = avconv%get_averaged_energy() + gnorm = avconv%get_averaged_gradient() + if (pr) then + write(env%unit,'("av. E:",1x,f14.7,1x,"->",1x,f14.7)') & + avconv%elog(avconv%nlog), energy + write(env%unit,'("av. G:",1x,f14.7,1x,"->",1x,f14.7)') & + avconv%glog(avconv%nlog), gnorm + end if + end if + ! adapt SCC acuracy if(gnorm .lt.0.004)then acc=acc_in @@ -964,7 +1005,7 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr) !> Calculation environment type(TEnvironment), intent(inout) :: env - + !> Calculator class(TCalculator), intent(inout) :: calc @@ -1091,5 +1132,175 @@ subroutine read_hessian(ihess, ndim, hessian, error) endif enddo end subroutine read_hessian -end module xtb_optimizer +subroutine geoconvav(nc, e, g, val, deriv) + implicit none + !> total number of E/G points + integer nc + !> total energy in Eh + real*8 e(*) + !> norm of Cartesian gradient (in TM: |dE/dxyz|) + real*8 g(*) + !> av. energy in Eh to be used further + real*8 val + !> av. gradient + real*8 deriv + + integer nav, low + integer i, j + parameter (nav=5) ! average over last nav + real*8 eav,gav + + ! only apply it if sufficient number of points i.e. a "tail" can exist + ! with the censo blockl = 8 default, this can first be effective in the second + if(nc.lt.3*nav)then + val = e(nc) + deriv = g(nc) + return + endif + + low = max(1,nc-nav+1) + j = 0 + eav = 0 + do i=nc,low,-1 + j = j + 1 + eav = eav + e(i) + gav = gav + g(i) + enddo + val = eav / float(j) + + low = max(1,nc-nav+1) + j = 0 + gav = 0 + do i=nc,low,-1 + j = j + 1 + gav = gav + g(i) + enddo + ! adjust the gradient norm to xtb "conventions" because e.g. a noisy + ! DCOSMO-RS gradient for large cases can never (even on average) + ! become lower than the "-opt normal" thresholds + deriv=gav / float(j) / 2.d0 +end subroutine geoconvav + +pure function new_convergence_log(nmax) result(self) + integer, intent(in) :: nmax + type(convergence_log) :: self + self%nlog = 0 + allocate(self%elog(nmax)) + allocate(self%glog(nmax)) +end function new_convergence_log + +function load_turbomole_log(nmax) result(self) + use xtb_mctc_resize, only : resize + use xtb_mctc_systools, only : getline + integer, intent(in) :: nmax + type(convergence_log) :: self + integer :: nlog, ilog + real(wp) :: etmp, gtmp + real(wp), allocatable :: elog(:), glog(:) + character(len=:), allocatable :: line + logical :: exist + integer, parameter :: initial_size = 50 + integer :: unit, stat + nlog = 0 + inquire(file="gradient", exist=exist) + if (exist) then + open(file="gradient", newunit=unit) + allocate(elog(initial_size), glog(initial_size)) + call getline(unit, line, stat) + do while(stat == 0) + if (index(line, "cycle") > 0) then + if (tokenize_cycle(line, etmp, gtmp)) then + if (nlog >= size(elog)) then + call resize(elog, size(elog)+size(elog)/2+1) + call resize(glog, size(glog)+size(glog)/2+1) + end if + nlog = nlog + 1 + elog(nlog) = etmp + glog(nlog) = gtmp + end if + end if + call getline(unit, line, stat) + end do + close(unit) + end if + self = convergence_log(nmax + nlog) + do ilog = 1, nlog + call self%set_eg_log(elog(ilog), glog(ilog)) + end do +end function load_turbomole_log + +function tokenize_cycle(line, e, g) result(stat) + character(len=*), intent(in) :: line + real(wp), intent(out) :: e + real(wp), intent(out) :: g + logical :: stat + integer :: stat1, stat2 + read(line(33:51), *, iostat=stat1) e + read(line(66:), *, iostat=stat2) g + stat = stat1 == 0 .and. stat2 == 0 +end function tokenize_cycle + +pure function get_averaged_energy(self) result(val) + class(convergence_log), intent(in) :: self + real(wp) :: eav, val + integer :: i, j, low + integer, parameter :: nav = 5 + + ! only apply it if sufficient number of points i.e. a "tail" can exist + ! with the censo blockl = 8 default, this can first be effective in the second + if(self%nlog.lt.3*nav)then + val = self%elog(self%nlog) + else + low = max(1, self%nlog-nav+1) + j = 0 + eav = 0 + do i = self%nlog, low, -1 + j = j + 1 + eav = eav + self%elog(i) + enddo + val = eav / float(j) + end if + +end function get_averaged_energy + +pure function get_averaged_gradient(self) result(deriv) + class(convergence_log), intent(in) :: self + real(wp) :: gav, deriv + integer :: i, j, low + integer, parameter :: nav = 5 + + ! only apply it if sufficient number of points i.e. a "tail" can exist + ! with the censo blockl = 8 default, this can first be effective in the second + if(self%nlog.lt.3*nav) then + deriv = self%glog(self%nlog) + else + low = max(1, self%nlog-nav+1) + j = 0 + gav = 0 + do i = self%nlog, low, -1 + j = j + 1 + gav = gav + self%glog(i) + enddo + ! adjust the gradient norm to xtb "conventions" because e.g. a noisy + ! DCOSMO-RS gradient for large cases can never (even on average) + ! become lower than the "-opt normal" thresholds + deriv = gav / float(j) / 2.d0 + end if + +end function get_averaged_gradient + +pure subroutine set_eg_log(self, e, g) + use xtb_mctc_resize, only : resize + class(convergence_log), intent(inout) :: self + real(wp), intent(in) :: e, g + if (self%nlog >= size(self%elog)) then + call resize(self%elog, size(self%elog)+size(self%elog)/2+1) + call resize(self%glog, size(self%glog)+size(self%glog)/2+1) + end if + self%nlog = self%nlog + 1 + self%elog(self%nlog) = e + self%glog(self%nlog) = g +end subroutine set_eg_log + +end module xtb_optimizer diff --git a/src/relaxation_engine.f90 b/src/relaxation_engine.f90 index b993efc53..cb8f9bbb1 100644 --- a/src/relaxation_engine.f90 +++ b/src/relaxation_engine.f90 @@ -23,6 +23,7 @@ module xtb_relaxation_engine use xtb_mctc_convert, only : fstoau, amutoau use xtb_mctc_fileTypes, only : fileType use xtb_type_environment, only : TEnvironment + use xtb_optimizer, only : convergence_log, load_turbomole_log use xtb_bfgs use xtb_david2 implicit none @@ -389,6 +390,7 @@ subroutine l_ancopt & use xtb_type_calculator use xtb_xtb_calculator use xtb_gfnff_calculator + use xtb_type_dummycalc use xtb_type_data use xtb_type_timer @@ -469,6 +471,7 @@ subroutine l_ancopt & real(wp), allocatable :: hdiag(:) real(wp), allocatable :: anc(:) real(wp), allocatable :: xyz0(:,:) + type(convergence_log), allocatable :: avconv ! settings for the LBFGS optimizer type(lbfgs_options) :: opt @@ -498,6 +501,20 @@ subroutine l_ancopt & ! check for user input regarding the maximum number of cycles if (maxcycle_in > 0) maxcycle = maxcycle_in + ! Activate averaged convergence criterium + if (optset%average_conv) then + select type(calc) + class is(TDummyCalculator) + avconv = load_turbomole_log(maxcycle) + if (avconv%nlog > 0 .and. pr) then + write(env%unit, '(a, 1x, i0, 1x, a)') & + "Convergence averaging initialized with", avconv%nlog, "entries" + end if + class default + avconv = convergence_log(maxcycle) + end select + end if + ! open the logfile, the log is bound to unit 942, so we cannot use newunit ! and have to hope that nobody else is currently occupying this identifier opt%ilog = ilog @@ -685,7 +702,7 @@ subroutine l_ancopt & call lbfgs_relax & & (env,iter,thiscycle,opt,molopt, & & chk,calc,energy,egap,gradient,sigma,nvar,hdiag,trafo,anc,xyz0, & - & converged,fail,timer) + & converged,fail,timer,avconv) thiscycle = min(ceiling(thiscycle*opt%cycle_inc),2*opt%micro_cycle) @@ -819,7 +836,7 @@ end subroutine lbfgs_step subroutine lbfgs_relax & & (env,iter,maxcycle,opt,mol, & & chk,calc,energy,egap,g_xyz,sigma,nvar,hdiag,trafo,anc,xyz0, & - & converged,fail,timer) + & converged,fail,timer,avconv) use xtb_type_molecule use xtb_type_restart @@ -876,6 +893,8 @@ subroutine lbfgs_relax & real(wp), intent(in) :: xyz0(:,:) !> timer for profiling the relaxation procedure type(tb_timer), intent(inout) :: timer + !> Average convergence memory + type(convergence_log), intent(inout), optional :: avconv type(scc_results) :: res @@ -1042,6 +1061,18 @@ subroutine lbfgs_relax & !call wrlog(mol%n,xyz,attyp,energy,gnorm,.false.) if (profile) call timer%measure(7) + if (present(avconv)) then + call avconv%set_eg_log(energy, gnorm) + energy = avconv%get_averaged_energy() + gnorm = avconv%get_averaged_gradient() + if (pr) then + write(env%unit,'("av. E:",1x,f14.7,1x,"->",1x,f14.7)') & + avconv%elog(avconv%nlog), energy + write(env%unit,'("av. G:",1x,f14.7,1x,"->",1x,f14.7)') & + avconv%glog(avconv%nlog), gnorm + end if + end if + ! check for convergence of the energy change ediff = energy - elast econverged = ediff < 0.0_wp .and. abs(ediff) < opt%e_thr diff --git a/src/set_module.f90 b/src/set_module.f90 index 0dafb8ac5..226e80244 100644 --- a/src/set_module.f90 +++ b/src/set_module.f90 @@ -225,6 +225,7 @@ subroutine write_set_opt(ictrl) write(ictrl,'(3x,"ts=",i0)') bool2int(tsopt) write(ictrl,'(3x,"tsroot=",i0)') tsroot write(ictrl,'(3x,"exact rf=",g0)') bool2string(optset%exact_rf) + write(ictrl,'(3x,"average conv=",g0)') bool2string(optset%average_conv) end subroutine write_set_opt subroutine write_set_thermo(ictrl) @@ -1525,6 +1526,9 @@ subroutine set_opt(env,key,val) if (.not.allocated(opt_outfile)) opt_outfile = val case('logfile') if (.not.allocated(opt_logfile)) opt_logfile = val + case('average conv') + if (getValue(env,val,ldum).and.set18) optset%average_conv = ldum + set18 = .false. end select end subroutine set_opt diff --git a/src/setparam.f90 b/src/setparam.f90 index a0cbe404b..495224ca4 100644 --- a/src/setparam.f90 +++ b/src/setparam.f90 @@ -112,7 +112,8 @@ module xtb_setparam ! maximum coordinate displacement in ancopt maxdispl_opt = 1.000_wp, & ! lowest force constant in ANC generation (should be > 0.005) - hlow_opt = 0.010_wp ) + hlow_opt = 0.010_wp, & + average_conv = .false.) integer, parameter :: p_modh_read = -2 integer, parameter :: p_modh_unit = -1 integer, parameter :: p_modh_lindh = 1 diff --git a/src/type/setvar.f90 b/src/type/setvar.f90 index 29c4c2b2e..4b8501847 100644 --- a/src/type/setvar.f90 +++ b/src/type/setvar.f90 @@ -87,6 +87,7 @@ module xtb_type_setvar ! lowest force constant in ANC generation (should be > 0.005) real(wp) :: hlow_opt = 0.0_wp logical :: exact_rf = .false. + logical :: average_conv = .false. end type ancopt_setvar type modhess_setvar