Fixing implicit solvation in tblite and safer MKL handling (#331)

Implicit solvation in tblite should now be consistent with xtb

.gitmodules

@@ -19,5 +19,5 @@
 	url = https://github.com/fortran-lang/test-drive.git
 [submodule "subprojects/tblite"]
 	path = subprojects/tblite
-	url = https://github.com/tblite/tblite.git
-    branch = main
+	url = https://github.com/pprcht/tblite.git
+    branch = xtb_solvation

CMakeLists.txt

@@ -54,6 +54,16 @@ if(NOT TARGET "OpenMP::OpenMP_Fortran" AND WITH_OpenMP)
   find_package("OpenMP" REQUIRED)
 endif()
 
+# Check if we are using OpenBLAS (need a precompiler definition if yes)
+if(LAPACK_LIBRARIES)
+  string(FIND "${LAPACK_LIBRARIES}" "openblas" _openblas_in_lapack)
+
+  if(NOT _openblas_in_lapack EQUAL -1)
+    message(STATUS "libopenblas was found as part of LAPACK")
+    add_compile_definitions(WITH_OPENBLAS)
+  endif()
+endif()
+
 # Fortran unit test interface (also used by other subprojects)
 if(NOT TARGET "test-drive::test-drive" AND WITH_TESTS)
    find_package("test-drive" REQUIRED)

README.md

@@ -57,10 +57,8 @@ Working and tested builds of CREST (mostly on Ubuntu 20.04 LTS):
 
 | Build System | Compiler | Linear Algebra Backend | Build type     | Status     | Note |
 |--------------|----------|------------------------|:--------------:|:----------:|:----:|
-| CMake 3.28.3 | GNU (gcc 10.3.0)  | [OpenBLAS](https://github.com/xianyi/OpenBLAS) (with OpenMP) | dynamic | ✅ ||
-| CMake 3.28.3 | GNU (gcc 10.3.0)  |  [MKL shared (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ✅ ||
-| CMake 3.28.3 | GNU (gcc 9.5.0)  |  [MKL shared (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ✅ ||
-| CMake 3.28.3 | [Intel (`ifort`/`icc` 2021.9.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html)   | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ⚠️  | OpenMP problem ([#285](https://github.com/crest-lab/crest/issues/285)) |
+| CMake 3.30.2 | GNU (gcc 14.1.0)  | [libopenblas 0.3.27](https://anaconda.org/conda-forge/libopenblas) | dynamic | ✅ ||
+| CMake 3.28.3 | [Intel (`ifort`/`icc` 2021.9.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html)   | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ⚠️  | OpenMP/MKL problem ([#285](https://github.com/crest-lab/crest/issues/285)) |
 | Meson 1.2.0 | [Intel (`ifort`/`icc` 2021.9.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html)   | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | static  | ✅ | [continuous release build](https://github.com/crest-lab/crest/releases/tag/latest) |
 | Meson 1.2.0 | [Intel (`ifort` 2021.9.0/`icx` 2023.1.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html)   | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | static  | ✅ ||
 

src/CMakeLists.txt

@@ -37,6 +37,7 @@ list(APPEND srcs
   "${dir}/axis_module.f90"
   "${dir}/biasmerge.f90"
   "${dir}/bondconstraint.f90"
+  "${dir}/canonical.f90"
   "${dir}/ccegen.f90"
   "${dir}/choose_settings.f90"
   "${dir}/classes.f90"
@@ -79,7 +80,6 @@ list(APPEND srcs
   "${dir}/strucreader.f90"
   "${dir}/symmetry2.f90"
   "${dir}/symmetry_i.c"
-  "${dir}/testmol.f90"
   "${dir}/timer.f90"
   "${dir}/trackorigin.f90"
   "${dir}/utilmod.f90"

src/algos/parallel.f90

@@ -99,21 +99,27 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread)
   integer :: i,j,k,l,io,ich,ich2,c,z,job_id,zcopy
   logical :: pr,wr,ex
   type(calcdata),allocatable :: calculations(:)
-  integer :: T  !> number of parallel running instances
   real(wp) :: energy,gnorm
   real(wp),allocatable :: grad(:,:),grads(:,:,:)
   integer :: thread_id,vz,job
   character(len=80) :: atmp
   real(wp) :: percent,runtime
 
   type(timer) :: profiler
+  integer :: T,Tn  !> threads and threads per core
+  logical :: nested
 
 !>--- check if we have any calculation settings allocated
   if (env%calc%ncalculations < 1) then
     write (stdout,*) 'no calculations allocated'
     return
   end if
 
+!>--- prepare calculation objects for parallelization (one per thread)
+  call new_ompautoset(env,'auto_nested',nall,T,Tn)
+  nested = env%omp_allow_nested
+
+
 !>--- prepare objects for parallelization
   T = env%threads
   allocate (calculations(T),source=env%calc)
@@ -153,7 +159,7 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread)
 !>--- loop over ensemble
   !$omp parallel &
   !$omp shared(env,calculations,nat,nall,at,xyz,eread,grads,c,k,z,pr,wr) &
-  !$omp shared(ich,ich2,mols)
+  !$omp shared(ich,ich2,mols, nested,Tn)
   !$omp single
   do i = 1,nall
 
@@ -162,6 +168,9 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread)
     !$omp task firstprivate( vz ) private(i,j,job,energy,io,thread_id,zcopy)
     call initsignal()
 
+    !>--- OpenMP nested region threads
+    if (nested) call ompmklset(Tn)
+
     thread_id = OMP_GET_THREAD_NUM()
     job = thread_id+1
     !>--- modify calculation spaces
@@ -259,14 +268,15 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
   integer :: i,j,k,l,io,ich,ich2,c,z,job_id,zcopy
   logical :: pr,wr,ex
   type(calcdata),allocatable :: calculations(:)
-  integer :: T  !> number of parallel running instances
   real(wp) :: energy,gnorm
   real(wp),allocatable :: grads(:,:,:)
   integer :: thread_id,vz,job
   character(len=80) :: atmp
   real(wp) :: percent,runtime
   type(calcdata),pointer :: mycalc
   type(timer) :: profiler
+  integer :: T,Tn  !> threads and threads per core
+  logical :: nested
 
 !>--- decide wether to skip this call
   if (trackrestart(env)) then
@@ -287,8 +297,11 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
     return
   end if
 
+!>--- prepare calculation objects for parallelization (one per thread)
+  call new_ompautoset(env,'auto_nested',nall,T,Tn)
+  nested = env%omp_allow_nested
+
 !>--- prepare objects for parallelization
-  T = env%threads
   allocate (calculations(T),source=mycalc)
   allocate (mols(T),molsnew(T))
   do i = 1,T
@@ -299,6 +312,9 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
       if (.not.ex) then
         io = makedir(trim(mycalc%calcs(j)%calcspace))
       end if
+      if(calculations(i)%calcs(j)%id == jobtype%tblite)then
+         calculations(i)%optnewinit=.true.
+      endif
       write (atmp,'(a,"_",i0)') sep,i
       calculations(i)%calcs(j)%calcspace = mycalc%calcs(j)%calcspace//trim(atmp)
       call calculations(i)%calcs(j)%printid(i,j)
@@ -331,7 +347,7 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
 !>--- loop over ensemble
   !$omp parallel &
   !$omp shared(env,calculations,nat,nall,at,xyz,eread,grads,c,k,z,pr,wr,dump) &
-  !$omp shared(ich,ich2,mols,molsnew)
+  !$omp shared(ich,ich2,mols,molsnew, nested,Tn)
   !$omp single
   do i = 1,nall
 
@@ -340,6 +356,9 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
     !$omp task firstprivate( vz ) private(j,job,energy,io,atmp,gnorm,thread_id,zcopy)
     call initsignal()
 
+    !>--- OpenMP nested region threads
+    if (nested) call ompmklset(Tn)
+
     thread_id = OMP_GET_THREAD_NUM()
     job = thread_id+1
     !>--- modify calculation spaces
@@ -1036,7 +1055,7 @@ subroutine parallel_md_finish_printout(MD,vz,io,profiler)
   if (io == 0) then
     write (btmp,'(a,1x,i3,a)') trim(atmp),vz,' completed successfully'
   else
-    write (btmp,'(a,1x,i3,a)') trim(atmp),vz,' terminated with early'
+    write (btmp,'(a,1x,i3,a)') trim(atmp),vz,' terminated EARLY'
   end if
   call profiler%write_timing(stdout,vz,trim(btmp))
 

src/algos/playground.f90

@@ -31,15 +31,7 @@ subroutine crest_playground(env,tim)
   use crest_data
   use crest_calculator
   use strucrd 
-  use optimize_module
-  use tblite_api
-  use wiberg_mayer, only: write_wbo 
-  use adjacency
-  use zdata
-  use probabilities_module
-  use parse_csv
-  use cregen_interface
-  use dynamics_module
+  use canonical_mod
   implicit none
   type(systemdata),intent(inout) :: env
   type(timer),intent(inout)      :: tim
@@ -55,21 +47,11 @@ subroutine crest_playground(env,tim)
   logical,allocatable :: rings(:,:)
   integer,allocatable :: tmp(:)
   logical :: connected,fail
-  type(zmolecule) :: zmol
 
   real(wp) :: energy
   real(wp),allocatable :: grad(:,:),geo(:,:),csv(:,:)
-  integer,allocatable :: na(:),nb(:),nc(:),at2(:)
-  integer :: nat2
-  real(wp),allocatable :: mu(:)
-  real(wp) :: kappa,rrad
-
-  integer :: nsim
-  character(len=:),allocatable :: ensnam
-  type(mddata) :: mddat
-  type(mddata),allocatable :: mddats(:)
-
 
+  type(canonical_sorter) ::  can
 !========================================================================================!
   call tim%start(14,'Test implementation') 
 !========================================================================================!
@@ -88,23 +70,20 @@ subroutine crest_playground(env,tim)
   write(*,*) 
 !========================================================================================!
 
-!>--- sets the MD length according to a flexibility measure
-  call md_length_setup(env)
-!>--- create the MD calculator saved to env
-  call env_to_mddat(env)
-
-
-    nsim = 3
-    call crest_search_multimd_init(env,mol,mddat,nsim)
-    allocate (mddats(nsim), source=mddat)
-    call crest_search_multimd_init2(env,mddats,nsim)
+  allocate(grad(3,mol%nat), source=0.0_wp)
+  call env2calc(env,calc,mol)
+  calc%calcs(1)%rdwbo=.true.
+  call calc%info(stdout)
 
-    call crest_search_multimd(env,mol,mddats,nsim)
-!>--- a file called crest_dynamics.trj should have been written
-    ensnam = 'crest_dynamics.trj'
-!>--- deallocate for next iteration
-    if(allocated(mddats))deallocate(mddats)
+  call engrad(mol,calc,energy,grad,io)
+  call calculation_summary(calc,mol,energy,grad)
+
 
+  write(stdout,*)
+  write(stdout,*) 'CANGEN algorithm' 
+  call can%init(mol,calc%calcs(1)%wbo,'apsp+')
+  call can%stereo(mol)
+  call can%rankprint(mol) 
 
 !========================================================================================!
   call tim%stop(14)