diff --git a/src/dynamic.f90 b/src/dynamic.f90
index 98bfb4338..170ceb4a2 100644
--- a/src/dynamic.f90
+++ b/src/dynamic.f90
@@ -147,7 +147,9 @@ subroutine md(env,mol,chk,calc, &
    use xtb_type_calculator
    use xtb_type_restart
    use xtb_type_data
+   use xtb_gfnff_calculator
    use xtb_shake, only: do_shake,ncons
+   use xtb_gfnff_shake, only: gff_do_shake => do_shake, gff_ncons => ncons
    use xtb_setparam
    use xtb_fixparam
    use xtb_scanparam
@@ -310,11 +312,20 @@ subroutine md(env,mol,chk,calc, &
       call rdmdrestart(mol%n,mol%xyz,velo)
    endif
 
-   if(set%shake_md) then
-      write(*,*) 'SHAKE on. # bonds  :',ncons,' all:',.not.set%xhonly
-   else
-      write(*,*) 'SHAKE off'
-   endif
+   select type(calc)
+      class default
+         if(set%shake_md) then
+            write(*,*) 'GFN-x SHAKE on. # bonds  :',ncons,' all:',.not.set%xhonly
+         else
+            write(*,*) 'GFN-x SHAKE off'
+         endif
+      type is(TGFFCalculator)
+         if (set%shake_md) then
+            write(*,*) 'GFNFF SHAKE on. # bonds  :',gff_ncons,' all:',.not.set%xhonly
+         else
+            write(*,*) 'GFNFF SHAKE off'
+         endif
+      end select
    if(thermostat)   write(*,*) 'Berendsen THERMOSTAT on'
    if(equi)       write(*,*) 'EQUILIBRATION mode'
    if( gmd)       write(*,*) 'GMD mode'
@@ -615,7 +626,14 @@ subroutine md(env,mol,chk,calc, &
       ! SHAKE (apply constraint at t+dt)
       !ccccccccccccccccccc
 
-      if(set%shake_md) call do_shake(mol%n,xyzo,mol%xyz,vel,acc,mass,tstep)
+      select type(calc)
+      class default
+         if(set%shake_md) call do_shake(mol%n,xyzo,mol%xyz,vel,acc,mass,tstep)
+      type is(TGFFCalculator)
+         if (set%shake_md) then
+            call gff_do_shake(mol%n,xyzo,mol%xyz,vel,acc,mass,tstep)
+         endif
+      end select
 
       ! update velocities
       velo = vel
diff --git a/src/gfnff/shake_module.f90 b/src/gfnff/shake_module.f90
index 89c12e338..9adb6e43c 100644
--- a/src/gfnff/shake_module.f90
+++ b/src/gfnff/shake_module.f90
@@ -16,20 +16,33 @@
 ! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
 
        module xtb_gfnff_shake
+       use xtb_param_atomicrad, only : atomicRad
+       use xtb_setparam, only: set
        implicit none
-       integer :: shake_mode
+       private :: atomicRad
+       integer, parameter :: ndim = 100000
+       logical :: pair_exists
        integer :: ncons = 0
        integer :: nconsu = 0
+       integer :: conslistu(2,ndim)
        integer, allocatable :: conslist(:,:)
        real*8,  allocatable :: distcons(:)
        real*8,  allocatable :: dro(:,:)
        real*8,  allocatable :: dr (:,:)
        integer, parameter :: maxcyc = 5000
        real*8, parameter :: tolshake = 1.d-6
+       integer, private, parameter :: metal(1:86) = [&
+      & 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, &
+      & 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, &
+      & 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, &
+      & 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
+      & 1, 0, 0, 0, 0, 0]
 
        contains
 
        subroutine init_shake(nat,at,xyz,topo)
+       use xtb_mctc_convert, only : autoaa
+       use xtb_fixparam, only : shakeset
        use xtb_gfnff_topology, only : TGFFTopology
        implicit none
        type(TGFFTopology), intent(in) :: topo
@@ -39,23 +52,92 @@ subroutine init_shake(nat,at,xyz,topo)
        integer :: iat,jat,i,j,jmin
        real*8  :: minrij,rij,wthr
        real*8  :: drij(3)
-
+       
        integer list(nat*(nat+1)/2),lin,ij
-
-      ncons  = topo%nbond
-      allocate(conslist(2,ncons),distcons(ncons),&
-     &         dro(3,ncons),dr(4,ncons))
-
-      conslist(1:2,1:ncons)=topo%blist(1:2,1:ncons)
-      do i = 1, ncons
-         iat = conslist(1,i)
-         jat = conslist(2,i)
-         drij(1:3) = xyz(1:3,iat) - xyz(1:3,jat)
-         distcons(i) = drij(1)**2 + drij(2)**2 + drij(3)**2
-      enddo
-
-      return
-      end subroutine init_shake
+       
+       real*8 rco
+       logical metalbond,rcut
+       
+       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       nconsu=0
+
+       ! constrain X-H only
+       if(set%shake_mode.eq.1)then
+          do i = 1, nat
+                if(at(i).eq.1) then
+                   minrij=1000.d0
+                   do j = 1, nat
+                      if(j.ne.i)then
+                            rij=(xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2
+                            if(rij.lt.minrij) then
+                               minrij=rij
+                               jmin=j
+                            endif
+                      endif
+                   enddo
+                   if(at(jmin).eq.1) then
+                      if(jmin.gt.i) then
+                            nconsu = nconsu + 1
+                            conslistu(1,nconsu) = i
+                            conslistu(2,nconsu) = jmin
+                      endif
+                   else
+                      nconsu = nconsu + 1
+                      conslistu(1,nconsu) = i
+                      conslistu(2,nconsu) = jmin
+                   endif
+                endif
+          enddo
+       endif
+       ! all bonds
+       if(set%shake_mode.eq.2)then
+          nconsu  = topo%nbond
+          conslistu(1:2,1:nconsu)=topo%blist(1:2,1:nconsu)
+       endif
+       print *, "SHAKE is ", set%shake_md, " ncons = ", nconsu, "shake_mode = ", set%shake_mode
+       
+       if(shakeset%n.gt.0)then
+             do i = 1, shakeset%n, 2
+                pair_exists = .false.
+                do j = 1, nconsu
+                   if ((conslistu(1,j) == shakeset%atoms(i) .and. conslistu(2,j) == shakeset%atoms(i+1)) .or. &
+                      (conslistu(1,j) == shakeset%atoms(i+1) .and. conslistu(2,j) == shakeset%atoms(i))) then
+                      pair_exists = .true.
+                      if (pair_exists) write(*,*) 'input pair',  shakeset%atoms(i), '-', shakeset%atoms(i+1), 'for shake is already considered.'
+                      exit
+                   endif
+                enddo
+                
+                if (.not. pair_exists) then
+                   nconsu = nconsu + 1
+                   conslistu(1,nconsu) = shakeset%atoms(i)
+                   conslistu(2,nconsu) = shakeset%atoms(i+1)
+                   write(*,*) 'SHAKE user input constraining bond ', &
+                      &      conslistu(1:2,nconsu)
+                endif
+             enddo
+       endif
+       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+99     ncons  = nconsu
+       if(nconsu.lt.1) then
+          set%shake_md = .false.
+          print *, "gfnff shake_md turned off"
+          return
+       endif
+     print *, "SHAKE is ", set%shake_md, " ncons = ", ncons, "shake_mode = ", set%shake_mode
+       allocate(conslist(2,ncons),distcons(ncons),&
+       &         dro(3,ncons),dr(4,ncons))
+    
+       conslist(1:2,1:ncons)=conslistu(1:2,1:ncons) 
+       do i = 1, ncons
+          iat = conslist(1,i)
+          jat = conslist(2,i)
+          drij(1:3) = xyz(1:3,iat) - xyz(1:3,jat)
+          distcons(i) = drij(1)**2 + drij(2)**2 + drij(3)**2
+       enddo
+       
+       return
+       end subroutine init_shake
 
        subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep)
        implicit none
@@ -87,8 +169,7 @@ subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep)
          dro(1:3,i) = xyzo(1:3,iat) - xyzo(1:3,jat)
         enddo
 
-!100     continue
-        do  while (.not.conv.and.icyc.le.maxcyc) !goto 100
+100     continue
 
           maxdev = 0.d0
 
@@ -132,8 +213,8 @@ subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep)
 
           icyc = icyc + 1
 
-        end do
-        !if(.not.conv.and.icyc.le.maxcyc) goto 100
+        
+        if(.not.conv.and.icyc.le.maxcyc) goto 100
 
         if(conv) then
          velo = velo + (xyzt-xyz)*tau1
diff --git a/src/setparam.f90 b/src/setparam.f90
index f88c9feeb..39d2087c6 100644
--- a/src/setparam.f90
+++ b/src/setparam.f90
@@ -312,8 +312,8 @@ module xtb_setparam
    integer  :: md_hmass = 4
 !  shake on (=0: off which is default) for X-H bonds only (=1),
 !  or all bonds (=2) or user defined bonds (=3)
-   integer  :: shake_mode = 2
-   logical  :: shake_md = .true.
+   integer  :: shake_mode = 0
+   logical  :: shake_md = .false.
    logical  :: xhonly = .true.
    logical  :: honly = .false.
    logical :: forcewrrestart = .false.