[pull] master from crest-lab:master

src/calculator/calc_type.f90

@@ -97,6 +97,8 @@ module calc_type
     character(len=:),allocatable :: shortflag   !> shorter job description
 
 !>--- gradient format specifications
+    logical :: numgrad = .false.      !> run numerical gradient (expensive!)
+    real(wp) :: gradstep = 0.0005_wp  !> displacement for numerical gradient
     logical :: rdgrad = .true.
     integer :: gradtype = 0
     integer :: gradfmt = 0
@@ -1091,6 +1093,9 @@ subroutine calculation_settings_info(self,iunit)
     character(len=*),parameter :: fmt3 = '(" :",2x,a20," : ",a)'
     character(len=*),parameter :: fmt4 = '(" :",1x,a)'
     character(len=20) :: atmp
+    logical :: gxtbwarn
+
+    gxtbwarn=.false.
 
     if (allocated(self%description)) then
       write (iunit,'(" :",1x,a)') trim(self%description)
@@ -1111,7 +1116,12 @@ subroutine calculation_settings_info(self,iunit)
     end if
     if (any((/jobtype%orca,jobtype%xtbsys,jobtype%turbomole, &
     &  jobtype%generic,jobtype%terachem/) == self%id)) then
-      write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary)
+      if(index(self%binary,'gxtb').ne.0)then
+        write(iunit,fmt4) 'g-xTB (development version)'
+        gxtbwarn = .true.
+      else  
+        write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary)
+      endif
     end if
     if (self%refine_lvl > 0) then
       write (atmp,*) 'refinement stage'
@@ -1169,6 +1179,11 @@ subroutine calculation_settings_info(self,iunit)
       endif
     end if
 
+    if(gxtbwarn)then
+       write(iunit,fmt4) 'WARNING: This currently is the development version of g-xTB.'
+       write(iunit,fmt4) 'WARNING: Gradients are NUMERICAL (i.e., expensive and noisy!)' 
+    endif
+
   end subroutine calculation_settings_info
 
 !=========================================================================================!
@@ -1199,10 +1214,15 @@ subroutine create_calclevel_shortcut(self,levelstring)
       self%id = jobtype%turbomole
       self%rdgrad = .false.
       self%binary = 'gp3'
-    case ('gxtb')
+    case ('gxtb','gxtb_dev')
       self%id = jobtype%turbomole 
       self%rdgrad = .false.       
-      self%binary = 'gxtb'         
+      self%binary = 'gxtb'
+      self%rdwbo = .false.
+      if(index(levelstring,'_dev').ne.0)then
+        self%other = '-grad'
+        self%rdgrad=.true.
+      endif  
     case ('orca')
       self%id = jobtype%orca
 

src/calculator/calculator.F90

@@ -179,6 +179,9 @@ subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
         !==========================================!
         call potential_core(molptr,calc,i,iostatus)
         !==========================================!
+        !> and numerical gradient, if selected
+        !==========================================!
+        call numgrad_core(molptr,calc,i,iostatus)
         !==========================================!
 
         if (iostatus /= 0) then
@@ -382,6 +385,68 @@ subroutine potential_core(molptr,calc,id,iostatus)
 
   end subroutine potential_core
 
