#!/usr/bin/env python

"""
Run sats and generate some configurations. With self-interstitial atoms.

"""

import sys

from quippy.potential import Potential, Minim

from sats.core.dynamics import molecular_dynamics
from sats.core.io import castep_write
from sats.configurations.interstitial import interstitial_dumbbell
from sats.configurations.interstitial import interstitial_crowdion
from sats.configurations.interstitial import interstitial_tetrahedral
from sats.configurations.interstitial import interstitial_octahedral

fs_lattice_parameter = 3.18
fcc_lattice_parameter = 3.95
hcp_lattice_parameter = 2.79
species = 74

fs_potential = Potential('IP FS')
fs_potential.set_calc_args({'E_scale': 0.99519, 'r_scale': 0.99302})

lattice_parameter = fs_lattice_parameter
supercell = (3, 3, 3)

interstitial = sys.argv[1]

if 'db' in interstitial:
    lattice, _db, direction = interstitial.split('_')
    bulk = interstitial_dumbbell(lattice, lattice_parameter, species='W',
                                 direction=direction, supercell=supercell)
elif 'crw' in interstitial:
    lattice, _crw = interstitial.split('_')
    bulk = interstitial_crowdion(lattice, lattice_parameter, species='W',
                                 supercell=supercell)
elif 'tet' in interstitial:
    lattice, _tet = interstitial.split('_')
    bulk = interstitial_tetrahedral(lattice, lattice_parameter, species='W',
                                    supercell=supercell)
elif 'oct' in interstitial:
    lattice, _oct = interstitial.split('_')
    bulk = interstitial_octahedral(lattice, lattice_parameter, species='W',
                                   supercell=supercell)
else:
    print("Unknown interstitial {0}".format(interstitial))
    raise SystemExit

bulk.name = interstitial

castep_write(bulk, interstitial, optimise=True, fix_lattice=False)

bulk.set_cutoff(fs_potential.cutoff() + 2.0)
bulk.set_calculator(fs_potential)

minimiser = Minim(bulk, relax_positions=True, relax_cell=True)
minimiser.run()

temperature = int(sys.argv[2])
# From 10.1007/BF01184339, thermal expansion of tungsten is small, but for
# different temperatures we can expand a bit:
# 3000 K V/V0 = 1.11;
# 5000 K V/V0 = 1.29
#
if temperature == 1000:
    bulk.set_lattice(bulk.lattice*1.003, scale_positions=True)
elif temperature == 3000:
    bulk.set_lattice(bulk.lattice*(1.11**(1.0/3)), scale_positions=True)
elif temperature == 5000:
    bulk.set_lattice(bulk.lattice*(1.29**(1.0/3)), scale_positions=True)

molecular_dynamics(bulk, fs_potential, temperature)
