Repository: materialsproject/custodian
Pull request #184 (2021-11-15 05:00:36)

[VASP] More robust zbrent fix
`zbrent`: Currently, Custodian switches to IBRION = 1 and copies CONTCAR to POSCAR. For materials that are large and have flexible structures, this is rarely sufficient.* I do not change the original approach here, but rather I have added a _second_ fallback option of tightening EDIFF (per the recommendation in the stdout when this error occurs, shown below) and increasing NELMIN to 6 (i.e. between 4 and 8 as noted in the [VASP manual](https://www.vasp.at/wiki/index.php/IBRION#IBRION.3D1:_ionic_relaxation_.28RMM-DIIS.29.)) in this PR. These ensure that the energy/forces are of higher quality. See the two quotes from the stdout and VASP manual below for context:

```
 curvature:   0.00 expect dE= 0.000E+00 dE for cont linesearch  0.000E+00
 ZBRENT: fatal error in bracketing
     please rerun with smaller EDIFF, or copy CONTCAR
     to POSCAR and continue 
```


> Since the algorithm [IBRION=1] builds up an approximation of the Hessian matrix it requires very accurate forces, otherwise it will fail to converge. An efficient way to achieve this is to set NELMIN to a value between 4 and 8 (for simple bulk materials 4 is usually adequate, whereas 8 might be required for complex surfaces where the charge density converges very slowly). This forces a minimum of 4 to 8 electronic steps between each ionic step, and guarantees that the forces are well converged at each step.

If the new `vtst_fixes` argument is set to `True` (default: `False`), it will also switch to the [FIRE](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.97.170201) optimizer when correcting a second failed `zbrent` run. I have run into this error hundreds of times and had a near-100% success rate with FIRE, but it's not enabled by default for compatibility purposes and in case others disagree with this decision for some reason.

Custodian didn't have any tests at all for `zbrent`, so I have added those too.

*This is based on personal experience with thousands of calculations on porous materials (a large portion of structures at https://materialsproject.org/mofs run into this issue and aren't fixed by a simple CONTCAR --> POSCAR swap). There is still too much numerical noise in the energy when this error occurs near the local minimum in the PES.


Patch:
diff --git a/custodian/vasp/handlers.py b/custodian/vasp/handlers.py
index 48fc9b80..de450d3c 100644
--- a/custodian/vasp/handlers.py
+++ b/custodian/vasp/handlers.py
@@ -104,6 +104,7 @@ class VaspErrorHandler(ErrorHandler):
         output_filename="vasp.out",
         natoms_large_cell=100,
         errors_subset_to_catch=None,
+        vtst_fixes=False,
     ):
         """
         Initializes the handler with the output file to check.
@@ -123,6 +124,8 @@ class VaspErrorHandler(ErrorHandler):
                 errors, and not others. If you wish to only excluded one or
                 two of the errors, you can create this list by the following
                 lines:
+            vtst_fixes (bool): Whether to consider VTST optimizers. Defaults to
+                False for compatibility purposes.
 
                 ```
                 subset = list(VaspErrorHandler.error_msgs.keys())
@@ -136,7 +139,8 @@ class VaspErrorHandler(ErrorHandler):
         self.error_count = Counter()
         # threshold of number of atoms to treat the cell as large.
         self.natoms_large_cell = natoms_large_cell
-        self.errors_subset_to_catch = errors_subset_to_catch or list(self.error_msgs.keys())
+        self.errors_subset_to_catch = errors_subset_to_catch or list(VaspErrorHandler.error_msgs.keys())
+        self.vtst_fixes = vtst_fixes
         self.logger = logging.getLogger(self.__class__.__name__)
 
     def check(self):
@@ -366,8 +370,26 @@ class VaspErrorHandler(ErrorHandler):
             self.error_count["brions"] += 1
 
         if "zbrent" in self.errors:
-            actions.append({"dict": "INCAR", "action": {"_set": {"IBRION": 1}}})
+            # ZBRENT is caused by numerical noise in the forces, often near the PES minimum
+            # This is often a severe problem for systems with many atoms and flexible
+            # structures (e.g. zeolites, MOFs)
             actions.append({"file": "CONTCAR", "action": {"_file_copy": {"dest": "POSCAR"}}})
+            if self.error_count["zbrent"] == 0:
+                # First try changing IBRION to 1 and continuing
+                actions.append({"dict": "INCAR", "action": {"_set": {"IBRION": 1}}})
+            elif self.error_count["zbrent"] == 1:
+                # If that fails, tighten the energy convergence criteria and raise minimum number of SCF steps
+                ediff = vi["INCAR"].get("EDIFF", 1e-4)
+                if ediff > 1e-6:
+                    actions.append({"dict": "INCAR", "action": {"_set": {"EDIFF": 1e-6}}})
+                else:
+                    actions.append({"dict": "INCAR", "action": {"_set": {"EDIFF": ediff / 10}}})
+                if vi["INCAR"].get("NELMIN", 2) < 6:
+                    actions.append({"dict": "INCAR", "action": {"_set": {"NELMIN": 6}}})
+                if self.vtst_fixes is True:
+                    # FIRE almost always resolves this issue but requires VTST to be installed
+                    actions.append({"dict": "INCAR", "action": {"_set": {"IOPT": 7, "IBRION": 3, "POTIM": 0}}})
+            self.error_count["zbrent"] += 1
 
         if "too_few_bands" in self.errors:
             if "NBANDS" in vi["INCAR"]:


Test patch:
diff --git a/custodian/vasp/tests/test_handlers.py b/custodian/vasp/tests/test_handlers.py
index 6c78f65a..245ece0d 100644
--- a/custodian/vasp/tests/test_handlers.py
+++ b/custodian/vasp/tests/test_handlers.py
@@ -164,6 +164,57 @@ class VaspErrorHandlerTest(unittest.TestCase):
             [{"action": {"_set": {"ISMEAR": 0, "SIGMA": 0.05}}, "dict": "INCAR"}],
         )
 
+    def test_zbrent(self):
+        h = VaspErrorHandler("vasp.zbrent")
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["IBRION"], 1)
+
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["EDIFF"], 1e-6)
+        self.assertEqual(i["NELMIN"], 6)
+
+        shutil.copy("INCAR.orig", "INCAR")
+        h = VaspErrorHandler("vasp.zbrent")
+        h.vtst_fixes = True
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["IBRION"], 1)
+        self.assertEqual(i["EDIFF"], 0.0004)
+
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["EDIFF"], 1e-6)
+        self.assertEqual(i["NELMIN"], 6)
+        self.assertEqual(i["IBRION"], 3)
+        self.assertEqual(i["IOPT"], 7)
+        self.assertEqual(i["POTIM"], 0)
+
+        shutil.copy("INCAR.ediff", "INCAR")
+        h = VaspErrorHandler("vasp.zbrent")
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["IBRION"], 1)
+        self.assertEqual(i["EDIFF"], 1e-6)
+
+        h.check()
+        d = h.correct()
+        self.assertEqual(d["errors"], ["zbrent"])
+        i = Incar.from_file("INCAR")
+        self.assertEqual(i["EDIFF"], 1e-7)
+        self.assertEqual(i["NELMIN"], 6)
+
     def test_brmix(self):
         h = VaspErrorHandler("vasp.brmix")
         self.assertEqual(h.check(), True)
diff --git a/test_files/INCAR.ediff b/test_files/INCAR.ediff
new file mode 100644
index 00000000..28f50c60
--- /dev/null
+++ b/test_files/INCAR.ediff
@@ -0,0 +1,24 @@
+ALGO = Fast
+EDIFF = 1e-6
+ENCUT = 520
+IBRION = 2
+ICHARG = 1
+ISIF = 3
+ISMEAR = -5
+ISPIN = 2
+LDAU = True
+LDAUJ = 0 0 0 0
+LDAUL = 0 2 2 0
+LDAUTYPE = 2
+LDAUU = 0 4.0 6 0
+LMAXMIX = 4
+LORBIT = 11
+LREAL = Auto
+LWAVE = False
+MAGMOM = 2*0.6 1*5 5*0.6
+NELM = 100
+NELMIN = 6
+NPAR = 1
+NSW = 99
+PREC = Accurate
+SIGMA = 0.05
diff --git a/test_files/vasp.zbrent b/test_files/vasp.zbrent
new file mode 100644
index 00000000..cc2683b9
--- /dev/null
+++ b/test_files/vasp.zbrent
@@ -0,0 +1,42 @@
+       N       E                     dE             d eps       ncg     rms          rms(c)
+DAV:   1    -0.808978395638E+03   -0.10606E-05   -0.49230E-05  2080   0.160E-02    0.333E-02
+RMM:   2    -0.808978711288E+03   -0.31565E-03   -0.12893E-04  1630   0.152E-02    0.298E-01
+RMM:   3    -0.808978163978E+03    0.54731E-03   -0.10917E-04  1576   0.167E-02    0.806E-02
+RMM:   4    -0.808978368321E+03   -0.20434E-03   -0.33170E-05  1336   0.843E-03    0.444E-02
+RMM:   5    -0.808978372128E+03   -0.38079E-05   -0.72285E-06   961   0.454E-03
+  59 F= -.80897837E+03 E0= -.80897978E+03  d E =-.204339E-02  mag=    -0.3409
+ curvature:   0.00 expect dE=-0.247E-08 dE for cont linesearch -0.117E-08
+ ZBRENT: interpolating
+ opt :   0.6785  next Energy=  -808.979501 (dE=-0.317E-02)
+ bond charge predicted
+       N       E                     dE             d eps       ncg     rms          rms(c)
+DAV:   1    -0.808978358854E+03    0.94668E-05   -0.10023E-05  2032   0.697E-03    0.182E-02
+RMM:   2    -0.808978312836E+03    0.46018E-04   -0.97341E-06  1244   0.486E-03
+  60 F= -.80897831E+03 E0= -.80898022E+03  d E =-.198409E-02  mag=    -0.3378
+ curvature:   0.00 expect dE=-0.405E-09 dE for cont linesearch -0.143E-09
+ ZBRENT: interpolating
+ opt :   0.6785  next Energy=  -808.979501 (dE=-0.317E-02)
+ bond charge predicted
+       N       E                     dE             d eps       ncg     rms          rms(c)
+DAV:   1    -0.808978272998E+03    0.85855E-04   -0.12172E-05  2032   0.996E-03    0.246E-02
+RMM:   2    -0.808978269206E+03    0.37918E-05   -0.29084E-05  1345   0.590E-03
+  61 F= -.80897827E+03 E0= -.80898068E+03  d E =-.194046E-02  mag=    -0.3340
+ curvature:   0.00 expect dE=-0.338E-09 dE for cont linesearch -0.426E-10
+ ZBRENT:  can not reach accuracy
+ opt :   0.0000  next Energy=     0.000000 (dE= 0.809E+03)
+ bond charge predicted
+       N       E                     dE             d eps       ncg     rms          rms(c)
+DAV:   1    -0.808977102563E+03    0.11704E-02   -0.10614E-01  3392   0.465E-01    0.338E-01
+RMM:   2    -0.809064633783E+03   -0.87531E-01   -0.17818E-02  3195   0.104E-01    0.746E+00
+RMM:   3    -0.808928314960E+03    0.13632E+00   -0.17136E-02  2844   0.122E-01    0.621E-01
+RMM:   4    -0.808973971408E+03   -0.45656E-01   -0.36113E-03  2723   0.868E-02    0.422E-01
+RMM:   5    -0.808976894323E+03   -0.29229E-02   -0.60123E-04  2128   0.386E-02    0.241E-01
+RMM:   6    -0.808976679739E+03    0.21458E-03   -0.11445E-04  1662   0.160E-02    0.243E-01
+RMM:   7    -0.808976527328E+03    0.15241E-03   -0.10336E-04  1656   0.154E-02    0.236E-01
+RMM:   8    -0.808976418689E+03    0.10864E-03   -0.86241E-05  1593   0.155E-02    0.179E-01
+RMM:   9    -0.808976445212E+03   -0.26523E-04   -0.54527E-05  1600   0.139E-02
+  62 F= -.80897645E+03 E0= -.80897894E+03  d E =-.116469E-03  mag=    -0.3177
+ curvature:   0.00 expect dE= 0.000E+00 dE for cont linesearch  0.000E+00
+ ZBRENT: fatal error in bracketing
+     please rerun with smaller EDIFF, or copy CONTCAR
+     to POSCAR and continue
\ No newline at end of file