ignore the SCC error in dftb

pyxtal/db.py

@@ -8,6 +8,7 @@
 from pyxtal import pyxtal
 import pymatgen.analysis.structure_matcher as sm
 from pyxtal.util import ase2pymatgen
+from ase.calculators.calculator import CalculationFailed
 
 def dftb_opt_par(ids, xtals, skf_dir, steps, folder, symmetrize, criteria):
     """
@@ -29,7 +30,7 @@ def dftb_opt_par(ids, xtals, skf_dir, steps, folder, symmetrize, criteria):
     os.chdir(cwd)
     return results
 
-def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.04):
+def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.05):
     """
     Single DFTB optimization for a given atomic xtal
 
@@ -44,21 +45,33 @@ def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.04)
     cwd = os.getcwd()
     atoms = xtal.to_ase(resort=False)
     eng = None
+    stress = None
     try:
         if symmetrize:
-            s = DFTB_relax(atoms, skf_dir, True, 250,
+            s = DFTB_relax(atoms, skf_dir, True, int(steps/2),
                            kresol=kresol*1.2, folder='.',
-                           scc_error = 0.1,
+                           scc_iter = 100,
                            logfile='ase.log')
-            s = DFTB_relax(atoms, skf_dir, True, steps,
+            s = DFTB_relax(atoms, skf_dir, True, int(steps/2),
                            kresol=kresol, folder='.',
                            scc_error = 1e-5,
+                           scc_iter = 100,
                            logfile='ase.log')
+            stress = np.sum(s.get_stress()[:3])/0.006241509125883258/3
         else:
-            my = DFTB(atoms, skf_dir, kresol=kresol*1.2, folder='.', scc_error=0.1)
-            s, eng = my.run(mode='vc-relax', step=250)
-            my = DFTB(s, skf_dir, kresol=kresol, folder='.', scc_error=1e-4)
-            s, eng = my.run(mode='vc-relax', step=steps)
+            my = DFTB(atoms, skf_dir, kresol=kresol*1.5, folder='.',
+                      scc_error=0.1, scc_iter=100)
+            s, eng = my.run(mode='vc-relax', step=int(steps/2))
+            my = DFTB(s, skf_dir, kresol=kresol, folder='.',
+                      scc_error=1e-4, scc_iter=100)
+            s, eng = my.run(mode='vc-relax', step=int(steps/2))
+            s = my.struc
+    except CalculationFailed:
+        # This is due to covergence error in geometry optimization
+        # Here we simply read the last energy
+        my.calc.read_results()
+        eng = my.calc.results['energy']
+        s = my.struc
     except:
         s = None
         print("Problem in DFTB Geometry optimization", id)
@@ -71,15 +84,15 @@ def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.04)
             eng = s.get_potential_energy()/len(s)
         else:
             eng /= len(s)
-        stress = np.sum(s.get_stress()[:3])/0.006241509125883258/3
 
         if criteria is not None:
             status = xtal.check_validity(criteria)
         else:
             status = True
 
         header = "{:4d}".format(id)
-        dicts = {'validity': status, 'energy': eng, 'stress': stress}
+        dicts = {'validity': status, 'energy': eng}
+        if stress is not None: dicts['stress'] = stress
         print(xtal.get_xtal_string(header=header, dicts=dicts))
 
         return xtal, eng, status

pyxtal/interface/dftb.py

@@ -1,4 +1,4 @@
-import os
+import os, re
 from ase.io import read
 from ase.optimize.fire import FIRE
 from ase.constraints import ExpCellFilter
@@ -10,7 +10,8 @@
 import numpy as np
 
 
-def make_Hamiltonian(skf_dir, atom_types, disp, kpts, scc_error=1e-06, write_band=False, use_omp=False):
+def make_Hamiltonian(skf_dir, atom_types, disp, kpts, scc_error=1e-06,
+                     scc_iter=500, write_band=False, use_omp=False):
     """
     Generate the DFTB Hamiltonian for DFTB+
     """
@@ -64,7 +65,7 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, scc_error=1e-06, write_ban
 
     kwargs = {'Hamiltonian_SCC': 'yes',
               'Hamiltonian_SCCTolerance': scc_error, #1e-06,
-              'Hamiltonian_MaxSCCIterations': 1000,
+              'Hamiltonian_MaxSCCIterations': scc_iter, #1000,
               #'Hamiltonian_Mixer': 'DIIS{}', #Default is Broyden
               #'Hamiltonian_Dispersion': dispersion,
               'slako_dir': skf_dir,
@@ -159,7 +160,7 @@ def make_Hamiltonian(skf_dir, atom_types, disp, kpts, scc_error=1e-06, write_ban
 def DFTB_relax(struc, skf_dir, opt_cell=False, step=500, \
                fmax=0.1, kresol=0.10, folder='tmp', disp='D3', \
                mask=None, symmetrize=True, logfile=None, \
-               scc_error=1e-6, use_omp=False):
+               scc_error=1e-6, scc_iter=500, use_omp=False):
     """
     DFTB optimizer based on ASE
 
@@ -269,6 +270,7 @@ def __init__(self, struc, skf_dir,
                  prefix = 'geo_final',
                  use_omp = False,
                  scc_error = 1e-6,
+                 scc_iter = 500,
                 ):
 
         self.struc = struc
@@ -281,6 +283,7 @@ def __init__(self, struc, skf_dir,
         self.prefix = prefix
         self.use_omp = use_omp
         self.scc_error = scc_error
+        self.scc_iter = scc_iter
         if not os.path.exists(self.folder):
            os.makedirs(self.folder)
 
@@ -302,6 +305,7 @@ def get_calculator(self, mode, step=500, ftol=1e-3, FixAngles=False, eVperA=True
         atom_types = set(self.struc.get_chemical_symbols())
         kwargs = make_Hamiltonian(self.skf_dir, atom_types, self.disp, self.kpts,
                                   scc_error=self.scc_error,
+                                  scc_iter=self.scc_iter,
                                   use_omp=self.use_omp)
 
         if mode in ['relax', 'vc-relax']:
@@ -369,11 +373,11 @@ def run(self, mode, step=500, ftol=1e-3, FixAngles=False, md_params={}):
         cwd = os.getcwd()
         os.chdir(self.folder)
 
-        calc = self.get_calculator(mode, step, ftol, FixAngles, md_params=md_params)
-        self.struc.set_calculator(calc)
+        self.calc = self.get_calculator(mode, step, ftol, FixAngles, md_params=md_params)
+        self.struc.set_calculator(self.calc)
         # self.struc.write('geo_o.gen', format='dftb')
         # execute the simulation
-        calc.calculate(self.struc)
+        self.calc.calculate(self.struc)
         if mode in ['relax', 'vc-relax']:
             final = read(self.prefix+'.gen')
         else:
@@ -382,7 +386,7 @@ def run(self, mode, step=500, ftol=1e-3, FixAngles=False, md_params={}):
         # get the final energy
         energy = self.struc.get_potential_energy()
 
-        with open(self.label+'.out') as f:
+        with open(self.label + '.out') as f:
             l = f.readlines()
             self.version = l[2]
         os.chdir(cwd)
@@ -626,34 +630,53 @@ def read_results(self):
         """ all results are read from results.tag file
             It will be destroyed after it is read to avoid
             reading it once again after some runtime error """
-
         with open(os.path.join(self.directory, 'results.tag'), 'r') as fd:
             self.lines = fd.readlines()
+        if len(self.lines) == 0:
+            #print("READ RESULTS from test.out")
+            self.results['energy'] = self.read_energy()
+            self.results['forces'] = None
+            self.results['stress'] = None
+        else:
+            self.atoms = self.atoms_input
+            self.results['energy'] = float(self.lines[1])*Hartree
+            forces = self.read_forces()
+            self.results['forces'] = forces
+
+            # stress stuff begins
+            sstring = 'stress'
+            have_stress = False
+            stress = list()
+            for iline, line in enumerate(self.lines):
+                if sstring in line:
+                    have_stress = True
+                    start = iline + 1
+                    end = start + 3
+                    for i in range(start, end):
+                        cell = [float(x) for x in self.lines[i].split()]
+                        stress.append(cell)
+            if have_stress:
+                stress = -np.array(stress) * Hartree / Bohr**3
+                self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
+            # stress stuff ends
+
+            # calculation was carried out with atoms written in write_input
+            os.remove(os.path.join(self.directory, 'results.tag'))
+
+    def read_energy(self):
+        """
+        If SCC is not converged, read the last step energy from test.out
+        """
+        outfile = self.label + '.out'
+        energies = []
+        with open(os.path.join(self.directory, outfile), 'r') as fd:
+            lines = fd.readlines()
+        for line in lines:
+            m = re.match(r'Total Energy:\s+[-\d.]+ H\s+([-.\d]+) eV', line)
+            if m:
+                energies.append(float(m.group(1)))
+        return energies[-1]
 
-        self.atoms = self.atoms_input
-        self.results['energy'] = float(self.lines[1])*Hartree
-        forces = self.read_forces()
-        self.results['forces'] = forces
-
-        # stress stuff begins
-        sstring = 'stress'
-        have_stress = False
-        stress = list()
-        for iline, line in enumerate(self.lines):
-            if sstring in line:
-                have_stress = True
-                start = iline + 1
-                end = start + 3
-                for i in range(start, end):
-                    cell = [float(x) for x in self.lines[i].split()]
-                    stress.append(cell)
-        if have_stress:
-            stress = -np.array(stress) * Hartree / Bohr**3
-            self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
-        # stress stuff ends
-
-        # calculation was carried out with atoms written in write_input
-        os.remove(os.path.join(self.directory, 'results.tag'))
 
     def read_forces(self):
         """Read Forces from dftb output file (results.tag)."""