# Set up an alpha-quartz SiO2-crystal. lattice = Hexagonal(4.916*Angstrom, 5.4054*Angstrom) # Define elements elements = [Silicon, Silicon, Silicon, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen] # Define coordinates fractional_coordinates = [[ 0.4697 , 0. , 0. ], [ 0. , 0.4697 , 0.666666666667], [ 0.5303 , 0.5303 , 0.333333333333], [ 0.4135 , 0.2669 , 0.1191 ], [ 0.2669 , 0.4135 , 0.547567 ], [ 0.7331 , 0.1466 , 0.785767 ], [ 0.5865 , 0.8534 , 0.214233 ], [ 0.8534 , 0.5865 , 0.452433 ], [ 0.1466 , 0.7331 , 0.8809 ]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- potential_set = TremoloXPotentialSet('Trinastic_JCP2013') potential_set.addParticleType(ParticleType(symbol='O', mass=15.9994*atomic_mass_unit, charge=-1.2)) potential_set.addParticleType(ParticleType(symbol='Ti', mass=47.867*atomic_mass_unit, charge=2.4)) potential_set.addParticleType(ParticleType(symbol='Si', mass=28.0855*atomic_mass_unit, charge=2.4)) potential_set.addParticleType(ParticleType(symbol='Ta', mass=180.948*atomic_mass_unit, charge=3.0)) potential_set.addParticleType(ParticleType(symbol='Hf', mass=178.49*atomic_mass_unit, charge=2.4)) # O-O interactions potential_set.addPotential(BuckinghamPotential(particleType1='O', particleType2='O', A=1388.77*eV, rho=0.3623*Angstrom, r_i=6.0*Angstrom, r_cut=7.5*Angstrom)) potential_set.addPotential(LennardJonesMNPotential(particleType1='O', particleType2='O', A=0.0*eV*Angstrom**12, B=175.0*eV*Angstrom**6, m=12, n=6, r_cut=10.0*Angstrom)) # Si-O interactions potential_set.addPotential(BuckinghamPotential(particleType1='O', particleType2='Si', A=18003.76*eV, rho=0.2052*Angstrom, r_i=6.0*Angstrom, r_cut=7.5*Angstrom)) potential_set.addPotential(LennardJonesMNPotential(particleType1='O', particleType2='Si', A=0.0*eV*Angstrom**12, B=133.54*eV*Angstrom**6, m=12, n=6, r_cut=10.0*Angstrom)) # Ta-O (weak) interactions potential_set.addPotential(BuckinghamPotential(particleType1='O', particleType2='Ta', A=100067.01*eV, rho=0.1319*Angstrom, r_i=6.0*Angstrom, r_cut=7.5*Angstrom)) potential_set.addPotential(LennardJonesMNPotential(particleType1='O', particleType2='Ta', A=0.0*eV*Angstrom**12, B=6.05*eV*Angstrom**6, m=12, n=6, r_cut=10.0*Angstrom)) potential_set.addPotential(MorsePotential(particleType1='O', particleType2='Ta', E_0=0.3789*eV, k=1.6254*Angstrom**-1, r_0=2.5445*Angstrom, r_i=6.5*Angstrom, r_cut=7.0*Angstrom)) # Hf-O (weak) interactions potential_set.addPotential(BuckinghamPotential(particleType1='O', particleType2='Hf', A=12372.16*eV, rho=0.2286*Angstrom, r_i=6.0*Angstrom, r_cut=7.5*Angstrom)) potential_set.addPotential(LennardJonesMNPotential(particleType1='O', particleType2='Hf', A=0.0*eV*Angstrom**12, B=81.36*eV*Angstrom**6, m=12, n=6, r_cut=10.0*Angstrom)) potential_set.addPotential(MorsePotential(particleType1='O', particleType2='Hf', E_0=0.3474*eV, k=1.6230*Angstrom**-1, r_0=2.048*Angstrom, r_i=6.5*Angstrom, r_cut=7.0*Angstrom)) # Ti-O (weak) interactions potential_set.addPotential(BuckinghamPotential(particleType1='O', particleType2='Ti', A=5105.12*eV, rho=0.2253*Angstrom, r_i=6.0*Angstrom, r_cut=7.5*Angstrom)) potential_set.addPotential(LennardJonesMNPotential(particleType1='O', particleType2='Ti', A=0.0*eV*Angstrom**12, B=20.0*eV*Angstrom**6, m=12, n=6, r_cut=10.0*Angstrom)) potential_set.addPotential(MorsePotential(particleType1='O', particleType2='Ti', E_0=0.3478*eV, k=1.9*Angstrom**-1, r_0=1.96*Angstrom, r_i=6.5*Angstrom, r_cut=7.0*Angstrom)) # # Ti-O (strong) interactions # potential_set.addPotential(BuckinghamPotential(particleType1='O', # particleType2='Ti', # A=5505.12*eV, # rho=0.2253*Angstrom, # r_i=6.0*Angstrom, # r_cut=7.5*Angstrom)) # potential_set.addPotential(LennardJonesMNPotential(particleType1='O', # particleType2='Ti', # A=0.0*eV*Angstrom**12, # B=20.0*eV*Angstrom**6, # m=12, # n=6, # r_cut=10.0*Angstrom)) # potential_set.addPotential(MorsePotential(particleType1='O', # particleType2='Ti', # E_0=0.5478*eV, # k=1.9*Angstrom**-1, # r_0=1.96*Angstrom, # r_i=6.5*Angstrom, # r_cut=7.0*Angstrom)) potential_set.setCoulombSolver(CoulombSPME(r_cut=9.0*Angstrom)) calculator = TremoloXCalculator(parameters=potential_set) bulk_configuration.setCalculator(calculator) bulk_configuration.update() bulk_configuration = OptimizeGeometry( bulk_configuration, max_forces=0.001*eV/Ang, max_stress=0.001*eV/Ang**3, max_steps=200, max_step_length=0.2*Ang, trajectory_filename=None, optimizer_method=LBFGS(), ) nlprint(bulk_configuration) # ------------------------------------------------------------- # Elastic constants # ------------------------------------------------------------- elastic_constants = ElasticConstants( bulk_configuration, optimizer=LBFGS(), max_forces=0.001*eV/Angstrom, max_steps=200, max_step_length=0.5*Angstrom, scf_restart_step_length=0.1*Angstrom, eta_max=0.002, n_eta=3, enable_symmetry=True, fit_order=1 ) nlprint(elastic_constants)