Skip to content

Commit

Permalink
adjusted optimization script for new setup
Browse files Browse the repository at this point in the history
  • Loading branch information
@awvwgk
awvwgk committed Oct 9, 2019
1 parent 6270484 commit 650f412
Showing 1 changed file with 173 additions and 85 deletions.
258 changes: 173 additions & 85 deletions scripts/xtb_opt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,73 +16,145 @@
# You should have received a copy of the GNU Lesser General Public License
# along with xtb. If not, see <https://www.gnu.org/licenses/>.

"""Main wrapper script to perform geometry optimizations with xtb."""
"""Wrapper script to perform geometry optimizations with xtb Python wrapper."""

if __name__ == "__main__":
import argparse

def get_cmd_args():
"""Commandline arguments for xtb wrapper in script mode"""
parser = argparse.ArgumentParser(prog='xtb',
description='Wrapper for xTB calculation.')

parser.add_argument("-f", "--format", dest="filetype", action="store",
default=None, help="Format of input geometry (automatic)")

parser.add_argument(dest="filename", action="store",
help="Input geometry")

parser.add_argument("-x", "--method", dest="method", action="store",
default='gfn2', help="SQM method for calculation (gfn2)")

parser.add_argument("-c", "--charge", dest="charge", action="store",
type=int, default=0, help="Total charge (0)")

parser.add_argument("--etemp", dest="etemp", action="store",
type=float, default=300.0,
help="Electronic temperature (300K)")

parser.add_argument("--accuracy", dest="acc", action="store",
type=float, default=1.0,
help="Calculation accuracy (1.0)")

parser.add_argument("--maxiter", dest="maxiter", action="store",
type=int, default=250,
help="Maximum number of SCC iterations (250)")

parser.add_argument("--lcovb", dest="rnn", action="store",
type=float, default=5.0,
help="Longest covalent bond for preconditioning (5.0)")

parser.add_argument("--econv", dest="econv", action="store",
type=float, default=5.4423e-4,
help="Convergence threshold for energy")

parser.add_argument("--gconv", dest="fconv", action="store",
type=float, default=3.88938e-05,
help="Convergence threshold for forces")

parser.add_argument("--linesearch", dest="armijo", action="store_true",
help="Use line search")

parser.add_argument("--precon", dest="precon", action="store_true",
help="Use preconditioning")

parser.add_argument("--optcell", dest="optcell", action="store_true",
help="Relax cellparameters")

parser.add_argument("--logfile", dest="logfile", action="store",
default='-', help="File for logging information (STDOUT)")

parser.add_argument("--trajectory", dest="trajectory", action="store",
default='xtbopt.traj',
help="File for trajectory output (xtbopt.traj)")

parser.add_argument("--solvent", dest="solvent", action="store",
default="none", help="Solvent for GBSA (none)")

parser.add_argument("--library", dest="library", action="store",
default=None, help="Path to libxtb.so (None)")
parser = argparse.ArgumentParser(
prog="xtb", description="Wrapper for xTB calculation."
)

parser.add_argument(
"-f",
"--format",
dest="filetype",
action="store",
default=None,
help="Format of input geometry (automatic)",
)

parser.add_argument(dest="filename", action="store", help="Input geometry")

parser.add_argument(
"-x",
"--method",
dest="method",
action="store",
default="gfn2",
help="SQM method for calculation (gfn2)",
)

parser.add_argument(
"-c",
"--charge",
dest="charge",
action="store",
type=int,
default=0,
help="Total charge (0)",
)

parser.add_argument(
"--etemp",
dest="etemp",
action="store",
type=float,
default=300.0,
help="Electronic temperature (300K)",
)

parser.add_argument(
"--accuracy",
dest="acc",
action="store",
type=float,
default=1.0,
help="Calculation accuracy (1.0)",
)

parser.add_argument(
"--maxiter",
dest="maxiter",
action="store",
type=int,
default=250,
help="Maximum number of SCC iterations (250)",
)

parser.add_argument(
"--lcovb",
dest="rnn",
action="store",
type=float,
default=5.0,
help="Longest covalent bond for preconditioning (5.0)",
)

parser.add_argument(
"--econv",
dest="econv",
action="store",
type=float,
default=5.4423e-4,
help="Convergence threshold for energy",
)

parser.add_argument(
"--gconv",
dest="fconv",
action="store",
type=float,
default=3.88938e-05,
help="Convergence threshold for forces",
)

parser.add_argument(
"--linesearch",
dest="armijo",
action="store_true",
help="Use line search"
)

parser.add_argument(
"--precon",
dest="precon",
action="store_true",
help="Use preconditioning"
)

parser.add_argument(
"--optcell",
dest="optcell",
action="store_true",
help="Relax cellparameters",
)

