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.
- create the geometry file:
from sknano.generators import MWNTBundleGenerator
from ase.io import read, write
from ase.visualize import view
mwntbundle = MWNTBundleGenerator(Nwalls=3, min_wall_diameter=5, Lz=5, bundle_geometry='hexagon')
mwntbundle.save('mwcnt_boundle.xyz')
mwcnt_ase = read('mwcnt_boundle.xyz')
view(mwcnt_ase)
- 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('Masses\n')
f.write(' 1 12.010700 # C\n')
f.write('Atoms\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]) )
f.write('\n')
n_atoms, coords = read_xyz("mwcnt_boundle.xyz")
write_lammps(coords)
- 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
Marco
# 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}
#---------------STRESS-------------------------------------
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
# ----------------------------------------------
# NPT
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
displace_atoms G1 move 0 0 10
run ${Nruns}
write_data restart.npt
undump Dnpt
unfix mynpt
unfix fnpt
print '---- DONE ----'