+  subroutine numgrad_core(molptr,calc,id,iostatus)
+!*******************************************************
+!* subroutine numgrad
+!* routine to perform a numerical gradient calculation
+!*******************************************************
+    implicit none
+    type(coord),intent(in)       :: molptr
+    type(calcdata),intent(inout) :: calc
+    integer,intent(in)           :: id
+    integer,intent(out)          :: iostatus
+
+    integer :: i,j,k,l,ich,och,io,pnat
+    type(coord),allocatable :: moltmp
+    real(wp) :: energy,el,er, step,step2
+    real(wp),allocatable :: ngrd(:,:)
+    !real(wp),parameter :: step = 0.0005_wp
+    !real(wp),parameter :: step2 = 0.5_wp/step
+
+    if (id > calc%ncalculations) return
+    if (.not.calc%calcs(id)%numgrad) return
+
+    pnat = molptr%nat
+    step = calc%calcs(id)%gradstep
+    step2 = 0.5_wp/step
+
+    !> back up energy 
+    energy = calc%etmp(id)
+
+    !> allocate temprorary gradient space
+    !$omp critical
+    allocate(ngrd(3,pnat), source=0.0_wp)
+    allocate(moltmp, source=molptr)
+    !$omp end critical
+
+    do i = 1,molptr%nat
+      do j = 1,3
+        moltmp%xyz(j,i) = moltmp%xyz(j,i)+step
+        call potential_core(moltmp,calc,id,iostatus)
+        er = calc%etmp(id)
+
+        moltmp%xyz(j,i) = moltmp%xyz(j,i)-2*step
+        call potential_core(moltmp,calc,id,iostatus)
+        el = calc%etmp(id)
+
+        moltmp%xyz(j,i) = moltmp%xyz(j,i)+step
+        ngrd(j,i) = step2*(er-el)
+      end do
+    end do
+
+    !> transfer tmp gradient to the calc object
+    calc%grdtmp(:,1:pnat,id) = ngrd(:,1:pnat)
+    !$omp critical
+    deallocate(moltmp)
+    deallocate(ngrd)
+    !$omp end critical
+
+    !> restore the energy
+    calc%etmp(id) = energy
+
+    return
+  end subroutine numgrad_core
+
 !========================================================================================!
 !========================================================================================!
 !========================================================================================!

src/calculator/gradreader.f90

@@ -176,27 +176,30 @@ subroutine rd_grad_tm(iunit,nat,energy,grad,iostatus)
     integer,intent(out) :: iostatus
     integer :: c,io,n,i,j
     character(len=128) :: atmp
-    character(len=20) :: btmp(8)
+    character(len=20) :: btmp(10)
     real(wp) :: dum
     logical :: readblock
 
     iostatus = 0
     energy = 0.0_wp
     grad(:,:) = 0.0_wp
 
-    c = 0
+    c = 1
     readblock = .false.
     do
       read (iunit,'(a)',iostat=io) atmp
       if (io < 0) exit !> EOF exit
       atmp = adjustl(atmp)
       if (atmp(1:4) == '$end') readblock = .false.
       if( readblock ) then      
+
         if(index(atmp,'cycle').ne.0)then
-          read(atmp,*) btmp(1:2),j,btmp(3:6),energy,btmp(7:8),dum
+          read(atmp,*) btmp(1:8)
+          read(btmp(7),*) energy
         elseif(c < nat)then !> skip coords
           c = c + 1  
         else !> read grad
+          !backspace(iunit)
           call rd_grad_n3(iunit,nat,grad,iostatus)
           exit
         endif
@@ -270,7 +273,7 @@ subroutine rd_grad_3n(iunit,nat,grad,iostatus)
     grad(:,:) = 0.0_wp
 
     c = 0
-    do i = 1,n
+    do i = 1,nat
       do j = 1,3
         read (iunit,*,iostat=io) dum
         if (io < 0) then
@@ -301,7 +304,7 @@ subroutine rd_grad_n3(iunit,nat,grad,iostatus)
     grad(:,:) = 0.0_wp
 
     c = 0
-    do i = 1,n
+    do i = 1,nat
       read (iunit,*,iostat=io) dum(1:3)
       if (io < 0) then
         iostatus = 3

src/calculator/printouts.F90

@@ -131,7 +131,8 @@ subroutine calculation_summary(calc,mol,energy,grad,molnew,iounit,print)
     end if
 
 !>--- gradients
-    if (all(calc%calcs(:)%rdgrad.eqv..false.)) then
+    if (all(calc%calcs(:)%rdgrad.eqv..false.) .and. & 
+     & all(calc%calcs(:)%numgrad.eqv..false.) ) then
       write (iunit,*)
       write (iunit,'(a)') '> No gradients calculated'
     else if (present(grad)) then