Problem with using VASP with ase to predict the properties of materials

ASE
1

I am a newbie using VASP, currently I am trying to calculate the formation energy per atom and pre bandgap using dft with VSAP.

I wrote script for calculating the bandgap

# -*- coding: utf-8 -*-
from ase.io import read
from ase.calculators.vasp import Vasp
from ase.optimize import BFGS
import numpy as np
import os

os.environ['VASP_PP_PATH'] = '/scratch'

# Load the structure from a CIF file
structure = read('mp-1218989.cif')

# Set up VASP calculator for relaxation
calc_relax = Vasp(xc='PBE',  # Exchange-correlation functional
                  encut=400,  # Plane-wave cutoff
                  kpts=(3, 3, 3),  # k-point grid for relaxation
                  ibrion=2,  # Ion relaxation: CG algorithm
                  nsw=50,  # Number of steps for ionic relaxation
                  ismear=0,  # Gaussian smearing
                  sigma=0.05,  # Width of smearing in eV
                  lreal='Auto',  # Use real-space projection for GPU acceleration
                  command='mpirun -np 1 vasp_gpu')  # VASP run command

# Attach the calculator to the structure for relaxation
structure.set_calculator(calc_relax)

# Optimize the structure
opt = BFGS(structure)
opt.run(fmax=0.01)

# Now setting up the NSCF calculation for band structure
calc_nscf = Vasp(xc='PBE',  # Exchange-correlation functional
                 encut=400,  # Plane-wave cutoff
                 kpts=(11, 11, 11),  # denser k-point mesh for band structure
                 ismear=0,  # Gaussian smearing (must use ISMEAR < -1 or 0)
                 sigma=0.05,
                 ibrion=-1,  # No ionic relaxation
                 nsw=0,  # No ionic steps
                 icharg=11,  # Non-self-consistent field run
                 command='mpirun -np 1 vasp_gpu')

# Attach the non-self-consistent calculator
structure.set_calculator(calc_nscf)

# Calculate the electronic structure (no need to run optimize)
structure.calc.calculate(structure)

# Extract eigenvalues
eigenvalues = structure.calc.get_eigenvalues(kpt=0, spin=0)  # For the first k-point and spin

# Calculate bandgap
occupied = np.max(eigenvalues[eigenvalues < 0])
unoccupied = np.min(eigenvalues[eigenvalues > 0])
bandgap = unoccupied - occupied

print('Bandgap (PBE):', bandgap, 'eV')

I still didn’t figured out how to calculate the formation energy

The above code is not utilising the gpu properly, I have tried to increase the value here mpirun -np 1 vasp_gpu to more than 1 but its just showing my system is not capable, but cleary I have 2 gpus and 20 cpus in my system.

  1. I am planning to perfom this on 100 molecules so any better method ?

please help if any one knows simple method to automate the calculation of the properties using vasp

2

No idea about the GPU side of things, but I notice that you have some redundancy in the geometry optimisation:

  • The VASP parameters include ibrion/nsw options that will cause each call to VASP to use the internal geometry optimiser
  • In the Python script you use an ASE BFGS object to optimise the geometry. This will repeatedly modify the positions and request VASP for forces

You won’t get sensible results using both of these at the same time, I recommend either to

  • use the internal VASP options, calculate the total energy or forces with e.g. structure.get_forces() and then read the CONTCAR to get a new optimised structure, or
  • set IBRION=-1, NSW=1 in the VASP parameters and use the ASE optimiser.

To get the final DFT energy with use structure.get_total_energy(). To get a formation energy you will need to compare this with the calculated energies of reference materials.

1 Like