# -*- coding: utf-8 -*- # ------------------------------------------------------------- # Bulk Configuration # ------------------------------------------------------------- # Set up lattice vector_a = [16.291800000000002, 0.0, 0.0]*Angstrom vector_b = [0.0, 16.291800000000002, 0.0]*Angstrom vector_c = [0.0, 0.0, 16.291800000000002]*Angstrom lattice = UnitCell(vector_a, vector_b, vector_c) # Define elements elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon] # Define coordinates fractional_coordinates = [[ 0.041666666667, 0.041666666667, 0.041666666667], [ 0.041666666667, 0.041666666667, 0.375 ], [ 0.041666666667, 0.041666666667, 0.708333333333], [ 0.041666666667, 0.375 , 0.041666666667], [ 0.041666666667, 0.375 , 0.375 ], [ 0.041666666667, 0.375 , 0.708333333333], [ 0.041666666667, 0.708333333333, 0.041666666667], [ 0.041666666667, 0.708333333333, 0.375 ], [ 0.041666666667, 0.708333333333, 0.708333333333], [ 0.375 , 0.041666666667, 0.041666666667], [ 0.375 , 0.041666666667, 0.375 ], [ 0.375 , 0.041666666667, 0.708333333333], [ 0.375 , 0.375 , 0.041666666667], [ 0.375 , 0.375 , 0.375 ], [ 0.375 , 0.375 , 0.708333333333], [ 0.375 , 0.708333333333, 0.041666666667], [ 0.375 , 0.708333333333, 0.375 ], [ 0.375 , 0.708333333333, 0.708333333333], [ 0.708333333333, 0.041666666667, 0.041666666667], [ 0.708333333333, 0.041666666667, 0.375 ], [ 0.708333333333, 0.041666666667, 0.708333333333], [ 0.708333333333, 0.375 , 0.041666666667], [ 0.708333333333, 0.375 , 0.375 ], [ 0.708333333333, 0.375 , 0.708333333333], [ 0.708333333333, 0.708333333333, 0.041666666667], [ 0.708333333333, 0.708333333333, 0.375 ], [ 0.708333333333, 0.708333333333, 0.708333333333], [ 0.125 , 0.125 , 0.125 ], [ 0.125 , 0.125 , 0.458333333333], [ 0.125 , 0.125 , 0.791666666667], [ 0.125 , 0.458333333333, 0.125 ], [ 0.125 , 0.458333333333, 0.458333333333], [ 0.125 , 0.458333333333, 0.791666666667], [ 0.125 , 0.791666666667, 0.125 ], [ 0.125 , 0.791666666667, 0.458333333333], [ 0.125 , 0.791666666667, 0.791666666667], [ 0.458333333333, 0.125 , 0.125 ], [ 0.458333333333, 0.125 , 0.458333333333], [ 0.458333333333, 0.125 , 0.791666666667], [ 0.458333333333, 0.458333333333, 0.125 ], [ 0.458333333333, 0.458333333333, 0.458333333333], [ 0.458333333333, 0.458333333333, 0.791666666667], [ 0.458333333333, 0.791666666667, 0.125 ], [ 0.458333333333, 0.791666666667, 0.458333333333], [ 0.458333333333, 0.791666666667, 0.791666666667], [ 0.791666666667, 0.125 , 0.125 ], [ 0.791666666667, 0.125 , 0.458333333333], [ 0.791666666667, 0.125 , 0.791666666667], [ 0.791666666667, 0.458333333333, 0.125 ], [ 0.791666666667, 0.458333333333, 0.458333333333], [ 0.791666666667, 0.458333333333, 0.791666666667], [ 0.791666666667, 0.791666666667, 0.125 ], [ 0.791666666667, 0.791666666667, 0.458333333333], [ 0.791666666667, 0.791666666667, 0.791666666667], [ 0.208333333333, 0.208333333333, 0.041666666667], [ 0.208333333333, 0.208333333333, 0.375 ], [ 0.208333333333, 0.208333333333, 0.708333333333], [ 0.208333333333, 0.541666666667, 0.041666666667], [ 0.208333333333, 0.541666666667, 0.375 ], [ 0.208333333333, 0.541666666667, 0.708333333333], [ 0.208333333333, 0.875 , 0.041666666667], [ 0.208333333333, 0.875 , 0.375 ], [ 0.208333333333, 0.875 , 0.708333333333], [ 0.541666666667, 0.208333333333, 0.041666666667], [ 0.541666666667, 0.208333333333, 0.375 ], [ 0.541666666667, 0.208333333333, 0.708333333333], [ 0.541666666667, 0.541666666667, 0.041666666667], [ 0.541666666667, 0.541666666667, 0.375 ], [ 0.541666666667, 0.541666666667, 0.708333333333], [ 0.541666666667, 0.875 , 0.041666666667], [ 0.541666666667, 0.875 , 0.375 ], [ 0.541666666667, 0.875 , 0.708333333333], [ 0.875 , 0.208333333333, 0.041666666667], [ 0.875 , 0.208333333333, 0.375 ], [ 0.875 , 0.208333333333, 0.708333333333], [ 0.875 , 0.541666666667, 0.041666666667], [ 0.875 , 0.541666666667, 0.375 ], [ 0.875 , 0.541666666667, 0.708333333333], [ 0.875 , 0.875 , 0.041666666667], [ 0.875 , 0.875 , 0.375 ], [ 0.875 , 0.875 , 0.708333333333], [ 0.291666666667, 0.291666666667, 0.125 ], [ 0.291666666667, 0.291666666667, 0.458333333333], [ 0.291666666667, 0.291666666667, 0.791666666667], [ 0.291666666667, 0.625 , 0.125 ], [ 0.291666666667, 0.625 , 0.458333333333], [ 0.291666666667, 0.625 , 0.791666666667], [ 0.291666666667, 0.958333333333, 0.125 ], [ 0.291666666667, 0.958333333333, 0.458333333333], [ 0.291666666667, 0.958333333333, 0.791666666667], [ 0.625 , 0.291666666667, 0.125 ], [ 0.625 , 0.291666666667, 0.458333333333], [ 0.625 , 0.291666666667, 0.791666666667], [ 0.625 , 0.625 , 0.125 ], [ 0.625 , 0.625 , 0.458333333333], [ 0.625 , 0.625 , 0.791666666667], [ 0.625 , 0.958333333333, 0.125 ], [ 0.625 , 0.958333333333, 0.458333333333], [ 0.625 , 0.958333333333, 0.791666666667], [ 0.958333333333, 0.291666666667, 0.125 ], [ 0.958333333333, 0.291666666667, 0.458333333333], [ 0.958333333333, 0.291666666667, 0.791666666667], [ 0.958333333333, 0.625 , 0.125 ], [ 0.958333333333, 0.625 , 0.458333333333], [ 0.958333333333, 0.625 , 0.791666666667], [ 0.958333333333, 0.958333333333, 0.125 ], [ 0.958333333333, 0.958333333333, 0.458333333333], [ 0.958333333333, 0.958333333333, 0.791666666667], [ 0.208333333333, 0.041666666667, 0.208333333333], [ 0.208333333333, 0.041666666667, 0.541666666667], [ 0.208333333333, 0.041666666667, 0.875 ], [ 0.208333333333, 0.375 , 0.208333333333], [ 0.208333333333, 0.375 , 0.541666666667], [ 0.208333333333, 0.375 , 0.875 ], [ 0.208333333333, 0.708333333333, 0.208333333333], [ 0.208333333333, 0.708333333333, 0.541666666667], [ 0.208333333333, 0.708333333333, 0.875 ], [ 0.541666666667, 0.041666666667, 0.208333333333], [ 0.541666666667, 0.041666666667, 0.541666666667], [ 0.541666666667, 0.041666666667, 0.875 ], [ 0.541666666667, 0.375 , 0.208333333333], [ 0.541666666667, 0.375 , 0.541666666667], [ 0.541666666667, 0.375 , 0.875 ], [ 0.541666666667, 0.708333333333, 0.208333333333], [ 0.541666666667, 0.708333333333, 0.541666666667], [ 0.541666666667, 0.708333333333, 0.875 ], [ 0.875 , 0.041666666667, 0.208333333333], [ 0.875 , 0.041666666667, 0.541666666667], [ 0.875 , 0.041666666667, 0.875 ], [ 0.875 , 0.375 , 0.208333333333], [ 0.875 , 0.375 , 0.541666666667], [ 0.875 , 0.375 , 0.875 ], [ 0.875 , 0.708333333333, 0.208333333333], [ 0.875 , 0.708333333333, 0.541666666667], [ 0.875 , 0.708333333333, 0.875 ], [ 0.291666666667, 0.125 , 0.291666666667], [ 0.291666666667, 0.125 , 0.625 ], [ 0.291666666667, 0.125 , 0.958333333333], [ 0.291666666667, 0.458333333333, 0.291666666667], [ 0.291666666667, 0.458333333333, 0.625 ], [ 0.291666666667, 0.458333333333, 0.958333333333], [ 0.291666666667, 0.791666666667, 0.291666666667], [ 0.291666666667, 0.791666666667, 0.625 ], [ 0.291666666667, 0.791666666667, 0.958333333333], [ 0.625 , 0.125 , 0.291666666667], [ 0.625 , 0.125 , 0.625 ], [ 0.625 , 0.125 , 0.958333333333], [ 0.625 , 0.458333333333, 0.291666666667], [ 0.625 , 0.458333333333, 0.625 ], [ 0.625 , 0.458333333333, 0.958333333333], [ 0.625 , 0.791666666667, 0.291666666667], [ 0.625 , 0.791666666667, 0.625 ], [ 0.625 , 0.791666666667, 0.958333333333], [ 0.958333333333, 0.125 , 0.291666666667], [ 0.958333333333, 0.125 , 0.625 ], [ 0.958333333333, 0.125 , 0.958333333333], [ 0.958333333333, 0.458333333333, 0.291666666667], [ 0.958333333333, 0.458333333333, 0.625 ], [ 0.958333333333, 0.458333333333, 0.958333333333], [ 0.958333333333, 0.791666666667, 0.291666666667], [ 0.958333333333, 0.791666666667, 0.625 ], [ 0.958333333333, 0.791666666667, 0.958333333333], [ 0.041666666667, 0.208333333333, 0.208333333333], [ 0.041666666667, 0.208333333333, 0.541666666667], [ 0.041666666667, 0.208333333333, 0.875 ], [ 0.041666666667, 0.541666666667, 0.208333333333], [ 0.041666666667, 0.541666666667, 0.541666666667], [ 0.041666666667, 0.541666666667, 0.875 ], [ 0.041666666667, 0.875 , 0.208333333333], [ 0.041666666667, 0.875 , 0.541666666667], [ 0.041666666667, 0.875 , 0.875 ], [ 0.375 , 0.208333333333, 0.208333333333], [ 0.375 , 0.208333333333, 0.541666666667], [ 0.375 , 0.208333333333, 0.875 ], [ 0.375 , 0.541666666667, 0.208333333333], [ 0.375 , 0.541666666667, 0.541666666667], [ 0.375 , 0.541666666667, 0.875 ], [ 0.375 , 0.875 , 0.208333333333], [ 0.375 , 0.875 , 0.541666666667], [ 0.375 , 0.875 , 0.875 ], [ 0.708333333333, 0.208333333333, 0.208333333333], [ 0.708333333333, 0.208333333333, 0.541666666667], [ 0.708333333333, 0.208333333333, 0.875 ], [ 0.708333333333, 0.541666666667, 0.208333333333], [ 0.708333333333, 0.541666666667, 0.541666666667], [ 0.708333333333, 0.541666666667, 0.875 ], [ 0.708333333333, 0.875 , 0.208333333333], [ 0.708333333333, 0.875 , 0.541666666667], [ 0.708333333333, 0.875 , 0.875 ], [ 0.125 , 0.291666666667, 0.291666666667], [ 0.125 , 0.291666666667, 0.625 ], [ 0.125 , 0.291666666667, 0.958333333333], [ 0.125 , 0.625 , 0.291666666667], [ 0.125 , 0.625 , 0.625 ], [ 0.125 , 0.625 , 0.958333333333], [ 0.125 , 0.958333333333, 0.291666666667], [ 0.125 , 0.958333333333, 0.625 ], [ 0.125 , 0.958333333333, 0.958333333333], [ 0.458333333333, 0.291666666667, 0.291666666667], [ 0.458333333333, 0.291666666667, 0.625 ], [ 0.458333333333, 0.291666666667, 0.958333333333], [ 0.458333333333, 0.625 , 0.291666666667], [ 0.458333333333, 0.625 , 0.625 ], [ 0.458333333333, 0.625 , 0.958333333333], [ 0.458333333333, 0.958333333333, 0.291666666667], [ 0.458333333333, 0.958333333333, 0.625 ], [ 0.458333333333, 0.958333333333, 0.958333333333], [ 0.791666666667, 0.291666666667, 0.291666666667], [ 0.791666666667, 0.291666666667, 0.625 ], [ 0.791666666667, 0.291666666667, 0.958333333333], [ 0.791666666667, 0.625 , 0.291666666667], [ 0.791666666667, 0.625 , 0.625 ], [ 0.791666666667, 0.625 , 0.958333333333], [ 0.791666666667, 0.958333333333, 0.291666666667], [ 0.791666666667, 0.958333333333, 0.625 ], [ 0.791666666667, 0.958333333333, 0.958333333333]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- potentialSet = Tersoff_Si_1988() calculator = TremoloXCalculator(parameters=potentialSet) calculator.setVerletListsDelta(0.25*Angstrom) bulk_configuration.setCalculator(calculator) nlprint(bulk_configuration) bulk_configuration.update() nlsave('si_bulk_specific_heat.hdf5', bulk_configuration) # ------------------------------------------------------------- # Molecular Dynamics # ------------------------------------------------------------- initial_velocity = MaxwellBoltzmannDistribution( temperature=300.0*Kelvin, remove_center_of_mass_momentum=True ) method = NVTNoseHoover( time_step=1*femtoSecond, reservoir_temperature=300*Kelvin, thermostat_timescale=100*femtoSecond, heating_rate=0*Kelvin/picoSecond, chain_length=3, initial_velocity=initial_velocity, ) constraints = [FixCenterOfMass()] md_trajectory = MolecularDynamics( bulk_configuration, constraints=constraints, trajectory_filename='si_bulk_specific_heat.hdf5', steps=5000, log_interval=1, method=method ) bulk_configuration = md_trajectory.lastImage() # Calculate specific heat capacity. specific_heat_capacity = SpecificHeatCapacity( md_trajectory, min_temperature=0.0*Kelvin, max_temperature=1000.0*Kelvin, ) # Plot specific heat capacity. import pylab temperatures = specific_heat_capacity.temperatures().inUnitsOf(Kelvin) specific_heats = specific_heat_capacity.specificHeatCapacity().inUnitsOf(Joule/gram/Kelvin) pylab.plot(temperatures, specific_heats) pylab.xlabel('Temperature (K)') pylab.ylabel('Specific Heat Capacity (J/g/K)') pylab.savefig('si_bulk_specific_heat.png')