Commit cd89ac5d authored by Christian Greff's avatar Christian Greff
Browse files

Added scripts to produce data for the carbon nanotube simulations.

parent c499280b

deform_cnt_npt.py

0 → 100644
+123 −0
Original line number Diff line number Diff line 


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 BrennerPotential

from ase.build import nanotube

from ase.io import Trajectory, read, write, gromacs

from ase.optimize import FIRE

from ase.md.langevin import Langevin

from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary

from ase import units

import os

from pathlib import Path

import pandas as pd

#from tqdm import tqdm



plot_dir = 'data_CNT_sim/nvt_plots/'

directory = 'data_CNT_sim/'

data_directory = 'data_CNT_sim/nvt_deform/'

cnt = nanotube(5, 5, symbol='C', length=3, bond=1.42)

cnt.calc = BrennerPotential()

#adjust cell and pbc (for visualization)

cnt.cell[:] = cnt.cell[:] + np.array([[10, 0, 0],[0, 10, 0],[0, 0, 0]])

cnt.set_pbc([0,0,1])



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

start_traj.write(cnt)

gromacs.write_gromacs(f'{directory}cnt_initial.dump', cnt) # dump to lammps dump so ovito can show it



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



def deform_in_z(cnt, eps):

    cnt_array = cnt.cell[:]

    cnt_array_new = cnt_array + np.dot(cnt_array_original, (np.array([[0,0,0],[0,0,0],[0,0,eps]])))

    cnt.set_cell(cnt_array_new, scale_atoms=True)



def printenergy(atoms=cnt):

     epot = atoms.get_potential_energy()

     ekin = atoms.get_kinetic_energy()     

     temp = atoms.get_temperature()

     print(f'epot: {epot}, ekin: {ekin}, temperature: {temp}')



def run_simulation(loop_number, step_number, timestep, eps, start_deform=0):

    MaxwellBoltzmannDistribution(cnt, temperature_K=600, force_temp=True)

    #dyn.run(100)

    #Stationary(cnt)

    #eqiliberate befor running

    for i in range(loop_number):

            loop_time = i*step_number

            for step in range(start_deform):

                 dyn.run(1)

                 if (loop_time+step)%1 == 0:

                     data['time'].append((loop_time+step)*timestep)

                     data['temperature'].append(cnt.get_temperature())

                     data['epot'].append(cnt.get_potential_energy())

                     data['ekin'].append(cnt.get_kinetic_energy())

                     data['etot'].append(cnt.get_total_energy())

            deform_in_z(cnt, eps)

            #dyn.attach(printenergy)

            for step in range(start_deform, step_number): 

                Stationary(cnt)

                dyn.run(1)

                if (loop_time+step)%1 == 0:

                    data['time'].append((loop_time+step)*timestep)

                    data['temperature'].append(cnt.get_temperature())

                    data['epot'].append(cnt.get_potential_energy())

                    data['ekin'].append(cnt.get_kinetic_energy())

                    data['etot'].append(cnt.get_total_energy())

                    printenergy(cnt)

                    if i%1 == 0:

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

                        gromacs.write_gromacs(f'{data_directory}cnt_{i}.gro', cnt) # dump to lammps dump so ovito can show it

                

            print(f'Rep {i+1} of {loop_number}')

            if i%10 == 0:

                gromacs.write_gromacs(f'{data_directory}cnt_{i}.dump', cnt) # dump to lammps dump so ovito can show it

    with open(f'{data_directory}data.json', 'w') as fp:

        json.dump(data, fp)



if __name__ == '__main__':

    for f in Path(f'{data_directory}').glob('*.dump' and '*.gro'): #clean up output path for dump files

        try:

              f.unlink()

        except:

             print('no files in data dir')



    strain_rate = 3*10**12 #s^-1 

    #strain_rate = eps/(timestep * step_number*10**-15)

    start_deform = 0

    timestep = 0.1 #in fs

    step_number = 30

    loop_number = 100

    eps = strain_rate * timestep * (step_number-start_deform)*10**-15

    #eps = 0.001

    total_strain = eps*loop_number



    # to save time we get the relaxed stucture from the nve sim

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

    cnt = traj_relax_read[-1]

    cnt.calc = BrennerPotential() #when reading from trajectory the potential gets lost... why? no idea 



    dyn = Langevin(cnt, timestep=timestep*units.fs,temperature_K=300, friction=0.1/units.fs)

    data = {'time':[],'temperature':[], 'epot':[],'ekin':[], 'etot': []}

    cnt_array_original = cnt.cell[:]

    run_simulation(loop_number, step_number, timestep, eps, start_deform) 

    with open(f'{data_directory}data.json', 'r') as fp:

        data = json.load(fp)

    df = pd.DataFrame.from_dict(data)

    df.to_csv(f'{data_directory}cnt_npt_data.csv')

    print(f'Saving data to {data_directory}cnt_npt_data.csv')

    print('done')

deform_cnt_nve.py

0 → 100644
+132 −0
Original line number Diff line number Diff line 


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 BrennerPotential

from ase.build import nanotube

from ase.io import Trajectory, read, write, gromacs, extxyz

from ase.optimize import FIRE

from ase.md.verlet import VelocityVerlet

from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary

from ase import units

import os

from pathlib import Path

import pandas as pd

