Initial commit

checklist.txt

test.py

ase_slab.png

test_lammps_file

crystal_structure*

config_out

*.traj

{directory}*

data_CNT_sim/*

data_CNT_sim/nve_deform/*

data_min_energy/*

*.ss

*.csv

!data_min_energy/

import numpy as np

import matplotlib.pyplot as plt

import json

from scipy.optimize import curve_fit, root_scalar

#from ase.calculators.emt import EMT

from asap3 import EMT, EMTRasmussenParameters, EMT2013

from ase.lattice.cubic import FaceCenteredCubic, BodyCenteredCubic, SimpleCubic

from ase.io import Trajectory, read

from ase.optimize import FIRE

import csv 

from tqdm import tqdm

import os

#print(cell.cell[:]) # cell.cell only gives a 1-dim array type output :/ ?

strain_high = 0.001

strain_low = -0.001

directory = 'data_elastic_constants/'

#this object stores a strain tensor for the 3 different types we investigate

class StrainTensor:

    def __init__(self):

        self.tensor = np.zeros((3,3))

    def make_tensor(self, d_type, eps):

        if d_type == 'uniax':

            self.uniax(eps)

        elif d_type == 'biax':

            self.biax(eps)

        elif d_type == 'shear':

            self.shear(eps)

        else:

             raise ValueError(f"No deform type {d_type}")

    def uniax(self, eps):

        self.tensor[0,0] = eps

    def biax(self, eps):

        self.tensor[0,0] = eps

        self.tensor[1,1] = eps

    def shear(self, eps):

        self.tensor[0,1] = eps

        self.tensor[1,0] = eps

def deform_cell(cell_array, strain_tensor):

    # don't use * for multiplying matrices

    atoms_positions_new = np.dot(atoms_positions,(np.identity(3)+strain_tensor))

    cell_array_new = np.dot(cell_array,(np.identity(3)+strain_tensor))

    cell.cell[:] = cell_array_new

    cell.positions = atoms_positions_new

    return cell

def fit_func(eps, p_0, p_1, p_2, p_3):

        return p_0 + p_1*eps + p_2*eps**2 + p_3*eps**3

def derivative_fit_func(eps, p_1, p_2, p_3):

    return p_1 + 2*p_2*eps + 3*p_3*eps**2

def second_derivative_fit_func(eps,p_2,p_3):

    return 2*p_2 + 6*p_3*eps

def write_to_csv(d_type, strains, energy_densities):

    os.makedirs(os.path.dirname(directory), exist_ok=True)

    with open(f'{directory}{d_type}_strain_energydensity.csv', 'w', newline='') as new_csv:

        newarray = csv.writer(new_csv, delimiter=",")

        newarray.writerow(["strain", "energy density in eV/Angstrom^3"])

        data =  []

        for strain, energy_density in zip(strains, energy_densities):

            data.append([strain, energy_density])

        newarray.writerows(data)

    print(f"Wrote strain energy dentity data to {directory}{d_type}_strain_energydensity.csv")

if __name__ == '__main__':

    os.makedirs(os.path.dirname(directory), exist_ok=True)

    #relax initial structur

    lattice_const = float(input("Lattice constant for fcc: "))

    #with open("data_min_energy/lattice_consts", "r") as file: #load lattice constant

    #    lattice_const_dict = json.load(file) 

    #lattice_const = lattice_const_dict['W_2nd_ord']['fcc']

    strain_high = float(input("Enter value for highest strain (e.g. 0.02): ")) 

    strain_low = float(input("Enter value for lowest (negative) strain (e.g. -0.02)")) 

    N = int(input("How many data points should be evaluated (e.g. 50): "))

    cell = FaceCenteredCubic('Cu', latticeconstant=lattice_const, size=(3,3,3))

    cell.calc = EMT()

    # cell.calc = EMT(EMTRasmussenParameters()) in case other parameters need to be used

    #relax cell  

    traj_relaxed_cell = Trajectory(f'{directory}relaxed_cell.log', 'w')

    relax_cell =  FIRE(cell, trajectory=f'{directory}relaxed_cell.traj')

    relax_cell.run(fmax=10**-8)

    traj_relaxed_cell.write(cell)

    traj_relax_read = Trajectory(f'{directory}relaxed_cell.traj', 'r')

    cell_relaxed = traj_relax_read[-1]

    #get initial state of relaxed cell to work with

    cell_array = cell_relaxed.cell[:]

    atoms_positions = cell_relaxed.positions

    volume = cell_relaxed.get_volume()

    epsilons = np.linspace(strain_low, strain_high, N)

    deform_types = ['uniax', 'biax', 'shear']

    #deform the cell 

    for d_type in deform_types:

        cell = read(f'{directory}relaxed_cell.traj@0', 'r')

        #traj = Trajectory(f'{directory}{d_type}.traj', 'r')

        #cell = traj[0]

        tensor = StrainTensor()

        traj = Trajectory(f'{directory}{d_type}.traj', 'w')

        print(f'\n\n\n Calculating potential energy for {d_type}-deformation.')

        print(f'Calculating {N} for strains from {strain_low} to {strain_high}')

        for eps in epsilons:

            tensor.make_tensor(d_type, eps)

            cell = deform_cell(cell_array, tensor.tensor)

            relax_cell.run(fmax=10**-8)

            cell.calc = EMT()

            cell.get_potential_energy()

            traj.write(cell)

   # fig1, ax1 = plt.subplots()

   # fig2, (ax2, ax3) = plt.subplots(1,2)

   # ax1.set_xlabel('Strain $\epsilon$')

   # ax1.set_ylabel('Strain Energy density eV/$Angstrom^{3}$')

   # ax2.set_ylabel('$dW/d\epsilon$ eV/$Angstrom^{3}$')

   # ax3.set_ylabel('$d^{2}W/d\epsilon^{2}$ eV')

   # ax2.set_xlabel('Strain $\epsilon$')

   # ax3.set_xlabel('Strain $\epsilon$')

   # elastic_const_dict = {} #the calculated elastic constants go here

   # min_eps_dict = {} #here we save the root of the first derivative -> the epsilon at which the energy is minimal 

    #fit and find min

    for d_type in deform_types:

        configs = read(f'{directory}{d_type}.traj@0:', 'r')

        energies = [config.get_potential_energy() for config in configs]

        #[config.get_volume() for config in configs]

        energy_per_volumes = []

        for energy in energies:

            energy_per_volumes.append(energy/volume)

        write_to_csv(d_type, epsilons, energy_per_volumes) 

# the following lines include all the imports we need

import numpy as np

import matplotlib.pyplot as plt 

import json

import csv

from scipy.optimize import curve_fit, root_scalar

from ase.calculators.emt import EMT

#from asap3 import EMT

from ase.lattice.cubic import FaceCenteredCubic, BodyCenteredCubic, SimpleCubic

from ase.io import Trajectory, read

import os

from tqdm import tqdm

'''

List of crystal lattice types

sc: simple cubic

fcc: face centered cubic

bcc: body centered cubic

'''

lattice_types = ['sc', 'fcc', 'bcc']

#dictionary of initial values for lattice constants (all in Angstrom)

#in task 3.1 you estimated values for the lattice constants, insert in the dictionary the below

L_init_values = {'sc' : 2.5, 

                 'fcc' : 3.6,

                 'bcc': 2.7}

def create_cubic_lattice(lattice_type, L):

    switch_dict = {'sc' : SimpleCubic('Cu', latticeconstant=L, size=(1,1,1)),

                   'bcc' : BodyCenteredCubic('Cu', latticeconstant=L, size=(1,1,1)),

                   'fcc' : FaceCenteredCubic('Cu', latticeconstant=L, size=(1,1,1))}  

    return switch_dict[lattice_type]

def calculate_latticeconst_epot_traj(lattice_type, L_init):

    os.makedirs(os.path.dirname('data_min_energy/'), exist_ok=True)

    traj_file = Trajectory(f'data_min_energy/{lattice_type}.traj', 'w')

    lattice_constants_list = []

    L_high = L_init + 0.3 

    L_low = L_init - 0.3 

    N = 300

    iterator_values = np.linspace(L_low, L_high, N)

    print(f'\n\n\nLattice type : {lattice_type}')

    print(f"Calculating {N} data points for lattice constants {L_low:.2f} to {L_high:.2f}")

    for L in tqdm(iterator_values):

        lattice_constants_list.append(L) # append current lattice constant

        data = create_cubic_lattice(lattice_type, L)

        data.calc = EMT()

        data.get_potential_energy()

        traj_file.write(data)

    return lattice_constants_list

def write_data_dict(lattice_type ,lattice_constants_list):

    data_dict = {}

    L_init = L_init_values[lattice_type]

    os.makedirs(os.path.dirname('data_min_energy/'), exist_ok=True)

    configs = read(f'data_min_energy/{lattice_type}.traj@0:', 'r') #@0:  go though all the saved configurations

    epots_list = [data.get_potential_energy() for data in configs] #list potential energies

    data_dict[lattice_type] = [lattice_constants_list, epots_list] #array pairing epot w/ corresponding latticeconst

    return data_dict[lattice_type]

#save data from dicts to permanent file

def save_data_to_json(data_dict, filename):

    with open(filename, 'w') as data:

        json.dump(data_dict, data)

def write_to_csv(lattice_type, strains, energy_densities):

    os.makedirs(os.path.dirname('/data_min_energy'), exist_ok=True)

    with open(f'data_min_energy/{lattice_type}_latticeconst_energydensity.csv', 'w', newline='') as new_csv:

        newarray = csv.writer(new_csv, delimiter=",")

        newarray.writerow(["lattice constant in angstrom", "energy density in eV/Angstrom^3"])

        data =  []

        for strain, energy_density in zip(strains, energy_densities):

            data.append([strain, energy_density])

        newarray.writerows(data)

def fit_L_min_epot(func, dfunc, data_dict):

    fitting_params = {}

    func_name = func.__name__

    for lattice_type in lattice_types:

        x,y = data_dict[lattice_type]#[:60][:60] could adjust over which values to fit

        popt, pcov = curve_fit(func, x, y)

        ax1.plot(x, [func(x_i, *popt) for x_i in x],

                '--', label = f'{func_name} {lattice_type} lattice')

        fitting_params[lattice_type] = popt

    ax1.legend()

    fig1.savefig(f'data_min_energy/{func_name}_epot_lattice_plot')

    roots_dict = {}

    for lattice_type in  lattice_types:

        p_array = fitting_params[lattice_type][1:] # not from 0 since p_0 sits at 0

        L_init = L_init_values[lattice_type]

        dfunc_fit = lambda x: dfunc(x, p_array)

        roots_dict[lattice_type] = root_scalar(dfunc_fit, x0=L_init).root #finds min

    return roots_dict, func_name 

if __name__ =='__main__':

    print("Give initial values for lattice constants in Angstom (e.g. 1.9) \n")

    L_init_values["sc"] = float(input("Simple Cubic: "))

    L_init_values["fcc"] = float(input("Face Centered Cubic: "))

    L_init_values["bcc"] = float(input("Body Centered Cubic: "))

    data_dict = {}

    for lattice_type in lattice_types:

        L_init = L_init_values[lattice_type]

        lattice_constants_list = calculate_latticeconst_epot_traj(lattice_type, L_init)

        data_dict[lattice_type] = write_data_dict(lattice_type, lattice_constants_list)

    os.makedirs(os.path.dirname('/data_min_energy'), exist_ok=True)

    save_data_to_json(data_dict, 'data_min_energy/epot_latticeconst')

    with open('data_min_energy/epot_latticeconst', 'r') as data:

       data_dict = json.load(data)

    print("\n\n\n")

    for lattice_type in lattice_types:

        strains, energy_densities = data_dict[lattice_type]

        write_to_csv(lattice_type, strains, energy_densities)

    #plot the data

        print(f"Strain vs energy density written to csv-files to file 'data_min_energy/{lattice_type}_latticeconst_energydensity.csv'")

 No newline at end of file

# MD for Nanotech quick 'n clean

## Clone this git repo    

To work with this code make sure to clone this exact git repo:    

`git clone <url of repo>`

## Using Python Packages

Packages include functions, objects, etc. that are not natively included in python.    

In the course of this practical we will use a number of packages i.e. `ase` for the actual simulation, `numpy` for matrix calculations, `scipy` for curve fitting. To use functions and objects we need to import them:    

`import <insert_package_name_here> as <name_of_alias>`     

e.g `import numpy as np`    

if we want to use just one particular function we can use:    

`from <package> import <function_or_object_from_package>`   

e.g. `from ase.optimize import FIRE`    

## Setting up the environment    

Before any of the scripts will run, make sure to initialize a virtual environment with all the required packages.     

To do this open the directory of this code in a terminal. If you cloned it to your home directory you should acces it via:    

`cd <name of clone dir>`     

Then set up the environment here (this one sets up an environment called env).         

`python -m venv env`    

Now activate the environment:     

`source env/bin/activate`    

After this your command line should look like thise `(env) <yourusername>@<nameofyourmachine>:/`     

To now include all the packages for the python scripts run:    

`pip install -r requirements.txt`    

## Running the code for determining the lattice constants     

Run the script:    

`python3 minimize_energy_script.py`     

The script will ask you to provide initial guesses for the lattice constants and will calculate 300 data points in a range of -0.3 Angstrom and +0.3 Angstrom of your initial guess.      

The strain energydentsity csv files are written to the folder `data_min_energy/`     

Plot the data and find the minima.

## Running the code for calculating the elastic constants   

Run:    

`python3 elastic_constants_script.py`    

The script will ask you to provide a lattice constant (for the fcc system), a highest and lowest strain (also negative strain i.e. compression) and how many datapoints you would like to evaluate.    

The data is written to `data_elastic_constants`. Plot the data and fit a 2nd order polynomial to calccalcualte the elastic constants.