%matplotlib inline
import json
import numpy as np
import operator
import os
import pandas
import pylab as plt
import random
from ase.data import reference_states, atomic_numbers
from sklearn.neighbors import KernelDensity
from pyiron_atomistics import Project
import matplotlib as mpl
mpl.rc('font', family='Times New Roman')
plt.rcParams["mathtext.fontset"] = "stix"
from pyiron_atomistics.thermodynamics.interfacemethod import freeze_one_half, remove_selective_dynamics, set_server, create_job_template, fix_iso, fix_z_dir, half_velocity, minimize_pos, minimize_vol, next_calc, npt_solid, npt_liquid, next_step_funct, round_temperature_next, strain_circle, analyse_minimized_structure, get_press, get_center_point, get_strain_lst, get_nve_job_name, plot_solid_liquid_ratio, ratio_selection, plot_equilibration, plot_melting_point_prediction, calc_temp_iteration, get_initial_melting_temperature_guess, validate_convergence, initialise_iterators, get_voronoi_volume, check_for_holes
pr = Project('melting')
input_file = os.path.join(pr.path, 'input.json')
output_file = os.path.join(pr.path, 'output.json')
pr.job_table()
if os.path.exists(input_file) and os.stat(input_file).st_size != 0:
with open(input_file, 'r') as f:
input_dict = json.load(f)
else:
input_dict = {
"config": [
"pair_style morse 9.97749\n",
"pair_coeff * * 0.4174 1.3885 2.845\n"
],
"species": ["Al"],
"element": "Al"
}
input_dict
pot_dict = input_dict.copy()
if 'model' not in pot_dict.keys():
pot_dict['model'] = 'Lammps'
if 'name' not in pot_dict.keys():
pot_dict['name'] = 'CustomPotential'
if 'filename' not in pot_dict.keys():
pot_path = []
else:
pot_path = [os.path.abspath(pot_dict['filename'])]
potential = pandas.DataFrame({'Config': [pot_dict['config']],
'Filename': [pot_path],
'Model': [pot_dict['model']],
'Name': [pot_dict['name']],
'Species': [pot_dict['species']]
})
project_parameter = {
'project': pr,
'run_time_steps': 50000,
'nvt_run_time_steps': 10000,
'nve_run_time_steps': 10000,
'temperature_left': 0,
'temperature_right': 1000,
'strain_run_time_steps': 1000,
'convergence_criterion': 1,
'potential': potential,
'cpu_cores': 1,
'job_type': pr.job_type.Lammps,
'enable_h5md': False,
'points': 21,
'boundary_value': 0.25,
'ratio_boundary': 0.25,
'timestep_lst': [2, 2, 1],
'fit_range_lst': [0.05, 0.01, 0.01],
'nve_run_time_steps_lst': [25000, 20000, 50000],
'number_of_atoms': 8000,
}
for k in input_dict.keys():
project_parameter[k] = input_dict[k]
if 'crystalstructure' not in project_parameter.keys():
project_parameter['crystalstructure'] = reference_states[atomic_numbers[project_parameter['element']]]['symmetry']
if 'seed' not in project_parameter.keys():
project_parameter['seed'] = random.randint(0,99999)
project_parameter['seed']
# Values from a previous calculation can be inserted here to reproduce the results
step_dict = {}
if os.path.exists(output_file):
with open(output_file, 'r') as f:
step_dict_str = json.load(f)
for k,v in step_dict_str.items():
step_dict[int(k)] = v
step_dict
From here on the notebook is automated - no change required !#
step_count = 0
temperature_next = None
enable_iteration = True
convergence_goal_achieved = False
pr = project_parameter['project']
if 'lattice_constant' in project_parameter.keys():
a = project_parameter['lattice_constant']
else:
a = None
if project_parameter['crystalstructure'] == 'hcp':
basis = pr.create_ase_bulk(name=project_parameter['element'], crystalstructure=project_parameter['crystalstructure'].lower(), a=a, orthorhombic=True)
else:
basis = pr.create_ase_bulk(name=project_parameter['element'], crystalstructure=project_parameter['crystalstructure'].lower(), a=a, cubic=True)
basis_lst = [basis.repeat([i, i, i]) for i in range(5,30)]
basis = basis_lst[np.argmin([np.abs(len(b)-project_parameter['number_of_atoms']/2) for b in basis_lst])]
basis.analyse_ovito_cna_adaptive()
basis.plot3d()
# pr.remove_jobs(recursive=True)
timestep_iter, fit_range_iter, nve_run_time_steps_iter = initialise_iterators(project_parameter)
timestep = next(timestep_iter)
fit_range = next(fit_range_iter)
nve_run_time_steps = next(nve_run_time_steps_iter)
pr.job_table()
Step 1: set up the solid sample and roughly estimate a melting point#
ham_minimize_pos= minimize_pos(
structure=basis,
project_parameter=project_parameter
)
ham_minimize_vol = minimize_vol(
structure=ham_minimize_pos.get_structure(),
project_parameter=project_parameter
)
temperature_next, structure_left = get_initial_melting_temperature_guess(
project_parameter=project_parameter,
ham_minimize_vol=ham_minimize_vol,
temperature_next=temperature_next
)
temperature_next # +/- 100K
if step_count in step_dict.keys():
temperature_next = step_dict[step_count]['temperature_next']
else:
step_dict[step_count]= {'temperature_next': temperature_next}
step_dict
Step 2: set up the interface structure#
temperature_next, temperature_mean, temperature_left, temperature_right, strain_result_lst, pressure_result_lst = calc_temp_iteration(
basis=ham_minimize_vol.get_structure().repeat([1,1,2]),
temperature_next=temperature_next,
project_parameter=project_parameter,
timestep = timestep,
nve_run_time_steps=nve_run_time_steps,
fit_range=fit_range,
center=None,
debug_plot=True)
temperature_next, temperature_mean, temperature_left, temperature_right
Step 3: run NVT and NVE MD#
output = validate_convergence(
pr=pr,
temperature_left=temperature_left,
temperature_next=temperature_next,
temperature_right=temperature_right,
enable_iteration=enable_iteration,
timestep_iter=timestep_iter,
timestep_lst=project_parameter['timestep_lst'],
timestep=timestep,
fit_range_iter=fit_range_iter,
fit_range_lst=project_parameter['fit_range_lst'],
fit_range=fit_range,
nve_run_time_steps_iter=nve_run_time_steps_iter,
nve_run_time_steps_lst=project_parameter['nve_run_time_steps_lst'],
nve_run_time_steps=nve_run_time_steps,
strain_result_lst=strain_result_lst,
pressure_result_lst=pressure_result_lst,
step_count=step_count,
step_dict=step_dict,
boundary_value=project_parameter['boundary_value'],
ratio_boundary=project_parameter['ratio_boundary'],
convergence_goal=project_parameter['convergence_criterion'],
output_file=output_file
)
convergence_goal_achieved, enable_iteration, step_count, step_dict, timestep, fit_range, nve_run_time_steps, boundary_value, ratio_boundary, temperature_next, center = output
step_dict[step_count]
convergence_goal_achieved
Step 4: full cycle, predict the final melting point#
temperature_estimate_lst, temperature_calculated_lst = [], []
while len(temperature_calculated_lst) < 3 or not convergence_goal_achieved and len(temperature_calculated_lst) < 10:
temperature_estimate_lst.append(temperature_next)
temperature_next, temperature_mean, temperature_left, temperature_right, strain_result_lst, pressure_result_lst = calc_temp_iteration(
basis=ham_minimize_vol.get_structure().repeat([1,1,2]),
temperature_next=temperature_next,
project_parameter=project_parameter,
timestep=timestep,
nve_run_time_steps=nve_run_time_steps,
fit_range=fit_range,
center=center,
debug_plot=True)
print(temperature_next, temperature_mean, temperature_left, temperature_right)
output = validate_convergence(pr=pr,
temperature_left=temperature_left,
temperature_next=temperature_next,
temperature_right=temperature_right,
enable_iteration=enable_iteration,
timestep_iter=timestep_iter,
timestep_lst=project_parameter['timestep_lst'],
timestep=timestep,
fit_range_iter=fit_range_iter,
fit_range_lst=project_parameter['fit_range_lst'],
fit_range=fit_range,
nve_run_time_steps_iter=nve_run_time_steps_iter,
nve_run_time_steps_lst=project_parameter['nve_run_time_steps_lst'],
nve_run_time_steps=nve_run_time_steps,
strain_result_lst=strain_result_lst,
pressure_result_lst=pressure_result_lst,
step_count=step_count,
step_dict=step_dict,
boundary_value=project_parameter['boundary_value'],
ratio_boundary=project_parameter['ratio_boundary'],
convergence_goal=project_parameter['convergence_criterion'],
output_file=output_file)
convergence_goal_achieved, enable_iteration, step_count, step_dict, timestep, fit_range, nve_run_time_steps, boundary_value, ratio_boundary, temperature_next, center = output
print(step_dict[step_count])
temperature_calculated_lst.append(temperature_next)
if not os.path.exists(output_file):
with open(output_file, 'w') as f:
json.dump(step_dict, f)
step_dict
#plot the convergence of loop calculations
first_key = list(step_dict.keys())[-1]
second_key = 'temperature_next'
allt = [step_dict[key][second_key] for key in list(step_dict.keys())]
plt.plot(np.arange(1, len(allt)), allt[0:-1], 'ro-', label=r"$T^e$")
plt.plot(np.arange(1, len(allt)), allt[1:], 'bo-', label=r"$T^p$")
plt.legend(fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.xlabel('Number of loops', fontsize=14)
plt.ylabel('Temperature (K)', fontsize=14)
plt.show()
df = pr.job_table()
df