Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def testtriclinic_relaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'triclinic',
optimizer=FIRE, fmax=self.fmax,
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)
def testtriclinic_unrelaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'triclinic',
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.2)
def test_Grochola(self):
a = FaceCenteredCubic('Au', size=[2,2,2])
calc = EAM('Au-Grochola-JCP05.eam.alloy')
a.set_calculator(calc)
FIRE(StrainFilter(a, mask=[1,1,1,0,0,0]), logfile=None).run(fmax=0.001)
a0 = a.cell.diagonal().mean()/2
self.assertTrue(abs(a0-4.0701)<2e-5)
self.assertTrue(abs(a.get_potential_energy()/len(a)+3.924)<0.0003)
C, C_err = fit_elastic_constants(a, symmetry='cubic', verbose=False)
C11, C12, C44 = Voigt_6x6_to_cubic(C)
self.assertTrue(abs((C11-C12)/GPa-32.07)<0.7)
self.assertTrue(abs(C44/GPa-45.94)<0.5)
def testmonoclinic_unrelaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'monoclinic',
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.2)
def testorthorhombic_relaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'orthorhombic',
optimizer=FIRE, fmax=self.fmax,
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.1)
def testcubic_relaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'cubic',
optimizer=FIRE, fmax=self.fmax,
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=1e-2)
if abs(C_err).max() > 0.:
self.assertArrayAlmostEqual(C_err/units.GPa, self.C_err_ref_relaxed, tol=1e-2)
def testmonoclinic_relaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'monoclinic',
optimizer=FIRE, fmax=self.fmax,
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)
def testorthorhombic_unrelaxed(self):
C, C_err = fit_elastic_constants(self.at0, 'orthorhombic',
verbose=False, graphics=False)
self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.1)
cryst = parameter('cryst').copy()
cryst.set_pbc(True)
# Double check elastic constants. We're just assuming this is really a periodic
# system. (True if it comes out of the cluster routines.)
compute_elastic_constants = parameter('compute_elastic_constants', False)
elastic_fmax = parameter('elastic_fmax', 0.01)
elastic_symmetry = parameter('elastic_symmetry', 'triclinic')
elastic_optimizer = parameter('elastic_optimizer', ase.optimize.FIRE)
if compute_elastic_constants:
cryst.set_calculator(calc)
log_file = open('elastic_constants.log', 'w')
C, C_err = fit_elastic_constants(cryst, verbose=False,
symmetry=elastic_symmetry,
optimizer=elastic_optimizer,
logfile=log_file,
fmax=elastic_fmax)
log_file.close()
logger.pr('Measured elastic constants (in GPa):')
logger.pr(np.round(C*10/GPa)/10)
crk = crack.CubicCrystalCrack(parameter('crack_surface'),
parameter('crack_front'),
Crot=C/GPa)
else:
if has_parameter('C'):
crk = crack.CubicCrystalCrack(parameter('crack_surface'),
parameter('crack_front'),
C=parameter('C'))
atom_types={'W': 1}, keep_alive=True)
calculator = lammps
unit_cell.set_calculator(calculator)
# simple calculation to get the lattice constant and cohesive energy
# alat0 = W.cell[0][1] - W.cell[0][0]
sf = StrainFilter(unit_cell) # or UnitCellFilter(W) -> to minimise wrt pos, cell
opt = FIRE(sf)
opt.run(fmax=1e-4) # max force in eV/A
alat = unit_cell.cell[0][1] - unit_cell.cell[0][0]
# print("a0 relaxation %.4f --> %.4f" % (a0, a))
# e_coh = W.get_potential_energy()
# print("Cohesive energy %.4f" % e_coh)
Cij, Cij_err = fit_elastic_constants(unit_cell,
symmetry="cubic",
delta=delta)
Cij = Cij/GPa # unit conversion to GPa
elasticMatrix3x3 = Cij[:3, :3]
# average of diagonal elements: C11, C22, C33
C11 = elasticMatrix3x3.diagonal().mean()
# make mask to extract non diagonal elements
mask = np.ones((3, 3), dtype=bool)
np.fill_diagonal(mask, False)
# average of all non diagonal elements from 1 to 3
C12 = elasticMatrix3x3[mask].mean()
# average of diagonal elements from 4 till 6: C44, C55, C66,