# ******************************************************************************
# pysimm.apps.random_walk module
# ******************************************************************************
#
# psuedo random walk algorithm written using pysimm tools
#
# ******************************************************************************
# License
# ******************************************************************************
# The MIT License (MIT)
#
# Copyright (c) 2016 Michael E. Fortunato, Coray M. Colina
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
from time import strftime
from itertools import permutations, izip
import numpy as np
from pysimm import system, lmps, forcefield, calc
from pysimm import error_print
[docs]def find_last_backbone_vector(s, m):
"""pysimm.apps.random_walk.find_last_backbone_vector
Finds vector between backbone atoms in terminal monomer. Requires current system s, and reference monomer m.
Args:
s: :class:`~pysimm.system.System` object
m: :class:`~pysimm.system.System` object
Returns:
list of vector components
"""
head_pos = [0, 0, 0]
tail_pos = [0, 0, 0]
for p in s.particles[-1*m.particles.count:]:
if p.linker == 'head':
head_pos = [p.x, p.y, p.z]
elif p.linker == 'tail':
tail_pos = [p.x, p.y, p.z]
return [head_pos[0] - tail_pos[0], head_pos[1] - tail_pos[1], head_pos[2] - tail_pos[2]]
[docs]def copolymer(m, nmon, s_=None, **kwargs):
"""pysimm.apps.random_walk.copolymer
Builds copolymer using random walk methodology using pattern
Args:
m: list of reference monomer :class:`~pysimm.system.System`s
nmon: total number of monomers to add to chain
s_: :class:`~pysimm.system.System` in which to build polymer chain (None)
settings: dictionary of simulation settings
density: density at which to build polymer (0.3)
forcefield: :class:`~pysimm.forcefield.Forcefield` object to acquire new force field parameters
capped: True/False if monomers are capped
unwrap: True to unwrap final system
traj: True to build xyz trajectory of polymer growth (True)
pattern: list of pattern for monomer repeat units, should match length of m ([1 for _ in range(len(m))])
limit: during MD, limit atomic displacement by this max value (LAMMPS ONLY)
sim: :class:`~pysimm.lmps.Simulation` object for relaxation between polymer growth
Returns:
new copolymer :class:`~pysimm.system.System`
"""
m = [x.copy() for x in m]
settings = kwargs.get('settings', {})
density = kwargs.get('density', 0.3)
f = kwargs.get('forcefield')
capped = kwargs.get('capped')
unwrap = kwargs.get('unwrap')
traj = kwargs.get('traj', True)
pattern = kwargs.get('pattern', [1 for _ in range(len(m))])
limit = kwargs.get('limit', 0.1)
sim = kwargs.get('sim')
for m_ in m:
m_.add_particle_bonding()
for p in m_.particles:
if p.type.name.find('@') >= 0 and p.type.name.split('@')[0].find('H'):
p.linker = 'head'
elif p.type.name.find('@') >= 0 and p.type.name.split('@')[0].find('T'):
p.linker = 'tail'
m_.remove_linker_types()
if s_ is None:
s = system.replicate(m[0], 1, density=density/nmon)
else:
s = system.replicate(m[0], 1, s_=s_, density=density/nmon)
print('%s: %s/%s monomers added' % (strftime('%H:%M:%S'), 1, nmon))
for p in s.particles:
if p.linker == 'head':
last_head = p
elif p.linker == 'tail':
last_tail = p
for m_ in m:
if capped:
m_.particles.remove(1)
m_.remove_spare_bonding()
m_.add_particle_bonding()
s.add_particle_bonding()
if traj:
s.write_xyz('random_walk.xyz')
temp_nmon = 1
while True:
m_ = m.pop(0)
m.append(m_)
p_ = pattern.pop(0)
pattern.append(p_)
if temp_nmon == 1 and p_ == 1:
m_ = m.pop(0)
m.append(m_)
p_ = pattern.pop(0)
pattern.append(p_)
elif temp_nmon == 1:
p_ -= 1
for insert in range(p_):
head = None
tail = None
backbone_vector = np.array([last_head.x - last_tail.x,
last_head.y - last_tail.y,
last_head.z - last_tail.z])
ref_head = None
ref_tail = None
for p in m_.particles:
if p.linker == 'head':
ref_head = p
elif p.linker == 'tail':
ref_tail = p
if ref_head and ref_tail:
ref_backbone_vector = np.array([ref_head.x - ref_tail.x,
ref_head.y - ref_tail.y,
ref_head.z - ref_tail.z])
rot_matrix = calc.find_rotation(ref_backbone_vector, backbone_vector)
m_.rotate(around=ref_tail, rot_matrix=rot_matrix)
translation_vector = [last_tail.x - ref_tail.x,
last_tail.y - ref_tail.y,
last_tail.z - ref_tail.z]
for p in m_.particles:
p.x = p.x + translation_vector[0] + 3*backbone_vector[0]
p.y = p.y + translation_vector[1] + 3*backbone_vector[1]
p.z = p.z + translation_vector[2] + 3*backbone_vector[2]
else:
print('reference molecule has no head or tail')
n = m_.copy()
if capped:
s.particles.remove(s.particles.count)
s.remove_spare_bonding()
s.add_particle_bonding()
s.add(n, change_dim=False)
s.add_particle_bonding()
head = last_head
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'tail':
tail = p
s.make_new_bonds(head, tail, f)
temp_nmon += 1
print('%s: %s/%s monomers added' % (strftime('%H:%M:%S'), temp_nmon, nmon))
if unwrap:
s.unwrap()
if sim is None:
sim = lmps.Simulation(s, name='relax_%03d' % (temp_nmon), log='relax.log', **settings)
sim.add_md(ensemble='nve', limit=limit, **settings)
sim.add_min(**settings)
if isinstance(sim, lmps.Simulation):
sim.system = s
sim.name = 'relax_%03d' % (temp_nmon)
sim.run(np=settings.get('np'))
if unwrap:
s.unwrap()
if unwrap:
s.wrap()
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'head':
last_head = p
elif p.linker == 'tail':
last_tail = p
if temp_nmon >= nmon:
break
if unwrap:
if not s.unwrap():
error_print('something went wrong')
return s
if traj:
s.write_xyz('random_walk.xyz', append=True)
if unwrap:
s.wrap()
for p in s.particles:
if p not in s.molecules[p.molecule.tag].particles:
s.molecules[p.molecule.tag].particles.add(p)
s.write_lammps('polymer.lmps')
s.unwrap()
s.write_xyz('polymer.xyz')
return s
[docs]def random_walk(m, nmon, s_=None, **kwargs):
"""pysimm.apps.random_walk.random_walk
Builds homopolymer using random walk methodology
Args:
m: reference monomer :class:`~pysimm.system.System`
nmon: total number of monomers to add to chain
s_: :class:`~pysimm.system.System` in which to build polymer chain (None)
extra_bonds: EXPERMINTAL, True if making ladder backbone polymer
settings: dictionary of simulation settings
density: density at which to build polymer (0.3)
forcefield: :class:`~pysimm.forcefield.Forcefield` object to acquire new force field parameters
capped: True/False if monomers are capped
unwrap: True to unwrap final system
traj: True to build xyz trajectory of polymer growth (True)
limit: during MD, limit atomic displacement by this max value (LAMMPS ONLY)
sim: :class:`~pysimm.lmps.Simulation` object for relaxation between polymer growth
Returns:
new polymer :class:`~pysimm.system.System`
"""
m = m.copy()
extra_bonds = kwargs.get('extra_bonds', False)
settings = kwargs.get('settings', {})
density = kwargs.get('density', 0.3)
f = kwargs.get('forcefield')
capped = kwargs.get('capped')
unwrap = kwargs.get('unwrap')
traj = kwargs.get('traj', True)
limit = kwargs.get('limit', 0.1)
sim = kwargs.get('sim')
m.add_particle_bonding()
for p in m.particles:
if p.type.name.find('@') >= 0 and p.type.name.split('@')[0].find('H'):
p.linker = 'head'
elif p.type.name.find('@') >= 0 and p.type.name.split('@')[0].find('T'):
p.linker = 'tail'
m.remove_linker_types()
if s_ is None:
s = system.replicate(m, 1, density=density/nmon)
else:
s = system.replicate(m, 1, s_=s_, density=None)
print('%s: %s/%s monomers added' % (strftime('%H:%M:%S'), 1, nmon))
if traj:
s.write_xyz('random_walk.xyz')
if capped:
m.particles.remove(1)
m.remove_spare_bonding()
m.add_particle_bonding()
for insertion in range(nmon - 1):
head = None
tail = None
backbone_vector = np.array(find_last_backbone_vector(s, m))
for p, p_ in izip(s.particles[-1*m.particles.count:], m.particles):
p_.x = p.x + 3*backbone_vector[0]
p_.y = p.y + 3*backbone_vector[1]
p_.z = p.z + 3*backbone_vector[2]
n = m.copy()
if capped:
s.particles.remove(s.particles.count)
s.remove_spare_bonding()
s.add_particle_bonding()
if extra_bonds:
heads = []
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'head':
heads.append(p)
else:
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'head':
head = p
s.add(n, change_dim=False)
s.add_particle_bonding()
if extra_bonds:
tails = []
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'tail':
tails.append(p)
else:
for p in s.particles[-1*n.particles.count:]:
if p.linker == 'tail':
tail = p
for p in s.particles:
if not p.bonded_to:
print(p.tag)
if head and tail:
s.make_new_bonds(head, tail, f)
print('%s: %s/%s monomers added' % (strftime('%H:%M:%S'), insertion+2, nmon))
elif extra_bonds and len(heads) == len(tails):
for h, t in izip(heads, tails):
s.make_new_bonds(h, t, f)
print('%s: %s/%s monomers added' % (strftime('%H:%M:%S'), insertion+2, nmon))
else:
print('cannot find head and tail')
if sim is None:
sim = lmps.Simulation(s, name='relax_%03d' % (insertion+2), log='relax.log', **settings)
sim.add_md(ensemble='nve', limit=limit, **settings)
sim.add_min(**settings)
if isinstance(sim, lmps.Simulation):
sim.system = s
sim.name = 'relax_%03d' % (insertion+2)
sim.run(np=settings.get('np'))
s.unwrap()
if traj:
s.write_xyz('random_walk.xyz', append=True)
if unwrap:
s.wrap()
for p in s.particles:
if p not in s.molecules[p.molecule.tag].particles:
s.molecules[p.molecule.tag].particles.add(p)
s.write_lammps('polymer.lmps')
s.unwrap()
s.write_xyz('polymer.xyz')
return s