modifications for safer thread handling (#289)

.github/workflows/build.yml

@@ -146,7 +146,7 @@ jobs:
       run: >-
         github-release ${{ env.RELEASE_COMMAND }}
         --tag ${{ env.RELEASE_TAG }}
-        --name "Continous release version"
+        --name "Continuous release version"
         --description "$DESCRIPTION"
         --pre-release
       env:

src/calculator/calculator.F90

@@ -138,8 +138,10 @@ subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
     iostatus = 0
     dum1 = 1.0_wp
     dum2 = 1.0_wp
+    if(n>0)then
     calc%etmp = 0.0_wp
     calc%grdtmp = 0.0_wp
+    endif
     !$omp end critical
 
 !**************************************
@@ -635,17 +637,33 @@ subroutine constrhess(nat,at,xyz,calc,phess)
     real(wp) :: energy,el,er
     real(wp),allocatable :: hess(:,:)
     logical :: consgeo
-    integer :: i,j,k,n3,io
+    integer :: i,j,k,n3,io,ncons
 
-    if (calc%nconstraints <= 0) return
+    ncons = calc%nconstraints
+    if (ncons <= 0) return
     !>--- skip if only nonadiabatic constraints
 
-    dummycalc = calc !> new dummy calculation
-    dummycalc%id = -1000  !> set to something arbitrary so that only constraints are considered
+    !> new dummy calculation
+    !> set to something arbitrary so that only constraints are considered
+    !> and copy the neccesities
+    dummycalc%id = -1000
     dummycalc%ncalculations = 0
     dummycalc%pr_energies = .false.
+
+    !> copy the constraints
+    dummycalc%nfreeze = calc%nfreeze 
+    dummycalc%nconstraints = ncons
+    !$omp critical
+    allocate(dummycalc%cons( ncons ))
+    do i=1,ncons
+      dummycalc%cons(i) = calc%cons(i)   
+    enddo
+    if(calc%nfreeze > 0)then
+     dummycalc%freezelist = calc%freezelist
+    endif
     n3 = nat*3
     allocate (hess(n3,n3),source=0.0_wp)
+    !$omp end critical
 
     call numhess1(nat,at,xyz,dummycalc,hess,io)
 
@@ -656,7 +674,7 @@ subroutine constrhess(nat,at,xyz,calc,phess)
         phess(k) = phess(k)+0.5_wp*(hess(j,i)+hess(i,j))
       end do
     end do
-
+    
     deallocate (hess)
     return
   end subroutine constrhess

src/ompmklset.F90

@@ -17,9 +17,9 @@
 ! along with crest.  If not, see <https://www.gnu.org/licenses/>.
 !================================================================================!
 
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 !c OMP and MKL parallelization settings
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 
 subroutine ompmklset(threads)
   use omp_lib
@@ -32,9 +32,9 @@ subroutine ompmklset(threads)
 #endif
 end subroutine ompmklset
 
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 !c OMP and MKL parallelization settings (short routine)
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 
 subroutine ompenvset(omp)
   use iomod
@@ -47,9 +47,10 @@ subroutine ompenvset(omp)
 
 end subroutine ompenvset
 
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 !c OMP and MKL autoset switchcase routine
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
+
 subroutine new_ompautoset(env,modus,maxjobs,parallel_jobs,cores_per_job)
   use omp_lib
   use crest_data
@@ -82,12 +83,10 @@ subroutine new_ompautoset(env,modus,maxjobs,parallel_jobs,cores_per_job)
       Tfloor = floor(Tfrac)
       cores_per_job = max(nint(Tfloor),1)
     end if
-    if (index(modus,'_nested') .ne. 0) then
+    if (index(modus,'_nested') .ne. 0 .and. cores_per_job > 1) then
       if (env%omp_allow_nested) then
         !> We should never need more than two active nested layers
         call omp_set_max_active_levels(2)
-        !else
-        !  parallel_jobs = T
       end if
     end if
 
@@ -115,9 +114,10 @@ subroutine new_ompautoset(env,modus,maxjobs,parallel_jobs,cores_per_job)
 
 end subroutine new_ompautoset
 
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 !c get omp/mkl automatically from the global variables
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
+
 subroutine ompgetauto(threads,omp,maxrun)
   use omp_lib
   use iomod
@@ -138,9 +138,9 @@ subroutine ompgetauto(threads,omp,maxrun)
 
 end subroutine ompgetauto
 
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 !c print omp/mkl threads that are used at the moment
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!
 subroutine ompprint_intern()
   use omp_lib
   implicit none