fix bug in scaling of KUNDROT volume derivatives for GK

mmff/0README

@@ -13,16 +13,16 @@
                                Usage Instructions
                                ##################
 
-     If MMFF executables are available, you may use the mmff.prm parameter
-     file just as you would any other TINKER parameter file. The only caveat
-     is that all MMFF "bond type = 1" bonds must be listed in the Tinker
-     keyfile via the MMFF-PIBOND keyword. The keyword takes as arguments the
-     atom numbers of atoms in type=1 bonds. Up to ten bonded atom pairs (ie,
-     20 atom numbers) are listed following the MMFF-PIBOND keyword. Multiple
-     MMFF-PIBOND keywords can be present. The MMFF bond type=1 is used for
-     conjugated bonds that are "single" bonds in reasonable resonance forms,
-     but could otherwise be double bonds in other structures. For example,
-     the middle C-C bond in 1,3-butadiene has an MMFF bond type of 1.
+     The mmff.prm parameter file may be used just as you would any other
+     Tinker parameter file. The only caveat is that all MMFF "bond type = 1"
+     bonds must be listed in the Tinker keyfile via the MMFF-PIBOND keyword.
+     This keyword takes as arguments the atom numbers of atoms in type=1
+     bonds. Up to ten bonded atom pairs (ie, 20 atom numbers) are listed
+     following the MMFF-PIBOND keyword. Multiple MMFF-PIBOND keywords can
+     be present. The MMFF bond type=1 is used for conjugated bonds that are
+     "single" bonds in reasonable resonance forms, but could otherwise be
+     double bonds in other structures. For example, the middle C-C bond in
+     1,3-butadiene has an MMFF bond type of 1.
 
                              #####################
                              Validation Test Cases

source/surface.f

@@ -52,21 +52,28 @@ subroutine surface (rad,weight,probe,surf,asurf)
       integer nsize,nfudge
       integer nredundant
       integer, allocatable :: redlist(:)
-      real*8 surf,usurf
+      real*8 surf,usurf,eps
       real*8 probe,alpha
       real*8 rad(*)
       real*8 weight(*)
       real*8 asurf(*)
       real*8, allocatable :: radii(:)
       real*8, allocatable :: asurfx(:)
       real*8, allocatable :: coords(:,:)
+      logical dowiggle
       character*6 symmtyp
 c
 c
-c     use Richmond method for small symmetric structures
+c     check coordinates for linearity, planarity and symmetry
 c
+      symmtyp = 'NONE'
       call chksymm (symmtyp)
-      if (n.le.65 .and. symmtyp.ne.'NONE') then
+      dowiggle = .false.
+      if (n.le.200 .and. symmtyp.ne.'NONE')  dowiggle = .true.
+c
+c     use Richmond method for small symmetric structures
+c
+      if (dowiggle) then
          call richmond (n,x,y,z,rad,weight,probe,surf,asurf)
          return
       end if
@@ -91,6 +98,13 @@ subroutine surface (rad,weight,probe,surf,asurf)
          if (rad(i) .ne. 0.0d0)  radii(i) = rad(i) + probe
       end do
 c
+c     random coordinate perturbation to avoid numerical issues
+c
+      if (dowiggle) then
+         eps = 0.001d0
+         call wiggle (n,coords,eps)
+      end if
+c
 c     transfer coordinates, complete to minimum of four spheres
 c     if needed, set Delaunay and alpha complex arrays
 c
@@ -179,7 +193,7 @@ subroutine surface1 (rad,weight,probe,surf,asurf,dsurf)
       integer nsize,nfudge
       integer nredundant
       integer, allocatable :: redlist(:)
-      real*8 surf,usurf
+      real*8 surf,usurf,eps
       real*8 probe,alpha
       real*8 rad(*)
       real*8 weight(*)
@@ -189,13 +203,20 @@ subroutine surface1 (rad,weight,probe,surf,asurf,dsurf)
       real*8, allocatable :: asurfx(:)
       real*8, allocatable :: coords(:,:)
       real*8, allocatable :: dsurfx(:,:)
+      logical dowiggle
       character*6 symmtyp
 c
 c
-c     use Richmond method for small symmetric structures
+c     check coordinates for linearity, planarity and symmetry
 c
+      symmtyp = 'NONE'
       call chksymm (symmtyp)
-      if (n.le.65 .and. symmtyp.ne.'NONE') then
+      dowiggle = .false.
+      if (n.le.200 .and. symmtyp.ne.'NONE')  dowiggle = .true.
+c
+c     use Richmond method for small symmetric structures
+c
+      if (dowiggle) then
          call richmond1 (n,x,y,z,rad,weight,probe,surf,asurf,dsurf)
          return
       end if
@@ -221,6 +242,13 @@ subroutine surface1 (rad,weight,probe,surf,asurf,dsurf)
          if (rad(i) .ne. 0.0d0)  radii(i) = rad(i) + probe
       end do
 c
+c     random coordinate perturbation to avoid numerical issues
+c
+      if (dowiggle) then
+         eps = 0.001d0
+         call wiggle (n,coords,eps)
+      end if
+c
 c     transfer coordinates, complete to minimum of four spheres
 c     if needed, set Delaunay and alpha complex arrays
 c

source/unionball.f

