{"id":806,"date":"2023-04-18T09:17:17","date_gmt":"2023-04-18T09:17:17","guid":{"rendered":"https:\/\/doublelayer.eu\/vilab\/?p=806"},"modified":"2023-04-18T09:19:52","modified_gmt":"2023-04-18T09:19:52","slug":"dft-geometry-optimizers","status":"publish","type":"post","link":"https:\/\/doublelayer.eu\/vilab\/2023\/04\/18\/dft-geometry-optimizers\/","title":{"rendered":"DFT geometry optimizers"},"content":{"rendered":"\n

These are undeservedly less attention to optimizers than density functionals (concerning Jacob’s ladder). It is not even covered in the recent review: Best-Practice DFT Protocols for Basic Molecular Computational Chemistry<\/a>. At the same time, in my current projects, the most resource-demanding was geometry optimization \u2013 the time spent on optimizing structures was much longer than a single-point calculation. Papers that introduce new (AI-based) optimizers promise significant speed-up. However, there are always some problems:<\/p>\n\n\n\n

    \n
  1. The tested systems are different from my electrochemical interfaces.<\/li>\n\n\n\n
  2. The code is not available or difficult to install.<\/li>\n\n\n\n
  3. The code is outdated and contains bugs.<\/li>\n\n\n\n
  4. Optimizers perform worse than the common ones, like QuasiNewton in ASE.<\/li>\n<\/ol>\n\n\n\n

    ASE wiki lists all internal and some external optimizers<\/a> and provides their comparison<\/a>. I have checked the most promising on a high-entropy alloy slab.<\/p>\n\n\n\n

    \"\"<\/figure>\n\n\n\n

    Observation 1.<\/strong> QuasiNewton outperforms all other optimizers. Point. I have run a standard GPAW\/DFT\/PBE\/PW optimization with various optimizers:<\/p>\n\n\n\n

    Observation 2.<\/strong> Pre-optimizing the slab with a cheaper method does not reduce the number of optimization steps. I have preoptimized the geometry with TBLITE\/DFTB\/GFN1-xTB to continue with GPAW\/DFT\/PBE\/PW. Preoptimization takes just some minutes and the obtained geometry looks similar to the DFT one but that does not reduce the number of DFT optimization steps. <\/p>\n\n\n\n

    Optimizer<\/td>N steps*<\/td>Time$<\/td>N steps*#<\/td>Total time#<\/td><\/tr>
    BFGS<\/td>16<\/td>02:44:27<\/td>17<\/td>03:01:26<\/td><\/tr>
    LBFGS<\/td>15<\/td>02:30:35<\/td>16<\/td>02:55:04<\/td><\/tr>
    BondMin<\/td>12<\/td>02:46:27<\/td>13<\/td>02:45:07<\/td><\/tr>
    GPMin<\/td>12<\/td>05:26:23<\/td>31<\/td>08:14:22<\/td><\/tr>
    MLMin<\/td>38<\/td>verylong<\/td>28<\/td>12:31:29<\/td><\/tr>
    FIRE<\/td>38<\/td>05:06:56<\/td>44<\/td>05:56:54<\/td><\/tr>
    QuasiNewton<\/td>8<\/td>01:36:23<\/td>9<\/td>02:00:10<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

    Note * \u2013 the printed number of steps might different from the actuall number of calculations because each calculator has a different way of reporting that number.<\/p>\n\n\n\n

    Note $ \u2013 the time between the end of the first and last steps. <\/p>\n\n\n\n

    Note # \u2013 started from the TBLITE\/DFTB\/GFN1-xTB preoptimized geometry.<\/p>\n\n\n\n

    N.B! I have done my test only once in two runs: starting with slab.xyz and preoptized geometry. Runs were on similar nodes and all optimizations were done on the same node.<\/p>\n\n\n\n

    Conclusion.<\/strong> Do not believe in claims in articles advertizing new optimizers \u2013 Run your tests before using them.<\/p>\n\n\n\n

    A practical finding.<\/strong> The usual problem with calculations that require many optimization steps is that they need to fit into HPC time limits. On the restart, ASE usually rewrites the trajectory. Some optimizers (GPMin and AI-based) could benefit from reading the full trajectory. So, I started writing two trajectories and a restart file like this.<\/p>\n\n\n\n

    # Restarting\nif os.path.exists(f'{name}_last.gpw') == True and os.stat(f'{name}_last.gpw').st_size > 0:\n    atoms,calc = restart(f'{name}_last.gpw', txt=None)\n    parprint(f'Restart from the gpw geometry.')\nelif os.path.exists(f'{name}_full.traj') == True and os.stat(f'{name}_full.traj').st_size > 0:\n    atoms = read(f'{name}_full.traj',-1)\n    parprint(f'Restart with the traj geometry.')\nelse:\n    atoms = read(f'{name}_init.xyz')\n    parprint(f'Start with the initial xyz geometry.')\n\n# Optimizing\nopt = QuasiNewton(atoms, trajectory=f'{name}.traj', logfile=f'{name}.log')\ntraj= Trajectory(f'{name}_full.traj', 'a', atoms)\nopt.attach(traj.write, interval=1)\ndef writegpw():\n    calc.write(f'{name}_last.gpw')\nopt.attach(writegpw, interval=1)\nopt.run(fmax=0.05, steps=42)<\/code><\/pre>\n\n\n\n
    Here are some details on the tests.<\/p>\n\n\n\n
    My gpaw_opt.py for DFT calculations on 24 cores:<\/p>\n\n\n\n
    # Load modules\nfrom ase import Atom, Atoms\nfrom ase.build import add_adsorbate, fcc100, fcc110, fcc111, fcc211, molecule\nfrom ase.calculators.mixing import SumCalculator\nfrom ase.constraints import FixAtoms, FixedPlane, FixInternals\nfrom ase.data.vdw_alvarez import vdw_radii\nfrom ase.db import connect\nfrom ase.io import write, read\nfrom ase.optimize import BFGS, GPMin, LBFGS, FIRE, QuasiNewton\nfrom ase.parallel import parprint\nfrom ase.units import Bohr\nfrom bondmin import BondMin\nfrom catlearn.optimize.mlmin import MLMin\nfrom dftd4.ase import DFTD4\nfrom gpaw import GPAW, PW, FermiDirac, PoissonSolver, Mixer, restart\nfrom gpaw.dipole_correction import DipoleCorrection\nfrom gpaw.external import ConstantElectricField\nfrom gpaw.utilities import h2gpts\nimport numpy as np\nimport os\n\natoms = read('slab.xyz')\natoms.set_constraint([FixAtoms(indices=[atom.index for atom in atoms if atom.tag in [1,2]])])\n\n# Set calculator\nkwargs = dict(poissonsolver={'dipolelayer':'xy'},\n              xc='RPBE',\n              kpts=(4,4,1),\n              gpts=h2gpts(0.18, atoms.get_cell(), idiv=4),\n              mode=PW(400),\n              basis='dzp',\n              parallel={'augment_grids':True,'sl_auto':True,'use_elpa':True},\n             )\ncalc = GPAW(**kwargs)\n\n#atoms.calc = SumCalculator([DFTD4(method='RPBE'), calc])\n#atoms.calc = calc\n\n# Optimization paramters\nmaxf = 0.05\n\n# Run optimization\n###############################################################################\n\n# 2.A. Optimize structure using MLMin (CatLearn).\ninitial_mlmin = atoms.copy()\ninitial_mlmin.set_calculator(calc)\nmlmin_opt = MLMin(initial_mlmin, trajectory='results_mlmin.traj')\nmlmin_opt.run(fmax=maxf, kernel='SQE', full_output=True)\n\n# 2.B Optimize using GPMin.\ninitial_gpmin = atoms.copy()\ninitial_gpmin.set_calculator(calc)\ngpmin_opt = GPMin(initial_gpmin, trajectory='results_gpmin.traj', logfile='results_gpmin.log', update_hyperparams=True)\ngpmin_opt.run(fmax=maxf)\n\n# 2.C Optimize using LBFGS.\ninitial_lbfgs = atoms.copy()\ninitial_lbfgs.set_calculator(calc)\nlbfgs_opt = LBFGS(initial_lbfgs, trajectory='results_lbfgs.traj', logfile='results_lbfgs.log')\nlbfgs_opt.run(fmax=maxf)\n\n# 2.D Optimize using FIRE.\ninitial_fire = atoms.copy()\ninitial_fire.set_calculator(calc)\nfire_opt = FIRE(initial_fire, trajectory='results_fire.traj', logfile='results_fire.log')\nfire_opt.run(fmax=maxf)\n\n# 2.E Optimize using QuasiNewton.\ninitial_qn = atoms.copy()\ninitial_qn.set_calculator(calc)\nqn_opt = QuasiNewton(initial_qn, trajectory='results_qn.traj', logfile='results_qn.log')\nqn_opt.run(fmax=maxf)\n\n# 2.F Optimize using BFGS.\ninitial_bfgs = atoms.copy()\ninitial_bfgs.set_calculator(calc)\nbfgs_opt = LBFGS(initial_bfgs, trajectory='results_bfgs.traj', logfile='results_bfgs.log')\nbfgs_opt.run(fmax=maxf)\n\n# 2.G. Optimize structure using BondMin.\ninitial_bondmin = atoms.copy()\ninitial_bondmin.set_calculator(calc)\nbondmin_opt = BondMin(initial_bondmin, trajectory='results_bondmin.traj',logfile='results_bondmin.log')\nbondmin_opt.run(fmax=maxf)\n\n# Summary of the results\n###############################################################################\n\nfire_results = read('results_fire.traj', ':')\nparprint('Number of function evaluations using FIRE:',\n         len(fire_results))\n\nlbfgs_results = read('results_lbfgs.traj', ':')\nparprint('Number of function evaluations using LBFGS:',\n         len(lbfgs_results))\n\ngpmin_results = read('results_gpmin.traj', ':')\nparprint('Number of function evaluations using GPMin:',\n         gpmin_opt.function_calls)\n\nbfgs_results = read('results_bfgs.traj', ':')\nparprint('Number of function evaluations using BFGS:',\n         len(bfgs_results))\n\nqn_results = read('results_qn.traj', ':')\nparprint('Number of function evaluations using QN:',\n         len(qn_results))\n\ncatlearn_results = read('results_mlmin.traj', ':')\nparprint('Number of function evaluations using MLMin:',\n         len(catlearn_results))\n\nbondmin_results = read('results_bondmin.traj', ':')\nparprint('Number of function evaluations using BondMin:',\n         len(bondmin_results))<\/code><\/pre>\n\n\n\n
    Initial slab.xyz file: <\/p>\n\n\n\n
    45\nLattice=\"8.529357696932532 0.0 0.0 4.264678848466266 7.386640443507905 0.0 0.0 0.0 29.190908217261956\" Properties=species:S:1:pos:R:3:tags:I:1 pbc=\"T T F\"\nIr       0.00000000       1.62473838      10.00000000        5\nRu       2.81412943       1.62473838      10.00000000        5\nPt       5.62825885       1.62473838      10.00000000        5\nPd       1.40706471       4.06184595      10.00000000        5\nAg       4.22119414       4.06184595      10.00000000        5\nAg       7.03532356       4.06184595      10.00000000        5\nAg       2.81412943       6.49895353      10.00000000        5\nRu       5.62825885       6.49895353      10.00000000        5\nPt       8.44238828       6.49895353      10.00000000        5\nPt       0.00000000       0.00000000      12.29772705        4\nAg       2.81412943       0.00000000      12.29772705        4\nRu       5.62825885       0.00000000      12.29772705        4\nRu       1.40706471       2.43710757      12.29772705        4\nIr       4.22119414       2.43710757      12.29772705        4\nAg       7.03532356       2.43710757      12.29772705        4\nAg       2.81412943       4.87421514      12.29772705        4\nIr       5.62825885       4.87421514      12.29772705        4\nPd       8.44238828       4.87421514      12.29772705        4\nPd       1.40706471       0.81236919      14.59545411        3\nIr       4.22119414       0.81236919      14.59545411        3\nPt       7.03532356       0.81236919      14.59545411        3\nAg       2.81412943       3.24947676      14.59545411        3\nIr       5.62825885       3.24947676      14.59545411        3\nIr       8.44238828       3.24947676      14.59545411        3\nPd       4.22119414       5.68658433      14.59545411        3\nPt       7.03532356       5.68658433      14.59545411        3\nAg       9.84945299       5.68658433      14.59545411        3\nPd       0.00000000       1.62473838      16.89318116        2\nPd       2.81412943       1.62473838      16.89318116        2\nAg       5.62825885       1.62473838      16.89318116        2\nPt       1.40706471       4.06184595      16.89318116        2\nAg       4.22119414       4.06184595      16.89318116        2\nAg       7.03532356       4.06184595      16.89318116        2\nRu       2.81412943       6.49895353      16.89318116        2\nRu       5.62825885       6.49895353      16.89318116        2\nRu       8.44238828       6.49895353      16.89318116        2\nIr       0.00000000       0.00000000      19.19090822        1\nAg       2.81412943       0.00000000      19.19090822        1\nPt       5.62825885       0.00000000      19.19090822        1\nPd       1.40706471       2.43710757      19.19090822        1\nAg       4.22119414       2.43710757      19.19090822        1\nPd       7.03532356       2.43710757      19.19090822        1\nAg       2.81412943       4.87421514      19.19090822        1\nRu       5.62825885       4.87421514      19.19090822        1\nIr       8.44238828       4.87421514      19.19090822        1<\/code><\/pre>\n\n\n\n
    My tblite_opt.py for DFTB calcualation with just one core. It takes some minutes but eventually crashes \ud83d\ude41<\/p>\n\n\n\n
    # Load modules\nfrom ase import Atom, Atoms\nfrom ase.build import add_adsorbate, fcc100, fcc110, fcc111, fcc211, molecule\nfrom ase.calculators.mixing import SumCalculator\nfrom ase.constraints import FixAtoms, FixedPlane, FixInternals\nfrom ase.data.vdw_alvarez import vdw_radii\nfrom ase.db import connect\nfrom ase.io import write, read\nfrom ase.optimize import BFGS, GPMin, LBFGS, FIRE, QuasiNewton\nfrom ase.parallel import parprint\nfrom ase.units import Bohr\nfrom tblite.ase import TBLite\nimport numpy as np\nimport os\n\n# https:\/\/tblite.readthedocs.io\/en\/latest\/users\/ase.html\n\natoms = read('slab.xyz')\natoms.set_constraint([FixAtoms(indices=[atom.index for atom in atoms if atom.tag in [1,2]])])\n\n# Set calculator\ncalc = TBLite(method=\"GFN1-xTB\",accuracy=1000,electronic_temperature=300,max_iterations=300)\natoms.set_calculator(calc)\nqn_opt = QuasiNewton(atoms, trajectory='results_qn.traj', logfile='results_qn.log', maxstep=0.1)\nqn_opt.run(fmax=0.1)<\/code><\/pre>\n\n\n\n
    To compare structures I have used MDanalysis, which unfortunately does not work with ASE traj, so I prepared xyz-files with “ase convert -n -1 file.traj file.xyz”<\/p>\n\n\n\n
    import MDAnalysis as mda\nfrom MDAnalysis.analysis.rms import rmsd\nimport sys\n\ndef coord(file_name):\n    file  = mda.Universe(f\"{file_name}.xyz\")\n    atoms = file.select_atoms(\"index 1:9\")\n    return  atoms.positions.copy()\n\nprint(rmsd(coord(sys.argv[1]),coord(sys.argv[2])))<\/code><\/pre>\n\n\n\n
    An instruction on installation of GPAW<\/a>. TBLITE<\/a> can be installed as “conda install -c conda-forge tblite”.<\/p>\n","protected":false},"excerpt":{"rendered":"
    These are undeservedly less attention to optimizers than density functionals (concerning Jacob’s ladder). It is not even covered in the recent review: Best-Practice DFT Protocols for Basic Molecular Computational Chemistry. At the same time, in my current projects, the most resource-demanding was geometry optimization \u2013 the time spent on optimizing structures was much longer than a… Read More<\/a><\/p>\n","protected":false},"author":5,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[26,12,40],"tags":[48,53,31,34],"class_list":["post-806","post","type-post","status-publish","format-standard","hentry","category-collaborations","category-know-how","category-opensource","tag-ase","tag-gpaw","tag-hpc","tag-software"],"_links":{"self":[{"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/posts\/806","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/users\/5"}],"replies":[{"embeddable":true,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/comments?post=806"}],"version-history":[{"count":2,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/posts\/806\/revisions"}],"predecessor-version":[{"id":812,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/posts\/806\/revisions\/812"}],"wp:attachment":[{"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/media?parent=806"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/categories?post=806"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/doublelayer.eu\/vilab\/wp-json\/wp\/v2\/tags?post=806"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}