Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Allow adjusting of spring exponent #357

Merged
merged 2 commits into from
Oct 1, 2020
+36 −16
Merged

Allow adjusting of spring exponent #357

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
71ef89f
Allow adjusting of spring exponent
awvwgk Oct 1, 2020
fe06083
Fix expo writing
awvwgk Oct 1, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 26 additions & 4 deletions src/constrain_param.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/constrain_pot.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 0 additions & 2 deletions src/scanparam.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:", &
Expand Down
14 changes: 7 additions & 7 deletions src/type/setvar.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down