#!/usr/bin/env python

"""Bulk Al(fcc) equation of state

This script calculates the total energy of an fcc Al unit cell
at several different volumes for fitting to an equation of state

john kitchin May 20, 2005
"""
import os
from Dacapo import Dacapo
from ASE import Atom,ListOfAtoms
from ASE.Utilities.GeometricTransforms import SetUnitCellVolume

bulk = ListOfAtoms([Atom('Al', (0, 0, 0))] )

b = 4.05
bulk.SetUnitCell([(0,b/2,b/2),
                  (b/2,0,b/2),
                  (b/2,b/2,0)])

calc = Dacapo(kpts=(4,4,4),               
		  planewavecutoff=350,
		  nbands=6,                   
		  usesymm=False)
    
bulk.SetCalculator(calc)

initial_volume = bulk.GetUnitCellVolume()

energies = []
volumes = []

for fraction in [0.9, 0.95, 1.0, 1.05, 1.1]:

    ncfile = 'bulkAl_%1.2f.nc' % fraction
    #only run calculations if the files don't exist
    if not os.path.exists(ncfile):
        new_volume = fraction * initial_volume

        SetUnitCellVolume(bulk,new_volume)

        calc.SetNetCDFFile(ncfile)

        bulk.GetPotentialEnergy()

    atoms = Dacapo.ReadAtoms(ncfile)
    energies.append(atoms.GetPotentialEnergy())
    volumes.append(atoms.GetUnitCellVolume())

from ASE.Utilities.EquationOfState import EquationOfState
eos = EquationOfState('Murnaghan',volumes,energies)
print eos
eos.GetPlot()
eos.SavePlot('Al_murn.png')

V0 = eos.GetV0()

'''
the volume of a primitive fcc cell is equal to
V=1/4*a^3 where a is the lattice constant.
'''
a = (V0*4)**(1./3.)
print 'lattice constant = %1.2f angstroms' % a

'''
Typical output:

V0(A^3)       B(Gpa)   E0 (eV)    pressure(GPa) Murnaghan
 16.78        85.91     -56.4668       -0.00
chisq = 0.0000
lattice constant = 4.05 angstroms
'''