Merge pull request #1540 from peastman/gbnerror

Fixed nan with GBn and GBn2

Merge pull request #1540 from peastman/gbnerror
Fixed nan with GBn and GBn2
37787af9 · peastman · GitHub · ac8452c7 · 162f5cc7 · 37787af9
Commit 37787af9 authored Jul 13, 2016 by peastman Committed by GitHub Jul 13, 2016
Hide whitespace changes
Inline Side-by-side

Showing with 2 additions and 2 deletions

wrappers/python/simtk/openmm/app/internal/customgbforces.py wrappers/python/simtk/openmm/app/internal/customgbforces.py +2 -2

No files found.
--- a/wrappers/python/simtk/openmm/app/internal/customgbforces.py
+++ b/wrappers/python/simtk/openmm/app/internal/customgbforces.py
@@ -760,7 +760,7 @@ class GBSAGBnForce(CustomAmberGBForceBase):
        self.addComputedValue("I", "Ivdw+neckScale*Ineck;"
                                   "Ineck=step(radius1+radius2+neckCut-r)*getm0(radindex1,radindex2)/(1+100*(r-getd0(radindex1,radindex2))^2+"
                                   "0.3*1000000*(r-getd0(radindex1,radindex2))^6);"
-                                   "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                   "Ivdw=select(step(r+sr2-or1), 0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r), 0);"
                                   "U=r+sr2;"
                                   "L=max(or1, D);"
                                   "D=abs(r-sr2);"
@@ -851,7 +851,7 @@ class GBSAGBn2Force(GBSAGBnForce):
        self.addComputedValue("I", "Ivdw+neckScale*Ineck;"
                                   "Ineck=step(radius1+radius2+neckCut-r)*getm0(radindex1,radindex2)/(1+100*(r-getd0(radindex1,radindex2))^2+"
                                   "0.3*1000000*(r-getd0(radindex1,radindex2))^6);"
-                                   "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                   "Ivdw=select(step(r+sr2-or1), 0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r), 0);"
                                   "U=r+sr2;"
                                   "L=max(or1, D);"
                                   "D=abs(r-sr2);"