#from tqdm import tqdm



plot_dir = 'data_CNT_sim/nve_plots/'

directory = 'data_CNT_sim/'

data_directory = 'data_CNT_sim/nve_deform/'

cnt = nanotube(5, 5, symbol='C', length=3, bond=1.42)

cnt.calc = BrennerPotential()

#adjust cell and pbc (for visualization)

cnt.cell[:] = cnt.cell[:] + np.array([[10, 0, 0],[0, 10, 0],[0, 0, 0]])

cnt.set_pbc([0,0,1])



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

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

start_traj.write(cnt)

gromacs.write_gromacs(f'{directory}cnt_initial.gro', cnt) #.gro to lammps.gro so ovito can show it



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



def relax_structure(cell):

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

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

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



def deform_in_z(cnt, eps):

    cnt_rel_array = cnt.cell[:]

    cnt_array_new = cnt_rel_array + np.dot(cnt_array_original, (np.array([[0,0,0],[0,0,0],[0,0,eps]])))

    cnt_rel.set_cell(cnt_array_new, scale_atoms=True)



def printenergy(atoms=cnt):

     epot = atoms.get_potential_energy()

     ekin = atoms.get_kinetic_energy()     

     temp = atoms.get_temperature()

     print(f'epot: {epot}, ekin: {ekin}, temperature: {temp}')



def run_simulation(loop_number, step_number, timestep, eps,filename, start_deform=0):

    #Stationary(cnt)

    #dyn.run(500) #run before simulation (w/o recording data)

    #eqiliberate befor running

    for i in range(loop_number):

    

        loop_time = i*step_number

        for step in range(start_deform): #for if we want to save data before deforming

                dyn.run(1)

                if (loop_time+step)%1 == 0:

                    data['time'].append((loop_time+step)*timestep)

                    data['temperature'].append(cnt.get_temperature())

                    data['epot'].append(cnt_rel.get_potential_energy())

                    data['ekin'].append(cnt_rel.get_kinetic_energy())

                    data['etot'].append(cnt_rel.get_total_energy())

        deform_in_z(cnt_rel, eps)

        #dyn.attach(printenergy)

        for step in range(start_deform, step_number): 

            #Stationary(cnt)

            dyn.run(1)

            if (loop_time+step)%1 == 0:

                data['time'].append((loop_time+step)*timestep)

                data['temperature'].append(cnt_rel.get_temperature())

                data['epot'].append(cnt_rel.get_potential_energy())

                data['ekin'].append(cnt_rel.get_kinetic_energy())

                data['etot'].append(cnt_rel.get_total_energy())

                printenergy(cnt_rel)

                if i%1 == 0:

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

                    gromacs.write_gromacs(f'{data_directory}cnt_{i}.gro', cnt_rel) # dump to lammps.gro so ovito can show it

            

        print(f'Rep {i+1} of {loop_number}')

        if i%10 == 0:

            gromacs.write_gromacs(f'{data_directory}cnt_{i}.gro', cnt_rel) # dump to lammps.gro so ovito can show it

    with open(f'{data_directory}{filename}_data.json', 'w') as fp:

        json.dump(data, fp)



if __name__ == '__main__':

    for f in Path(f'{data_directory}').glob('*.dump' and '*.gro' and '*.xyz'): #clean up output path for dump files

        try:

              f.unlink()

        except:

             print('no files in data dir')



    #adjustable variables 

    strain_rate = 2*10**12#s^-1 

    #strain_rate = eps/(timestep * step_number*10**-15)

    start_deform = 0

    timestep = 0.1 #in fs

    step_number = 20

    loop_number = 200

    #implicite variables

    eps = strain_rate * timestep * (step_number-start_deform)*10**-15

    #eps = 0.001

    total_strain = eps*loop_number

    filename = 'sr{sr:2.2e}_etot{etot:2.2e}_steps{steps}_loops{loops}'.format(sr=strain_rate, etot=total_strain,

                                                                                  steps=step_number, loops=loop_number)

    relax_structure(cnt) 

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

    cnt_rel = traj_relax_read[-1]

    cnt_rel.calc = BrennerPotential()

    

    dyn = VelocityVerlet(cnt_rel, timestep=timestep*units.fs)

    data = {'time':[],'temperature':[], 'epot':[],'ekin':[], 'etot': []}

    cnt_array_original = cnt_rel.cell[:]

    run_simulation(loop_number, step_number, timestep, eps, filename, start_deform) 

    with open(f'{data_directory}{filename}_data.json', 'r') as fp:

        data = json.load(fp)

    df = pd.DataFrame.from_dict(data)

    df.to_csv(f'{data_directory}cnt_nve_data.csv')

    print(f'Saving data to {data_directory}cnt_nve_data.csv')

    print('done')

requirements.txt

+3 −0
Original line number Diff line number Diff line
@@ -7,9 +7,12 @@ kiwisolver==1.4.8 
matplotlib==3.10.3

numpy==2.2.5

packaging==25.0

pandas==2.2.3

pillow==11.2.1

pyparsing==3.2.3

python-dateutil==2.9.0.post0

pytz==2025.2

scipy==1.15.3

six==1.17.0

tqdm==4.67.1

tzdata==2025.2
Impressum   Privacy Policy   WW8 Homepage   Explore   Help