@@ -110,21 +110,19 @@ subroutine unionball (n,x,y,z,rad,weight,probe,doderiv,dovol,
          if (rad(i) .ne. 0.0d0)  radii(i) = rad(i) + probe
       end do
 c
-c     test for symmetrical system with possible numerical issues
+c     check coordinates for linearity, planarity and symmetry
 c
       symmtyp = 'NONE'
       call chksymm (symmtyp)
-      if (n.le.65 .and. symmtyp.ne.'NONE') then
-         write (iout,10)  symmtyp
-   10    format (/,' UNIONBALL  --  Warning, ',a6,' Symmetry;'
-     &              ' Wiggling Coordinates')
-      end if
+      dowiggle = .false.
+      if (n.le.200 .and. symmtyp.ne.'NONE')  dowiggle = .true.
 c
 c     random coordinate perturbation to avoid numerical issues
 c
-      dowiggle = .false.
-      if (n.le.65 .and. symmtyp.ne.'NONE')  dowiggle = .true.
       if (dowiggle) then
+         write (iout,10)  symmtyp
+   10    format (/,' UNIONBALL  --  Warning, ',a6,' Symmetry;'
+     &              ' Wiggling Coordinates')
          eps = 0.001d0
          call wiggle (n,coords,eps)
       end if

source/volume.f

@@ -55,7 +55,7 @@ subroutine volume (rad,weight,probe,surf,vol,asurf,avol)
       integer nsize,nfudge
       integer nredundant
       integer, allocatable :: redlist(:)
-      real*8 surf,usurf
+      real*8 surf,usurf,eps
       real*8 vol,uvol,voln
       real*8 reentrant
       real*8 probe,alpha
@@ -67,13 +67,20 @@ subroutine volume (rad,weight,probe,surf,vol,asurf,avol)
       real*8, allocatable :: asurfx(:)
       real*8, allocatable :: avolx(:)
       real*8, allocatable :: coords(:,:)
+      logical dowiggle
       character*6 symmtyp
 c
 c
-c     use Connolly method for small symmetric structures
+c     check coordinates for linearity, planarity and symmetry
 c
+      symmtyp = 'NONE'
       call chksymm (symmtyp)
-      if (n.le.65 .and. symmtyp.ne.'NONE') then
+      dowiggle = .false.
+      if (n.le.200 .and. symmtyp.ne.'NONE')  dowiggle = .true.
+c
+c     use Connolly method for small symmetric structures
+c
+      if (dowiggle) then
          reentrant = 0.0d0
          call connolly (n,x,y,z,rad,probe,reentrant,surf,vol)
          vol = vol * weight(1)
@@ -105,6 +112,13 @@ subroutine volume (rad,weight,probe,surf,vol,asurf,avol)
          if (rad(i) .ne. 0.0d0)  radii(i) = rad(i) + probe
       end do
 c
+c     random coordinate perturbation to avoid numerical issues
+c
+      if (dowiggle) then
+         eps = 0.001d0
+         call wiggle (n,coords,eps)
+      end if
+c
 c     transfer coordinates, complete to minimum of four spheres
 c     if needed, set Delaunay and alpha complex arrays
 c
@@ -202,7 +216,7 @@ subroutine volume1 (rad,weight,probe,surf,vol,asurf,avol,
       integer nsize,nfudge
       integer nredundant
       integer, allocatable :: redlist(:)
-      real*8 surf,usurf
+      real*8 surf,usurf,eps
       real*8 vol,uvol,voln
       real*8 reentrant
       real*8 probe,alpha
@@ -218,20 +232,30 @@ subroutine volume1 (rad,weight,probe,surf,vol,asurf,avol,
       real*8, allocatable :: coords(:,:)
       real*8, allocatable :: dsurfx(:,:)
       real*8, allocatable :: dvolx(:,:)
+      logical dowiggle
       character*6 symmtyp
 c
 c
-c     use other methods for small symmetric structures
+c     check coordinates for linearity, planarity and symmetry
 c
+      symmtyp = 'NONE'
       call chksymm (symmtyp)
-      if (n.le.65 .and. symmtyp.ne.'NONE') then
+      dowiggle = .false.
+      if (n.le.200 .and. symmtyp.ne.'NONE')  dowiggle = .true.
+c
+c     use arc-based methods for small symmetric structures
+c
+      if (dowiggle) then
          reentrant = 0.0d0
          call connolly (n,x,y,z,rad,probe,reentrant,surf,vol)
          call kundrot1 (n,x,y,z,rad,probe,dvol)
          call richmond1 (n,x,y,z,rad,weight,probe,surf,asurf,dsurf)
          vol = vol * weight(1)
          voln = vol / dble(n)
          do i = 1, n
+            dvol(1,i) = dvol(1,i) * weight(i)
+            dvol(2,i) = dvol(2,i) * weight(i)
+            dvol(3,i) = dvol(3,i) * weight(i)
             avol(i) = voln
          end do
          return
@@ -260,6 +284,13 @@ subroutine volume1 (rad,weight,probe,surf,vol,asurf,avol,
          if (rad(i) .ne. 0.0d0)  radii(i) = rad(i) + probe
       end do
 c
+c     random coordinate perturbation to avoid numerical issues
+c
+      if (dowiggle) then
+         eps = 0.001d0
+         call wiggle (n,coords,eps)
+      end if
+c
 c     transfer coordinates, complete to minimum of four spheres
 c     if needed, set Delaunay and alpha complex arrays
 c