Dear Lammps User
I built the nanofiber (a crystal structure) in Material Studio Version 8 and tried to use msi2lmps.exe to transform it, but came across this error: ‘cannot read symmetry file’. I am sending the email to request a possible solution.
Has anyone used the moltemplate to build the crystal structure before?
best regards
cheng
Hi Cheng
At the same time you posted your question on the LAMMPS mailing
list, I received an email from someone else asking whether moltemplate
could be used to build a nanofiber of crystalline cellulose crystal Iβ
(Bcc). I don't know if you are also trying to simulate celluose Iβ,
or something else. Either way, the answer to both of your questions
is yes.
Whether moltemplate is the most convenient tool to use is another story.
Here is an extremely simple example which builds a crystal of
water ice (type "1h") using a small 8-molecule unit cell:
https://raw.githubusercontent.com/jewettaij/moltemplate/master/examples/all_atom/force_field_explicit_parameters/ice_crystal_SPCE/images/ice_rect8_crystal_3x2x2_LR.jpg
https://raw.githubusercontent.com/jewettaij/moltemplate/master/examples/all_atom/force_field_explicit_parameters/ice_crystal_SPCE/images/ice_rect8_unitcell.png
For input files, see:
https://github.com/jewettaij/moltemplate/tree/master/examples/all_atom/force_field_explicit_parameters/ice_crystal_SPCE
Moltemplate can build any crystalline system, but it has the
following limitations:
1) You can have molecules containing bonds located inside a crystal
unit cell, however bonds connecting atoms in different unit cells are
not generated automatically. This includes bonds connecting atoms
across periodic boundaries of the simulation box. (That later
limitation is currently shared with many other molecule builders,
including topotools I believe.) However it is not too hard add these
bonds yourself manually, because the atom naming system it uses makes
it easy to identify atoms within the crystal. (But you do have to add
these bonds manually.) For the case of a bundle of polymers (such as
celluose Iβ), this is particularly easy to build in moltemplate,
because there are only bonds along the polymer axis. (See details
below.)
2) Moltemplate does not know how to interpret atom types in PDB files
or infer bonds between atoms. (It can read atom coordinates, but you
still have to define the molecules, their atom types, and their bond
connectivity yourself manually.) Therefore I don't recommend that you
build the system in another program and import the PDB file into
moltemplate. Instead if you decide to use moltemplate, I recommend
that you define the unit cell yourself, and then use moltemplate
commands to replicate the unit cell to build a crystal, and then
define bonds to link atoms in different unit cells together (as
mentioned above).
3) For the reason listed above (2), moltempate is not a good tool for
building all-atom simulations containing complex biomolecules like
proteins. If you anticipate that you will eventually need to add
enzymes or and other proteins, don't use moltemplate.
4) Currently moltemplate only includes force field parameters for a
few all-atom force-fields: OPLSAA, GAFF, GAFF2, and COMPASS. Your
life will be much easier if you happen to be using one of these force
fields. For the COMPASS force field, moltemplate only includes the
limited set of parameters which have been published publicly. However
if you have obtained an .FRC file with a wider list of COMPASS
parameters, you can convert that force-field file into moltemplate
format using the "msifrc2lt.py" script (available on github within the
moltemplate download. Unfortunately, this script does not yet work
with PCFF or CVFF).
If I'm not mistaken, the monclinic unit cell for cellulose Iβ
contains 4 sugar rings (two different types of rings, both
duplicated), and has 46 atoms. This is a small enough system that I
can imagine the number of angles, dihedrals, and impropers for each
sugar ring might be small enough that you could manually enumerate all
of them and copy the parameters from the files you generated with
materials studio (or some other molecule builder), assuming they are
in a human readable format. Then you could create a new .LT file
containing an abbreviated force-field file which contains only the
parameters for the atoms that you care about in your small system,
similar to the "oplsaa_simple.lt" example located here:
http://moltemplate.org/examples/force_fields/alkane_chain_single/oplsaa_simple.lt
In this way, you could prepare a system using one of the
force-fields that is not yet included with moltemplate.
---- manually adding bonds between atoms in different unit cells ----
Generally:
Suppose you create a crystal using the following commands:
cells = new UnitCell
[10].move(17.0,0,0)
[5].move(0, 10, 6)
[4].move(0, 0, 14)
You could bond all of the "C1" atoms with from one cell with the "O"
atom from the next cell (in the +x direction) using this notation:
write("Data Bond List") {
$bond:b1 $atom:cells[0][0][0]/C1 $atom:cells[1][0][0]/O
$bond:b2 $atom:cells[1][0][0]/C1 $atom:cells[2][0][0]/O
$bond:b3 $atom:cells[2][0][0]/C1 $atom:cells[3][0][0]/O
# :
$bond:b9 $atom:cells[8][0][0]/C1 $atom:cells[9][0][0]/O
#bind the end of the polymer to the beginning (periodic boundaries):
$bond:b10 $atom:cells[9][0][0]/C1 $atom:cells[0][0][0]/O
THEN YOU HAVE TO REPEAT THIS FOR j,k in cells[i][j][k]
$bond:b10 $atom:cells[0][1][0]/C1 $atom:cells[1][1][0]/O
$bond:b11 $atom:cells[1][1][0]/C1 $atom:cells[2][1][0]/O
$bond:b12 $atom:cells[2][1][0]/C1 $atom:cells[3][1][0]/O
$bond:b13 $atom:cells[3][1][0]/C1 $atom:cells[9][1][0]/O
# :
#bind the end of the polymer to the beginning (periodic boundaries):
$bond:b19 $atom:cells[8][1][0]/C1 $atom:cells[9][1][0]/O
$bond:b20 $atom:cells[9][1][0]/C1 $atom:cells[0][1][0]/O
:
}
There are a lot of these bonds, and I typically write a short python
script with a for-loop to generate a text file containing these
commands which I copy into the .LT file manually. (Incredibly,
moltemplate does not support for-loops yet.)
You don't have to define angles and dihedrals which cross the unit
cell boundaries. Moltemplate will use whatever force-field rules you
have supplied to create angles and dihedrals for these bonds and look
up their parameters.
Note: The lines above only adds bonds along the X direction. If your
crystal has bonds between unit cells along the Y and Z directions, you
have to add bonds in those directions as well.
------ cellulose Iβ ------
This is not the easiest way to build cellulose Iβ. Moltemplate is
hierarchical. (http://moltemplate.org/examples/menger.html). Once
you define a larger object (like a polymer) you can replicate it,
without needing to redefine all of the bonds inside of it.
It's not clear to me whether, you need to use triclinic boundary
conditions for monoclinic cellulose Iβ. For simple hexagonal unit
cells, you can use rectangular boundary conditions. Life is always
easier if you can use rectangular boundary conditions, but you can
simulate triclinic periodic systems in LAMMPS if you need to. In the
examples below, I assume you are using rectangular boundary
conditions.
Probably the easiest way to build cellulose Iβ would be to define a
moltemplate object containing a single polymer chain that runs along
the x axis (for example) which runs the entire length of the unit cell
(and includes the bond connecting one end of the polymer to the
opposite end).
Then replicated it many times in the y and z directions (for
example), arranged in a 2D (hexagonal?) crystal. You can add offsets
in the X direction which alternate to give the staggered appearance
shown in the figure below:
https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcQxhSVDwV_DlzxtWLg-N-GgidoWPeSRUU1G0MONBQJlvF3GdCSsYw
(The procedure to do that is outlined at the very end of this email post.)
The most similar moltemplate example is probably the hexadecane
example. In that example, a monomer subunit (CH2):
https://raw.githubusercontent.com/jewettaij/moltemplate/master/examples/all_atom/force_field_OPLSAA/hexadecane/images/ch2_ry60_LR.jpg
..is replicated many times to create a polymer
https://raw.githubusercontent.com/jewettaij/moltemplate/master/examples/all_atom/force_field_OPLSAA/hexadecane/images/hexadecane_LR.jpg
..and this polymer is copied many times:
https://raw.githubusercontent.com/jewettaij/moltemplate/master/examples/all_atom/force_field_OPLSAA/hexadecane/images/hexadecane_12x12x2_t%3D0_LR.jpg
---- Details for cellulose Iβ ------
However this differs from cellulose Iβ, because in that case, you
would want the polymer to run the entire length of the simulation box,
and you would not have terminal groups at the ends of the polymer (the
CH3 groups at either end). Instead, you would add a bond connecting
one end of the polymer to the opposite end. Once you have defined
such a polymer, you could use the "new" command create a larger object
consisting of two instances of the same polymer, one of which is
shifted in the X and Y directions (for example) relative to the other.
(If I'm not mistaken, you need 2 polymers, because one of them is
shifted with respect to the other along the polymer axis.) Then
duplicate this 2-polymer object in the YZ plane to fill the volume of
the simulation box.
I don't know if this long post helped at all.
Good luck.
Andrew