From 7dc64294a0c4b3ab920db7c4521a7742d0409fc2 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Wed, 15 May 2024 22:45:24 -0400 Subject: [PATCH] ignore the SCC error in dftb --- pyxtal/db.py | 33 ++++++++++----- pyxtal/interface/dftb.py | 89 +++++++++++++++++++++++++--------------- 2 files changed, 79 insertions(+), 43 deletions(-) diff --git a/pyxtal/db.py b/pyxtal/db.py index 9e34ef07..39f7635e 100644 --- a/pyxtal/db.py +++ b/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,7 +84,6 @@ 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) @@ -79,7 +91,8 @@ def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.04) 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 diff --git a/pyxtal/interface/dftb.py b/pyxtal/interface/dftb.py index dbd3fc5c..3e3ce06a 100644 --- a/pyxtal/interface/dftb.py +++ b/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)."""