Update metaopt algorithm

man/xcontrol.7.adoc

@@ -257,7 +257,7 @@ $metadyn (6.1 only)
 *scale factor*='int','real',...::
     scales the factor 'int' with 'real'
 *alp*='real'::
-    width of the Gaussian potential used in the rmsd criteria
+    width of the Gaussian potential used in the rmsd criteria or width of the Gaussian for random displacement for --metaopt
 *coord*='file'::
     external structures to initialize the rmsd criteria (xmol format required)
 *atoms*: 'list',...::

src/CMakeLists.txt

@@ -124,6 +124,7 @@ list(APPEND srcs
   "${dir}/pseudodiag.f90"
   "${dir}/qpot.f90"
   "${dir}/qsort.f90"
+  "${dir}/random_generators.f90"
   "${dir}/rdcoord2.f90"
   "${dir}/read_gfn_param.f90"
   "${dir}/readin.f90"

src/meson.build

@@ -122,6 +122,7 @@ srcs += files(
   'pseudodiag.f90',
   'qpot.f90',
   'qsort.f90',
+  'random_generators.f90',
   'rdcoord2.f90',
   'read_gfn_param.f90',
   'readin.f90',

src/optimizer.f90

@@ -304,7 +304,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, &
    character(len=9):: hessfmt
 
    ! print ANCopt header !
-   call ancopt_header(env%unit,set%veryverbose)
+   if (pr) call ancopt_header(env%unit,set%veryverbose)
 
    if(mol%n.eq.1) return ! skip optimization for 1 atom 
 

src/prog/main.F90

@@ -98,6 +98,7 @@ module xtb_prog_main
    use xtb_ptb_calculator, only: TPTBCalculator
    use xtb_solv_cpx, only: TCpcmx
    use xtb_dipro, only: get_jab, jab_input
+   use random_generators, only: normal_distribution
    !> PTB related modules
    use xtb_main_json, only: main_ptb_json
 
@@ -1156,25 +1157,44 @@ subroutine xtbMain(env, argParser)
             metaset%xyz(:, :, metaset%nstruc) = mol%xyz
             ! randomize structure to avoid zero RMSD
             do i = 1, mol%n
-               do j = 1, 3
-                  call random_number(er)
-                  mol%xyz(j, i) = mol%xyz(j, i) + 1.0e-6_wp * er
-               end do
+               if (.not. fixset%is_fixed(i)) then
+                  do j = 1, 3
+                     mol%xyz(j, i) = mol%xyz(j, i) + normal_distribution(metaset%global_width, 0.0_wp)
+                  end do
+               end if
             end do
+            block
+               real(wp) :: Xc(3), Yc(3) !< centroids of structures
+               real(wp) :: U(3,3) !< rotation matrix, not used
+               real(wp) :: dummy(1,1) !< gradient of rmsd; if .true. must have a proper dimensions
+               real(wp) :: rmsd_val !< computed RMSD
+               call rmsd(mol%n, mol%xyz, metaset%xyz(:, :, metaset%nstruc), &
+                 &       0, U, Xc, Yc, rmsd_val, .false., dummy)
+               write (env%unit, '("RMSD of distorted structure: ", F10.6)') rmsd_val
+            end block
             call geometry_optimization &
                &     (env, mol, chk, calc, &
                &      egap, set%etemp, set%maxscciter, set%optset%maxoptcycle, etot, g, sigma, &
                &      set%optset%optlev, set%verbose, .true., murks)
+            block
+               real(wp) :: Xc(3), Yc(3) !< centroids of structures
+               real(wp) :: U(3,3) !< rotation matrix, not used
+               real(wp) :: dummy(1,1) !< gradient of rmsd; if .true. must have a proper dimensions
+               real(wp) :: rmsd_val !< computed RMSD
+               call rmsd(mol%n, mol%xyz, metaset%xyz(:, :, metaset%nstruc), &
+                 &       0, U, Xc, Yc, rmsd_val, .false., dummy)
+               write (env%unit, '("RMSD between previously optimized and newly optimized structures: ", F10.6)') rmsd_val
+            end block
             if (.not. set%verbose) then
                write (env%unit, '("current energy:",1x,f20.8)') etot
             end if
+            call writeMolecule(mol, ich, fileType%xyz, energy=etot, gnorm=norm2(g))
             if (murks) then
-               call close_file(ich)
                write (env%unit, '(/,3x,"***",1x,a,1x,"***",/)') &
                   "FAILED TO CONVERGE GEOMETRY OPTIMIZATION"
-               call touch_file('NOT_CONVERGED')
+               mol%xyz = metaset%xyz(:, :, metaset%nstruc)
             end if
-            call writeMolecule(mol, ich, fileType%xyz, energy=etot, gnorm=norm2(g))
+            flush (env%unit)
          end do
          call close_file(ich)
          call stop_timing(6)

src/random_generators.f90

@@ -0,0 +1,78 @@
+! This file is part of xtb.
+!
+! Copyright (C) 2025 Igor S. Gerasimov
+!
+! 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 random_generators
+   use xtb_mctc_accuracy, only: sp, dp
+
+   implicit none
+
+   private
+
+   interface normal_distribution
+     module procedure normal_distribution_sp
+     module procedure normal_distribution_dp
+   end interface normal_distribution
+
+   public normal_distribution
+
+contains
+
+   !>
+   !> @brief   Generates random numbers according to the Normal (or Gaussian) random number distribution
+   !>
+   !> @details Implements Box-Muller transform for converting uniform distribution to normal
+   !>
+   !> @param[in] sigma   standard deviation
+   !> @param[in] mu      mean of distribution
+   !> @return            pseudorandom value
+   real(sp) function normal_distribution_sp(sigma, mu) result(randval)
+      real(sp), intent(in) :: sigma, mu
+      real(sp), parameter :: two_pi = 2.0_sp * 4.0_sp * atan(1.0_sp)
+      real(sp) :: u(2), mag
+
+      u = 0.0_sp
+      do while (u(1) == 0.0_sp)
+         call random_number(u)
+      end do
+      mag = sigma * sqrt(-2.0_sp * log(u(1)))
+      randval = mag * cos(two_pi * u(2)) + mu
+
+   end function normal_distribution_sp
+
+   !>
+   !> @brief   Generates random numbers according to the Normal (or Gaussian) random number distribution
+   !>
+   !> @details Implements Box-Muller transform for converting uniform distribution to normal
+   !>
+   !> @param[in] sigma   standard deviation
+   !> @param[in] mu      mean of distribution
+   !> @return            pseudorandom value
+   real(dp) function normal_distribution_dp(sigma, mu) result(randval)
+      real(dp), intent(in) :: sigma, mu
+      real(dp), parameter :: two_pi = 2.0_dp * 4.0_dp * atan(1.0_dp)
+      real(dp) :: u(2), mag
+
+      u = 0.0_dp
+      do while (u(1) == 0.0_dp)
+         call random_number(u)
+      end do
+      mag = sigma * sqrt(-2.0_dp * log(u(1)))
+      randval = mag * cos(two_pi * u(2)) + mu
+
+   end function normal_distribution_dp
+
+end module random_generators

src/type/setvar.f90

@@ -130,6 +130,7 @@ module xtb_type_setvar
    contains
       procedure :: allocate => allocate_fix
       procedure :: deallocate => deallocate_fix
+      procedure :: is_fixed => is_fixed_fix
    end type fix_setvar
 
    type :: constr_setvar
@@ -389,6 +390,18 @@ subroutine deallocate_fix(self)
    if(allocated(self%val))   deallocate( self%val )
 end subroutine deallocate_fix
 
+!>
+!> @brief check if atom is fixed
+!>
+logical function is_fixed_fix(self, idx) result(is)
+   class(fix_setvar), intent(in) :: self
+   integer, intent(in) :: idx
+   is = .false.
+   if (allocated(self%atoms)) then
+      if (any(self%atoms == idx)) is = .true.
+   end if
+end function is_fixed_fix
+
 subroutine allocate_constr(self,nat,nval,fc,expo)
    implicit none
    class(constr_setvar) :: self