Allow adjusting of spring exponent (#357)

src/constrain_param.f90

@@ -462,14 +462,14 @@ subroutine set_constr(env,key,val,nat,at,idMap,xyz)
       potset%pos%n = size(list)
 
    case('DISTANCE')
-      if (narg.ne.3) then
+      if (narg.lt.3 .or. narg.gt.4) then
          call env%error('not enough arguments to constrain a distance',source)
          return
       endif
       ioffset = 2*potset%dist%n
       potset%dist%n = potset%dist%n+1
    !  part 1: get the constrained atoms
-      do i = 1, narg-1
+      do i = 1, 2
          if (getValue(env,trim(argv(i)),idum)) then
             potset%dist%atoms(ioffset+i) = idum
          else
@@ -482,17 +482,39 @@ subroutine set_constr(env,key,val,nat,at,idMap,xyz)
       i = potset%dist%atoms(ioffset+1)
       j = potset%dist%atoms(ioffset+2)
       dist = sqrt(sum((xyz(:,i)-xyz(:,j))**2))
-      if (trim(argv(narg)).eq.'auto') then
+      if (trim(argv(3)).eq.'auto') then
          potset%dist%val(potset%dist%n) = dist
       else
-         if (getValue(env,trim(argv(narg)),ddum)) then
+         if (getValue(env,trim(argv(3)),ddum)) then
             potset%dist%val(potset%dist%n) = ddum * aatoau
          else
             call env%warning("Something went wrong in set_constr_ 'distance'. (2)",source)
             potset%dist%n = potset%dist%n-1 ! remove invalid contrain
             return
          endif
       endif
+      if (narg.eq.4) then
+         if (getValue(env,trim(argv(4)),idum)) then
+            if (idum < 2 .or. mod(idum, 2).ne.0) then
+               call env%warning("Invalid spring exponent given", source)
+               potset%dist%n = potset%dist%n-1 ! remove invalid contrain
+               return
+            end if
+            potset%dist%expo(potset%dist%n) = real(idum, wp)
+         else
+            call env%warning("Something went wrong in set_constr_ 'distance'. (3)",source)
+            potset%dist%n = potset%dist%n-1 ! remove invalid contrain
+            return
+         end if
+         write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//&
+            '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å",1x,"with expo",1x,i0)') &
+            i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa, &
+            nint(potset%dist%expo(potset%dist%n))
+      else
+         write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//&
+            '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å")') &
+            i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa
+      end if
 
    case('ANGLE')
       if (narg.ne.4) then

src/constrain_pot.f90

@@ -193,10 +193,10 @@ subroutine constrain_dist(fix,n,at,xyz,g,e)
       d=r-r0
 !     e=e+fix%fc*d*d
 !     ff=fix%fc*2.0d0*d
-      dum= d**fix%expo
-      dum2=d**(fix%expo-1.0_wp)
+      dum= d**fix%expo(m)
+      dum2=d**(fix%expo(m)-1.0_wp)
       e=e+fix%fc*dum
-      ff=fix%fc*fix%expo*dum2
+      ff=fix%fc*fix%expo(m)*dum2
       dum=ff/r
       g(:,j)=g(:,j)+dum*rij
       g(:,i)=g(:,i)-dum*rij

src/scanparam.f90

@@ -228,8 +228,6 @@ subroutine pot_info(iunit,n,at,xyz)
             val0*autoaa, val*autoaa
       enddo
       write(iunit,'(a)')
-      write(iunit,'(5x,a,f14.7)') "  constraining potential exponent:", &
-         potset%dist%expo
       write(iunit,'(5x,a,f14.7)') " applied force constant per dist.:", &
          potset%dist%fc
       write(iunit,'(5x,a,f14.7)') "effective force constant per atom:", &

src/type/setvar.f90

@@ -111,7 +111,7 @@ module xtb_type_setvar
       integer  :: n    = 0
       integer  :: nval = 0
       real(wp) :: fc   = 0.0_wp
-      real(wp) :: expo = 0.0_wp
+      real(wp),allocatable :: expo(:)
       integer, allocatable :: atoms(:)
       real(wp),allocatable :: val(:)
    contains
@@ -356,7 +356,10 @@ subroutine allocate_fix(self,nat,nval,fc,expo)
    call self%deallocate
    if (present(nval)) self%nval = nval
    if (present(fc))   self%fc = fc
-   if (present(expo)) self%expo = expo
+   if (present(nval) .and. present(expo)) then
+      allocate(self%expo(nval))
+      self%expo(:) = expo
+   end if
    allocate( self%atoms(nat), source = 0 )
    if (present(nval)) allocate( self%val(nval),  source = 0.0_wp )
 end subroutine allocate_fix
@@ -365,7 +368,7 @@ subroutine deallocate_fix(self)
    class(fix_setvar) :: self
    self%n = 0
    self%fc = 0.0_wp
-   self%expo = 0.0_wp
+   if(allocated(self%expo)) deallocate( self%expo )
    if(allocated(self%atoms)) deallocate( self%atoms )
    if(allocated(self%val))   deallocate( self%val )
 end subroutine deallocate_fix
@@ -383,7 +386,7 @@ subroutine allocate_constr(self,nat,nval,fc,expo)
    allocate( self%lookup(nval), source = 0 )
    allocate( self%typeid(nval), source = 0 )
    call self%pos     %allocate(nat,nat*(nat+1)/2)
-   call self%dist    %allocate(nval*2,nval)
+   call self%dist    %allocate(nval*2,nval,expo=expo)
    call self%angle   %allocate(nval*3,nval)
    call self%dihedral%allocate(nval*4,nval)
    if (present(fc)) then
@@ -392,9 +395,6 @@ subroutine allocate_constr(self,nat,nval,fc,expo)
       self%angle   %fc = fc
       self%dihedral%fc = fc
    endif
-   if (present(expo)) then
-      self%dist%expo = expo
-   endif
 end subroutine allocate_constr
 
 subroutine deallocate_constr(self)