Page 1 of 1

bugs in neb code

Posted: Mon Apr 01, 2019 2:00 pm
by ezanardi
I've found two bugs in neb code

1) In tsase/neb/util.py vproj(v1,v2) it has to divide by the square of the modulus of v2 ,
but currently it divides by the modulus of v2. A patch to fix that bug (just changing vmag to vmag2) is:

[code]
diff a/tsase/neb/util.py b/tsase/neb/util.py
--- a/tsase/neb/util.py
+++ b/tsase/neb/util.py
@@ -61,7 +61,7 @@ def vproj(v1, v2):
Parameters:
v1, v2: numpy vectors
"""
- mag2 = vmag(v2)
+ mag2 = vmag2(v2)
if mag2 == 0:
printf("Can't project onto a zero vector", ERR)
return v1
[/code]

2) In both tsase/neb/pssneb.py and tsase/neb/ssneb.py forces calculation, when using dneb modification,
it has to multiply fsdneb by 2/pi * atan( |Fperp|²/|Fsperp|² ) but currently it multiplies by
2/pi * atan( |Fperp|/|Fsperp| ) (no squares). A patch to fix that bug (just changing vmag to vmag2) is:

[code]
diff a/tsase/neb/pssneb.py b/tsase/neb/pssneb.py
--- a/tsase/neb/pssneb.py
+++ b/tsase/neb/pssneb.py
@@ -334,8 +334,8 @@ class pssneb:

# New dneb where dneb force converges with (What?!)
if not self.dnebOrg:
- FperpSQ = vmag(self.path[i].fPerp)
- FsperpSQ = vmag(self.path[i].fsperp)
+ FperpSQ = vmag2(self.path[i].fPerp)
+ FsperpSQ = vmag2(self.path[i].fsperp)
if FsperpSQ > 0:
self.path[i].fsdneb *= 2.0 / pi * atan(FperpSQ / \
FsperpSQ)
diff a/tsase/neb/ssneb.py b/tsase/neb/ssneb.py
--- a/tsase/neb/ssneb.py
+++ b/tsase/neb/ssneb.py
@@ -409,8 +409,8 @@ class ssneb:

# dneb modification so that it will converge
if not self.dnebOrg:
- FperpSQ = vmag(self.path[i].fPerp)
- FsperpSQ = vmag(self.path[i].fsperp)
+ FperpSQ = vmag2(self.path[i].fPerp)
+ FsperpSQ = vmag2(self.path[i].fsperp)
if FsperpSQ > 0:
self.path[i].fsdneb *= 2.0/pi*atan(FperpSQ/FsperpSQ)

[/code]

Re: bugs in neb code

Posted: Mon Apr 01, 2019 6:08 pm
by graeme
Wow, thank you very much! The error in the projection algorithm is particularly worrisome. Anyway, again, thank you for the fix. I've updated the code.

Re: bugs in neb code

Posted: Fri Jul 26, 2019 10:30 pm
by xph
I was worried in the beginning, but later found out that all the "vproj" calls are projecting forces on to a unit vector (tangent direction), except in dneb. So dneb was totally messed up, but the regular (ss)neb was good.