#!/usr/bin/env python
"""
Lattice constant Test Driver for Hexagonal Structure

Computes the lattice constant for:
1. Any element,
2. HCP crystal,
3. At 0 K, 0 GPa,
by simplex minimization.

Structure Overview:
1. main() deals with I/O and get information from getLatticeConstant()
2. getLatticeConstant() uses an initial guess at the lattice constant 'a' and calls searchLatticeConstants()
3. searchLatticeConstants() minimizes getEnergy() with simplex algorithm
4. When lower energy is obtained with separation > init_upper_guess, collapseHandler() is called

Date: 2015/08/25
Author: Junhao Li <streaver91@gmail.com>

Last Update: 2016/02/26 Junhao Li

             2017/08/28 Daniel S. Karls
"""
# Import External Modules
# ASE Modules
from ase.build import bulk
from ase import Atoms

# KIM Modules
from ase.calculators.kim.kim import KIM

# Other Python Modules
import sys
import json
import math
import jinja2
import os
import copy
import numpy as np
from scipy.optimize import fmin
from collections import OrderedDict

# Constants for result.edn
SPACE_GROUP = 'P63/mmc'
WCYKOFF_CODE = '2a'
KEY_SOURCE_VALUE = 'source-value'
KEY_SOURCE_UNIT = 'source-unit'
KEY_SOURCE_UNCERT = 'source-std-uncert-value'
UNIT_LENGTH = 'angstrom'
UNIT_TEMPERATURE = 'K'
UNIT_PRESSURE = 'GPa'
UNIT_ENERGY = 'eV'
STRUCTURE_PROP_ID = 'tag:staff@noreply.openkim.org,2014-04-15:property/structure-hexagonal-crystal-npt'
ENERGY_PROP_ID = 'tag:staff@noreply.openkim.org,2014-04-15:property/cohesive-free-energy-hexagonal-crystal'

# Constants frequently used in calculation
SQRT3 = math.sqrt(3.0)
PERFECT_CA = math.sqrt(8.0 / 3.0)
HCP_CUBIC_CELL = np.array([1.0, SQRT3, 1.0])
HCP_CUBIC_POSITIONS = np.array([
    [0.0, 0.0, 0.0],
    [0.5, 0.5 * SQRT3, 0.0],
    [0.5, 0.5 / 3.0 * SQRT3, 0.5],
    [0.0, (0.5 + 0.5 / 3.0) * SQRT3, 0.5]
])
HCP_CUBIC_NATOMS = 4

# Minimization Convergence Criteria
FMIN_FTOL = 1e-10
FMIN_XTOL = 1e-10

def getEnergy(cellVector, meta, cache):
    # cellVector: 1d numpy array of size 1(2), containing a (and c)
    # meta: various information...
    # cache: cache previously created atoms

    # Determine a and c based on the length of the cellVector
    a = cellVector[0]
    if len(cellVector) == 2:
        c = cellVector[1]
    else:
        c = a * PERFECT_CA # Default to perfect c/a ratio

    cellSize = HCP_CUBIC_CELL * [a, a, c]

    repeat = (1, 1, 1)
    nAtoms = HCP_CUBIC_NATOMS

    # Get atoms from cache or create a new one if haven't been created
    sizes = cache['sizes']
    if repeat in sizes:
        sizeId = sizes.index(repeat)
        atoms = cache['atoms'][sizeId]
    else:
        # create a new atoms and save to cache
        print 'Creating new atoms:', repeat
        sizes.append(repeat)
        atoms = meta['prototype'].copy()
        atoms *= repeat
        atoms.set_calculator(KIM(meta['model']))

        cache['atoms'].append(atoms)

    # Scale atoms with the scaled cellSize
    atoms.set_cell(cellSize, scale_atoms = True)
    energy = atoms.get_potential_energy() / nAtoms
    return energy

def searchLatticeConstants(latticeConstantsGuess, meta, cache):
    # Simplex searching lattice constants that minimize potential energy of atoms
    # Searching starts from latticeConstantsGuess
    print 'Simplex search starting from:', latticeConstantsGuess
    tmpLatticeConstants, tmpEnergy = fmin(
        getEnergy,
        latticeConstantsGuess,
        args = (meta, cache),
        full_output = True,
        ftol = FMIN_FTOL,
        xtol = FMIN_XTOL,
    )[:2]
    return tmpLatticeConstants, tmpEnergy

def getInterlayerDist(latticeConstants):
    # Calculate the distance between atom 0 and 1 (refer to HCP_CUBIC_POSITIONS)
    hrDist = latticeConstants[0] / SQRT3
    vDist = latticeConstants[1] / 2
    return math.sqrt(hrDist * hrDist + vDist * vDist)

def collapseHandler():
    # Called when lower energy is obtained by separating atoms beyond cutoff
    print 'ERR: Lower energy is obtained by separating atoms beyond cutoff.'
    print 'ERR: System Collapsed. Structure Maybe Unstable.'
    sys.exit(0)

