Define stress-strain for a multi-wall carbon nanotubes bundle

Dear all,

I want to simulate a MWCNT bundle in which the central MWCNT is pulled out of the bundle.
How to define the stress/strain curve of such system?

I am building my system with scikit-nano and using LAMMPS with Airebo-M force field for simulation.

  1. create the geometry file:
from sknano.generators import MWNTBundleGenerator
from import read, write
from ase.visualize import view

mwntbundle = MWNTBundleGenerator(Nwalls=3, min_wall_diameter=5, Lz=5, bundle_geometry='hexagon')'')
mwcnt_ase = read('')
  1. Transform the so-produced xyz file into a lammps geometry file:
import numpy as np

def read_xyz(xyz_file):
    with open(xyz_file, 'r') as f:
        xyz_lines = f.readlines()
    n_atoms = int(xyz_lines[0])
    coords_str = np.array([ ll.split()[1:] for ll in xyz_lines[2:] ])
    coords_float = coords_str.astype(np.float)
    return n_atoms, coords_float

def write_lammps(coords):
    nat_per_mwcnt = len(coords)/7
    with open('mwcnt_boundle.lammps', "w+") as f:
        f.write('# LAMMPS data file \n\n')
        f.write(str(len(coords))+' atoms\n\n')
        f.write('1 atom types\n\n')
        f.write('-100 100 xlo xhi\n')
        f.write('-100 100 ylo yhi\n')
        f.write('-100 100 zlo zhi\n\n')
        f.write(' 1 12.010700 # C\n')
        for cc in range(len(coords)):
            crd = coords[cc]
            ## creating 7 different molecules, each one corresponding to a MWCNT
            if cc%nat_per_mwcnt == 0:
                nmol = int(cc//nat_per_mwcnt)
            f.write( '{:4d} {:2d} 1 0 {:6f} {:6f} {:6f}\n'.format(cc+1, nmol, crd[0], crd[1], crd[2]) )

n_atoms, coords = read_xyz("")
  1. At this point I can use LAMMPS to simulate the central mwcnt being pulled out of the bundle. The mwcnt correspond to residue 1. Please find below my input file. In a nutshell: I define a cilynder enclosing molecule G1 and pull it out with the displace_atoms command from LAMMPS. This will go into an external loop making it more effective.
    I have taken inspiration from Ref1 and Ref2.

I guess I can give a “global” definition of stress/strain, ie. the strain being the $\Delta$Z of the full system and the stress being the total sum of the stress of all atoms, a “surface” definition where I can explicitely impose a force on the surface of the central mwcnt, or a “local” one where only the atoms in the central mwcnt are considered. The final goal is to have an estimation of the resistance of the whole system to external longitudinal strain.

Here I am defining the strain as the $\Delta$Z between the top layer of the shifted mwcnt and its initial position, while the stress is defined with “compute stress/atom”.

Unsurprisingly, even with a larger system than the one produced here, I get a not very clear result, with strain/stress oscillating between -1e10 and 1e10 GPa.

Can anyone suggest a way to improve this simulation?

Thank you


# Initialisation

units           metal
dimension           3
boundary        p p p
atom_style      full 
newton          on

read_data       mwcnt_boundle.lammps 
pair_style      airebo/morse 3.0 1 1
pair_coeff      * * /data/mdi0316/WORK/MWCNT/src/Inputfiles/CH.airebo-m C

# time settings
variable        dt      equal   0.001
variable        tprod   equal   5
variable        Nruns    equal floor(${tprod}/${dt})
variable        Nevery   equal 10
variable        Nrepeat  equal 10
variable        Nfreq    equal 100
timestep        ${dt} #ps
print           "--- Time settings ---"
print           "    dt       = ${dt}"
print           "    tprod    = ${tprod}"
print           "    Nruns    = ${Nruns}"
print           "    Nevery   = ${Nevery}"
print           "    Nrepeat  = ${Nrepeat}"
print           "    Nfreq    = ${Nfreq}"

# computes and output variables settings
compute         cKe all ke
compute         cPe all pe
variable        totE equal c_cKe+c_cPe
compute         Tall all temp
compute         cP all pressure Tall
thermo          100
thermo_style    custom time temp pe v_totE press vol density

# ----------------------------------------------
# groups settings
group           G1 molecule 1
variable        xmax equal bound(G5,xmax)+0.5
variable        xmin equal bound(G5,xmin)-0.5
variable        zmax equal bound(G5,zmax)+0.5
variable        zmin equal bound(G5,zmin)-0.5
variable        Cdiam equal ${xmax}-${xmin}
variable        Cradi equal ${Cdiam}/2
variable        Cleng equal ${zmax}-${zmin}

variable        x_com_inner equal xcm(G1,x)
variable        y_com_inner equal xcm(G1,y)

region          innreg cylinder z ${x_com_inner} ${y_com_inner} ${Cradi} ${zmin} ${zmax}
variable        innvol equal ${Cleng}*3.14*${Cradi}*${Cradi}

compute         1 all stress/atom NULL
compute         2 all reduce sum c_1[3]
fix             1 all ave/time 1 1 100 c_2
variable        sigmazz equal f_1

## -------------- OUTPUTS ----------------------
fix             3 all ave/time 1 1 100 v_zmin v_zmax v_sigmazz file displace.dat

# ----------------------------------------------

fix             bar all temp/berendsen 300 300 0.1 #Tstart Tstop Tdamp = 0.1ps
fix             mynpt all npt temp 300 300 0.1 aniso 1 1 1
dump            Dnpt all atom 100 dump_npt.lammpstrj
fix             fnpt all ave/time ${Nevery} ${Nrepeat} ${Nfreq} c_Tall v_totE c_cP v_innvol file npt.dat
print           "--- Beginning of NPT run (${Nruns} steps) ---"
displace_atoms  G1 move 0 0 10
run             ${Nruns}
write_data      restart.npt
undump          Dnpt
unfix           mynpt
unfix           fnpt

print           '---- DONE ----'

Why should there be much resistance to pulling the NT out from the bundle?
Why should it be dependent on the distance you have pulled?
I would expect that pulled tube will glide out from the bundle without much resistance, even if that would be a very tight configuration and that there may be small “bumps” when the atoms have to pass each other, but those would be rather small since that surface-surface interaction is very smooth.
Given the typical size of fluctuation of data from atomic scale simulations, particularly properties that are very sensitive like stress, I would expect that you need some extreme amount of averaging to get meaningful data that consistently becomes visible above the noise.

How to define the stress/strain curve of such system?

Stress-strain curves are a macroscopic material mechanics concept. Strain is likely to be defined for a field of coordinate with regard to an initial configuration of the simulation box. If you want to evaluate the force needed to pull out your central NT out of the bundle, you are more likely to be interested in the Potential of Mean Force and the COLVARS module similarly to this work:

But if you want to compute a stress/strain curve of a whole NT bundle, you could stretch a periodic MWCNT along the periodic direction and compute the average stress components of the system, eventually until it breaks. But as Axel already mentioned, this can require a high amount of averaging, especially at finite temperature.