#!/usr/bin/env python
"""
This seems to work for fcc, bcc, sc and diamond
Things to figure out:
* why multiply strain on right
* why we can't use changing volumes
Last Update: 2019/04/13 Ellad Tadmor (added ability to do diamond)
"""
from __future__ import print_function
from ase.build import bulk
from ase.units import GPa
import os
import numpy as np
import numdifftools as ndt
from ase.calculators.kim.kim import KIM
from scipy.optimize import fmin
from scipy.optimize import minimize
import sys
import json
import jinja2
class ElasticConstants(object):
"""Determine the cubic elastic constants by numerically determining the Hessian"""
def __init__(self, calc, element, potentialname, crystalstructure, latticeconst):
self.calculator = calc
self.element = element
self.potentialname = potentialname
self.crystalstructure = crystalstructure
self.latticeconst = latticeconst
self.slab = self.create_slab()
self.o_cell = self.slab.get_cell()
self.slab.set_calculator(self.calculator)
self.o_volume = self.slab.get_volume()
def create_slab(self):
slab = bulk(self.element,a=self.latticeconst,crystalstructure=self.crystalstructure,cubic=True)
return slab
def voigt_to_matrix(self,voigt_vec):
"""Convert a voigt notation vector to a matrix """
matrix = np.zeros((3,3))
matrix[0,0] = voigt_vec[0]
matrix[1,1] = voigt_vec[1]
matrix[2,2] = voigt_vec[2]
matrix[ tuple([[1,2], [2,1]]) ] = voigt_vec[3]
matrix[ tuple([[0,2], [2,0]]) ] = voigt_vec[4]
matrix[ tuple([[0,1], [1,0]]) ] = voigt_vec[5]
return matrix
def get_energy_from_positions(self,pos):
natoms = len(pos)/3
self.slab.set_positions(np.reshape(pos,(natoms,3)))
energy = self.slab.get_potential_energy()
return energy
def get_gradient_from_positions(self,pos):
natoms = len(pos)/3
self.slab.set_positions(np.reshape(pos,(natoms,3)))
forces = self.slab.get_forces()
return -forces.flatten()
def energy_from_strain(self,strain_vec):
""" Apply a strain according to the strain_vec """
# self.slab = self.o_slab.copy()
# print strain_vec
strain_mat = self.voigt_to_matrix(strain_vec)
old_cell = self.o_cell
new_cell = old_cell + np.dot(old_cell, strain_mat)
# new_cell = old_cell + np.einsum('ij,aj->ai',strain_mat,old_cell)
self.slab.set_cell(new_cell,scale_atoms=True)
if self.crystalstructure == "diamond":
pos0 = self.slab.get_positions().flatten()
res = minimize(self.get_energy_from_positions,pos0,jac=self.get_gradient_from_positions)
energy = res["fun"]
else:
energy = self.slab.get_potential_energy()
# volume = self.slab.get_volume()
# n_of_atoms = self.slab.get_number_of_atoms()
return (energy/self.o_volume)/GPa
def energy_from_scale(self,scale):
# strain_mat = np.eye(3) * (scale)
old_cell = self.o_cell
# new_cell = old_cell + np.dot( strain_mat, old_cell)
new_cell = old_cell * (1 + scale)
self.slab.set_cell(new_cell,scale_atoms=True)
energy = self.slab.get_potential_energy()
# volume = self.slab.get_volume()
# n_of_atoms = self.slab.get_number_of_atoms()
return energy
def results(self):
""" Return the cubic elastic constants """
#get the minimum
self.minscale = float(fmin(self.energy_from_scale,0,xtol=0,ftol=1e-7))
self.oo_cell = self.o_cell.copy()
self.o_cell = self.o_cell * (1+self.minscale)
func = self.energy_from_strain
hess = ndt.Hessian(func, step=0.001, full_output=True)
elastic_constants, info = hess(np.zeros(6, dtype=float))
error_estimate = info.error_estimate
inds11 = tuple([[0,1,2],[0,1,2]])
C11 = elastic_constants[ inds11 ].mean()
C11sig = np.sqrt( ((1./3*error_estimate[ inds11 ])**2).sum() )
inds12 = tuple([[1,2,2,0,0,1],[0,0,1,1,2,2]])
C12 = ( elastic_constants[inds12].mean() )
C12sig = np.sqrt( ((1./6*error_estimate[inds12])**2).sum() )
inds44 = tuple([[3,4,5],[3,4,5]])
C44 = 1./4 * elastic_constants[inds44].mean()
C44sig = 1./4 * np.sqrt( ((1./3*error_estimate[inds44])**2).sum() )
B = 1./3 * ( C11 + 2 * C12 )
Bsig = np.sqrt( ( 1./3 * C11sig )**2 + ( 2./3 * C12sig )**2 )
excessinds = tuple([[3,4,5,3,4,5,3,4,5,0,1,2,4,5,0,1,2,3,5,0,1,2,3,4],
[0,0,0,1,1,1,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5]])
excess =(np.abs(elastic_constants[excessinds]).mean() )
excess_sig = np.sqrt( ((1./24*error_estimate[excessinds])**2).sum() )
results_dict = { 'C11': C11,
'C11_sig' : C11sig,
'C12' : C12,
'C12_sig' : C12sig,
'C44': C44,
'C44_sig' : C44sig,
'B' : B,
'B_sig' : Bsig,
'excess': excess,
'excess_sig' : excess_sig,
'units' : 'GPa',
"element": symbol,
"crystal_structure": lattice,
"space_group": space_groups[lattice],
"wyckoff_code": wyckoff_codes[lattice],
"lattice_constant": self.latticeconst,
"scale_discrepency": self.minscale,
}
return results_dict
symbol = raw_input()
lattice = raw_input()
model = raw_input()
latticeconst_result = raw_input()
space_groups = {"fcc": "Fm-3m", "bcc": "Im-3m", "sc": "Pm-3m", "diamond": "Fd-3m"}
wyckoff_codes = {"fcc": "4a", "bcc": "2a", "sc": "1a", "diamond": "8a"}
normed_basis = {
lattice: json.dumps(bulk(symbol, lattice, a=1, cubic=True).positions.tolist(), separators=(' ', ' '))
for lattice in space_groups.keys()
}
count = bulk(symbol, lattice, a=1, cubic=True).positions.shape[0] # For species
latticeconst = float(latticeconst_result) * 1e10
print(symbol, lattice, model, latticeconst)
calc = KIM(model)
bulkmodulus = ElasticConstants(calc, symbol, model, lattice, latticeconst)
results = bulkmodulus.results()
results.update({"basis_coordinates": normed_basis[lattice]})
# Repeat symbols to match normed basis
results.update({"species": '" "'.join([symbol] * count)})
# Echo results
print()
print("C11 = {:12.6f} +/- {:12.6f} GPa".format(results["C11"],results["C11_sig"]))
print("C12 = {:12.6f} +/- {:12.6f} GPa".format(results["C12"],results["C12_sig"]))
print("C44 = {:12.6f} +/- {:12.6f} GPa".format(results["C44"],results["C44_sig"]))
print("B = {:12.6f} +/- {:12.6f} GPa".format(results["B"],results["B_sig"]))
print()
template_environment = jinja2.Environment(
loader=jinja2.FileSystemLoader('/'),
block_start_string='@[',
block_end_string=']@',
variable_start_string='@<',
variable_end_string='>@',
comment_start_string='@#',
comment_end_string='#@',
undefined=jinja2.StrictUndefined,
)
#template the EDN output
with open(os.path.abspath("output/results.edn"), "w") as f:
template = template_environment.get_template(os.path.abspath("results.edn.tpl"))
f.write(template.render(**results))