Help with CPA: Setup and Convergence

[spattamatta]

Hello All,

I am trying to setup a CPA calculation of L12 structure (of type AB3) with A site occupied by a mix of "Ti, Al and Fe" and the three B sites each occupied by a mix of "Ni, Co and Fe".

Following the examples and using ASE MuST, I setup the system. I know how to setup the system if all CPA atoms are of the same kind, but not sure if I am doing it correctly in case of two kinds of CPA atoms. Can someone take a look at the python code below and help me.

Secondly, when I run the code, it runs for a while and fails with the following error:

ERROR: rho0(ir) < 0
VALUE: -0.1586309D-01
STOP AT constructChargeDensity
For id = 1, ir = 977, r(ir) = 1.89201989
ERROR: rho0(ir) < 0
VALUE: -0.1586309D-01
STOP AT constructChargeDensity
For id = 1, ir = 977, r(ir) = 1.89201989
ERROR: rho0(ir) < 0
VALUE: -0.1586309D-01
STOP AT constructChargeDensity
For id = 1, ir = 977, r(ir) = 1.89201989
ERROR: rho0(ir) < 0
VALUE: -0.1586309D-01
STOP AT constructChargeDensity

Following another issue on this GitHub page (#31) and the suggestions by Professor Yang Wang, I tried changing k-points, mixing types and values etc, but it did not help.

Thank you,
Subrahmanyam Pattamatta
The University of Hong Kong

---

CODE to the setup

import os
import ase
from ase.build import bulk
from ase.units import kJ, Ry, Bohr
from ase.lattice.compounds import L1_2
from must import MuST, generate_starting_potentials

if name == 'main':

# Setup MuST calculator
calc = MuST(
    #####################################
    # Position and Potential Data Files #
    #####################################
    default_in_pot='Ni_ss_pot Co_ss_pot Ti_ss_pot Fe_ss_pot Al_ss_pot', 
    pot_in_form=0, 
    default_out_pot='NiCoTiFeAl',
    pot_out_form=1,
    ##########################
    # SCF-related Parameters #
    ##########################
    stop_rout_name='main', 
    nscf=200,
    method=3, 
    out_to_scr='n', 
    temperature=0.0, 
    ####################################
    # LDA Potential-related Parameters #
    ####################################
    potential_type=0,
    xc=0, 
    uniform_grid=(64,64,64),
    ####################################
    # Energy (Ryd.) Contour Parameters #
    ####################################
    read_mesh=0, 
    n_egrids=30,
    erbot=-0.4,
    real_axis_method=0,
    real_axis_points=300,
    ################################
    # Magnetism-related Parameters #
    ################################
    spin=1,
    #int_espin=1,
    ########################################
    # Scattering Theory-related Parameters #
    ########################################
    lmax_T=3, 
    ndivin=1001, 
    #########################################
    # R-space or K-space Related Parameters #
    #########################################
    liz_cutoff=8.5 * Bohr,
    k_scheme=0, 
    kpts=(12, 12, 12),
    bzsym=1,
    ###################################
    # Mixing and Tolerance Parameters #
    ###################################
    mix_algo=2,
    mix_quantity=1,
    mix_param=0.1,
    etol=5.0E-6 * Ry,
    ptol=1.0E-7 * Ry,
    ########################################################
    # Mixing and Tolerance Parameters for Effective Medium #
    ########################################################
    em_iter=80,
    em_scheme=2,
    em_mix_param=(0.1, 0.1),
    em_eswitch=0.003,
    em_tmatrix_tol=1E-7,
    #######################
    # No. of cores to run #
    #######################
    ntasks=16,
    )

# Experimental value of L12 lattice parameter in Angstrom
a = 3.6

# (Ti,Al,Fe) (Ni,Co,Fe)3 of type AB3 i.e L12 structure
# Create TiNi3 Alloy first
atoms = L1_2(directions=[[1,0,0], [0,1,0], [0,0,1]],
             size=(1,1,1), 
             symbol=('Ti','Ni'),
             pbc = (1,1,1),
             latticeconstant=a)
atoms = ase.build.sort(atoms) # Reorders atoms alphabetically -> Ni3 Ti

# Add CPA sites to make it (Ni43.3 Co23.7 Fe8)3 (Ti14.4 Al8.6 Fe2) 
atoms.info = {'CPA': [{'index': 0, 'Ni': 43.3/75, 'Co': 23.7/75, 'Fe': 8/75},
                      {'index': 1, 'Ni': 43.3/75, 'Co': 23.7/75, 'Fe': 8/75},
                      {'index': 2, 'Ni': 43.3/75, 'Co': 23.7/75, 'Fe': 8/75},
                      {'index': 3, 'Ti': 14.4/25, 'Al':  8.6/25, 'Fe': 2/25}]}    

# Move to new work directory
work_dir = 'work'
if not os.path.exists(work_dir):
    os.makedirs(work_dir)
os.chdir(work_dir)

generate_starting_potentials(atoms, crystal_type=1, a=a, moment=0.0, nspins=1, cpa=True)

# Attach calculator to the atoms object
atoms.set_calculator(calc)

# Trigger KKR calculation using .get_potential_energy() in eV
energy = atoms.get_potential_energy()
print(energy)

[wangy2014]

Hi, the problem you have seen is likely caused by divergence of the DFT potential. That is the self-consistent calculation is not convergent. The potential divergence could be caused by divergence of the CPA calculation. You may look into the output log file (o_n00*) to check the CPA calculation is indeed convergent. If not, you may need to increase the number of k-points and/or modify the mixing scheme and mixing parameter for the CPA calculation.