Skip to content

Instantly share code, notes, and snippets.

@tabedzki

tabedzki/lammps_init_v2.py

Last active July 13, 2021 15:32
Show Gist options
  • Save tabedzki/be1513d3db0df4584d273eb23da1eb66 to your computer and use it in GitHub Desktop.
Save tabedzki/be1513d3db0df4584d273eb23da1eb66 to your computer and use it in GitHub Desktop.
Download ZIP
Polymer creation for LAMMPS
Raw
lammps_init_v2.py
import numpy as np
import pbc_utils as pbc
import sys
import os
import math
class input_config:
def __init__(self, xbox, ybox, zbox):
self.natoms = 0
self.nbonds = 0
self.nmasses = 0
self.ndihedrals = 0
self.nimpropers = 0
self.masses = []
self.ang_types = []
self.bond_types = []
self.bonds = np.array([], dtype=np.int64).reshape(0,4)
self.nbond_types = 0
self.nangles = 0
self.nang_types = 0
self.x = np.array([], dtype=np.int64).reshape(0,6)
# self.x = np.zeros((0, 6), 'd')
self.RESID = np.zeros((0, 3), 'd')
self.L = np.zeros(3, 'd')
self.L[0] = float(xbox)
self.L[1] = float(ybox)
self.L[2] = float(zbox)
self.lo = -(self.L)/2
self.hi = (self.L)/2
self.xlo = self.lo[0]
self.ylo = self.lo[1]
self.zlo = self.lo[2]
self.xhi = self.hi[0]
self.yhi = self.hi[1]
self.zhi = self.hi[2]
self.typ = np.zeros((0, 3), 'd')
self.num_chns = 0
self.periodic = False
def add_diblock_rho0(self, part1, part2, frac, chl, rho0, Lbond, bond_type):
num_chns = int(self.L[0] * self.L[1] * self.L[2] * rho0/chl)
print(num_chns)
self.add_diblock(part1, part2, frac, chl, num_chns, Lbond, bond_type)
def add_diblock(self, part1, part2, frac, chl, num_chns, Lbond,bond_type, rmin = 0.0):
if ( not part1 in self.masses):
self.masses.append(part1)
self.nmasses += 1
if ( not part2 in self.masses):
self.masses.append(part2)
self.nmasses += 1
if ( not bond_type in self.bond_types):
self.bond_types.append(bond_type)
self.nbond_types+= 1
# resid = self.natoms + 1
ns_loc = chl * num_chns
xloc = np.zeros((ns_loc, 6), 'd')
# bond_loc = np.zeros((0, 4), 'd')
# bond_loc = np.([], dtype=np.float).reshape(0,4)
nbonds_loc = num_chns * (chl - 1)
bond_loc = np.empty((nbonds_loc,4), int)
# self.nbonds
natoms = self.natoms
self.natoms += chl * num_chns
print(self.natoms)
chn_id = self.num_chns
self.num_chns += chl
bond_count = 0
print("STOP")
print( num_chns, chl)
for i_ch in range(num_chns):
for i_monomer in range(chl):
natoms += 1
# print(i_monomer, chl)
# xloc[i_ch]
# print(float(i_monomer)/float(chl))
if float(i_monomer)/float(chl) < frac:
xloc[i_ch*chl+i_monomer,2] = part1
else:
xloc[i_ch*chl+i_monomer,2] = part2
# print(i_monomer)
xloc[i_ch*chl+i_monomer,0] = natoms
xloc[i_ch*chl+i_monomer,1] = chn_id + i_ch
if i_monomer == 0:
# print("STOP")
xloc[i_ch*chl,3] = self.xlo + np.random.random_sample() * self.L[0]
xloc[i_ch*chl,4] = self.ylo + np.random.random_sample() * self.L[1]
xloc[i_ch*chl,5] = self.zlo + np.random.random_sample() * self.L[2]
# print(xloc[i_ch * chl])
else:
# print(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ]))
# bond_loc = np.concatenate((bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ])))
# bond_loc = np.concatenate(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1]), axis = 0)
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = bond_type
bond_loc[bond_count, 2] = natoms
bond_loc[bond_count, 3] = natoms - 1
bond_count += 1
self.nbonds += 1
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
xprev = xloc[i_ch*chl+i_monomer-1,3]
yprev = xloc[i_ch*chl+i_monomer-1,4]
zprev = xloc[i_ch*chl+i_monomer-1,5]
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_monomer == 1:
restriction = False
else:
xpp = xloc[i_ch*chl+i_monomer-2,3]
ypp = xloc[i_ch*chl+i_monomer-2,4]
zpp = xloc[i_ch*chl+i_monomer-2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
xloc[i_ch*chl+i_monomer,3] = xx
xloc[i_ch*chl+i_monomer,4] = yy
xloc[i_ch*chl+i_monomer,5] = zz
# print(self.x, xloc)
self.x = np.concatenate([self.x, xloc])
self.bonds = np.vstack([self.bonds, bond_loc])
def add_diblock_angle(self, part1, part2, frac, chl, num_chns, Lbond,bond_type,angle_type = None, rmin = 0.0):
if ( not part1 in self.masses):
self.masses.append(part1)
self.nmasses += 1
if ( not part2 in self.masses):
self.masses.append(part2)
self.nmasses += 1
if ( not bond_type in self.bond_types):
self.bond_types.append(bond_type)
self.nbond_types+= 1
resid = self.natoms + 1
ns_loc = chl * num_chns
xloc = np.zeros((ns_loc, 6), 'd')
# bond_loc = np.zeros((0, 4), 'd')
# bond_loc = np.([], dtype=np.float).reshape(0,4)
nbonds_loc = num_chns * (chl - 1)
bond_loc = np.empty((nbonds_loc,4), int)
nangles_loc = num_chns * (chl -2 )
bond_loc = np.empty((nangles_loc,4), int)
# self.nbonds
natoms = self.natoms
self.natoms += chl * num_chns
chn_id = self.num_chns
self.num_chns += chl
bond_count = 0
if not angle_type == None:
if ( not angle_type in self.ang_types):
self.ang_types.append(part2)
print("STOP")
print( num_chns, chl)
for i_ch in range(num_chns):
for i_monomer in range(chl):
natoms += 1
# print(i_monomer, chl)
# xloc[i_ch]
# print(float(i_monomer)/float(chl))
if float(i_monomer)/float(chl) < frac:
xloc[i_ch*chl+i_monomer,2] = part1
else:
xloc[i_ch*chl+i_monomer,2] = part2
# print(i_monomer)
xloc[i_ch*chl+i_monomer,0] = natoms
xloc[i_ch*chl+i_monomer,1] = chn_id + i_ch
if i_monomer == 0:
# print("STOP")
xloc[i_ch*chl,3] = self.xlo + np.random.random_sample() * self.L[0]
xloc[i_ch*chl,4] = self.ylo + np.random.random_sample() * self.L[1]
xloc[i_ch*chl,5] = self.zlo + np.random.random_sample() * self.L[2]
# print(xloc[i_ch * chl])
else:
# if not i_monomer + 1 == chl:
# print(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ]))
# bond_loc = np.concatenate((bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ])))
# bond_loc = np.concatenate(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1]), axis = 0)
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = bond_type
bond_loc[bond_count, 2] = natoms
bond_loc[bond_count, 3] = natoms - 1
bond_count += 1
self.nbonds += 1
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
xprev = xloc[i_ch*chl+i_monomer-1,3]
yprev = xloc[i_ch*chl+i_monomer-1,4]
zprev = xloc[i_ch*chl+i_monomer-1,5]
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_monomer == 1:
restriction = False
else:
xpp = xloc[i_ch*chl+i_monomer-2,3]
ypp = xloc[i_ch*chl+i_monomer-2,4]
zpp = xloc[i_ch*chl+i_monomer-2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
xloc[i_ch*chl+i_monomer,3] = xx
xloc[i_ch*chl+i_monomer,4] = yy
xloc[i_ch*chl+i_monomer,5] = zz
# print(self.x, xloc)
self.x = np.concatenate([self.x, xloc])
self.bonds = np.vstack([self.bonds, bond_loc])
def add_homopolyer(self, part, chl, num_chns, Lbond, bond_type):
self.add_diblock(part, part, 1.0, chl, num_chns, Lbond, bond_type)
def add_homopolyer_rho0(self, part, chl, rho0, Lbond, bond_type):
num_chns = int(self.L[0] * self.L[1] * self.L[2] * rho0/chl)
print(num_chns)
self.add_diblock(part, part, 1.0, chl, num_chns, Lbond, bond_type)
def add_ms(self, part, ang_type, num_chns):
if ( not part in self.masses):
self.masses.append(part)
self.nmasses += 1
if ( not ang_type in self.ang_types):
self.ang_types.append(part)
self.nang_types += 1
def write(self, output):
otp = open(output, 'w')
otp.write("Generated by Chris' code\n\n")
line = "%d atoms\n" % (self.natoms )
otp.write(line)
line = "%d bonds\n" % len(self.bonds)
# print(self.bonds)
otp.write(line)
line = "%d angles\n" % (self.nangles)
otp.write(line)
line = "%d dihedrals\n" % (self.ndihedrals)
otp.write(line)
line = "%d impropers\n\n" % (self.ndihedrals)
otp.write(line)
line = "%d atom types\n" % len(self.masses)
otp.write(line)
line = "%d bond types\n" % len(self.bond_types)
otp.write(line)
line = "%d angle types\n" % self.nang_types
otp.write(line)
line = "%d dihedral types\n" % self.ndihedrals
otp.write(line)
line = "%d improper types\n" % self.nimpropers
otp.write(line)
line = "\n"
otp.write(line)
line = '%f %f xlo xhi\n' % (self.lo[0], self.hi[0])
otp.write(line)
line = '%f %f ylo yhi\n' % (self.lo[1], self.hi[1])
otp.write(line)
line = '%f %f zlo zhi\n\n' % (self.lo[2], self.hi[2])
otp.write(line)
if len(self.masses) > 0 :
otp.write("Masses \n\n")
for i, val in enumerate(self.masses):
line = "%d 1.000\n" % (val)
otp.write(line)
otp.write("\nAtoms \n\n")
for i, val in enumerate(self.x):
line = "{:d} {:d} {:d} {:f} {:f} {:f}\n"
idx,m,t,x,y,z = val
# print(i)
# print(val)
# print(self.x)
# line = ' '.join(map(str, val))
otp.write(line.format(int(idx),int(m),int(t),x,y,z))
# otp.write(line + "\n")
if len(self.bonds) > 0 :
otp.write("\nBonds \n\n")
for i, val in enumerate(self.bonds):
# line = "%d %d %d %f %f %f\n" %
line = ' '.join(map(str, val))
otp.write(line + "\n")
def add_comb(self, bb_part1, Nb,Ns, num_chns, ss_pt1, back_bond, bb_part2=None, frac_bb=1, ss_pt2=None,
frac_side=1.0, Lbond=1.0, freq=1,
back_side_bond=None, side_bond=None, rmin = 0.0):
if ( not bb_part1 in self.masses and not bb_part1 == None ):
self.masses.append(bb_part1)
self.nmasses += 1
if ( not bb_part2 in self.masses and not bb_part2 == None):
self.masses.append(bb_part2)
self.nmasses += 1
if ( not ss_pt1 in self.masses and not ss_pt1 == None ):
self.masses.append(ss_pt1)
self.nmasses += 1
if ( not ss_pt2 in self.masses and not ss_pt2 == None):
self.masses.append(ss_pt2)
self.nmasses += 1
if ( not back_bond in self.bond_types and not back_bond == None):
self.bond_types.append(back_bond)
self.nbond_types+= 1
if ( not back_side_bond in self.bond_types and not back_side_bond == None):
self.bond_types.append(back_side_bond)
self.nbond_types+= 1
if ( not side_bond in self.bond_types and not side_bond == None):
self.bond_types.append(side_bond)
self.nbond_types+= 1
if side_bond == None:
side_bond = back_bond
if back_side_bond == None:
back_side_bond = back_bond
resid = self.natoms + 1
ns_loc = (Nb + Ns * Nb//freq) * num_chns
# ns_loc = chl * num_chns
xloc = np.zeros((ns_loc, 6), 'd')
# bond_loc = np.zeros((0, 4), 'd')
# bond_loc = np.([], dtype=np.float).reshape(0,4)
nbonds_loc = num_chns * ( (Nb - 1) + Nb//freq * (Ns) )
bond_loc = np.empty((nbonds_loc,4), int)
# self.nbonds
natoms = self.natoms
molecule_len = Nb + Ns * Nb//freq
self.natoms += (Nb + Ns * Nb//freq) * num_chns
chn_id = self.num_chns
self.num_chns += num_chns
bond_count = 0
# print("STOP")
# print( num_chns, chl)
# print((Nb + Ns * Nb//freq) * num_chns)
for i_ch in range(num_chns):
for i_monomer in range(Nb):
# natoms += 1
# print(i_monomer, chl)
# xloc[i_ch]
# print(float(i_monomer)/float(chl))
if float(i_monomer)/float(Nb) < frac_bb:
xloc[i_ch*molecule_len+i_monomer,2] = bb_part1
else:
xloc[i_ch*molecule_len+i_monomer,2] = bb_part2
# print(i_monomer)
xloc[i_ch*molecule_len+i_monomer,0] = natoms + i_ch * molecule_len + i_monomer
xloc[i_ch*molecule_len+i_monomer,1] = chn_id + i_ch # molecule id
if i_monomer == 0:
# print("STOP")
xloc[i_ch*molecule_len,3] = self.xlo + np.random.random_sample() * self.L[0]
xloc[i_ch*molecule_len,4] = self.ylo + np.random.random_sample() * self.L[1]
xloc[i_ch*molecule_len,5] = self.zlo + np.random.random_sample() * self.L[2]
# print(xloc[i_ch * chl])
else:
# print(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ]))
# bond_loc = np.concatenate((bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ])))
# bond_loc = np.concatenate(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1]), axis = 0)
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = back_bond
bond_loc[bond_count, 2] = natoms + i_ch*molecule_len + i_monomer -1
bond_loc[bond_count, 3] = natoms + i_ch*molecule_len + i_monomer
# bond_count += 1
# self.nbonds += 1
# theta = 2 * np.pi * np.random.random_sample()
# bond_loc[bond_count, 3] = natoms - 1
bond_count += 1
self.nbonds += 1
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
xprev = xloc[i_ch*molecule_len+i_monomer-1,3]
yprev = xloc[i_ch*molecule_len+i_monomer-1,4]
zprev = xloc[i_ch*molecule_len+i_monomer-1,5]
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_monomer == 1:
restriction = False
else:
xpp = xloc[i_ch*molecule_len+i_monomer-2,3]
ypp = xloc[i_ch*molecule_len+i_monomer-2,4]
zpp = xloc[i_ch*molecule_len+i_monomer-2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
xloc[i_ch*molecule_len+i_monomer,3] = xx
xloc[i_ch*molecule_len+i_monomer,4] = yy
xloc[i_ch*molecule_len+i_monomer,5] = zz
for i_side in range(Ns):
# natoms += 1
index = natoms + i_ch * molecule_len + Nb + i_monomer // freq * Ns + i_side
indbb = natoms + i_ch * molecule_len + i_monomer
xloc[index-natoms,0] = index
xloc[index-natoms,1] = chn_id + i_ch # molecule id
if float(i_side)/float(Ns) < frac_side:
xloc[index-natoms,2] = ss_pt1
else:
xloc[index -natoms,2] = ss_pt2
if i_side == 0:
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = back_side_bond
bond_loc[bond_count, 2] = indbb
bond_loc[bond_count, 3] = index
bond_count += 1
self.nbonds += 1
xprev = xloc[indbb - natoms,3]
yprev = xloc[indbb - natoms,4]
zprev = xloc[indbb - natoms,5]
else:
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = side_bond
bond_loc[bond_count, 2] = index - 1
bond_loc[bond_count, 3] = index
bond_count += 1
self.nbonds += 1
xprev = xloc[index - 1-natoms,3]
yprev = xloc[index - 1-natoms,4]
zprev = xloc[index - 1-natoms,5]
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_side == 0:
restriction = False
# xpp = xloc[indbb,3]
# ypp = xloc[indbb,4]
# zpp = xloc[indbb,5]
else:
xpp = xloc[index -natoms - 2,3]
ypp = xloc[index -natoms - 2,4]
zpp = xloc[index -natoms - 2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
# print(index, molecule_len * num_chns)
xloc[index-natoms,3] = xx
xloc[index-natoms,4] = yy
xloc[index-natoms,5] = zz
# print(self.x, xloc)
self.x = np.concatenate([self.x, xloc])
self.bonds = np.vstack([self.bonds, bond_loc])
def add_comb_rho0(self, bb_part1, Nb,Ns, rho0, ss_pt1, back_bond, bb_part2=None, frac_bb=1, ss_pt2=None,
frac_side=1.0, Lbond=1.0, freq=1,
back_side_bond=None, side_bond=None, rmin = 0.0):
if ( not bb_part1 in self.masses and not bb_part1 == None ):
self.masses.append(bb_part1)
self.nmasses += 1
if ( not bb_part2 in self.masses and not bb_part2 == None):
self.masses.append(bb_part2)
self.nmasses += 1
if ( not ss_pt1 in self.masses and not ss_pt1 == None ):
self.masses.append(ss_pt1)
self.nmasses += 1
if ( not ss_pt2 in self.masses and not ss_pt2 == None):
self.masses.append(ss_pt2)
self.nmasses += 1
if ( not back_bond in self.bond_types and not back_bond == None):
self.bond_types.append(back_bond)
self.nbond_types+= 1
if ( not back_side_bond in self.bond_types and not back_side_bond == None):
self.bond_types.append(back_side_bond)
self.nbond_types+= 1
if ( not side_bond in self.bond_types and not side_bond == None):
self.bond_types.append(side_bond)
self.nbond_types+= 1
if side_bond == None:
side_bond = back_bond
if back_side_bond == None:
back_side_bond = back_bond
num_chns = int(self.L[0] * self.L[1] * self.L[2] * rho0 / (Nb + math.ceil(float(Nb)/freq ) * Ns))
print(num_chns)
resid = self.natoms + 1
ns_loc = (Nb + Ns * Nb//freq) * num_chns
# ns_loc = chl * num_chns
xloc = np.zeros((ns_loc, 6), 'd')
# bond_loc = np.zeros((0, 4), 'd')
# bond_loc = np.([], dtype=np.float).reshape(0,4)
nbonds_loc = num_chns * ( (Nb - 1) + Nb//freq * (Ns) )
bond_loc = np.empty((nbonds_loc,4), int)
# self.nbonds
natoms = self.natoms
molecule_len = Nb + Ns * Nb//freq
self.natoms += (Nb + Ns * Nb//freq) * num_chns
chn_id = self.num_chns
self.num_chns += num_chns
bond_count = 0
# print("STOP")
# print( num_chns, chl)
# print((Nb + Ns * Nb//freq) * num_chns)
for i_ch in range(num_chns):
for i_monomer in range(Nb):
# natoms += 1
# print(i_monomer, chl)
# xloc[i_ch]
# print(float(i_monomer)/float(chl))
if float(i_monomer)/float(Nb) < frac_bb:
xloc[i_ch*molecule_len+i_monomer,2] = bb_part1
else:
xloc[i_ch*molecule_len+i_monomer,2] = bb_part2
# print(i_monomer)
xloc[i_ch*molecule_len+i_monomer,0] = natoms + i_ch * molecule_len + i_monomer
xloc[i_ch*molecule_len+i_monomer,1] = chn_id + i_ch # molecule id
if i_monomer == 0:
# print("STOP")
xloc[i_ch*molecule_len,3] = self.xlo + np.random.random_sample() * self.L[0]
xloc[i_ch*molecule_len,4] = self.ylo + np.random.random_sample() * self.L[1]
xloc[i_ch*molecule_len,5] = self.zlo + np.random.random_sample() * self.L[2]
# print(xloc[i_ch * chl])
else:
# print(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ]))
# bond_loc = np.concatenate((bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1 ])))
# bond_loc = np.concatenate(bond_loc, np.array([self.nbonds, bond_type, natoms, natoms -1]), axis = 0)
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = back_bond
bond_loc[bond_count, 2] = natoms + i_ch*molecule_len + i_monomer -1
bond_loc[bond_count, 3] = natoms + i_ch*molecule_len + i_monomer
# bond_count += 1
# self.nbonds += 1
# theta = 2 * np.pi * np.random.random_sample()
# bond_loc[bond_count, 3] = natoms - 1
bond_count += 1
self.nbonds += 1
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
xprev = xloc[i_ch*molecule_len+i_monomer-1,3]
yprev = xloc[i_ch*molecule_len+i_monomer-1,4]
zprev = xloc[i_ch*molecule_len+i_monomer-1,5]
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_monomer == 1:
restriction = False
else:
xpp = xloc[i_ch*molecule_len+i_monomer-2,3]
ypp = xloc[i_ch*molecule_len+i_monomer-2,4]
zpp = xloc[i_ch*molecule_len+i_monomer-2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
xloc[i_ch*molecule_len+i_monomer,3] = xx
xloc[i_ch*molecule_len+i_monomer,4] = yy
xloc[i_ch*molecule_len+i_monomer,5] = zz
for i_side in range(Ns):
# natoms += 1
index = natoms + i_ch * molecule_len + Nb + i_monomer // freq * Ns + i_side
indbb = natoms + i_ch * molecule_len + i_monomer
xloc[index-natoms,0] = index
xloc[index-natoms,1] = chn_id + i_ch # molecule id
if float(i_side)/float(Ns) < frac_side:
xloc[index-natoms,2] = ss_pt1
else:
xloc[index -natoms,2] = ss_pt2
if i_side == 0:
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = back_side_bond
bond_loc[bond_count, 2] = indbb
bond_loc[bond_count, 3] = index
bond_count += 1
self.nbonds += 1
xprev = xloc[indbb - natoms,3]
yprev = xloc[indbb - natoms,4]
zprev = xloc[indbb - natoms,5]
else:
bond_loc[bond_count, 0] = self.nbonds
bond_loc[bond_count, 1] = side_bond
bond_loc[bond_count, 2] = index - 1
bond_loc[bond_count, 3] = index
bond_count += 1
self.nbonds += 1
xprev = xloc[index - 1-natoms,3]
yprev = xloc[index - 1-natoms,4]
zprev = xloc[index - 1-natoms,5]
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_side == 0:
restriction = False
# xpp = xloc[indbb,3]
# ypp = xloc[indbb,4]
# zpp = xloc[indbb,5]
else:
xpp = xloc[index -natoms - 2,3]
ypp = xloc[index -natoms - 2,4]
zpp = xloc[index -natoms - 2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
# print(index, molecule_len * num_chns)
xloc[index-natoms,3] = xx
xloc[index-natoms,4] = yy
xloc[index-natoms,5] = zz
# print(self.x, xloc)
self.x = np.concatenate([self.x, xloc])
self.bonds = np.vstack([self.bonds, bond_loc])
def add_triblock(self, part1, part2, part3, frac1, frac2, chl, num_chns, Lbond,bond_type12, bond_type23, rmin = 0.0):
if ( not part1 in self.masses):
self.masses.append(part1)
self.nmasses += 1
if ( not part2 in self.masses):
self.masses.append(part2)
self.nmasses += 1
if ( not part3 in self.masses):
self.masses.append(part3)
self.nmasses += 1
if ( not bond_type12 in self.bond_types):
self.bond_types.append(bond_type12)
self.nbond_types+= 1
if ( not bond_type23 in self.bond_types):
self.bond_types.append(bond_type23)
self.nbond_types+= 1
# resid = self.natoms + 1
ns_loc = chl * num_chns
xloc = np.zeros((ns_loc, 6), 'd')
# bond_loc = np.zeros((0, 4), 'd')
# bond_loc = np.([], dtype=np.float).reshape(0,4)
nbonds_loc = num_chns * (chl - 1)
bond_loc = np.empty((nbonds_loc,4), int)
# self.nbonds
natoms = self.natoms
self.natoms += chl * num_chns
print(self.natoms)
chn_id = self.num_chns
self.num_chns += chl
bond_count = 0
print("STOP")
print( num_chns, chl)
for i_ch in range(num_chns):
for i_monomer in range(chl):
natoms += 1
# print(i_monomer, chl)
# xloc[i_ch]
# print(float(i_monomer)/float(chl))
f_along = float(i_monomer)/float(chl)
if f_along < frac1:
xloc[i_ch*chl+i_monomer,2] = part1
elif f_along < frac1+frac2:
xloc[i_ch*chl+i_monomer,2] = part2
else:
xloc[i_ch*chl+i_monomer,2] = part3
# print(i_monomer)
xloc[i_ch*chl+i_monomer,0] = natoms
xloc[i_ch*chl+i_monomer,1] = chn_id + i_ch
if i_monomer == 0:
# print("STOP")
xloc[i_ch*chl,3] = self.xlo + np.random.random_sample() * self.L[0]
xloc[i_ch*chl,4] = self.ylo + np.random.random_sample() * self.L[1]
xloc[i_ch*chl,5] = self.zlo + np.random.random_sample() * self.L[2]
# print(xloc[i_ch * chl])
else:
if f_along >= frac1+ frac2:
bndtyp = bond_type23
else:
bndtyp = bond_type12
bond_loc[bond_count, 0] = self.nbonds
# bond_loc[bond_count, 1] = bond_type
bond_loc[bond_count, 1] = bndtyp
bond_loc[bond_count, 2] = natoms
bond_loc[bond_count, 3] = natoms - 1
bond_count += 1
self.nbonds += 1
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(theta)
xprev = xloc[i_ch*chl+i_monomer-1,3]
yprev = xloc[i_ch*chl+i_monomer-1,4]
zprev = xloc[i_ch*chl+i_monomer-1,5]
restriction = True
while restriction:
theta = 2 * np.pi * np.random.random_sample()
phi = np.pi * np.random.random_sample()
dx = Lbond * np.sin(phi) * np.cos(theta)
dy = Lbond * np.sin(phi) * np.sin(theta)
dz = Lbond * np.cos(phi)
xx = xprev + dx
yy = yprev + dy
zz = zprev + dz
if np.abs(zz) < self.L[2]/2. :
if i_monomer == 1:
restriction = False
else:
xpp = xloc[i_ch*chl+i_monomer-2,3]
ypp = xloc[i_ch*chl+i_monomer-2,4]
zpp = xloc[i_ch*chl+i_monomer-2,5]
dxp = xx - xpp
dyp = yy - ypp
dzp = zz - zpp
rpsq = dxp*dxp+dyp*dyp+dzp*dzp
rp = np.sqrt(rpsq)
if rp > rmin:
restriction = False
if self.periodic == True:
if xx > self.xhi:
xx -= self.L[0]
if yy > self.yhi:
yy -= self.L[1]
if xx < self.xlo:
xx += self.L[0]
if yy < self.ylo:
yy += self.L[1]
xloc[i_ch*chl+i_monomer,3] = xx
xloc[i_ch*chl+i_monomer,4] = yy
xloc[i_ch*chl+i_monomer,5] = zz
# print(self.x, xloc)
self.x = np.concatenate([self.x, xloc])
self.bonds = np.vstack([self.bonds, bond_loc])
Raw
make_comb_phase_diagram.py
import lammps_init_v2 as lmp
import numpy as np
box = lmp.input_config(10,10,10)
box.periodic = False
box.add_homopolyer(1, 10, 10, 0.5, 1)
box.add_homopolyer(2, 10, 10, 0.5, 1)
box.add_homopolyer(3, 10, 10, 0.5, 1)
# box.add_homopolyer(2, 10, 20, 0.5, 1)
# box.add_diblock(1,1,0.2,10,8,1,1)
# box.add_triblock(1,2,1,0.2,0.4,10,8,1,2,1)
# box.add_diblock(2,2,0.2,10,8,1,1)
box.write("tripolymer_melt.data-debugging")
Raw
pbc_utils.py
import numpy as nmp
import math as m
def pbc_vdr(x1, x2, box, boxh):
dr = x1 - x2
while (dr[0] >= boxh[0]):
dr[0] -= box[0]
while (dr[0] < -boxh[0]):
dr[0] += box[0]
while (dr[1] >= boxh[1]):
dr[1] -= box[1]
while (dr[1] < -boxh[1]):
dr[1] += box[1]
while (dr[2] >= boxh[2]):
dr[2] -= box[2]
while (dr[2] < -boxh[2]):
dr[2] += box[2]
return dr
def pbc_mdr(x1, x2, box, boxh):
dr = pbc_vdr(x1,x2,box,boxh)
mdr = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2]
return m.sqrt(mdr)
def pbc_pos_inbox(x1, box):
if (x1[0] >= box[0]):
x1[0] -= box[0]
elif (x1[0] < 0.0):
x1[0] += box[0]
if (x1[1] >= box[1]):
x1[1] -= box[1]
elif (x1[1] < 0.0):
x1[1] += box[1]
if (x1[2] >= box[2]):
x1[2] -= box[2]
elif (x1[2] < 0.0):
x1[2] += box[2]
return x1
def pbc_inbox(x1, box, boxh):
if (x1[0] >= boxh[0]):
x1[0] -= box[0]
elif (x1[0] < -boxh[0]):
x1[0] += box[0]
if (x1[1] >= boxh[1]):
x1[1] -= box[1]
elif (x1[1] < -boxh[1]):
x1[1] += box[1]
if (x1[2] >= boxh[2]):
x1[2] -= box[2]
elif (x1[2] < -boxh[2]):
x1[2] += box[2]
return x1
def pbc_mdr2(x1, x2, box, boxh):
dr = pbc_vdr(x1,x2,box,boxh)
mdr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2]
return mdr2
def pbc_nojumps(x,bx):
ns = nmp.shape(x)[1]
fr = nmp.shape(x)[0]
bxh = bx * 0.5
for t in range(1,fr):
for i in range(0,ns):
dr = pbc_vdr(x[t-1,i], x[t,i], bx[t], bxh[t]) #x[t-1] - x[t]
x[t] = x[t-1] - dr
return x
def guess_box(x):
# If "x" contains a trajectory
if ( len(nmp.shape(x)) == 3):
fr = nmp.shape(x)[0]
bx = nmp.zeros((fr,3),'d')
for i in range(0,fr):
for j in range(0,3):
bx[i,j] = nmp.max(x[i,:,j]) - nmp.min(x[i,:,j])
# If "x" is a single frame
else:
bx = nmp.zeros(3,'d')
for j in range(0,3):
bx[j] = nmp.max(x[:,j]) - nmp.min(x[:,j])
return bx
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment