I'm assuming that most of the bonds, angles and dihedrals for your
system are in the polymer, not the silicon crystal. If that is the
case, then can you not simply manually add bonds connecting the
polymer fragments together after you use the replicate command? (You
could do this by saving a restart file after you replicated the
system, converting the restart file to a data file using restart2data,
and then editing the data file to add a couple bonds to the "Bonds"
section.)
I believe its possible to write a separate code to define
all these parameters, but my structure would have ~ 70,000
atoms and I think it would be impractical to do so.
Indeed, it took me a while. Incidentally, you're welcome to use my
script if you want. I tried to make it very general. I'll be
announcing it here soon.
Take a look at http://moltemplate.org
With "moltemplate", you can define a molecule named "UnitCell" and
instantiate 3x3x3 copies of it to make a crystal. You can also make a
polymer that spans the 3x3x3 crystal and insert it into the crystal,
deleting atoms which overlap with the polymer. Then you can replicate
the entire system, and add bonds between opposite ends of the polymer
fragments.
Here is a quick and sloppy version of what I mean:
UnitCell0 {
# atomID molID atomType charge x y z
write("Data Atoms") {
$atom:Si1 $mol:/SiCrystal @atom:Si 0.0 0.00 0.00 0.00
$atom:Si2 $mol:/SiCrystal @atom:Si 0.0 0.50 0.00 0.50
$atom:Si3 $mol:/SiCrystal @atom:Si 0.0 0.25 0.25 0.75
$atom:Si4 $mol:/SiCrystal @atom:Si 0.0 0.75 0.25 0.25
$atom:Si5 $mol:/SiCrystal @atom:Si 0.0 0.00 0.00 0.50
$atom:Si6 $mol:/SiCrystal @atom:Si 0.0 0.50 0.50 0.00
$atom:Si7 $mol:/SiCrystal @atom:Si 0.0 0.25 0.75 0.25
$atom:Si8 $mol:/SiCrystal @atom:Si 0.0 0.75 0.75 0.75
}
# Coordinates are in reduced units.
# I'll use the "scale()" command later.
# (I probably got these coordinates wrong. I was in a hurry.)
write_once("Data Masses") {
@atom:Si 28.0
}
write_once("In Settings") {
# i j
pair_coeff @atom:Si @atom:Si CHOSE YOUR PARAMETERS
# Optional: put all the Si atoms in the same group
group silicon type @atom:Si
}
} # UnitCell0
# Create a new version of the unit cell whose coordinates are scaled up by 5.43
# (5.43 Angstroms is the length of a unit cell of Silicon I think.)
UnitCell = UnitCell0.scale(5.43)
# You can make a 3x3x3 array of UnitCells this way:
crystal = new UnitCell [3].move(5.43, 0.0, 0.0)
[3].move(0.0, 5.43, 0.0)
[3].move(0.0, 0.0, 5.43)
The vectors after the 3 move commands are the lattice-vectors in each
of the 3 directions.