setVerbosity(MinimalLog) # Set up lattice vector_a = [7.680028171823331, 0.0, 0.0]*Angstrom vector_b = [0.0, 7.680028171823331, 0.0]*Angstrom vector_c = [0.0, 0.0, 25.1459]*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, Hydrogen, Hydrogen, Hydrogen, Hydrogen] # Define coordinates fractional_coordinates = [[ 0.740207911277, 0.379384437585, 0.053456100246], [ 0.240207911277, 0.379384437585, 0.053456100246], [ 0.740207911277, 0.879384437585, 0.053456100246], [ 0.240207911277, 0.879384437585, 0.053456100246], [ 0.740207911277, 0.129384437585, 0.107447009301], [ 0.240207911277, 0.129384437585, 0.107447009301], [ 0.740207911277, 0.629384437585, 0.107447009301], [ 0.240207911277, 0.629384437585, 0.107447009301], [ 0.990207911277, 0.129384437585, 0.161437918356], [ 0.490207911277, 0.129384437585, 0.161437918356], [ 0.990207911277, 0.629384437585, 0.161437918356], [ 0.490207911277, 0.629384437585, 0.161437918356], [ 0.987311593647, 0.379580856217, 0.216683545941], [ 0.46906536299 , 0.379452955575, 0.210519270452], [ 0.990828621535, 0.879818769389, 0.221097303397], [ 0.501313791857, 0.879763541592, 0.216386823534], [ 0.717963374338, 0.381086870693, 0.261379244465], [ 0.239789848254, 0.380731133835, 0.272203911213], [ 0.737484438109, 0.88173757462 , 0.279891994496], [ 0.255658159092, 0.881172160909, 0.272089395915], [ 0.64737994545 , 0.155081686987, 0.322038410091], [ 0.332633178258, 0.134157771767, 0.320726001181], [ 0.646725726026, 0.60849461158 , 0.321539687743], [ 0.332046635146, 0.628081509536, 0.320510616304], [ 0.742738765525, 0.383003779559, 0.38202844159 ], [ 0.938737930305, 0.383942599828, 0.379847387695], [ 0.235546312668, 0.249513267782, 0.37718917881 ], [ 0.260710284058, 0.134421914451, 0.37650808722 ], [ 0.689906018453, 0.383970107874, 0.439716648972]] # Set up configuration bulk_configuration = BulkConfiguration( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # Add tags bulk_configuration.addTags('H_Si', [25, 26, 27, 28]) bulk_configuration.addTags('Selection 0', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]) bulk_configuration.addTags('constrained', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]) # ------------------------------------------------------------- # Calculator # ------------------------------------------------------------- #---------------------------------------- # PAW Data Set #---------------------------------------- basis_set = [ PAWPBESuggested.Silicon, PAWPBESuggested.Hydrogen, ] #---------------------------------------- # Exchange-Correlation #---------------------------------------- exchange_correlation = GGA.PBE density_mesh_cutoff = OptimizedFFTGridSampling( energy_cutoff=1520.0*eV, maximum_average_prime_factor=5.0, maximum_grid_enlargement=0.1, ) k_point_sampling = MonkhorstPackGrid( na=6, nb=6, nc=1, symmetries=[ ([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]], [ 0., 0., 0.]), ([[-1., 0., 0.], [ 0.,-1., 0.], [ 0., 0.,-1.]], [ 0., 0., 0.]), ], force_timereversal=True, shift_to_gamma=[True, True, True], ) exx_grid_cutoff = OptimizedFFTGridSampling( energy_cutoff=760.0*eV, maximum_average_prime_factor=5.0, maximum_grid_enlargement=0.1, ) exact_exchange_parameters = ExactExchangeParameters( aux_basis_tolerance=0.001, number_of_waves=1024, maximum_wave_vector=50.0, integral_tolerance=0.0001, relative_screening_tolerance=0.01, ) compensation_charge_mesh_cutoff = OptimizedFFTGridSampling( energy_cutoff=1520.0*eV, maximum_average_prime_factor=5.0, maximum_grid_enlargement=0.1, ) numerical_accuracy_parameters = NumericalAccuracyParameters( density_mesh_cutoff=density_mesh_cutoff, k_point_sampling=k_point_sampling, radial_step_size=0.001*Bohr, density_cutoff=1e-06, interaction_max_range=20.0*Angstrom, number_of_reciprocal_points=1024, reciprocal_energy_cutoff=1250.0*Hartree, bands_per_electron=1.2, occupation_method=GaussianSmearing(1000.0*Kelvin*boltzmann_constant), exx_grid_cutoff=exx_grid_cutoff, exact_exchange_parameters=exact_exchange_parameters, compensation_charge_mesh_cutoff=compensation_charge_mesh_cutoff, ) iteration_control_parameters = IterationControlParameters( tolerance=0.0001, max_steps=100, algorithm=PulayMixer(), damping_factor=0.2, number_of_history_steps=5, start_mixing_after_step=0, mixing_variable=HamiltonianVariable, preconditioner=Preconditioner.Off, linear_dependence_threshold=0.0, max_exx_updates=50, non_convergence_behavior=ContinueCalculation(), enable_scf_stop_file=True, ) poisson_solver = FastFourierSolver() initialization_method = BasisSetInitialization() density_matrix_method = PPCGSolver( absolute_tolerance=1e-08, relative_tolerance=1e-99, maximum_number_of_iterations=2, block_size=5, rr_period=3, buffer_states=0.05, initialization_method=initialization_method, ) algorithm_parameters = AlgorithmParameters( density_matrix_method=density_matrix_method, store_grids=True, store_basis_on_grid=Automatic, store_energy_density_matrix=Automatic, scf_restart_step_length=0.1*Angstrom, use_symmetries=Automatic, ) parallel_parameters = ParallelParameters( processes_per_neb_image=None, processes_per_individual=None, processes_per_bias_point=None, processes_per_saddle_search=1, ) checkpoint_handler = CheckpointHandler( time_interval=0.5*Hour, ) calculator = PlaneWaveCalculator( wave_function_cutoff=380.0*eV, basis_set=basis_set, exchange_correlation=exchange_correlation, numerical_accuracy_parameters=numerical_accuracy_parameters, iteration_control_parameters=iteration_control_parameters, poisson_solver=poisson_solver, algorithm_parameters=algorithm_parameters, parallel_parameters=parallel_parameters, checkpoint_handler=checkpoint_handler, charge=0.0, correction_extension=None, fixed_spin_moment=False, store_wave_functions=False, processes_per_kpoint=Automatic, ) bulk_configuration.setCalculator(calculator) nlprint(bulk_configuration) bulk_configuration.update() nlsave('SiH3-SiH.hdf5', bulk_configuration) # ------------------------------------------------------------- # Optimize Geometry # ------------------------------------------------------------- fix_atom_indices_0 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] constraints = [FixAtomConstraints(fix_atom_indices_0)] bulk_configuration = OptimizeGeometry( bulk_configuration, max_forces=0.02*eV/Ang, max_steps=200, max_step_length=0.2*Ang, constraints=constraints, trajectory_filename='SiH2-H2.hdf5', trajectory_interval=1, restart_strategy=RestartFromTrajectory(), disable_stress=True, optimizer_method=LBFGS(), enable_optimization_stop_file=True, ) nlsave('SiH2-H2.hdf5', bulk_configuration) nlprint(bulk_configuration) # ------------------------------------------------------------- # Total Energy # ------------------------------------------------------------- total_energy = TotalEnergy(bulk_configuration) nlsave('SiH2-H2.hdf5', total_energy) nlprint(total_energy)