parser.add_argument(
"--logfile",
dest="logfile",
action="store",
default="-",
help="File for logging information (STDOUT)",
)

parser.add_argument(
"--trajectory",
dest="trajectory",
action="store",
default="xtbopt.traj",
help="File for trajectory output (xtbopt.traj)",
)

parser.add_argument(
"--solvent",
dest="solvent",
action="store",
default="none",
help="Solvent for GBSA (none)",
)

return parser.parse_args()

Expand All @@ -91,11 +163,12 @@ def get_cmd_args():
from ase.optimize.precon import Exp, PreconFIRE
from ase.constraints import ExpCellFilter
from ase.units import Hartree
from xtb.calculators import GFN1, GFN2, GFN0
from xtb import GFN1, GFN2, GFN0

# overwrite convergence thresholds of PreconFIRE optimizer
class PatchedOptimizer(PreconFIRE):
"""Patches the ASE optimizer to use a different convergence threshold"""

econv = None
fconv = None
elast = None
Expand All @@ -104,7 +177,9 @@ def initialize(self):
PreconFIRE.initialize(self)
self.elast = None

def run(self, steps=100000000, econv=0.00054423, fconv=3.88938e-05): # pylint: disable=arguments-differ
def run(
self, steps=100000000, econv=0.00054423, fconv=3.88938e-05
): # pylint: disable=arguments-differ
self.econv = econv
self.fconv = fconv
PreconFIRE.run(self, 0.05, steps)
Expand All @@ -124,11 +199,11 @@ def converged(self, forces=None):
# check for convergence in forces
if forces is None:
forces = self.atoms.get_forces()
fnorm = sqrt((forces**2).sum())
fnorm = sqrt((forces ** 2).sum())
fconverged = fnorm < self.fconv

#print("--> norm(F):", fnorm, "converged?", fconverged)
#print("--> energy :", ecurr, "converged?", econverged)
# print("--> norm(F):", fnorm, "converged?", fconverged)
# print("--> energy :", ecurr, "converged?", econverged)

return econverged and fconverged

Expand All @@ -140,18 +215,21 @@ def converged(self, forces=None):
else:
mol = read(args.filename)

parameters = {'accuracy': args.acc,
'temperature': args.etemp,
'max_iterations': args.maxiter,
'solvent': args.solvent,
'libxtb_path': args.library,
'print_level': 2}
parameters = {
"accuracy": args.acc,
"electronic_temperature": args.etemp,
"max_iterations": args.maxiter,
"solvent": args.solvent,
"print_level": 2,
}

if args.method == 'gfn0':
if args.method == "gfn0":
calc = GFN0(**parameters)
elif args.method == 'gfn1':
elif args.method == "gfn1":
if mol.pbc.any():
raise Exception("GFN1-xTB is not available with PBC")
calc = GFN1(**parameters)
elif args.method == 'gfn2':
elif args.method == "gfn2":
if mol.pbc.any():
raise Exception("GFN2-xTB is not available with PBC")
calc = GFN2(**parameters)
Expand All @@ -161,7 +239,7 @@ def converged(self, forces=None):

e = mol.get_potential_energy()

print('Initial energy: eV, Eh', e, e/Hartree)
print("Initial energy: eV, Eh", e, e / Hartree)

if args.precon:
precon = Exp(A=3, r_NN=args.rnn + 0.1, r_cut=args.rnn + 0.6)
Expand All @@ -170,23 +248,33 @@ def converged(self, forces=None):

if args.optcell:
sf = ExpCellFilter(mol)
relax = PatchedOptimizer(sf, precon=precon, trajectory=args.trajectory,
logfile=args.logfile, use_armijo=args.armijo)
relax = PatchedOptimizer(
sf,
precon=precon,
trajectory=args.trajectory,
logfile=args.logfile,
use_armijo=args.armijo,
)
else:
relax = PatchedOptimizer(mol, precon=precon, trajectory=args.trajectory,
logfile=args.logfile, use_armijo=args.armijo)
relax = PatchedOptimizer(
mol,
precon=precon,
trajectory=args.trajectory,
logfile=args.logfile,
use_armijo=args.armijo,
)

try:
relax.run(econv=args.econv, fconv=args.fconv)
except KeyboardInterrupt:
print('User got impatient')
print("User got impatient")
except RuntimeError:
print('Optimization terminated due to internal error')
print("Optimization terminated due to internal error")

e = mol.get_potential_energy()
print('Final energy: eV, Eh', e, e/Hartree)
print("Final energy: eV, Eh", e, e / Hartree)

if mol.pbc.any():
write("xtbopt.POSCAR", mol, format='vasp')
write("xtbopt.POSCAR", mol, format="vasp")
else:
write("xtbopt.xyz", mol, format='xyz')
write("xtbopt.xyz", mol, format="xyz")

0 comments on commit 650f412

Please sign in to comment.