Constrained Optimization: Use VASP+ASE to Fix Bond Length, Angle, & Dihedral Angle during Relaxation

Fixing bond length, bond angles, and dihedral angles during structural optimization, also referred to as constrained optimization, is a crucial technique in molecular modeling. It enables the incorporation of experimental data or prior knowledge into the optimization process, resulting in more precise and realistic estimates of target properties. Additionally, it allows for the sampling of desired structures on the potential energy surface while keeping certain degrees of freedom fixed. This is particularly useful for gaining insights and searching for transition state structures. Despite the importance of constrained optimization, many commercial molecular modeling packages, such as VASP, lack support for this technique. In this tutorial, we will provide a script that utilizes the combination of the Atomic Simulation Environment (ASE) and VASP to enable constraint optimization. The approach involves running the optimization using one of ASE's built-in optimizers, while VASP acts as a single-point calculator for each ionic step.

Installation of ASE: https://wiki.fysik.dtu.dk/ase/install.html

Connect VASP with ASE & set VASP keywords in ASE: https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#environment-variables

from ase.io import read, write
from ase.calculators.vasp import Vasp
from ase.constraints import FixInternals
from ase.optimize.sciopt import SciPyFminCG
​
s = read('POSCAR')
#s.set_initial_magnetic_moments([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0])
s.pbc = True
​
calc = Vasp(xc='PBE',
           encut=600,
           kpts=[3,1,1],
           prec='Normal',
           sigma=0.1,
           ismear=0,
           isym=0,
           isif=2,
           ivdw=12,
           setups='recommended',
           lorbit = 11,
           lreal=True,
           ediff=1E-5,
           nelm=250,
           ispin=2)
           #lhfcalc = True,
           #hfscreen = 0.2,
           #aexx = 0.25,
           #aggax = 0.75,
           #aggac = 1.0,
           #aldac = 1.0)
​
dihedral_indices1 = [151, 37, 145, 103] # Refer to atom #152, 38, 146, 104, respectively
dihedral1 = [s.get_dihedral(*dihedral_indices1), dihedral_indices1]
c = FixInternals(dihedrals_deg=[dihedral1])
s.set_constraint(c)
s.set_calculator(calc)
s.get_potential_energy()
​
dyn=SciPyFminCG(s,trajectory='opt.traj',logfile='opt.log')
dyn.run(fmax=0.05)

Monitor the energy and maximum atomic force of each ionic step:

(phonopy) [xwb7910@quser32 IS]$ cat opt.log 
             Step     Time          Energy         fmax
*Force-consistent energies used in optimization.
SciPyFminCG:    0 02:39:33     -920.617143*       4.9303
SciPyFminCG:    1 03:05:29     -923.576702*       1.1004
SciPyFminCG:    2 03:44:46     -924.256351*       1.6563
SciPyFminCG:    3 04:04:20     -924.650707*       0.5656
SciPyFminCG:    4 04:29:52     -924.833918*       0.5231
SciPyFminCG:    5 05:02:50     -925.063383*       0.9259
...
SciPyFminCG:   32 16:16:37     -925.998534*       0.0680
SciPyFminCG:   33 16:34:16     -926.000330*       0.0663
SciPyFminCG:   34 17:05:12     -926.002758*       0.1025
SciPyFminCG:   35 17:23:46     -926.005156*       0.1038
SciPyFminCG:   36 17:40:14     -926.006960*       0.0491 #< 0.05, converged

---

If you found this tutorial helpful, consider joining our community by making a donation. Your support will enable me to continue creating valuable content while also supporting my baby's care. Together, we can make a meaningful impact. Thank you!