def getLatticeConstant(elem, model):
    # Initialize meta and cache
    meta = {
        'model': model
    }
    cache = {
        'sizes': [],
        'atoms': []
    }

    # Create prototype atoms
    atoms = Atoms(elem + '4',
        positions = HCP_CUBIC_POSITIONS,
        cell = HCP_CUBIC_CELL,
        pbc = [1, 1, 1]
    )

    meta['prototype'] = atoms

    # Relaxation with c/a Ratio Fixed to PERFECT_CA
    print 'Relaxation with c/a ratio fixed'
    init_upper_guess = 5.0
    minEnergy = 0
    for i in [x * 0.05 + 0.5 for x in range(-2, 3)]:
        print 'Simplex search starting from: %r * %i' % (init_upper_guess, i)
        tmpLatticeConstants, tmpEnergy = searchLatticeConstants(
            [init_upper_guess * i],
            meta,
            cache
        )
        print 'Tmp Lattice Constants:', tmpLatticeConstants
        print 'Tmp Energy:', tmpEnergy
        if tmpEnergy < minEnergy and tmpLatticeConstants[0] < init_upper_guess:
            # Make sure there are interatomic forces holding the structure
            minLatticeConstants = tmpLatticeConstants
            minEnergy = tmpEnergy
        print '--------'

    if minEnergy == 0:
        collapseHandler()

    # Relaxation With c/a Ratio Relaxed
    print 'Relaxation with c/a ratio relaxed'
    tmpA = minLatticeConstants[0]
    minEnergy = 0
    for i in [x * 0.05 + 1 for x in range(-4, 5)]:
        print 'Simplex search starting from c/a ratio: 1.633 *', i
        tmpLatticeConstants, tmpEnergy = searchLatticeConstants(
            [tmpA, tmpA * PERFECT_CA * i],
            meta,
            cache
        )
        print 'Tmp Lattice Constants:', tmpLatticeConstants
        print 'Tmp Energy:', tmpEnergy
        if tmpEnergy < minEnergy and \
            tmpLatticeConstants[0] < init_upper_guess and \
            getInterlayerDist(tmpLatticeConstants) < init_upper_guess:
            # Make sure there is both intralayer and interlayer forces
            minLatticeConstants = tmpLatticeConstants
            minEnergy = tmpEnergy
        print '--------'

    if minEnergy == 0:
        collapseHandler()

    print 'Lattice Constants:', minLatticeConstants
    print 'Energy:', minEnergy

    # Return Results
    return minLatticeConstants[0], minLatticeConstants[1], minEnergy

def V(value, unit = '', uncert = ''):
    # Generate OrderedDict for KIM JSON Output
    res = OrderedDict([
        (KEY_SOURCE_VALUE, value),
    ])
    if unit != '':
        res.update(OrderedDict([
            (KEY_SOURCE_UNIT, unit),
        ]))
    if uncert != '':
        res.update(OrderedDict([
            (KEY_SOURCE_UNCERT, uncert)
        ]))
    return res

def main():
    # Input Parameters
    elem = raw_input('Element = ')
    lattice = raw_input('Lattice = ')
    model = raw_input('Model = ')

    # Parameters for Debugging
    # ========
    # elem = 'Co'
    # lattice = 'hcp'
    # model = 'EAM_Dynamo_PurjaPun_Mishin_Co__MO_885079680379_001'

    # elem = 'Si'
    # lattice = 'hcp'
    # model = 'EDIP_BOP_Bazant_Kaxiras_Si__MO_958932894036_001'

    # elem = 'Al'
    # lattice = 'hcp'
    # model = 'EAM_Dynamo_Mishin_Farkas_Al__MO_651801486679_001'

    # elem = 'C'
    # lattice = 'hcp'
    # model = 'EAM_Dynamo_Hepburn_Ackland_FeC__MO_143977152728_001'
    # model = 'Tersoff_LAMMPS_Erhart_Albe_CSi__MO_903987585848_001'
    # ========

    # Print Inputs
    print 'Element:', elem
    print 'Lattice:', lattice # Not used here
    print 'Model:', model

    # Obtain Lattice Constants and Cohesive Energy
    a, c, energy = getLatticeConstant(elem, model)
    print 'Lattice Constants:', a, c

    # Output Results
    structureResults = OrderedDict([
        ('property-id', STRUCTURE_PROP_ID),
        ('instance-id', 1),
        ('cauchy-stress', V([0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UNIT_PRESSURE)),
    ])
    cohesiveEnergyResults = OrderedDict([
        ('property-id', ENERGY_PROP_ID),
        ('instance-id', 2),
        ('cohesive-free-energy', V(-energy, UNIT_ENERGY)),
    ])
    commonResults = OrderedDict([
        ('short-name', V(['hcp'])),
        ('species', V([elem, elem])),
        ('a', V(a, UNIT_LENGTH)),
        ('c', V(c, UNIT_LENGTH)),
        ('basis-atom-coordinates', V([[0.0, 0.0, 0.0], [2.0/3, 1.0/3, 0.5]])),
        ('space-group', V(SPACE_GROUP)),
        ('temperature', V(0, UNIT_TEMPERATURE)),
    ])
    structureResults.update(commonResults)
    cohesiveEnergyResults.update(commonResults)
    results = [
        structureResults,
        cohesiveEnergyResults,
    ]
    resultsString = json.dumps(results, separators = (' ', ' '), indent = 4)
    print resultsString
    with open(os.path.abspath('output/results.edn'), 'w') as f:
        f.write(resultsString)

if __name__ == '__main__':
    main()