mirror of
https://github.com/ParticulateFlow/LPP.git
synced 2025-12-08 06:37:46 +00:00
247 lines
7.0 KiB
Python
247 lines
7.0 KiB
Python
# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
|
|
# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
#
|
|
# Copyright (2005) Sandia Corporation. Under the terms of Contract
|
|
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
|
# certain rights in this software. This software is distributed under
|
|
# the GNU General Public License.
|
|
|
|
# chain tool
|
|
|
|
# Imports and external programs
|
|
|
|
from __future__ import print_function, division, absolute_import
|
|
import math
|
|
from data import data
|
|
|
|
oneline = "Create bead-spring chains for LAMMPS input"
|
|
|
|
docstr = """
|
|
c = chain(N,rho) setup box with N monomers at reduced density rho
|
|
c = chain(N,rho,1,1,2) x,y,z = aspect ratio of box (def = 1,1,1)
|
|
|
|
c.seed = 48379 set random # seed (def = 12345)
|
|
c.mtype = 2 set type of monomers (def = 1)
|
|
c.btype = 1 set type of bonds (def = 1)
|
|
c.blen = 0.97 set length of bonds (def = 0.97)
|
|
c.dmin = 1.02 set min dist from i-1 to i+1 site (def = 1.02)
|
|
|
|
c.id = "chain" set molecule ID to chain # (default)
|
|
c.id = "end1" set molecule ID to count from one end of chain
|
|
c.id = "end2" set molecule ID to count from either end of chain
|
|
|
|
c.build(100,10) create 100 chains, each of length 10
|
|
|
|
can be invoked multiple times interleaved with different settings
|
|
must fill box with total of N monomers
|
|
|
|
c.write("data.file") write out all built chains to LAMMPS data file
|
|
"""
|
|
|
|
# History
|
|
# 8/05, Steve Plimpton (SNL): original version
|
|
|
|
# ToDo list
|
|
|
|
# Variables
|
|
# n = number of monomers
|
|
# rhostar = reduced density
|
|
# seed = 12345
|
|
# mtype = type of monomers
|
|
# btype = type of bonds
|
|
# blen = length of bonds
|
|
# dmin = minimum distance from i-1 to i+1
|
|
# id = "chain","end1",or "end2"
|
|
# atoms = list of atoms
|
|
# bonds = list of bonds
|
|
# xprd,yprd,zprd = x,y,z box size
|
|
# xlo,ylo,zlo = -xyz prd / 2
|
|
# xhi,yhi,zhi = x,y,zprd /2
|
|
|
|
# Class definition
|
|
|
|
class chain:
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def __init__(self,n,rhostar,*list):
|
|
self.n = n
|
|
self.rhostar = rhostar
|
|
xaspect = yaspect = zaspect = 1.0
|
|
if len(list):
|
|
xaspect = list[0]
|
|
yaspect = list[1]
|
|
zaspect = list[2]
|
|
self.seed = 12345
|
|
self.mtype = 1
|
|
self.btype = 1
|
|
self.blen = 0.97
|
|
self.dmin = 1.02
|
|
self.id = "chain"
|
|
self.atoms = []
|
|
self.bonds = []
|
|
|
|
volume = n/rhostar
|
|
prd = pow(volume/xaspect/yaspect/zaspect,1.0/3.0)
|
|
self.xprd = xaspect * prd
|
|
self.xlo = -self.xprd/2.0
|
|
self.xhi = self.xprd/2.0
|
|
self.yprd = yaspect * prd
|
|
self.ylo = -self.yprd/2.0
|
|
self.yhi = self.yprd/2.0
|
|
self.zprd = zaspect * prd
|
|
self.zlo = -self.zprd/2.0
|
|
self.zhi = self.zprd/2.0
|
|
|
|
print("Simulation box: %g by %g by %g" % (self.xprd,self.yprd,self.zprd))
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def build(self,n,nper):
|
|
for ichain in range(n):
|
|
atoms = []
|
|
bonds = []
|
|
id_atom_prev = id_mol_prev = id_bond_prev = 0
|
|
if len(self.atoms):
|
|
id_atom_prev = self.atoms[-1][0]
|
|
id_mol_prev = self.atoms[-1][1]
|
|
if len(self.bonds):
|
|
id_bond_prev = self.bonds[-1][0]
|
|
|
|
for imonomer in range(nper):
|
|
if imonomer == 0:
|
|
x = self.xlo + self.random()*self.xprd
|
|
y = self.ylo + self.random()*self.yprd
|
|
z = self.zlo + self.random()*self.zprd
|
|
ix = iy = iz = 0
|
|
else:
|
|
restriction = True
|
|
while restriction:
|
|
rsq = 2.0
|
|
while rsq > 1.0:
|
|
dx = 2.0*self.random() - 1.0
|
|
dy = 2.0*self.random() - 1.0
|
|
dz = 2.0*self.random() - 1.0
|
|
rsq = dx*dx + dy*dy + dz*dz
|
|
r = math.sqrt(rsq)
|
|
dx,dy,dz = dx/r,dy/r,dz/r
|
|
x = atoms[-1][3] + dx*self.blen
|
|
y = atoms[-1][4] + dy*self.blen
|
|
z = atoms[-1][5] + dz*self.blen
|
|
restriction = False
|
|
if imonomer >= 2:
|
|
dx = x - atoms[-2][3]
|
|
dy = y - atoms[-2][4]
|
|
dz = z - atoms[-2][5]
|
|
if math.sqrt(dx*dx + dy*dy + dz*dz) <= self.dmin:
|
|
restriction = True
|
|
|
|
x,y,z,ix,iy,iz = self.pbc(x,y,z,ix,iy,iz)
|
|
idatom = id_atom_prev + imonomer + 1
|
|
if self.id == "chain":
|
|
idmol = id_mol_prev + 1
|
|
elif self.id == "end1":
|
|
idmol = imonomer + 1
|
|
elif self.id == "end2":
|
|
idmol = imonomer + 1
|
|
if idmol > nper//2:
|
|
idmol = nper - imonomer
|
|
else:
|
|
raise Exception("chain ID is not a valid value")
|
|
|
|
atoms.append([idatom,idmol,self.mtype,x,y,z,ix,iy,iz])
|
|
if imonomer:
|
|
bondid = id_bond_prev + imonomer
|
|
bonds.append([bondid,self.btype,idatom-1,idatom])
|
|
|
|
self.atoms += atoms
|
|
self.bonds += bonds
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def write(self,file):
|
|
if len(self.atoms) != self.n:
|
|
raise Exception("%d monomers instead of requested %d" % \
|
|
(len(self.atoms),self.n))
|
|
|
|
atlist = [atom[2] for atom in self.atoms]
|
|
atypes = max(atlist)
|
|
|
|
btypes = 0
|
|
if len(self.bonds):
|
|
btlist = [bond[1] for bond in self.bonds]
|
|
btypes = max(btlist)
|
|
|
|
# create the data file
|
|
|
|
d = data()
|
|
d.title = "LAMMPS FENE chain data file"
|
|
d.headers["atoms"] = len(self.atoms)
|
|
d.headers["bonds"] = len(self.bonds)
|
|
d.headers["atom types"] = atypes
|
|
d.headers["bond types"] = btypes
|
|
d.headers["xlo xhi"] = (self.xlo,self.xhi)
|
|
d.headers["ylo yhi"] = (self.ylo,self.yhi)
|
|
d.headers["zlo zhi"] = (self.zlo,self.zhi)
|
|
|
|
lines = []
|
|
for i in range(atypes): lines.append("%d 1.0\n" % (i+1))
|
|
d.sections["Masses"] = lines
|
|
|
|
lines = []
|
|
for atom in self.atoms:
|
|
line = "%d %d %d %g %g %g %d %d %d\n" % \
|
|
(atom[0], atom[1], atom[2], atom[3], atom[4], atom[5],
|
|
atom[6], atom[7], atom[8])
|
|
lines.append(line)
|
|
d.sections["Atoms"] = lines
|
|
|
|
lines = []
|
|
for bond in self.bonds:
|
|
line = "%d %d %d %d\n" % (bond[0], bond[1], bond[2], bond[3])
|
|
lines.append(line)
|
|
d.sections["Bonds"] = lines
|
|
|
|
d.write(file)
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def pbc(self,x,y,z,ix,iy,iz):
|
|
if x < self.xlo:
|
|
x += self.xprd
|
|
ix -= 1
|
|
elif x >= self.xhi:
|
|
x -= self.xprd
|
|
ix += 1
|
|
if y < self.ylo:
|
|
y += self.yprd
|
|
iy -= 1
|
|
elif y >= self.yhi:
|
|
y -= self.yprd
|
|
iy += 1
|
|
if z < self.zlo:
|
|
z += self.zprd
|
|
iz -= 1
|
|
elif z >= self.zhi:
|
|
z -= self.zprd
|
|
iz += 1
|
|
return x,y,z,ix,iy,iz
|
|
|
|
# --------------------------------------------------------------------
|
|
|
|
def random(self):
|
|
k = self.seed//IQ
|
|
self.seed = IA*(self.seed-k*IQ) - IR*k
|
|
if self.seed < 0:
|
|
self.seed += IM
|
|
return AM*self.seed
|
|
|
|
# --------------------------------------------------------------------
|
|
# random # generator constants
|
|
|
|
IM = 2147483647
|
|
AM = 1.0/IM
|
|
IA = 16807
|
|
IQ = 127773
|
|
IR = 2836
|