Using built-in ase optimization with vasp

ASE has some nice optimization tools built into it. We can use them in vasp too. This example is adapted from this test: https://trac.fysik.dtu.dk/projects/ase/browser/trunk/ase/test/vasp/vasp_Al_volrelax.py

First the VASP way.

from vasp import Vasp
from ase.build import bulk

Al = bulk('Al', 'fcc', a=4.5, cubic=True)
calc = Vasp(label='bulk/Al-lda-vasp',
            xc='LDA', isif=7, nsw=5,
            ibrion=1, ediffg=-1e-3,
            lwave=False, lcharg=False,
            atoms=Al)
print(calc.potential_energy)
print(calc)

-10.07430731
============================================================
VASP Calculation: /home/jovyan/dft-book/notebooks/04-bulk-systems/bulk/Al-lda-vasp
============================================================

Formula: Al4
Number of atoms: 4

Unit cell:
  Volume: 73.3685 ų
  a=4.1864 b=4.1864 c=4.1864 Å
  α=90.00° β=90.00° γ=90.00°

Energy: -10.074307 eV (-2.518577 eV/atom)
Max force: 0.000000 eV/Å
Pressure: -1.533 GPa
Stress (GPa): xx=1.533 yy=1.533 zz=1.533
              yz=0.000 xz=0.000 xy=0.000

Parameters:
  xc: LDA
  kpts: (1, 1, 1)
  ismear: 1
  sigma: 0.1
  ibrion: 1
  isif: 7
  nsw: 5
============================================================

Now, the ASE way using ExpCellFilter for cell optimization. This filter uses an exponential parametrization of the cell deformation which is more numerically stable than the older StrainFilter.

from vasp import Vasp
from ase.build import bulk
from ase.optimize import BFGS
from ase.constraints import ExpCellFilter  # More stable than StrainFilter
import os

Al = bulk('Al', 'fcc', a=4.5, cubic=True)

calc = Vasp(label='bulk/Al-lda-ase',
            xc='LDA',
            atoms=Al)

calc.calculate()

# Check if we have a completed calculation to work with
if os.path.exists(os.path.join(calc.directory, 'OUTCAR')):
    # Get initial energy to ensure calculator is properly attached
    Al.get_potential_energy()
    
    # Use ExpCellFilter instead of StrainFilter for better numerical stability
    ecf = ExpCellFilter(Al)
    opt = BFGS(ecf, logfile='relaxation.log')
    opt.run(fmax=0.1, steps=5)
    
    print('Stress:\n', calc.stress)
    print('Al post ASE volume relaxation\n', calc.load_atoms().get_cell())
    print(calc)
else:
    print('VASP calculation not yet complete. Run previous cell first.')
    print('Note: ASE cell optimization requires VASP to compute stress.')

Stress:
 [ 7.34461581  7.34461581  7.34461581 -0.         -0.         -0.        ]
Al post ASE volume relaxation
 Cell([4.5, 4.5, 4.5])
============================================================
VASP Calculation: /home/jovyan/dft-book/notebooks/04-bulk-systems/bulk/Al-lda-ase
============================================================

Formula: Al4
Number of atoms: 4

Unit cell:
  Volume: 91.1250 ų
  a=4.5000 b=4.5000 c=4.5000 Å
  α=90.00° β=90.00° γ=90.00°

Energy: -9.584487 eV (-2.396122 eV/atom)
Max force: 0.000000 eV/Å
Pressure: -1176.737 GPa
Stress (GPa): xx=1176.737 yy=1176.737 zz=1176.737
              yz=-0.000 xz=-0.000 xy=-0.000

Parameters:
  xc: LDA
  kpts: (1, 1, 1)
  ismear: 1
  sigma: 0.1
============================================================

Now for a detailed comparison:

from vasp import Vasp
import os
import numpy as np

calc1 = Vasp(label='bulk/Al-lda-vasp')
calc2 = Vasp(label='bulk/Al-lda-ase')

# Check if both calculations exist
if (os.path.exists(os.path.join(calc1.directory, 'OUTCAR')) and 
    os.path.exists(os.path.join(calc2.directory, 'OUTCAR'))):
    atoms1 = calc1.load_atoms()
    atoms2 = calc2.load_atoms()

    cellA = atoms1.get_cell()
    cellB = atoms2.get_cell()

    print('VASP-relaxed cell:')
    print(cellA)
    print(f'Volume: {atoms1.get_volume():.4f} ų\n')
    
    print('ASE-optimized cell:')
    print(cellB)
    print(f'Volume: {atoms2.get_volume():.4f} ų\n')

    match = (np.abs(cellA - cellB) < 0.01).all()
    print(f'Cells match within 0.01 Å: {match}')
    if not match:
        print('Note: If ASE optimization failed, the cell remains at initial geometry.')
else:
    print('One or both calculations not complete. Run previous cells first.')

VASP-relaxed cell:
Cell([4.18635928601072, 4.18635928601072, 4.18635928601072])
Volume: 73.3685 ų

ASE-optimized cell:
Cell([4.5, 4.5, 4.5])
Volume: 91.1250 ų

Cells match within 0.01 Å: False
Note: If ASE optimization failed, the cell remains at initial geometry.

This could be handy if you want to use any of the optimizers in mod:ase.optimize in conjunction with mod:ase.constraints, which are more advanced than what is in VASP.