mirror of https://github.com/lammps/lammps.git
import moltemplate 2.8.6 from andrew jewett
This commit is contained in:
Axel Kohlmeyer
parent
6bd5a3d69b
commit
7abc960d39
|
|
@ -0,0 +1,180 @@
|
|||
# from https://github.com/github/gitignore/blob/master/Python.gitignore
|
||||
|
||||
# Byte-compiled / optimized / DLL files
|
||||
__pycache__/
|
||||
*.py[cod]
|
||||
*$py.class
|
||||
|
||||
# C extensions
|
||||
*.so
|
||||
|
||||
# Distribution / packaging
|
||||
.Python
|
||||
env/
|
||||
build/
|
||||
develop-eggs/
|
||||
dist/
|
||||
downloads/
|
||||
eggs/
|
||||
.eggs/
|
||||
lib/
|
||||
lib64/
|
||||
parts/
|
||||
sdist/
|
||||
var/
|
||||
wheels/
|
||||
*.egg-info/
|
||||
.installed.cfg
|
||||
*.egg
|
||||
|
||||
# PyInstaller
|
||||
# Usually these files are written by a python script from a template
|
||||
# before PyInstaller builds the exe, so as to inject date/other infos into it.
|
||||
*.manifest
|
||||
*.spec
|
||||
|
||||
# Installer logs
|
||||
pip-log.txt
|
||||
pip-delete-this-directory.txt
|
||||
|
||||
# Unit test / coverage reports
|
||||
htmlcov/
|
||||
.tox/
|
||||
.coverage
|
||||
.coverage.*
|
||||
.cache
|
||||
nosetests.xml
|
||||
coverage.xml
|
||||
*,cover
|
||||
.hypothesis/
|
||||
|
||||
# Translations
|
||||
*.mo
|
||||
*.pot
|
||||
|
||||
# Django stuff:
|
||||
*.log
|
||||
local_settings.py
|
||||
|
||||
# Flask stuff:
|
||||
instance/
|
||||
.webassets-cache
|
||||
|
||||
# Scrapy stuff:
|
||||
.scrapy
|
||||
|
||||
# Sphinx documentation
|
||||
docs/_build/
|
||||
|
||||
# PyBuilder
|
||||
target/
|
||||
|
||||
# Jupyter Notebook
|
||||
.ipynb_checkpoints
|
||||
|
||||
# pyenv
|
||||
.python-version
|
||||
|
||||
# celery beat schedule file
|
||||
celerybeat-schedule
|
||||
|
||||
# dotenv
|
||||
.env
|
||||
|
||||
# virtualenv
|
||||
.venv/
|
||||
venv/
|
||||
ENV/
|
||||
|
||||
# Spyder project settings
|
||||
.spyderproject
|
||||
|
||||
# Rope project settings
|
||||
.ropeproject
|
||||
|
||||
# text-editor temporary files:
|
||||
*~
|
||||
|
||||
# misc rubbish
|
||||
deleteme*
|
||||
DELETEME*
|
||||
|
||||
######## files specific to moltemplate and lammps ########
|
||||
|
||||
# latex/bibtex temporary files for the moltemplate manual:
|
||||
moltemplate_manual*.aux
|
||||
moltemplate_manual*.bbl
|
||||
moltemplate_manual*.blg
|
||||
moltemplate_manual*.log
|
||||
moltemplate_manual*.out
|
||||
moltemplate_manual*.toc
|
||||
|
||||
######## files created by running LAMMPS: ########
|
||||
|
||||
log.lammps
|
||||
log.cite
|
||||
traj*.lammpstrj
|
||||
|
||||
######## files generated by running moltemplate: ########
|
||||
|
||||
system.data
|
||||
system.in
|
||||
system.in.init
|
||||
system.in.settings
|
||||
system.in.charges
|
||||
system.psf
|
||||
ttree_assignments.txt
|
||||
output_ttree/
|
||||
|
||||
# Sections from the LAMMPS data file generated by moltemplate.sh
|
||||
"Data Header"*
|
||||
"Data Atoms"*
|
||||
"Data Masses"*
|
||||
"Data Velocities"*
|
||||
"Data Bonds"*
|
||||
"Data Bond List"*
|
||||
"Data Bonds AtomId AtomId"*
|
||||
"Data Angles"*
|
||||
"Data Dihedrals"*
|
||||
"Data Impropers"*
|
||||
"Data Bond Coeffs"*
|
||||
"Data Angle Coeffs"*
|
||||
"Data Dihedral Coeffs"*
|
||||
"Data Improper Coeffs"*
|
||||
"Data Pair Coeffs"*
|
||||
"Data PairIJ Coeffs"*
|
||||
|
||||
# interactions-by-type (not id. This is not part of the LAMMPS standard.)
|
||||
"Data Charge By Bond"*
|
||||
"Data Bonds By Type"*
|
||||
"Data Angles By Type"*
|
||||
"Data Dihedrals By Type"*
|
||||
"Data Impropers By Type"*
|
||||
|
||||
# class2 data sections
|
||||
"Data BondBond Coeffs"*
|
||||
"Data BondAngle Coeffs"*
|
||||
"Data MiddleBondTorsion Coeffs"*
|
||||
"Data EndBondTorsion Coeffs"*
|
||||
"Data AngleTorsion Coeffs"*
|
||||
"Data AngleAngleTorsion Coeffs"*
|
||||
"Data BondBond13 Coeffs"*
|
||||
"Data AngleAngle Coeffs"*
|
||||
|
||||
# sections for non-point-like particles:
|
||||
"Data Ellipsoids"*
|
||||
"Data Lines"*
|
||||
"Data Triangles"*
|
||||
|
||||
# periodic boundary conditions
|
||||
"Data Boundary"*
|
||||
|
||||
# Sections from the LAMMPS input script(s) generated by moltemplate.sh
|
||||
|
||||
"In Init"*
|
||||
"In Settings"*
|
||||
"In Coords"*
|
||||
"In Charges"*
|
||||
#temporary file created by moltemplate.sh for storing coordinates
|
||||
tmp_atom_coords.dat
|
||||
|
||||
|
|
@ -1,9 +1,11 @@
|
|||
[](https://travis-ci.org/jewettaij/moltemplate.svg?branch=master)
|
||||
|
||||
Moltemplate
|
||||
===========
|
||||
|
||||
## Description
|
||||
|
||||
Moltemplate is a cross-platform text-based molecule builder for LAMMPS.
|
||||
Moltemplate is a *general* cross-platform text-based molecule builder for **LAMMPS** and **ESPResSo**. Moltemplate was intended for building custom coarse-grained molecular models, but it can be used to prepare realistic all-atom simulations as well. It currently supports the **OPLS**, **COMPASS**, **AMBER**(GAFF,GAFF2), **MARTINI**, **SDK**, **LOPLS**(2015), and **TraPPE**(1998) force fields, and includes approximately 40 examples. (New force fields and examples are added continually by users.)
|
||||
|
||||
## Typical usage
|
||||
|
||||
|
|
@ -45,7 +47,7 @@ Make sure that your default pip install bin directory is in your PATH. (This is
|
|||
|
||||
pip uninstall moltemplate
|
||||
|
||||
If you continue to run into difficulty, try installing moltemplate into a temporary virtual environment by installing "virtualenv", downloading moltemplate (to "~/moltemplate" in the example below), and running these commands:
|
||||
If you continue to run into difficulty, try installing moltemplate into a temporary virtual environment by installing "*virtualenv*", downloading moltemplate (to "~/moltemplate" in the example below), and running these commands:
|
||||
|
||||
cd ~/moltemplate
|
||||
virtualenv venv
|
||||
|
|
@ -53,7 +55,9 @@ If you continue to run into difficulty, try installing moltemplate into a tempor
|
|||
pip install .
|
||||
#(now do something useful with moltemplate...)
|
||||
|
||||
(You will have to "run source ~/moltemplate/venv/bin/activate" beforehand every time you want to run moltemplate.) If all this fails, then try installing moltemplate by manually updating your \$PATH environment variable. Instructions for doing that are included below.
|
||||
(You will have to "run source ~/moltemplate/venv/bin/activate" beforehand every time you want to run moltemplate.
|
||||
The *virtualenv* tool is
|
||||
[explained in detail here](http://docs.python-guide.org/en/latest/dev/virtualenvs/)) If all this fails, then try installing moltemplate by manually updating your \$PATH environment variable. Instructions for doing that are included below.
|
||||
|
||||
## Manual installation:
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
Before Width: | Height: | Size: 3.4 KiB After Width: | Height: | Size: 3.4 KiB |
|
Before Width: | Height: | Size: 15 KiB After Width: | Height: | Size: 15 KiB |
|
Before Width: | Height: | Size: 39 KiB After Width: | Height: | Size: 39 KiB |
|
Before Width: | Height: | Size: 24 KiB After Width: | Height: | Size: 24 KiB |
|
Before Width: | Height: | Size: 16 KiB After Width: | Height: | Size: 16 KiB |
|
Before Width: | Height: | Size: 12 KiB After Width: | Height: | Size: 12 KiB |
File diff suppressed because it is too large
Load Diff
|
|
@ -13,22 +13,22 @@ Explanation:
|
|||
|
||||
Usage:
|
||||
|
||||
genpoly_lt.py \\
|
||||
[-bond btype a1 a2] \\
|
||||
[-helix deltaphi] \\
|
||||
[-axis x,y,z] \\
|
||||
[-circular yes/no/connected] \\
|
||||
genpoly_lt.py \
|
||||
[-bond btype a1 a2] \
|
||||
[-helix deltaphi] \
|
||||
[-axis x,y,z] \
|
||||
[-circular yes/no/connected] \
|
||||
[-dir-indices ia ib] \
|
||||
[-angle atype a1 a2 a3 i1 i2 i3] \\
|
||||
[-dihedral dtype a1 a2 a3 a4 i1 i2 i3 i4] \\
|
||||
[-improper itype a1 a2 a3 a4 i1 i2 i3 i4] \\
|
||||
[-monomer-name mname] \\
|
||||
[-sequence sequence.txt] \\
|
||||
[-polymer-name pname] \\
|
||||
[-inherits ForceFieldObject] \\
|
||||
[-header "import \"monomer.lt\""] \\
|
||||
[-cuts cuts.txt] \\
|
||||
[-box paddingX,paddingY,paddingZ] \\
|
||||
[-angle atype a1 a2 a3 i1 i2 i3] \
|
||||
[-dihedral dtype a1 a2 a3 a4 i1 i2 i3 i4] \
|
||||
[-improper itype a1 a2 a3 a4 i1 i2 i3 i4] \
|
||||
[-monomer-name mname] \
|
||||
[-sequence sequence.txt] \
|
||||
[-polymer-name pname] \
|
||||
[-inherits ForceFieldObject] \
|
||||
[-header "import monomer.lt"] \
|
||||
[-cuts cuts.txt] \
|
||||
[-box paddingX,paddingY,paddingZ] \
|
||||
< coords.raw > polymer.lt
|
||||
|
||||
Arguments (optional):
|
||||
|
|
|
|||
|
|
@ -20,18 +20,18 @@ Benzene inherits GAFF {
|
|||
|
||||
# atomID molID atomType charge X Y Z
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:ca -0.115 5.274 1.999 -8.568
|
||||
$atom:C2 $mol @atom:ca -0.115 6.627 2.018 -8.209
|
||||
$atom:C3 $mol @atom:ca -0.115 7.366 0.829 -8.202
|
||||
$atom:C4 $mol @atom:ca -0.115 6.752 -0.379 -8.554
|
||||
$atom:C5 $mol @atom:ca -0.115 5.399 -0.398 -8.912
|
||||
$atom:C6 $mol @atom:ca -0.115 4.660 0.791 -8.919
|
||||
$atom:H11 $mol @atom:ha 0.115 4.704 2.916 -8.573
|
||||
$atom:H21 $mol @atom:ha 0.115 7.101 2.950 -7.938
|
||||
$atom:H31 $mol @atom:ha 0.115 8.410 0.844 -7.926
|
||||
$atom:H41 $mol @atom:ha 0.115 7.322 -1.296 -8.548
|
||||
$atom:H51 $mol @atom:ha 0.115 4.925 -1.330 -9.183
|
||||
$atom:H61 $mol @atom:ha 0.115 3.616 0.776 -9.196
|
||||
$atom:C1 $mol @atom:ca -0.115 -0.739 1.189 -0.00733
|
||||
$atom:C2 $mol @atom:ca -0.115 0.614 1.208 0.35167
|
||||
$atom:C3 $mol @atom:ca -0.115 1.353 0.019 0.35867
|
||||
$atom:C4 $mol @atom:ca -0.115 0.739 -1.189 0.00667
|
||||
$atom:C5 $mol @atom:ca -0.115 -0.614 -1.208 -0.35133
|
||||
$atom:C6 $mol @atom:ca -0.115 -1.353 -0.019 -0.35833
|
||||
$atom:H11 $mol @atom:ha 0.115 -1.309 2.106 -0.01233
|
||||
$atom:H21 $mol @atom:ha 0.115 1.088 2.14 0.62267
|
||||
$atom:H31 $mol @atom:ha 0.115 2.397 0.034 0.63467
|
||||
$atom:H41 $mol @atom:ha 0.115 1.309 -2.106 0.01267
|
||||
$atom:H51 $mol @atom:ha 0.115 -1.088 -2.14 -0.62233
|
||||
$atom:H61 $mol @atom:ha 0.115 -2.397 -0.034 -0.63533
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
|
|
|
|||
|
|
@ -60,8 +60,8 @@ d) Try entering these commands:
|
|||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc box -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc wrap -compound res -all -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
pbc box -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 36 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 127 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 51 KiB |
|
|
@ -15,19 +15,20 @@ import "gaff.lt"
|
|||
|
||||
Benzene inherits GAFF {
|
||||
|
||||
# atomID molID atomType charge X Y Z
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:ca 0.115 5.274 1.999 -8.568
|
||||
$atom:C2 $mol @atom:ca 0.115 6.627 2.018 -8.209
|
||||
$atom:C3 $mol @atom:ca 0.115 7.366 0.829 -8.202
|
||||
$atom:C4 $mol @atom:ca 0.115 6.752 -0.379 -8.554
|
||||
$atom:C5 $mol @atom:ca 0.115 5.399 -0.398 -8.912
|
||||
$atom:C6 $mol @atom:ca 0.115 4.660 0.791 -8.919
|
||||
$atom:H11 $mol @atom:ha -0.115 4.704 2.916 -8.573
|
||||
$atom:H21 $mol @atom:ha -0.115 7.101 2.950 -7.938
|
||||
$atom:H31 $mol @atom:ha -0.115 8.410 0.844 -7.926
|
||||
$atom:H41 $mol @atom:ha -0.115 7.322 -1.296 -8.548
|
||||
$atom:H51 $mol @atom:ha -0.115 4.925 -1.330 -9.183
|
||||
$atom:H61 $mol @atom:ha -0.115 3.616 0.776 -9.196
|
||||
$atom:C1 $mol @atom:ca -0.115 -0.739 1.189 -0.00733
|
||||
$atom:C2 $mol @atom:ca -0.115 0.614 1.208 0.35167
|
||||
$atom:C3 $mol @atom:ca -0.115 1.353 0.019 0.35867
|
||||
$atom:C4 $mol @atom:ca -0.115 0.739 -1.189 0.00667
|
||||
$atom:C5 $mol @atom:ca -0.115 -0.614 -1.208 -0.35133
|
||||
$atom:C6 $mol @atom:ca -0.115 -1.353 -0.019 -0.35833
|
||||
$atom:H11 $mol @atom:ha 0.115 -1.309 2.106 -0.01233
|
||||
$atom:H21 $mol @atom:ha 0.115 1.088 2.14 0.62267
|
||||
$atom:H31 $mol @atom:ha 0.115 2.397 0.034 0.63467
|
||||
$atom:H41 $mol @atom:ha 0.115 1.309 -2.106 0.01267
|
||||
$atom:H51 $mol @atom:ha 0.115 -1.088 -2.14 -0.62233
|
||||
$atom:H61 $mol @atom:ha 0.115 -2.397 -0.034 -0.63533
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
|
|
|
|||
|
|
@ -4,26 +4,25 @@ import "benzene.lt" # <- defines the "Benzene" molecule type.
|
|||
|
||||
# Periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0.0 80.00 xlo xhi
|
||||
0.0 80.00 ylo yhi
|
||||
0.0 80.00 zlo zhi
|
||||
0.0 48.00 xlo xhi
|
||||
0.0 48.00 ylo yhi
|
||||
0.0 48.00 zlo zhi
|
||||
}
|
||||
|
||||
# Create 1000 ethylenes and 500 benzenes
|
||||
# Create 216 ethylenes and 108 benzenes
|
||||
|
||||
ethylenes = new Ethylene[10].move(8.0, 0, 0)
|
||||
[10].move(0, 8.0, 0)
|
||||
[10].move(0, 0, 8.0)
|
||||
ethylenes = new Ethylene[6].move(8.0, 0, 0)
|
||||
[6].move(0, 8.0, 0)
|
||||
[6].move(0, 0, 8.0)
|
||||
|
||||
benzenes = new Benzene[10].move(8.0, 0, 0)
|
||||
[10].move(0, 8.0, 0)
|
||||
[5].move(0, 0, 16.0)
|
||||
benzenes = new Benzene[6].move(8.0, 0, 0)
|
||||
[6].move(0, 8.0, 0)
|
||||
[3].move(0, 0, 16.0)
|
||||
|
||||
# Now shift the positions of all of the benzene molecules,
|
||||
# to reduce the chance that they overlap with the ethylene molecules.
|
||||
|
||||
benzenes[*][*][*].move(4.0, 4.0, 4.0)
|
||||
|
||||
# Note: There is also an example in the OPLSAA directory which shows how to
|
||||
# generate the coordinates using PACKMOL. That allows us to omit all of
|
||||
# the coordinates and .move() commands. (This works with AMBER/GAFF too.)
|
||||
# Note: There is also an example which shows how to generate the coordinates
|
||||
# using PACKMOL. (That allows us to omit the coordinates and .move() commands.)
|
||||
|
|
|
|||
|
|
@ -35,7 +35,7 @@ include system.in.settings
|
|||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
timestep 2.0
|
||||
dump 1 all custom 5000 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
|
||||
thermo 500
|
||||
|
|
|
|||
|
|
@ -0,0 +1,38 @@
|
|||
# -------- WARNING: --------
|
||||
|
||||
This directory contains some examples of all-atom simulations using the COMPASS
|
||||
force field.
|
||||
|
||||
This software is experimental, and the force-fields and equilbration protocols
|
||||
have not been tested carefully by me. There is no gaurantee that simulations
|
||||
prepared using moltemplate will reproduce the behavior of other MD codes.
|
||||
|
||||
The moltemplate implementation of COMPASS currently relies on the same
|
||||
incomplete force-field file that "msi2lmp" uses ("compass_published.frc").
|
||||
Unfortunately this means that many force field parameters and some atom types
|
||||
(such as sp2-carbons) have not (yet) been publicly released and are not
|
||||
available.
|
||||
|
||||
Currently I recommend that users should run the "cleanup_moltemplate.sh"
|
||||
script after running "moltemplate.sh system.lt". Then manually check that
|
||||
the "system.in.settings" and "system.in.charges" files which remain
|
||||
make sense. Specifically, you must check that the angle_coeff,
|
||||
dihedral_coeff, bond_coeff commands are not full of zeros (in places
|
||||
where they should not be zero. This is another consequence of the
|
||||
fact that the .FRC files I mentioned above are incomplete.) It's a
|
||||
good idea to also check that the charges in the "system.in.charges"
|
||||
file seem reasonable (ie. not all zeros). (There is a list of
|
||||
warnings at the end of the "compass_published.lt" file. You can check
|
||||
to see if any of the bonds in your system are covered by these
|
||||
warnings.) Later on hopefully I'll add some automated way to warn
|
||||
users when these problems arise, but now you should check for them
|
||||
manually.
|
||||
|
||||
# -------- REQUEST FOR HELP: --------
|
||||
|
||||
If you notice a problem with these examples, please report it.
|
||||
Peer-review is the only way to improve this software (or any software).
|
||||
Other suggestions are also welcome!
|
||||
|
||||
(Contact jewett.aij@gmail.com, 2017-10-03)
|
||||
|
||||
|
|
@ -0,0 +1,37 @@
|
|||
This example is a simple simulation of a long alkane chain,
|
||||
in a vacuum at room temperature using the COMPASS force field.
|
||||
|
||||
1) Create the "system.data", "system.in.init", and "system.in.settings"
|
||||
files which LAMMPS will read by running:
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
|
||||
2) Run LAMMPS in this order:
|
||||
|
||||
lmp_mpi -i run.in.min # minimize the energy (to avoid atom overlap) before...
|
||||
lmp_mpi -i run.in.nvt # running the simulation at constant temperature
|
||||
|
||||
(The name of the LAMMPS executable, eg "lmp_mpi", may vary.)
|
||||
|
||||
---- Details ----
|
||||
|
||||
The "Alkane50" molecule, as well as the "CH2", and "CH3" monomers it contains
|
||||
use the COMPASS force-field. This means that when we define these molecules,
|
||||
we only specify the atom names, bond list, and coordinates.
|
||||
We do not have to list the atom charges, angles, dihedrals, or impropers.
|
||||
The rules for creating atomic charge and angle topology are contained in
|
||||
the "compass_published.lt" file created by step 3) above. The "ch2group.lt",
|
||||
"ch3group.lt", and "alkane50.lt" files all refer to "compass_published.lt",
|
||||
(as well as the "COMPASS" force-field object which it defines). Excerpt:
|
||||
|
||||
import "compass_published.lt"
|
||||
CH2 inherits COMPASS { ...
|
||||
CH3 inherits COMPASS { ...
|
||||
Alkane50 inherits COMPASS { ...
|
||||
|
||||
Alternatively, you can manually define a list of angles, dihedrals, and
|
||||
improper interactions in these files, instead of asking the force-field
|
||||
to generate them for you. You can also specify some of the angles and
|
||||
dihedrals explicitly, and let the force-field handle the rest.
|
||||
(Many of the examples which come with moltemplate do this.)
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
|
||||
# Note: By default, the system.data and system.in.settings files contain
|
||||
# extra information for atoms defined in OPLSAA which you are not using
|
||||
# in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
|
||||
cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
# --- Running LAMMPS ---
|
||||
#
|
||||
# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
|
||||
# input scripts which link to the input scripts and data files
|
||||
# you hopefully have created earlier with moltemplate.sh:
|
||||
# system.in.init, system.in.settings, system.data
|
||||
# If not, carry out the instructions in "README_setup.sh".
|
||||
#
|
||||
# -- Instructions: --
|
||||
# If "lmp_mpi" is the name of the command you use to invoke lammps,
|
||||
# then you would run lammps on these files this way:
|
||||
|
||||
|
||||
lmp_mpi -i run.in.min # minimization
|
||||
lmp_mpi -i run.in.nvt # minimization and simulation at constant volume
|
||||
|
||||
#(Note: The constant volume simulation lacks pressure equilibration. These are
|
||||
# completely separate simulations. The results of the constant pressure
|
||||
# simulation might be ignored when beginning the simulation at constant
|
||||
# volume. (This is because restart files in LAMMPS don't always work,
|
||||
# and I was spending a lot of time trying to convince people it was a
|
||||
# LAMMPS bug, instead of a moltemplate bug, so I disabled restart files.)
|
||||
# Read the "run.in.nvt" file to find out how to use the "read_restart"
|
||||
# command to load the results of the pressure-equilibration simulation,
|
||||
# before beginning a constant-volume run.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.npt
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# Optional:
|
||||
# To check for missing angle,dihedral params run moltemplate this way instead:
|
||||
# moltemplate.sh -checkff system.lt
|
||||
|
||||
|
||||
# Moltemplate generates various files with names ending in *.in* and *.data.
|
||||
# Move them to the directory where you plan to run LAMMPS (in this case "../")
|
||||
mv -f system.data system.in* ../
|
||||
|
||||
# Optional:
|
||||
# The "./output_ttree/" directory is full of temporary files generated by
|
||||
# moltemplate. They can be useful for debugging, but are usually thrown away.
|
||||
rm -rf output_ttree/
|
||||
|
||||
cd ../
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional:
|
||||
# Note: The system.data and system.in.settings files contain extra information
|
||||
# for atoms defined in OPLSAA which you are not using in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
#
|
||||
# cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
|
||||
------- To view a lammps trajectory in VMD --------
|
||||
|
||||
|
||||
1) Build a PSF file for use in viewing with VMD.
|
||||
|
||||
This step works with VMD 1.9 and topotools 1.2.
|
||||
(Older versions, like VMD 1.8.6, don't support this.)
|
||||
|
||||
|
||||
a) Start VMD
|
||||
b) Menu Extensions->Tk Console
|
||||
c) Enter:
|
||||
|
||||
(I assume that the the DATA file is called "system.data")
|
||||
|
||||
topo readlammpsdata system.data full
|
||||
animate write psf system.psf
|
||||
|
||||
2)
|
||||
|
||||
Later, to Load a trajectory in VMD:
|
||||
|
||||
Start VMD
|
||||
Select menu: File->New Molecule
|
||||
-Browse to select the PSF file you created above, and load it.
|
||||
(Don't close the window yet.)
|
||||
-Browse to select the trajectory file.
|
||||
If necessary, for "file type" select: "LAMMPS Trajectory"
|
||||
Load it.
|
||||
|
||||
---- A note on trajectory format: -----
|
||||
If the trajectory is a DUMP file, then make sure the it contains the
|
||||
information you need for pbctools (see below. I've been using this
|
||||
command in my LAMMPS scripts to create the trajectories:
|
||||
|
||||
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
It's a good idea to use an atom_style which supports molecule-ID numbers
|
||||
so that you can assign a molecule-ID number to each atom. (I think this
|
||||
is needed to wrap atom coordinates without breaking molecules in half.)
|
||||
|
||||
Of course, you don't have to save your trajectories in DUMP format,
|
||||
(other formats like DCD work fine) I just mention dump files
|
||||
because these are the files I'm familiar with.
|
||||
|
||||
3) ----- Wrap the coordinates to the unit cell
|
||||
(without cutting the molecules in half)
|
||||
|
||||
a) Start VMD
|
||||
b) Load the trajectory in VMD (see above)
|
||||
c) Menu Extensions->Tk Console
|
||||
d) Try entering these commands:
|
||||
|
||||
pbc wrap -compound res -all
|
||||
pbc box
|
||||
|
||||
----- Optional ----
|
||||
Sometimes the solvent or membrane obscures the view of the solute.
|
||||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc box -shiftcenterrel {0.0 0.15 0.0}
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
Alternately if you have a solute whose atoms are all of type 1,
|
||||
then you can also try this to center the box around it:
|
||||
|
||||
pbc wrap -sel type=1 -all -centersel type=2 -center com
|
||||
|
||||
4)
|
||||
You should check if your periodic boundary conditions are too small.
|
||||
To do that:
|
||||
select Graphics->Representations menu option
|
||||
click on the "Periodic" tab, and
|
||||
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
|
||||
|
||||
5) Optional: If you like, change the atom types in the PSF file so
|
||||
that VMD recognizes the atom types, use something like:
|
||||
|
||||
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
|
||||
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
|
||||
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
|
||||
|
||||
(If you do this, it might effect step 2 above.)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 4.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.8 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 50 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 59 KiB |
|
|
@ -0,0 +1,137 @@
|
|||
# This is a simple example showing how to build a long polymer
|
||||
# (in this case, an alkane chain). I split the
|
||||
# hexadecane molecule into individual CH2 and CH3 monomers.
|
||||
# I defined it this way so that you can easily modify
|
||||
# it to change the length of the alkane chain.
|
||||
|
||||
|
||||
import "ch2group.lt" # load the definition of the "CH2" object
|
||||
import "ch3group.lt" # load the definition of the "CH3" object
|
||||
|
||||
|
||||
|
||||
Alkane50 inherits COMPASS {
|
||||
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
|
||||
# This is a long polymer consisting of 48 CH2 groups and 2 CH3 end-caps.
|
||||
# Rather than create them one-by-one, I decided to create them all
|
||||
# using a single "new" command. Later, I can modify this array.
|
||||
|
||||
# Create an array of 50 "CH2" objects distributed along the X axis
|
||||
|
||||
monomers = new CH2 [50].rot(180,1,0,0).move(1.2533223,0,0)
|
||||
|
||||
# NOTE: the ".rot(180,1,0,0).move(1.2533223,0,0)" means that each
|
||||
# successive monomer is rotated 180 degrees (with respect to the previous
|
||||
# monomer), and then moved 1.2533223 Angstroms down the X axis.
|
||||
# Alternately, if you are reading the coordinates from a file, you don't have
|
||||
# to indicate the position & orientation of each monomer. In that case, use:
|
||||
# monomers = new CH2 [50]
|
||||
|
||||
|
||||
# ---- Now, modify the ends: ---
|
||||
# Delete the CH2 groups at the beginning and end, and replace them with CH3.
|
||||
|
||||
delete monomers[0]
|
||||
delete monomers[49]
|
||||
|
||||
monomers[0] = new CH3
|
||||
monomers[49] = new CH3
|
||||
|
||||
# Move the CH3 groups to the correct location at either end of the chain:
|
||||
|
||||
#monomers[0].move(0,0,0) # <--(this monomer is already in the correct place)
|
||||
monomers[49].rot(180.0,0,0,1).move(61.4127927,0,0) #61.4127927=49*1.2533223
|
||||
|
||||
|
||||
## NOTE: Alternately, you can define the polymer without deleting the ends:
|
||||
# monomers[0] = new CH3
|
||||
# monomers[1-48] = new CH2[48].rot(180,1,0,0).move(1.2533223,0,0)
|
||||
## Note: monomers[0] and monomers[1] overlap, so we move 1-48 to make room:
|
||||
# monomers[1-48].rot(180,1,0,0).move(1.2533223,0,0) # move many monomers
|
||||
## Now add the final monomer at the end:
|
||||
# monomers[49] = new CH3.rot(180.0,0,0,1).move(61.4127927,0,0)
|
||||
#
|
||||
## NOTE: Alternately, you can read the coordinates from a file.
|
||||
## In that case, you can use simpler commands:
|
||||
# monomers[0] = new CH3
|
||||
# monomers[1-48] = new CH2[48]
|
||||
# monomers[49] = new CH3
|
||||
|
||||
|
||||
|
||||
# Now add a list of bonds connecting the carbon atoms together:
|
||||
# (Angles, dihedrals, impropers will be automatically added later.)
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:b1 $atom:monomers[0]/C $atom:monomers[1]/C
|
||||
$bond:b2 $atom:monomers[1]/C $atom:monomers[2]/C
|
||||
$bond:b3 $atom:monomers[2]/C $atom:monomers[3]/C
|
||||
$bond:b4 $atom:monomers[3]/C $atom:monomers[4]/C
|
||||
$bond:b5 $atom:monomers[4]/C $atom:monomers[5]/C
|
||||
$bond:b6 $atom:monomers[5]/C $atom:monomers[6]/C
|
||||
$bond:b7 $atom:monomers[6]/C $atom:monomers[7]/C
|
||||
$bond:b8 $atom:monomers[7]/C $atom:monomers[8]/C
|
||||
$bond:b9 $atom:monomers[8]/C $atom:monomers[9]/C
|
||||
$bond:b10 $atom:monomers[9]/C $atom:monomers[10]/C
|
||||
$bond:b11 $atom:monomers[10]/C $atom:monomers[11]/C
|
||||
$bond:b12 $atom:monomers[11]/C $atom:monomers[12]/C
|
||||
$bond:b13 $atom:monomers[12]/C $atom:monomers[13]/C
|
||||
$bond:b14 $atom:monomers[13]/C $atom:monomers[14]/C
|
||||
$bond:b15 $atom:monomers[14]/C $atom:monomers[15]/C
|
||||
$bond:b16 $atom:monomers[15]/C $atom:monomers[16]/C
|
||||
$bond:b17 $atom:monomers[16]/C $atom:monomers[17]/C
|
||||
$bond:b18 $atom:monomers[17]/C $atom:monomers[18]/C
|
||||
$bond:b19 $atom:monomers[18]/C $atom:monomers[19]/C
|
||||
$bond:b20 $atom:monomers[19]/C $atom:monomers[20]/C
|
||||
$bond:b21 $atom:monomers[20]/C $atom:monomers[21]/C
|
||||
$bond:b22 $atom:monomers[21]/C $atom:monomers[22]/C
|
||||
$bond:b23 $atom:monomers[22]/C $atom:monomers[23]/C
|
||||
$bond:b24 $atom:monomers[23]/C $atom:monomers[24]/C
|
||||
$bond:b25 $atom:monomers[24]/C $atom:monomers[25]/C
|
||||
$bond:b26 $atom:monomers[25]/C $atom:monomers[26]/C
|
||||
$bond:b27 $atom:monomers[26]/C $atom:monomers[27]/C
|
||||
$bond:b28 $atom:monomers[27]/C $atom:monomers[28]/C
|
||||
$bond:b29 $atom:monomers[28]/C $atom:monomers[29]/C
|
||||
$bond:b30 $atom:monomers[29]/C $atom:monomers[30]/C
|
||||
$bond:b31 $atom:monomers[30]/C $atom:monomers[31]/C
|
||||
$bond:b32 $atom:monomers[31]/C $atom:monomers[32]/C
|
||||
$bond:b33 $atom:monomers[32]/C $atom:monomers[33]/C
|
||||
$bond:b34 $atom:monomers[33]/C $atom:monomers[34]/C
|
||||
$bond:b35 $atom:monomers[34]/C $atom:monomers[35]/C
|
||||
$bond:b36 $atom:monomers[35]/C $atom:monomers[36]/C
|
||||
$bond:b37 $atom:monomers[36]/C $atom:monomers[37]/C
|
||||
$bond:b38 $atom:monomers[37]/C $atom:monomers[38]/C
|
||||
$bond:b39 $atom:monomers[38]/C $atom:monomers[39]/C
|
||||
$bond:b40 $atom:monomers[39]/C $atom:monomers[40]/C
|
||||
$bond:b41 $atom:monomers[40]/C $atom:monomers[41]/C
|
||||
$bond:b42 $atom:monomers[41]/C $atom:monomers[42]/C
|
||||
$bond:b43 $atom:monomers[42]/C $atom:monomers[43]/C
|
||||
$bond:b44 $atom:monomers[43]/C $atom:monomers[44]/C
|
||||
$bond:b45 $atom:monomers[44]/C $atom:monomers[45]/C
|
||||
$bond:b46 $atom:monomers[45]/C $atom:monomers[46]/C
|
||||
$bond:b47 $atom:monomers[46]/C $atom:monomers[47]/C
|
||||
$bond:b48 $atom:monomers[47]/C $atom:monomers[48]/C
|
||||
$bond:b49 $atom:monomers[48]/C $atom:monomers[49]/C
|
||||
}
|
||||
|
||||
} # Alkane50
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163316030377
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.8924307629540046
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.6310438442242609
|
||||
|
|
@ -0,0 +1,76 @@
|
|||
# This file contains a definition for the "CH2" molecular subunit.
|
||||
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
|
||||
|
||||
# Then define "CH2":
|
||||
|
||||
|
||||
CH2 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
}
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
} # CH2
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH2.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
# This file contains a definition for the "CH3" molecular subunit.
|
||||
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
|
||||
|
||||
# Then define "CH3":
|
||||
|
||||
|
||||
CH3 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
$atom:H3 $mol:... @atom:h1 0.0 -0.892431 -0.631044 0.000000
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
$bond:CH3 $atom:C $atom:H3
|
||||
}
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
} # CH3
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH3.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
import "alkane50.lt" # Defines the "Alkane50" molecule
|
||||
|
||||
|
||||
polymer = new Alkane50
|
||||
|
||||
|
||||
# Specify the size of the world the polymer lives in:
|
||||
write_once("Data Boundary") {
|
||||
0.0 72.0 xlo xhi
|
||||
0.0 72.0 ylo yhi
|
||||
0.0 72.0 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
###############################################################################
|
||||
# Note: If you want to create multiple polymers, and/or mix them with other
|
||||
# molecules, just add more "new" commands, for example:
|
||||
# polymer1 = new Alkane50.move(0,0,10)
|
||||
# polymer2 = new Alkane50.move(0,0,20)
|
||||
# :
|
||||
# ...or use array notation, for example:
|
||||
# polymers = new Alkane50[20].move(0,0,10)
|
||||
#
|
||||
# Note: Multidimensional arrays can be used to fill a planar region or a volume
|
||||
# polymers = new Alkane50 [4].move(0, 0, 30.0)
|
||||
# [4].move(0, 30.0, 0)
|
||||
# [2].move(70.0, 0, 0)
|
||||
|
|
@ -0,0 +1,37 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
#
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data "system.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
include "system.in.charges"
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
|
||||
# -- minimization protocol --
|
||||
|
||||
# Note: The minimization step is not necessary in this example. However
|
||||
# in general, it's always a good idea to minimize the system beforehand.
|
||||
|
||||
thermo 50
|
||||
dump 1 all custom 50 traj_min.lammpstrj id mol type x y z ix iy iz
|
||||
minimize 1.0e-4 1.0e-6 100000 400000
|
||||
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also. I prefer "write_data" and "read_data".)
|
||||
|
||||
write_data system_after_min.data
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,38 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
# 2) You must minimize the system beforehand by using "run.in.min".
|
||||
# This will create the file "system_after_min.data" which this file reads.
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
# Read the coordinates generated by an earlier simulation
|
||||
|
||||
read_data "system_after_min.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
include "system.in.charges"
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 1000 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
|
||||
thermo 500
|
||||
#thermo_modify flush yes
|
||||
|
||||
run 1000000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -0,0 +1,37 @@
|
|||
This example is a simple simulation of many short alkane chains (butane) in
|
||||
a box near the boiling point at atmospheric pressure.(Please read "WARNING.TXT")
|
||||
|
||||
1) Create the "system.data", "system.in.init", and "system.in.settings"
|
||||
files which LAMMPS will read by running:
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
|
||||
2) Run LAMMPS in this order:
|
||||
|
||||
lmp_mpi -i run.in.npt # running the simulation at constant pressure
|
||||
lmp_mpi -i run.in.nvt # running the simulation at constant temperature
|
||||
|
||||
(The name of the LAMMPS executable, eg "lmp_mpi", may vary.)
|
||||
|
||||
---- Details ----
|
||||
|
||||
The "Butane" molecule, as well as the "CH2", and "CH3" monomers it contains
|
||||
use the COMPASS force-field. This means that when we define these molecules,
|
||||
we only specify the atom names, bond list, and coordinates.
|
||||
We do not have to list the atom charges, angles, dihedrals, or impropers.
|
||||
The rules for creating atomic charge and angle topology are contained in
|
||||
the "compass_published.lt" file created by step 3) above. The "ch2group.lt",
|
||||
"ch3group.lt", and "butane.lt" files all refer to "compass_published.lt",
|
||||
(as well as the "COMPASS" force-field object which it defines). Excerpt:
|
||||
|
||||
import "compass_published.lt"
|
||||
CH2 inherits COMPASS { ...
|
||||
CH3 inherits COMPASS { ...
|
||||
Butane inherits COMPASS { ...
|
||||
|
||||
Alternatively, you can manually define a list of angles, dihedrals, and
|
||||
improper interactions in these files, instead of asking the force-field
|
||||
to generate them for you. You can also specify some of the angles and
|
||||
dihedrals explicitly, and let the force-field handle the rest.
|
||||
(Many of the examples which come with moltemplate do this.)
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
|
||||
# Note: By default, the system.data and system.in.settings files contain
|
||||
# extra information for atoms defined in OPLSAA which you are not using
|
||||
# in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
|
||||
cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
# --- Running LAMMPS ---
|
||||
#
|
||||
# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
|
||||
# input scripts which link to the input scripts and data files
|
||||
# you hopefully have created earlier with moltemplate.sh:
|
||||
# system.in.init, system.in.settings, system.data
|
||||
# If not, carry out the instructions in "README_setup.sh".
|
||||
#
|
||||
# -- Instructions: --
|
||||
# If "lmp_mpi" is the name of the command you use to invoke lammps,
|
||||
# then you would run lammps on these files this way:
|
||||
|
||||
|
||||
lmp_mpi -i run.in.min # minimization
|
||||
lmp_mpi -i run.in.nvt # minimization and simulation at constant volume
|
||||
|
||||
#(Note: The constant volume simulation lacks pressure equilibration. These are
|
||||
# completely separate simulations. The results of the constant pressure
|
||||
# simulation might be ignored when beginning the simulation at constant
|
||||
# volume. (This is because restart files in LAMMPS don't always work,
|
||||
# and I was spending a lot of time trying to convince people it was a
|
||||
# LAMMPS bug, instead of a moltemplate bug, so I disabled restart files.)
|
||||
# Read the "run.in.nvt" file to find out how to use the "read_restart"
|
||||
# command to load the results of the pressure-equilibration simulation,
|
||||
# before beginning a constant-volume run.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.npt
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# Optional:
|
||||
# To check for missing angle,dihedral params run moltemplate this way instead:
|
||||
# moltemplate.sh -checkff system.lt
|
||||
|
||||
|
||||
# Moltemplate generates various files with names ending in *.in* and *.data.
|
||||
# Move them to the directory where you plan to run LAMMPS (in this case "../")
|
||||
mv -f system.data system.in* ../
|
||||
|
||||
# Optional:
|
||||
# The "./output_ttree/" directory is full of temporary files generated by
|
||||
# moltemplate. They can be useful for debugging, but are usually thrown away.
|
||||
rm -rf output_ttree/
|
||||
|
||||
cd ../
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional:
|
||||
# Note: The system.data and system.in.settings files contain extra information
|
||||
# for atoms defined in OPLSAA which you are not using in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
#
|
||||
# cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
|
||||
------- To view a lammps trajectory in VMD --------
|
||||
|
||||
|
||||
1) Build a PSF file for use in viewing with VMD.
|
||||
|
||||
This step works with VMD 1.9 and topotools 1.2.
|
||||
(Older versions, like VMD 1.8.6, don't support this.)
|
||||
|
||||
|
||||
a) Start VMD
|
||||
b) Menu Extensions->Tk Console
|
||||
c) Enter:
|
||||
|
||||
(I assume that the the DATA file is called "system.data")
|
||||
|
||||
topo readlammpsdata system.data full
|
||||
animate write psf system.psf
|
||||
|
||||
2)
|
||||
|
||||
Later, to Load a trajectory in VMD:
|
||||
|
||||
Start VMD
|
||||
Select menu: File->New Molecule
|
||||
-Browse to select the PSF file you created above, and load it.
|
||||
(Don't close the window yet.)
|
||||
-Browse to select the trajectory file.
|
||||
If necessary, for "file type" select: "LAMMPS Trajectory"
|
||||
Load it.
|
||||
|
||||
---- A note on trajectory format: -----
|
||||
If the trajectory is a DUMP file, then make sure the it contains the
|
||||
information you need for pbctools (see below. I've been using this
|
||||
command in my LAMMPS scripts to create the trajectories:
|
||||
|
||||
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
It's a good idea to use an atom_style which supports molecule-ID numbers
|
||||
so that you can assign a molecule-ID number to each atom. (I think this
|
||||
is needed to wrap atom coordinates without breaking molecules in half.)
|
||||
|
||||
Of course, you don't have to save your trajectories in DUMP format,
|
||||
(other formats like DCD work fine) I just mention dump files
|
||||
because these are the files I'm familiar with.
|
||||
|
||||
3) ----- Wrap the coordinates to the unit cell
|
||||
(without cutting the molecules in half)
|
||||
|
||||
a) Start VMD
|
||||
b) Load the trajectory in VMD (see above)
|
||||
c) Menu Extensions->Tk Console
|
||||
d) Try entering these commands:
|
||||
|
||||
pbc wrap -compound res -all
|
||||
pbc box
|
||||
|
||||
----- Optional ----
|
||||
Sometimes the solvent or membrane obscures the view of the solute.
|
||||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
pbc box -shiftcenterrel {-0.05 -0.05 -0.05} -width 0.9 -style tubes
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
Alternately if you have a solute whose atoms are all of type 1,
|
||||
then you can also try this to center the box around it:
|
||||
|
||||
pbc wrap -sel type=1 -all -centersel type=2 -center com
|
||||
|
||||
4)
|
||||
You should check if your periodic boundary conditions are too small.
|
||||
To do that:
|
||||
select Graphics->Representations menu option
|
||||
click on the "Periodic" tab, and
|
||||
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
|
||||
|
||||
5) Optional: If you like, change the atom types in the PSF file so
|
||||
that VMD recognizes the atom types, use something like:
|
||||
|
||||
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
|
||||
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
|
||||
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
|
||||
|
||||
(If you do this, it might effect step 2 above.)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 28 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 28 KiB |
|
|
@ -0,0 +1,39 @@
|
|||
import "ch2group.lt" # load the definition of the "CH2" object
|
||||
import "ch3group.lt" # load the definition of the "CH3" object
|
||||
|
||||
|
||||
Butane inherits COMPASS {
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
|
||||
# Create an array of 4 objects distributed along the X axis
|
||||
|
||||
monomer1 = new CH3
|
||||
monomer2 = new CH2.rot(180,1,0,0).move(1.2533223,0,0)
|
||||
monomer3 = new CH2.move(2.5066446,0,0)
|
||||
monomer4 = new CH3.rot(180,0,0,1).move(3.7599669,0,0)
|
||||
|
||||
# Now add a list of bonds connecting the carbon atoms together:
|
||||
# (Angles, dihedrals, impropers will be automatically added later.)
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:b1 $atom:monomer1/C $atom:monomer2/C
|
||||
$bond:b2 $atom:monomer2/C $atom:monomer3/C
|
||||
$bond:b3 $atom:monomer3/C $atom:monomer4/C
|
||||
}
|
||||
|
||||
} # Butane
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163316030377
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.8924307629540046
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.6310438442242609
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
# This file contains a definition for the "CH2" molecular subunit.
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
# Then define "CH2":
|
||||
|
||||
|
||||
|
||||
|
||||
CH2 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
}
|
||||
|
||||
} # CH2
|
||||
|
||||
|
||||
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names (eg. "c4", "h1") are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH2.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
# This file contains a definition for the "CH3" molecular subunit.
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
# Then define "CH3":
|
||||
|
||||
|
||||
|
||||
|
||||
CH3 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
$atom:H3 $mol:... @atom:h1 0.0 -0.892431 -0.631044 0.000000
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
$bond:CH3 $atom:C $atom:H3
|
||||
}
|
||||
|
||||
} # CH3
|
||||
|
||||
|
||||
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names (eg. "c4", "h1") are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH3.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,24 @@
|
|||
import "butane.lt" # <- defines the "Butane" molecule type.
|
||||
|
||||
# Periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0.0 62.4 xlo xhi
|
||||
0.0 62.4 ylo yhi
|
||||
0.0 62.4 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Generate an array of 864 = 12 x 12 x 6 Butane molecules
|
||||
# which (more or less) uniformly fills the simulation box:
|
||||
|
||||
molecules = new Butane [12].move(0, 0, 5.2)
|
||||
[12].move(0, 5.2, 0)
|
||||
[6].move(10.4, 0, 0)
|
||||
|
||||
|
||||
|
||||
# NOTE: The spacing between molecules is large. There should be extra room to
|
||||
# move during the initial stages of equilibration. However, you will have to
|
||||
# run the simulation at NPT conditions later to compress the system to a
|
||||
# more realistic density.
|
||||
|
|
@ -0,0 +1,103 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data system.data
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# To avoid explosions, I have a 4-step equilibraion process (expand, minimize,
|
||||
# reorient, compress). The system (as defined in the "system.data" file)
|
||||
# is already expanded. That means there are 3 steps left:
|
||||
|
||||
dump dumpeq1 all custom 50 traj_eq1_min.lammpstrj id mol type x y z ix iy iz
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed press vol
|
||||
thermo 50
|
||||
thermo_modify norm yes
|
||||
|
||||
# -- Equilibration: part 1: initial minimization --
|
||||
|
||||
# Note: In general, it's always a good idea to minimize the system at first.
|
||||
|
||||
minimize 1.0e-5 1.0e-7 100000 400000
|
||||
undump dumpeq1
|
||||
|
||||
write_data system_after_eq1_min.data
|
||||
|
||||
# -- Equilibration part 2: reorienting the molecules (NVT) --
|
||||
|
||||
timestep 1.0
|
||||
dump dumpeq2 all custom 200 traj_eq2_reorient.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
# Run the system at high temperature (at constant volume) to reorient the
|
||||
# the molecules (which would otherwise be pointing in the same direction).
|
||||
|
||||
# To speed it up, I randomize the atomic positions for a few thousand steps
|
||||
# using fix langevin (and fix nve). Then I switch to fix nvt (Nose-Hoover).
|
||||
# (If I start with fix nvt (Nose-Hoover), it seems to get "stuck" for a while.)
|
||||
|
||||
|
||||
fix fxlan all langevin 900.0 900.0 120 48279
|
||||
fix fxnve all nve
|
||||
|
||||
run 4000
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnve
|
||||
# Now continue the simulation at high temperature using fix nvt (Nose-Hoover).
|
||||
fix fxnvt all nvt temp 900.0 900.0 100.0
|
||||
|
||||
run 10000
|
||||
undump dumpeq2
|
||||
|
||||
|
||||
write_data system_after_eq2_reorient.data
|
||||
|
||||
unfix fxnvt
|
||||
|
||||
|
||||
# -- equilibration part 3: Equilibrating the density (NPT) --
|
||||
|
||||
# Originally, the simulation box (in "system.data" and "system.lt") was
|
||||
# unrealistically large. The spacing between the molecules was large also.
|
||||
# I did this to enable the molecules to move freely and reorient themselves.
|
||||
# After doing that, we should run the simulation under NPT conditions to
|
||||
# allow the simulation box to contract to it's natural size.
|
||||
# To help it collapse, we begin the simulation at a relatively high pressure
|
||||
# Later on, we will slowly decrease it to 1 bar.
|
||||
|
||||
# First cool the system. (Do this at high pressure to avoid bubble formation.)
|
||||
|
||||
dump dumpeq3 all custom 200 traj_eq3_npt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnpt all npt temp 900.0 260.0 100.0 iso 500.0 500.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 20000
|
||||
|
||||
# At the very end of the previous simulation, the temperature dropped below
|
||||
# the boiling point. Run the simulation for longer at these conditions to
|
||||
# give it a chance for the vapor -> liquid transition to complete.
|
||||
# We will also slowly decrease the pressure to 1 bar.
|
||||
|
||||
unfix fxnpt
|
||||
fix fxnpt all npt temp 260.0 260.0 100.0 iso 500.0 1.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 100000
|
||||
|
||||
write_data system_after_eq3_npt.data
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
# 2) You must equilibrate the system beforehand using "run.in.npt".
|
||||
# This will create the file "system_after_npt.data" which this file reads.
|
||||
# (Note: I have not verified that this equilibration protocol works well.)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
|
||||
# Read the coordinates generated by an earlier NPT simulation
|
||||
|
||||
read_data system_after_eq3_npt.data
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also. I prefer "write_data" and "read_data".)
|
||||
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 500 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 260.0 260.0 500.0 tchain 1
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed
|
||||
thermo 100
|
||||
thermo_modify norm yes
|
||||
|
||||
run 50000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
This example is a simple simulation of many long alkane chains (hexadecane) in a
|
||||
box near the boiling point at atmospheric pressure. Please read "WARNING.TXT".
|
||||
|
||||
NOTE: This particular example uses the COMPASS force-field
|
||||
However, moltemplate is not limited to COMPASS.
|
||||
|
||||
1) Create the "system.data", "system.in.init", and "system.in.settings"
|
||||
files which LAMMPS will read by running:
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
|
||||
2) Run LAMMPS in this order:
|
||||
|
||||
lmp_mpi -i run.in.npt # running the simulation at constant pressure
|
||||
lmp_mpi -i run.in.nvt # running the simulation at constant temperature
|
||||
|
||||
(The name of the LAMMPS executable, eg "lmp_mpi", may vary.)
|
||||
|
||||
---- Details ----
|
||||
|
||||
The "Hexadecane" molecule, as well as the "CH2", and "CH3" monomers it contains
|
||||
use the COMPASS force-field. This means that when we define these molecules,
|
||||
we only specify the atom names, bond list, and coordinates.
|
||||
We do not have to list the atom charges, angles, dihedrals, or impropers.
|
||||
The rules for creating atomic charge and angle topology are contained in
|
||||
the "compass_published.lt" file created by step 3) above. The "ch2group.lt",
|
||||
"ch3group.lt", and "hexadecane.lt" files all refer to "compass_published.lt",
|
||||
(as well as the "COMPASS" force-field object which it defines). Excerpt:
|
||||
|
||||
import "compass_published.lt"
|
||||
CH2 inherits COMPASS { ...
|
||||
CH3 inherits COMPASS { ...
|
||||
Hexadecane inherits COMPASS { ...
|
||||
|
||||
Alternatively, you can manually define a list of angles, dihedrals, and
|
||||
improper interactions in these files, instead of asking the force-field
|
||||
to generate them for you. You can also specify some of the angles and
|
||||
dihedrals explicitly, and let the force-field handle the rest.
|
||||
(Many of the examples which come with moltemplate do this.)
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
|
||||
# Note: By default, the system.data and system.in.settings files contain
|
||||
# extra information for atoms defined in OPLSAA which you are not using
|
||||
# in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
|
||||
cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
# --- Running LAMMPS ---
|
||||
#
|
||||
# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
|
||||
# input scripts which link to the input scripts and data files
|
||||
# you hopefully have created earlier with moltemplate.sh:
|
||||
# system.in.init, system.in.settings, system.data
|
||||
# If not, carry out the instructions in "README_setup.sh".
|
||||
#
|
||||
# -- Instructions: --
|
||||
# If "lmp_mpi" is the name of the command you use to invoke lammps,
|
||||
# then you would run lammps on these files this way:
|
||||
|
||||
|
||||
lmp_mpi -i run.in.npt # minimization and simulation at constant pressure
|
||||
lmp_mpi -i run.in.nvt # simulation at constant volume
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.npt
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# Optional:
|
||||
# To check for missing angle,dihedral params run moltemplate this way instead:
|
||||
# moltemplate.sh -checkff system.lt
|
||||
|
||||
|
||||
# Moltemplate generates various files with names ending in *.in* and *.data.
|
||||
# Move them to the directory where you plan to run LAMMPS (in this case "../")
|
||||
mv -f system.data system.in* ../
|
||||
|
||||
# Optional:
|
||||
# The "./output_ttree/" directory is full of temporary files generated by
|
||||
# moltemplate. They can be useful for debugging, but are usually thrown away.
|
||||
rm -rf output_ttree/
|
||||
|
||||
cd ../
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional:
|
||||
# Note: The system.data and system.in.settings files contain extra information
|
||||
# for atoms defined in OPLSAA which you are not using in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
#
|
||||
# cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
|
||||
------- To view a lammps trajectory in VMD --------
|
||||
|
||||
|
||||
1) Build a PSF file for use in viewing with VMD.
|
||||
|
||||
This step works with VMD 1.9 and topotools 1.2.
|
||||
(Older versions, like VMD 1.8.6, don't support this.)
|
||||
|
||||
|
||||
a) Start VMD
|
||||
b) Menu Extensions->Tk Console
|
||||
c) Enter:
|
||||
|
||||
(I assume that the the DATA file is called "system.data")
|
||||
|
||||
topo readlammpsdata system.data full
|
||||
animate write psf system.psf
|
||||
|
||||
2)
|
||||
|
||||
Later, to Load a trajectory in VMD:
|
||||
|
||||
Start VMD
|
||||
Select menu: File->New Molecule
|
||||
-Browse to select the PSF file you created above, and load it.
|
||||
(Don't close the window yet.)
|
||||
-Browse to select the trajectory file.
|
||||
If necessary, for "file type" select: "LAMMPS Trajectory"
|
||||
Load it.
|
||||
|
||||
---- A note on trajectory format: -----
|
||||
If the trajectory is a DUMP file, then make sure the it contains the
|
||||
information you need for pbctools (see below. I've been using this
|
||||
command in my LAMMPS scripts to create the trajectories:
|
||||
|
||||
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
It's a good idea to use an atom_style which supports molecule-ID numbers
|
||||
so that you can assign a molecule-ID number to each atom. (I think this
|
||||
is needed to wrap atom coordinates without breaking molecules in half.)
|
||||
|
||||
Of course, you don't have to save your trajectories in DUMP format,
|
||||
(other formats like DCD work fine) I just mention dump files
|
||||
because these are the files I'm familiar with.
|
||||
|
||||
3) ----- Wrap the coordinates to the unit cell
|
||||
(without cutting the molecules in half)
|
||||
|
||||
a) Start VMD
|
||||
b) Load the trajectory in VMD (see above)
|
||||
c) Menu Extensions->Tk Console
|
||||
d) Try entering these commands:
|
||||
|
||||
pbc wrap -compound res -all
|
||||
pbc box
|
||||
|
||||
----- Optional ----
|
||||
Sometimes the solvent or membrane obscures the view of the solute.
|
||||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc box -shiftcenterrel {0.0 0.15 0.0}
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
Alternately if you have a solute whose atoms are all of type 1,
|
||||
then you can also try this to center the box around it:
|
||||
|
||||
pbc wrap -sel type=1 -all -centersel type=2 -center com
|
||||
|
||||
4)
|
||||
You should check if your periodic boundary conditions are too small.
|
||||
To do that:
|
||||
select Graphics->Representations menu option
|
||||
click on the "Periodic" tab, and
|
||||
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
|
||||
|
||||
5) Optional: If you like, change the atom types in the PSF file so
|
||||
that VMD recognizes the atom types, use something like:
|
||||
|
||||
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
|
||||
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
|
||||
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
|
||||
|
||||
(If you do this, it might effect step 2 above.)
|
||||
|
|
@ -0,0 +1,13 @@
|
|||
# -------- WARNING: --------
|
||||
|
||||
This software is experimental, and the force-fields and equilbration protocols
|
||||
have not been tested carefully by me. There is no gaurantee that the simulation
|
||||
will reproduce the behavior of real hexadecane molecules,
|
||||
|
||||
# -------- REQUEST FOR HELP: --------
|
||||
|
||||
However, if you notice a problem with this example, please report it.
|
||||
Peer-review is the only way to improve this software (or any software).
|
||||
Other suggestions are also welcome!
|
||||
|
||||
(Contact jewett.aij@gmail.com, 2017-10-03)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 26 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 35 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.9 KiB |
|
|
@ -0,0 +1,70 @@
|
|||
# This file contains a definition for the "CH2" molecular subunit.
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
# Then define "CH2":
|
||||
|
||||
|
||||
|
||||
|
||||
CH2 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
}
|
||||
|
||||
} # CH2
|
||||
|
||||
|
||||
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names (eg. "c4", "h1") are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH2.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
# This file contains a definition for the "CH3" molecular subunit.
|
||||
# First, load the COMPASS force field parameters we will need:
|
||||
|
||||
import "compass_published.lt" # <-- defines the "COMPASS" force field
|
||||
|
||||
# (The "compass_published.lt" file is located in the "force_fields"
|
||||
# subdirectory distributed with moltemplate.)
|
||||
# Then define "CH3":
|
||||
|
||||
|
||||
|
||||
|
||||
CH3 inherits COMPASS {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:c4 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:h1 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:h1 0.0 0.000000 0.631044 -0.892431
|
||||
$atom:H3 $mol:... @atom:h1 0.0 -0.892431 -0.631044 0.000000
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
$bond:CH3 $atom:C $atom:H3
|
||||
}
|
||||
|
||||
} # CH3
|
||||
|
||||
|
||||
|
||||
|
||||
# COMMENTS:
|
||||
# 1) Atom type names (eg. "c4", "h1") are defined in "compass_published.lt".
|
||||
# 2) In this example, the atomic charge of an atom is calculated by summing
|
||||
# partial charge contributions from neighboring atoms bonded to this atom.
|
||||
# (according to the rules in "compass_published.lt"). For this reason,
|
||||
# we can ignore the "charge" column in the "Data Atoms" section. Just
|
||||
# leave theses charges as "0.0" for now. Moltemplate will recalculate them.
|
||||
# 3) The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH3.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,79 @@
|
|||
# This example looks complicated because I split the
|
||||
# hexadecane molecule into individual CH2 and CH3 monomers.
|
||||
#
|
||||
# I defined it this way so that you can easily modify
|
||||
# it to change the length of the alkane chain.
|
||||
|
||||
|
||||
import "compass_published.lt" # load the "COMPASS" force-field information
|
||||
import "ch2group.lt" # load the definition of the "CH2" object
|
||||
import "ch3group.lt" # load the definition of the "CH3" object
|
||||
|
||||
|
||||
|
||||
Hexadecane inherits COMPASS {
|
||||
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
|
||||
|
||||
# Now create an array of 16 "CH2" objects distributed along the X axis
|
||||
|
||||
monomers = new CH2 [16].rot(180,1,0,0).move(1.2533223,0,0)
|
||||
|
||||
# Each monomer is rotated 180 degrees with respect to the previous
|
||||
# monomer, and then moved 1.2533223 Angstroms down the X axis.
|
||||
|
||||
# ---- Now, modify the ends: ---
|
||||
# Delete the CH2 groups at the beginning and end, and replace them with CH3.
|
||||
|
||||
delete monomers[0]
|
||||
delete monomers[15]
|
||||
|
||||
monomers[0] = new CH3
|
||||
monomers[15] = new CH3
|
||||
|
||||
# Move the CH3 groups to the correct location at either end of the chain:
|
||||
|
||||
monomers[15].rot(180.0,0,0,1).move(18.7998345,0,0)
|
||||
# Note: 18.7998345 = (16-1) * 1.2533223
|
||||
|
||||
|
||||
# Now add a list of bonds connecting the carbon atoms together:
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:b1 $atom:monomers[0]/C $atom:monomers[1]/C
|
||||
$bond:b2 $atom:monomers[1]/C $atom:monomers[2]/C
|
||||
$bond:b3 $atom:monomers[2]/C $atom:monomers[3]/C
|
||||
$bond:b4 $atom:monomers[3]/C $atom:monomers[4]/C
|
||||
$bond:b5 $atom:monomers[4]/C $atom:monomers[5]/C
|
||||
$bond:b6 $atom:monomers[5]/C $atom:monomers[6]/C
|
||||
$bond:b7 $atom:monomers[6]/C $atom:monomers[7]/C
|
||||
$bond:b8 $atom:monomers[7]/C $atom:monomers[8]/C
|
||||
$bond:b9 $atom:monomers[8]/C $atom:monomers[9]/C
|
||||
$bond:b10 $atom:monomers[9]/C $atom:monomers[10]/C
|
||||
$bond:b11 $atom:monomers[10]/C $atom:monomers[11]/C
|
||||
$bond:b12 $atom:monomers[11]/C $atom:monomers[12]/C
|
||||
$bond:b13 $atom:monomers[12]/C $atom:monomers[13]/C
|
||||
$bond:b14 $atom:monomers[13]/C $atom:monomers[14]/C
|
||||
$bond:b15 $atom:monomers[14]/C $atom:monomers[15]/C
|
||||
}
|
||||
|
||||
} # Hexadecane
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163316030377
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.8924307629540046
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.6310438442242609
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
import "hexadecane.lt" # <- defines the "Hexadecane" molecule type.
|
||||
|
||||
# Periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0.0 62.4 xlo xhi
|
||||
0.0 62.4 ylo yhi
|
||||
0.0 62.4 zlo zhi
|
||||
}
|
||||
|
||||
molecules = new Hexadecane [12].move(0, 0, 5.2)
|
||||
[12].move(0, 5.2, 0)
|
||||
[2].move(31.2, 0, 0)
|
||||
|
||||
|
||||
|
||||
|
||||
# NOTE: The spacing between molecules is large. There should be extra room to
|
||||
# move during the initial stages of equilibration. However, you will have to
|
||||
# run the simulation at NPT conditions later to compress the system to a
|
||||
# more realistic density.
|
||||
|
|
@ -0,0 +1,102 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data system.data
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# To avoid explosions, I have a 4-step equilibraion process (expand, minimize,
|
||||
# reorient, compress). The system (as defined in the "system.data" file)
|
||||
# is already expanded. That means there are 3 steps left:
|
||||
|
||||
dump dumpeq1 all custom 50 traj_eq1_min.lammpstrj id mol type x y z ix iy iz
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed press vol
|
||||
thermo 50
|
||||
thermo_modify norm yes
|
||||
|
||||
# -- Equilibration: part 1: initial minimization --
|
||||
|
||||
# Note: In general, it's always a good idea to minimize the system at first.
|
||||
|
||||
minimize 1.0e-5 1.0e-7 100000 400000
|
||||
undump dumpeq1
|
||||
|
||||
write_data system_after_eq1_min.data
|
||||
|
||||
# -- Equilibration part 2: reorienting the molecules (NVT) --
|
||||
|
||||
timestep 1.0
|
||||
dump dumpeq2 all custom 200 traj_eq2_reorient.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
# Run the system at high temperature (at constant volume) to reorient the
|
||||
# the molecules (which would otherwise be pointing in the same direction).
|
||||
|
||||
# To speed it up, I randomize the atomic positions for a few thousand steps
|
||||
# using fix langevin (and fix nve). Then I switch to fix nvt (Nose-Hoover).
|
||||
# (If I start with fix nvt (Nose-Hoover), it seems to get "stuck" for a while.)
|
||||
|
||||
|
||||
fix fxlan all langevin 900.0 900.0 120 48279
|
||||
fix fxnve all nve
|
||||
|
||||
run 4000
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnve
|
||||
# Now continue the simulation at high temperature using fix nvt (Nose-Hoover).
|
||||
fix fxnvt all nvt temp 900.0 900.0 100.0
|
||||
|
||||
run 50000
|
||||
undump dumpeq2
|
||||
|
||||
|
||||
write_data system_after_eq2_reorient.data
|
||||
|
||||
unfix fxnvt
|
||||
|
||||
# -- equilibration part 3: Equilibrating the density (NPT) --
|
||||
|
||||
# Originally, the simulation box (in "system.data" and "system.lt") was
|
||||
# unrealistically large. The spacing between the molecules was large also.
|
||||
# I did this to enable the molecules to move freely and reorient themselves.
|
||||
# After doing that, we should run the simulation under NPT conditions to
|
||||
# allow the simulation box to contract to it's natural size. We do that here:
|
||||
# To help it collapse, we begin the simulation at a relatively high pressure
|
||||
# Later on, we will slowly decrease it to 1 bar.
|
||||
|
||||
# First cool the system. (Do this at high pressure to avoid bubble formation.)
|
||||
|
||||
dump dumpeq3 all custom 200 traj_eq3_npt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnpt all npt temp 900.0 260.0 100.0 iso 500.0 500.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 20000
|
||||
|
||||
# At the very end of the previous simulation, the temperature dropped below
|
||||
# the boiling point. Run the simulation for longer at these conditions to
|
||||
# give it a chance for the vapor -> liquid transition to complete.
|
||||
# We will also slowly decrease the pressure to 1 bar.
|
||||
|
||||
unfix fxnpt
|
||||
fix fxnpt all npt temp 260.0 260.0 100.0 iso 500.0 1.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 100000
|
||||
|
||||
write_data system_after_eq3_npt.data
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
# 2) You must equilibrate the system beforehand using "run.in.npt".
|
||||
# This will create the file "system_after_npt.data" which this file reads.
|
||||
# (Note: I have not verified that this equilibration protocol works well.)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
|
||||
# Read the coordinates generated by an earlier NPT simulation
|
||||
|
||||
read_data system_after_eq3_npt.data
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also. I prefer "write_data" and "read_data".)
|
||||
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 500 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed
|
||||
thermo 100
|
||||
thermo_modify norm yes
|
||||
|
||||
run 50000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -0,0 +1,37 @@
|
|||
This example is a simple simulation of many short alkane chains (butane) in a
|
||||
box near the boiling point at atmospheric pressure. Please read "WARNING.TXT".
|
||||
|
||||
1) Create the "system.data", "system.in.init", and "system.in.settings"
|
||||
files which LAMMPS will read by running:
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
|
||||
2) Run LAMMPS in this order:
|
||||
|
||||
lmp_mpi -i run.in.npt # running the simulation at constant pressure
|
||||
lmp_mpi -i run.in.nvt # running the simulation at constant temperature
|
||||
|
||||
(The name of the LAMMPS executable, eg "lmp_mpi", may vary.)
|
||||
|
||||
---- Details ----
|
||||
|
||||
The "Butane50" molecule, as well as the "CH2", and "CH3" monomers it contains
|
||||
use the OPLSAA force-field. This means that when we define these molecules,
|
||||
we only specify the atom names, bond list, and coordinates.
|
||||
We do not have to list the atom charges, angles, dihedrals, or impropers.
|
||||
The rules for creating atomic charge and angle topology are contained in
|
||||
the "oplsaa.lt" file created by step 3) above. The "ch2group.lt",
|
||||
"ch3group.lt", and "butane.lt" files all refer to "oplsaa.lt",
|
||||
(as well as the "OPLSAA" force-field object which it defines). Excerpt:
|
||||
|
||||
import "oplsaa.lt"
|
||||
CH2 inherits OPLSAA { ...
|
||||
CH3 inherits OPLSAA { ...
|
||||
Butane inherits OPLSAA { ...
|
||||
|
||||
Alternatively, you can manually define a list of angles, dihedrals, and
|
||||
improper interactions in these files, instead of asking the force-field
|
||||
to generate them for you. You can also specify some of the angles and
|
||||
dihedrals explicitly, and let the force-field handle the rest.
|
||||
(Many of the examples which come with moltemplate do this.)
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
|
||||
# Note: By default, the system.data and system.in.settings files contain
|
||||
# extra information for atoms defined in OPLSAA which you are not using
|
||||
# in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
|
||||
cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
# --- Running LAMMPS ---
|
||||
#
|
||||
# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
|
||||
# input scripts which link to the input scripts and data files
|
||||
# you hopefully have created earlier with moltemplate.sh:
|
||||
# system.in.init, system.in.settings, system.data
|
||||
# If not, carry out the instructions in "README_setup.sh".
|
||||
#
|
||||
# -- Instructions: --
|
||||
# If "lmp_mpi" is the name of the command you use to invoke lammps,
|
||||
# then you would run lammps on these files this way:
|
||||
|
||||
|
||||
lmp_mpi -i run.in.min # minimization
|
||||
lmp_mpi -i run.in.nvt # minimization and simulation at constant volume
|
||||
|
||||
#(Note: The constant volume simulation lacks pressure equilibration. These are
|
||||
# completely separate simulations. The results of the constant pressure
|
||||
# simulation might be ignored when beginning the simulation at constant
|
||||
# volume. (This is because restart files in LAMMPS don't always work,
|
||||
# and I was spending a lot of time trying to convince people it was a
|
||||
# LAMMPS bug, instead of a moltemplate bug, so I disabled restart files.)
|
||||
# Read the "run.in.nvt" file to find out how to use the "read_restart"
|
||||
# command to load the results of the pressure-equilibration simulation,
|
||||
# before beginning a constant-volume run.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.npt
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# Optional:
|
||||
# To check for missing angle,dihedral params run moltemplate this way instead:
|
||||
# moltemplate.sh -checkff system.lt
|
||||
|
||||
|
||||
# Moltemplate generates various files with names ending in *.in* and *.data.
|
||||
# Move them to the directory where you plan to run LAMMPS (in this case "../")
|
||||
mv -f system.data system.in* ../
|
||||
|
||||
# Optional:
|
||||
# The "./output_ttree/" directory is full of temporary files generated by
|
||||
# moltemplate. They can be useful for debugging, but are usually thrown away.
|
||||
rm -rf output_ttree/
|
||||
|
||||
cd ../
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional:
|
||||
# Note: The system.data and system.in.settings files contain extra information
|
||||
# for atoms defined in OPLSAA which you are not using in this simulation.
|
||||
# This is harmless, but if you to delete this information from your
|
||||
# system.in.settings and system.in.data files, run this script:
|
||||
#
|
||||
# cleanup_moltemplate.sh
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
|
||||
------- To view a lammps trajectory in VMD --------
|
||||
|
||||
|
||||
1) Build a PSF file for use in viewing with VMD.
|
||||
|
||||
This step works with VMD 1.9 and topotools 1.2.
|
||||
(Older versions, like VMD 1.8.6, don't support this.)
|
||||
|
||||
|
||||
a) Start VMD
|
||||
b) Menu Extensions->Tk Console
|
||||
c) Enter:
|
||||
|
||||
(I assume that the the DATA file is called "system.data")
|
||||
|
||||
topo readlammpsdata system.data full
|
||||
animate write psf system.psf
|
||||
|
||||
2)
|
||||
|
||||
Later, to Load a trajectory in VMD:
|
||||
|
||||
Start VMD
|
||||
Select menu: File->New Molecule
|
||||
-Browse to select the PSF file you created above, and load it.
|
||||
(Don't close the window yet.)
|
||||
-Browse to select the trajectory file.
|
||||
If necessary, for "file type" select: "LAMMPS Trajectory"
|
||||
Load it.
|
||||
|
||||
---- A note on trajectory format: -----
|
||||
If the trajectory is a DUMP file, then make sure the it contains the
|
||||
information you need for pbctools (see below. I've been using this
|
||||
command in my LAMMPS scripts to create the trajectories:
|
||||
|
||||
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
It's a good idea to use an atom_style which supports molecule-ID numbers
|
||||
so that you can assign a molecule-ID number to each atom. (I think this
|
||||
is needed to wrap atom coordinates without breaking molecules in half.)
|
||||
|
||||
Of course, you don't have to save your trajectories in DUMP format,
|
||||
(other formats like DCD work fine) I just mention dump files
|
||||
because these are the files I'm familiar with.
|
||||
|
||||
3) ----- Wrap the coordinates to the unit cell
|
||||
(without cutting the molecules in half)
|
||||
|
||||
a) Start VMD
|
||||
b) Load the trajectory in VMD (see above)
|
||||
c) Menu Extensions->Tk Console
|
||||
d) Try entering these commands:
|
||||
|
||||
pbc wrap -compound res -all
|
||||
pbc box
|
||||
|
||||
----- Optional ----
|
||||
Sometimes the solvent or membrane obscures the view of the solute.
|
||||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
pbc box -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
Alternately if you have a solute whose atoms are all of type 1,
|
||||
then you can also try this to center the box around it:
|
||||
|
||||
pbc wrap -sel type=1 -all -centersel type=2 -center com
|
||||
|
||||
4)
|
||||
You should check if your periodic boundary conditions are too small.
|
||||
To do that:
|
||||
select Graphics->Representations menu option
|
||||
click on the "Periodic" tab, and
|
||||
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
|
||||
|
||||
5) Optional: If you like, change the atom types in the PSF file so
|
||||
that VMD recognizes the atom types, use something like:
|
||||
|
||||
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
|
||||
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
|
||||
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
|
||||
|
||||
(If you do this, it might effect step 2 above.)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 28 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 5.0 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.8 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 28 KiB |
|
|
@ -0,0 +1,42 @@
|
|||
import "ch2group.lt" # load the definition of the "CH2" object
|
||||
import "ch3group.lt" # load the definition of the "CH3" object
|
||||
|
||||
|
||||
Butane inherits OPLSAA {
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
|
||||
# Create an array of 4 objects distributed along the X axis
|
||||
|
||||
monomer1 = new CH3
|
||||
monomer2 = new CH2.rot(180,1,0,0).move(1.2533223,0,0)
|
||||
monomer3 = new CH2.move(2.5066446,0,0)
|
||||
monomer4 = new CH3.rot(180,0,0,1).move(3.7599669,0,0)
|
||||
|
||||
# Now add a list of bonds connecting the carbon atoms together:
|
||||
# (Angles, dihedrals, impropers will be automatically added later.)
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:b1 $atom:monomer1/C $atom:monomer2/C
|
||||
$bond:b2 $atom:monomer2/C $atom:monomer3/C
|
||||
$bond:b3 $atom:monomer3/C $atom:monomer4/C
|
||||
}
|
||||
|
||||
} # Butane
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163316030377
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.8924307629540046
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.6310438442242609
|
||||
|
|
@ -0,0 +1,67 @@
|
|||
# This file contains a definition for the "CH2" molecular subunit.
|
||||
|
||||
# First, load the OPLS force field parameters we will need.
|
||||
# This file is located in the "force_fields" subdirectory
|
||||
# of the moltemplate distribution.
|
||||
|
||||
import "oplsaa.lt" # <-- defines the standard "OPLSAA" force field
|
||||
|
||||
|
||||
# Then define "CH2":
|
||||
|
||||
|
||||
CH2 inherits OPLSAA {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:81 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:85 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:85 0.0 0.000000 0.631044 -0.892431
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
}
|
||||
|
||||
} # CH2
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH2.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
# ---- NOTES: ----
|
||||
#
|
||||
# Atom type numbers (@atom:80,81,85) are defined in "oplsaa.lt",
|
||||
# @atom:80 "Alkane CH3-"
|
||||
# @atom:81 "Alkane -CH2-"
|
||||
# @atom:85 "Alkane H-C CH3"
|
||||
# @atom:85 "Alkane H-C CH2"
|
||||
# In this example, atomic charges are generated by atom type (according to the
|
||||
# rules in oplsaa.lt), and can be omitted. Just leave them as "0.00" for now.
|
||||
# The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
#
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,69 @@
|
|||
# This file contains a definition for the "CH3" molecular subunit.
|
||||
|
||||
# First, load the OPLS force field parameters we will need.
|
||||
# This file is located in the "force_fields" subdirectory
|
||||
# of the moltemplate distribution.
|
||||
|
||||
import "oplsaa.lt" # <-- defines the standard "OPLSAA" force field
|
||||
|
||||
|
||||
# Then define "CH3":
|
||||
|
||||
|
||||
CH3 inherits OPLSAA {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C $mol:... @atom:80 0.0 0.000000 0.000000 0.000000
|
||||
$atom:H1 $mol:... @atom:85 0.0 0.000000 0.631044 0.892431
|
||||
$atom:H2 $mol:... @atom:85 0.0 0.000000 0.631044 -0.892431
|
||||
$atom:H3 $mol:... @atom:85 0.0 -0.892431 -0.631044 0.000000
|
||||
}
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:CH1 $atom:C $atom:H1
|
||||
$bond:CH2 $atom:C $atom:H2
|
||||
$bond:CH3 $atom:C $atom:H3
|
||||
}
|
||||
|
||||
} # CH3
|
||||
|
||||
|
||||
|
||||
|
||||
# Optional: Shift all the coordinates in the +Y direction by 0.4431163.
|
||||
# This way, the carbon atom is no longer located at 0,0,0, but the
|
||||
# axis of an alkane chain containing this monomer is at 0,0,0.
|
||||
# (This makes it more convenient to construct a polymer later.
|
||||
# If this is confusing, then simply add 0.4431163 to the Y
|
||||
# coordinates in the "Data Atoms" section above.)
|
||||
|
||||
CH3.move(0,0.4431163,0)
|
||||
|
||||
|
||||
|
||||
|
||||
# ---- NOTES: ----
|
||||
#
|
||||
# Atom type numbers (@atom:80,81,85) are defined in "oplsaa.lt",
|
||||
# @atom:80 "Alkane CH3-"
|
||||
# @atom:81 "Alkane -CH2-"
|
||||
# @atom:85 "Alkane H-C CH3"
|
||||
# @atom:85 "Alkane H-C CH2"
|
||||
# In this example, atomic charges are generated by atom type (according to the
|
||||
# rules in oplsaa.lt), and can be omitted. Just leave them as "0.00" for now.
|
||||
# The "..." in "$mol:..." tells moltemplate that this molecule may be part
|
||||
# of a larger molecule, and (if so) to use the larger parent object's
|
||||
# molecule id number as it's own.
|
||||
#
|
||||
######### (scratchwork calculations for the atomic coordinates) #########
|
||||
#
|
||||
# Lcc = 1.5350 # length of the C-C bond (Sp3)
|
||||
# Lch = 1.0930 # length of the C-H bond
|
||||
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
|
||||
# DeltaXc = Lcc*sin(theta/2) # = 1.2533222517240594
|
||||
# DeltaYc = Lcc*cos(theta/2) # = 0.8862326632060754
|
||||
# # 0.5*DeltaYc = 0.4431163
|
||||
# DeltaZh = Lch*sin(theta/2) # = 0.892431
|
||||
# DeltaYh = Lch*cos(theta/2) # = 0.631044
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
import "butane.lt" # <- defines the "Butane" molecule type.
|
||||
|
||||
# Periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0.0 62.4 xlo xhi
|
||||
0.0 62.4 ylo yhi
|
||||
0.0 62.4 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Generate an array of 864 = 12 x 12 x 6 Butane molecules
|
||||
# which (more or less) uniformly fills the simulation box:
|
||||
|
||||
molecules = new Butane [12].move(0, 0, 5.2)
|
||||
[12].move(0, 5.2, 0)
|
||||
[6].move(10.4, 0, 0)
|
||||
|
||||
|
||||
|
||||
|
||||
# NOTE: The spacing between molecules is large. There should be extra room to
|
||||
# move during the initial stages of equilibration. However, you will have to
|
||||
# run the simulation at NPT conditions later to compress the system to a
|
||||
# more realistic density.
|
||||
|
|
@ -0,0 +1,103 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data system.data
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# To avoid explosions, I have a 4-step equilibraion process (expand, minimize,
|
||||
# reorient, compress). The system (as defined in the "system.data" file)
|
||||
# is already expanded. That means there are 3 steps left:
|
||||
|
||||
dump dumpeq1 all custom 50 traj_eq1_min.lammpstrj id mol type x y z ix iy iz
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed press vol
|
||||
thermo 50
|
||||
thermo_modify norm yes
|
||||
|
||||
# -- Equilibration: part 1: initial minimization --
|
||||
|
||||
# Note: In general, it's always a good idea to minimize the system at first.
|
||||
|
||||
minimize 1.0e-5 1.0e-7 100000 400000
|
||||
undump dumpeq1
|
||||
|
||||
write_data system_after_eq1_min.data
|
||||
|
||||
# -- Equilibration part 2: reorienting the molecules (NVT) --
|
||||
|
||||
timestep 1.0
|
||||
dump dumpeq2 all custom 200 traj_eq2_reorient.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
# Run the system at high temperature (at constant volume) to reorient the
|
||||
# the molecules (which would otherwise be pointing in the same direction).
|
||||
|
||||
# To speed it up, I randomize the atomic positions for a few thousand steps
|
||||
# using fix langevin (and fix nve). Then I switch to fix nvt (Nose-Hoover).
|
||||
# (If I start with fix nvt (Nose-Hoover), it seems to get "stuck" for a while.)
|
||||
|
||||
|
||||
fix fxlan all langevin 900.0 900.0 120 48279
|
||||
fix fxnve all nve
|
||||
|
||||
run 4000
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnve
|
||||
# Now continue the simulation at high temperature using fix nvt (Nose-Hoover).
|
||||
fix fxnvt all nvt temp 900.0 900.0 100.0
|
||||
|
||||
run 10000
|
||||
undump dumpeq2
|
||||
|
||||
|
||||
write_data system_after_eq2_reorient.data
|
||||
|
||||
unfix fxnvt
|
||||
|
||||
|
||||
# -- equilibration part 3: Equilibrating the density (NPT) --
|
||||
|
||||
# Originally, the simulation box (in "system.data" and "system.lt") was
|
||||
# unrealistically large. The spacing between the molecules was large also.
|
||||
# I did this to enable the molecules to move freely and reorient themselves.
|
||||
# After doing that, we should run the simulation under NPT conditions to
|
||||
# allow the simulation box to contract to it's natural size.
|
||||
# To help it collapse, we begin the simulation at a relatively high pressure
|
||||
# Later on, we will slowly decrease it to 1 bar.
|
||||
|
||||
# First cool the system. (Do this at high pressure to avoid bubble formation.)
|
||||
|
||||
dump dumpeq3 all custom 200 traj_eq3_npt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnpt all npt temp 900.0 260.0 100.0 iso 500.0 500.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 20000
|
||||
|
||||
# At the very end of the previous simulation, the temperature dropped below
|
||||
# the boiling point. Run the simulation for longer at these conditions to
|
||||
# give it a chance for the vapor -> liquid transition to complete.
|
||||
# We will also slowly decrease the pressure to 1 bar.
|
||||
|
||||
unfix fxnpt
|
||||
fix fxnpt all npt temp 260.0 260.0 100.0 iso 500.0 1.0 1000.0 drag 2.0
|
||||
|
||||
timestep 1.0
|
||||
run 100000
|
||||
|
||||
write_data system_after_eq3_npt.data
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (Follow the instructions in README_setup.sh,
|
||||
# or run the file as a script using ./README_setup.sh)
|
||||
# 2) You must equilibrate the system beforehand using "run.in.npt".
|
||||
# This will create the file "system_after_npt.data" which this file reads.
|
||||
# (Note: I have not verified that this equilibration protocol works well.)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
|
||||
# Read the coordinates generated by an earlier NPT simulation
|
||||
|
||||
read_data system_after_eq3_npt.data
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also. I prefer "write_data" and "read_data".)
|
||||
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
include system.in.charges
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 500 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 260.0 260.0 500.0 tchain 1
|
||||
thermo_style custom step temp pe etotal epair ebond eangle edihed
|
||||
thermo 100
|
||||
thermo_modify norm yes
|
||||
|
||||
run 50000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -1,8 +1,5 @@
|
|||
This is an example of how to use the OPLSAA force-field in LAMMPS
|
||||
|
||||
As of 2016-11-21, this code has not been tested for accuracy.
|
||||
(See the WARNING.TXT file.)
|
||||
|
||||
step 1)
|
||||
To build the files which LAMMPS needs, follow the instructions in:
|
||||
README_setup.sh
|
||||
|
|
|
|||
|
|
@ -60,8 +60,8 @@ d) Try entering these commands:
|
|||
It can help to shift the location of the periodic boundary box
|
||||
To shift the box in the y direction (for example) do this:
|
||||
|
||||
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc box -shiftcenterrel {0.0 0.15 0.0}
|
||||
pbc wrap -compound res -all -shiftcenterrel {-0.05 -0.05 -0.05}
|
||||
pbc box -shiftcenterrel {-0.05 -0.05 -0.05} -width 0.5
|
||||
|
||||
Distances are measured in units of box-length fractions, not Angstroms.
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 36 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 127 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 51 KiB |
|
|
@ -1,7 +1,7 @@
|
|||
import "oplsaa.lt"
|
||||
|
||||
# The "oplsaa.lt" file contains force-field definitions and masses for the
|
||||
# atoms in your system. See oplsaa_lt_generator/README.TXT for details.
|
||||
# The "oplsaa.lt" file contains force-field definitions,
|
||||
# partial charges and masses for the atoms in your system.
|
||||
|
||||
# Note:
|
||||
# Atom type @atom:90 corresponds to "Aromatic C"
|
||||
|
|
@ -9,25 +9,32 @@ import "oplsaa.lt"
|
|||
|
||||
Benzene inherits OPLSAA {
|
||||
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:90 0.00 5.274 1.999 -8.568 # "Aromatic C"
|
||||
$atom:C2 $mol @atom:90 0.00 6.627 2.018 -8.209 # "Aromatic C"
|
||||
$atom:C3 $mol @atom:90 0.00 7.366 0.829 -8.202 # "Aromatic C"
|
||||
$atom:C4 $mol @atom:90 0.00 6.752 -0.379 -8.554 # "Aromatic C"
|
||||
$atom:C5 $mol @atom:90 0.00 5.399 -0.398 -8.912 # "Aromatic C"
|
||||
$atom:C6 $mol @atom:90 0.00 4.660 0.791 -8.919 # "Aromatic C"
|
||||
$atom:H11 $mol @atom:91 0.00 4.704 2.916 -8.573 # "Aromatic H-C"
|
||||
$atom:H21 $mol @atom:91 0.00 7.101 2.950 -7.938 # "Aromatic H-C"
|
||||
$atom:H31 $mol @atom:91 0.00 8.410 0.844 -7.926 # "Aromatic H-C"
|
||||
$atom:H41 $mol @atom:91 0.00 7.322 -1.296 -8.548 # "Aromatic H-C"
|
||||
$atom:H51 $mol @atom:91 0.00 4.925 -1.330 -9.183 # "Aromatic H-C"
|
||||
$atom:H61 $mol @atom:91 0.00 3.616 0.776 -9.196 # "Aromatic H-C"
|
||||
# atom-id mol-id atom-type charge X Y Z # comment
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:C1 $mol @atom:90 0.00 -0.739 1.189 -0.00733 # 90 = "Aromatic C"
|
||||
$atom:C2 $mol @atom:90 0.00 0.614 1.208 0.35167 # 90 = "Aromatic C"
|
||||
$atom:C3 $mol @atom:90 0.00 1.353 0.019 0.35867 # 90 = "Aromatic C"
|
||||
$atom:C4 $mol @atom:90 0.00 0.739 -1.189 0.00667 # 90 = "Aromatic C"
|
||||
$atom:C5 $mol @atom:90 0.00 -0.614 -1.208 -0.35133 # 90 = "Aromatic C"
|
||||
$atom:C6 $mol @atom:90 0.00 -1.353 -0.019 -0.35833 # 90 = "Aromatic C"
|
||||
$atom:H11 $mol @atom:91 0.00 -1.309 2.106 -0.01233 # 91 = "Aromatic H-C"
|
||||
$atom:H21 $mol @atom:91 0.00 1.088 2.14 0.62267 # 91 = "Aromatic H-C"
|
||||
$atom:H31 $mol @atom:91 0.00 2.397 0.034 0.63467 # 91 = "Aromatic H-C"
|
||||
$atom:H41 $mol @atom:91 0.00 1.309 -2.106 0.01267 # 91 = "Aromatic H-C"
|
||||
$atom:H51 $mol @atom:91 0.00 -1.088 -2.14 -0.62233 # 91 = "Aromatic H-C"
|
||||
$atom:H61 $mol @atom:91 0.00 -2.397 -0.034 -0.63533 # 91 = "Aromatic H-C"
|
||||
}
|
||||
# Note: You don't have to specify the charge in this example because
|
||||
# we are using the OPLSAA force-field assigns this by atom-type.
|
||||
# Just leave these numbers as 0.00 for now.
|
||||
# Note: LAMMPS expects an integer in the 2nd column (the Molecule-ID number).
|
||||
# If we put "$mol" there, moltemplate will generate this integer for you
|
||||
|
||||
write('Data Bond List') {
|
||||
# A list of the bonds in the molecule
|
||||
# BondID AtomID1 AtomID2
|
||||
|
||||
write("Data Bond List") {
|
||||
$bond:C12 $atom:C1 $atom:C2
|
||||
$bond:C23 $atom:C2 $atom:C3
|
||||
$bond:C34 $atom:C3 $atom:C4
|
||||
|
|
@ -42,4 +49,7 @@ Benzene inherits OPLSAA {
|
|||
$bond:C6H6 $atom:C6 $atom:H61
|
||||
}
|
||||
|
||||
# In the "Data Bond List" section we don't have to specify the bond type.
|
||||
# The bond-type will be determined by the atom type (according to "oplsaa.lt")
|
||||
|
||||
} # Benzene
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
import "oplsaa.lt"
|
||||
|
||||
# The "oplsaa.lt" file contains force-field definitions and masses for the
|
||||
# atoms in your system. See oplsaa_lt_generator/README.TXT for details.
|
||||
# The "oplsaa.lt" file contains force-field definitions,
|
||||
# partial charges and masses for the atoms in your system.
|
||||
|
||||
# Note:
|
||||
# Atom type 88 corresponds to "Alkene H2-C="
|
||||
|
|
@ -11,19 +11,24 @@ import "oplsaa.lt"
|
|||
|
||||
Ethylene inherits OPLSAA {
|
||||
|
||||
# atom-id mol-id atom-type charge X Y Z
|
||||
# atom-id mol-id atom-type charge X Y Z # comment
|
||||
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:88 0.00 -0.6695 0.000000 0.000000
|
||||
$atom:C2 $mol @atom:88 0.00 0.6695 0.000000 0.000000
|
||||
$atom:H11 $mol @atom:89 0.00 -1.234217 -0.854458 0.000000
|
||||
$atom:H12 $mol @atom:89 0.00 -1.234217 0.854458 0.000000
|
||||
$atom:H21 $mol @atom:89 0.00 1.234217 -0.854458 0.000000
|
||||
$atom:H22 $mol @atom:89 0.00 1.234217 0.854458 0.000000
|
||||
$atom:C1 $mol @atom:88 0.00 -0.6695 0.00000 0.000 #88 = "Alkene H2-C="
|
||||
$atom:C2 $mol @atom:88 0.00 0.6695 0.00000 0.000 #88 = "Alkene H2-C="
|
||||
$atom:H11 $mol @atom:89 0.00 -1.23422 -0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H12 $mol @atom:89 0.00 -1.23422 0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H21 $mol @atom:89 0.00 1.23422 -0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H22 $mol @atom:89 0.00 1.23422 0.85446 0.000 #89 = "Alkene H-C="
|
||||
}
|
||||
# Note: You don't have to specify the charge in this example because
|
||||
# we are using the OPLSAA force-field assigns this by atom-type.
|
||||
# Just leave these numbers as 0.00 for now.
|
||||
# Note: LAMMPS expects an integer in the 2nd column (the Molecule-ID number).
|
||||
# If we put "$mol" there, moltemplate will generate this integer for you
|
||||
|
||||
# A list of the bonds in the molecule
|
||||
# BondID AtomID1 AtomID2
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:C12 $atom:C1 $atom:C2
|
||||
|
|
@ -33,11 +38,8 @@ Ethylene inherits OPLSAA {
|
|||
$bond:C2H2 $atom:C2 $atom:H22
|
||||
}
|
||||
|
||||
# In the "Data Bond List" section we don't have to specify the bond type.
|
||||
# The bond-type will be determined by the atom type (according to "oplsaa.lt")
|
||||
|
||||
} # Ethylene
|
||||
|
||||
|
||||
|
||||
# Note: You don't need to supply the partial partial charges of the atoms.
|
||||
# If you like, just fill the fourth column with zeros ("0.000").
|
||||
# Moltemplate and LAMMPS will automatically assign the charge later
|
||||
|
||||
|
|
|
|||
|
|
@ -4,20 +4,20 @@ import "benzene.lt" # <- defines the "Benzene" molecule type.
|
|||
|
||||
# Periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0.0 80.00 xlo xhi
|
||||
0.0 80.00 ylo yhi
|
||||
0.0 80.00 zlo zhi
|
||||
0.0 48.00 xlo xhi
|
||||
0.0 48.00 ylo yhi
|
||||
0.0 48.00 zlo zhi
|
||||
}
|
||||
|
||||
# Create 1000 ethylenes and 500 benzenes
|
||||
# Create 216 ethylenes and 108 benzenes
|
||||
|
||||
ethylenes = new Ethylene[10].move(8.0, 0, 0)
|
||||
[10].move(0, 8.0, 0)
|
||||
[10].move(0, 0, 8.0)
|
||||
ethylenes = new Ethylene[6].move(8.0, 0, 0)
|
||||
[6].move(0, 8.0, 0)
|
||||
[6].move(0, 0, 8.0)
|
||||
|
||||
benzenes = new Benzene[10].move(8.0, 0, 0)
|
||||
[10].move(0, 8.0, 0)
|
||||
[5].move(0, 0, 16.0)
|
||||
benzenes = new Benzene[6].move(8.0, 0, 0)
|
||||
[6].move(0, 8.0, 0)
|
||||
[3].move(0, 0, 16.0)
|
||||
|
||||
# Now shift the positions of all of the benzene molecules,
|
||||
# to reduce the chance that they overlap with the ethylene molecules.
|
||||
|
|
|
|||
|
|
@ -40,7 +40,7 @@ include system.in.settings
|
|||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
timestep 2.0
|
||||
dump 1 all custom 5000 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
|
||||
thermo 500
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
import "oplsaa.lt"
|
||||
|
||||
# The "oplsaa.lt" file contains force-field definitions and masses for the
|
||||
# atoms in your system. See oplsaa_lt_generator/README.TXT for details.
|
||||
# The "oplsaa.lt" file contains force-field definitions,
|
||||
# partial charges and masses for the atoms in your system.
|
||||
|
||||
# Note:
|
||||
# Atom type @atom:90 corresponds to "Aromatic C"
|
||||
|
|
@ -9,31 +9,32 @@ import "oplsaa.lt"
|
|||
|
||||
Benzene inherits OPLSAA {
|
||||
|
||||
# We just need a list of atom types and bonds.
|
||||
#
|
||||
# You don't have to specify the charge in this example because we are
|
||||
# using the OPLSAA force-field assigns this by atom-type.
|
||||
#
|
||||
# You also don't have to specify the coordinates, because
|
||||
# you are using PACKMOL to generate them for you.
|
||||
# Just leave these numbers as 0.00 for now..
|
||||
# atom-id mol-id atom-type charge X Y Z # comment
|
||||
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:C2 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:C3 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:C4 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:C5 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:C6 $mol @atom:90 0.00 0.00 0.00 0.00 # "Aromatic C"
|
||||
$atom:H11 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
$atom:H21 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
$atom:H31 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
$atom:H41 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
$atom:H51 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
$atom:H61 $mol @atom:91 0.00 0.00 0.00 0.00 # "Aromatic H-C"
|
||||
write("Data Atoms") {
|
||||
$atom:C1 $mol @atom:90 0.00 -0.739 1.189 -0.00733 # 90 = "Aromatic C"
|
||||
$atom:C2 $mol @atom:90 0.00 0.614 1.208 0.35167 # 90 = "Aromatic C"
|
||||
$atom:C3 $mol @atom:90 0.00 1.353 0.019 0.35867 # 90 = "Aromatic C"
|
||||
$atom:C4 $mol @atom:90 0.00 0.739 -1.189 0.00667 # 90 = "Aromatic C"
|
||||
$atom:C5 $mol @atom:90 0.00 -0.614 -1.208 -0.35133 # 90 = "Aromatic C"
|
||||
$atom:C6 $mol @atom:90 0.00 -1.353 -0.019 -0.35833 # 90 = "Aromatic C"
|
||||
$atom:H11 $mol @atom:91 0.00 -1.309 2.106 -0.01233 # 91 = "Aromatic H-C"
|
||||
$atom:H21 $mol @atom:91 0.00 1.088 2.14 0.62267 # 91 = "Aromatic H-C"
|
||||
$atom:H31 $mol @atom:91 0.00 2.397 0.034 0.63467 # 91 = "Aromatic H-C"
|
||||
$atom:H41 $mol @atom:91 0.00 1.309 -2.106 0.01267 # 91 = "Aromatic H-C"
|
||||
$atom:H51 $mol @atom:91 0.00 -1.088 -2.14 -0.62233 # 91 = "Aromatic H-C"
|
||||
$atom:H61 $mol @atom:91 0.00 -2.397 -0.034 -0.63533 # 91 = "Aromatic H-C"
|
||||
}
|
||||
# Note: You don't have to specify the charge in this example because
|
||||
# we are using the OPLSAA force-field assigns this by atom-type.
|
||||
# Just leave these numbers as 0.00 for now.
|
||||
# Note: LAMMPS expects an integer in the 2nd column (the Molecule-ID number).
|
||||
# If we put "$mol" there, moltemplate will generate this integer for you
|
||||
|
||||
write('Data Bond List') {
|
||||
# A list of the bonds in the molecule
|
||||
# BondID AtomID1 AtomID2
|
||||
|
||||
write("Data Bond List") {
|
||||
$bond:C12 $atom:C1 $atom:C2
|
||||
$bond:C23 $atom:C2 $atom:C3
|
||||
$bond:C34 $atom:C3 $atom:C4
|
||||
|
|
@ -48,4 +49,7 @@ Benzene inherits OPLSAA {
|
|||
$bond:C6H6 $atom:C6 $atom:H61
|
||||
}
|
||||
|
||||
# In the "Data Bond List" section we don't have to specify the bond type.
|
||||
# The bond-type will be determined by the atom type (according to "oplsaa.lt")
|
||||
|
||||
} # Benzene
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
import "oplsaa.lt"
|
||||
|
||||
# The "oplsaa.lt" file contains force-field definitions and masses for the
|
||||
# atoms in your system. See oplsaa_lt_generator/README.TXT for details.
|
||||
# The "oplsaa.lt" file contains force-field definitions,
|
||||
# partial charges and masses for the atoms in your system.
|
||||
|
||||
# Note:
|
||||
# Atom type 88 corresponds to "Alkene H2-C="
|
||||
|
|
@ -11,16 +11,24 @@ import "oplsaa.lt"
|
|||
|
||||
Ethylene inherits OPLSAA {
|
||||
|
||||
# atom-id mol-id atom-type charge X Y Z
|
||||
# atom-id mol-id atom-type charge X Y Z # comment
|
||||
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol @atom:88 0.000 -0.6695 0.000000 0.000000
|
||||
$atom:C2 $mol @atom:88 0.000 0.6695 0.000000 0.000000
|
||||
$atom:H11 $mol @atom:89 0.000 -1.234217 -0.854458 0.000000
|
||||
$atom:H12 $mol @atom:89 0.000 -1.234217 0.854458 0.000000
|
||||
$atom:H21 $mol @atom:89 0.000 1.234217 -0.854458 0.000000
|
||||
$atom:H22 $mol @atom:89 0.000 1.234217 0.854458 0.000000
|
||||
$atom:C1 $mol @atom:88 0.00 -0.6695 0.00000 0.000 #88 = "Alkene H2-C="
|
||||
$atom:C2 $mol @atom:88 0.00 0.6695 0.00000 0.000 #88 = "Alkene H2-C="
|
||||
$atom:H11 $mol @atom:89 0.00 -1.23422 -0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H12 $mol @atom:89 0.00 -1.23422 0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H21 $mol @atom:89 0.00 1.23422 -0.85446 0.000 #89 = "Alkene H-C="
|
||||
$atom:H22 $mol @atom:89 0.00 1.23422 0.85446 0.000 #89 = "Alkene H-C="
|
||||
}
|
||||
# Note: You don't have to specify the charge in this example because
|
||||
# we are using the OPLSAA force-field assigns this by atom-type.
|
||||
# Just leave these numbers as 0.00 for now.
|
||||
# Note: LAMMPS expects an integer in the 2nd column (the Molecule-ID number).
|
||||
# If we put "$mol" there, moltemplate will generate this integer for you
|
||||
|
||||
# A list of the bonds in the molecule
|
||||
# BondID AtomID1 AtomID2
|
||||
|
||||
write('Data Bond List') {
|
||||
$bond:C12 $atom:C1 $atom:C2
|
||||
|
|
@ -30,11 +38,8 @@ Ethylene inherits OPLSAA {
|
|||
$bond:C2H2 $atom:C2 $atom:H22
|
||||
}
|
||||
|
||||
# In the "Data Bond List" section we don't have to specify the bond type.
|
||||
# The bond-type will be determined by the atom type (according to "oplsaa.lt")
|
||||
|
||||
} # Ethylene
|
||||
|
||||
|
||||
|
||||
# Note: You don't need to supply the partial partial charges of the atoms.
|
||||
# If you like, just fill the fourth column with zeros ("0.000").
|
||||
# Moltemplate and LAMMPS will automatically assign the charge later
|
||||
|
||||
|
|
|
|||
|
|
@ -1,14 +1,14 @@
|
|||
12
|
||||
Benzene
|
||||
C1 5.274 1.999 -8.568
|
||||
C2 6.627 2.018 -8.209
|
||||
C3 7.366 0.829 -8.202
|
||||
C4 6.752 -0.379 -8.554
|
||||
C5 5.399 -0.398 -8.912
|
||||
C6 4.660 0.791 -8.919
|
||||
H11 4.704 2.916 -8.573
|
||||
H21 7.101 2.950 -7.938
|
||||
H31 8.410 0.844 -7.926
|
||||
H41 7.322 -1.296 -8.548
|
||||
H51 4.925 -1.330 -9.183
|
||||
H61 3.616 0.776 -9.196
|
||||
C1 -0.739 1.189 -0.00733
|
||||
C2 0.614 1.208 0.35167
|
||||
C3 1.353 0.019 0.35867
|
||||
C4 0.739 -1.189 0.00667
|
||||
C5 -0.614 -1.208 -0.35133
|
||||
C6 -1.353 -0.019 -0.35833
|
||||
H11 -1.309 2.106 -0.01233
|
||||
H21 1.088 2.14 0.62267
|
||||
H31 2.397 0.034 0.63467
|
||||
H41 1.309 -2.106 0.01267
|
||||
H51 -1.088 -2.14 -0.62233
|
||||
H61 -2.397 -0.034 -0.63533
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
This example is a simple simulation of many long alkane chains (hexadecane) in a
|
||||
box at room temperature and atmospheric pressure. Please read "WARNING.TXT".
|
||||
box near the boiling point atmospheric pressure. Please read "WARNING.TXT".
|
||||
|
||||
NOTE: This particular example uses the OPLSAA force-field
|
||||
However, moltemplate is not limited to OPLSAA.
|
||||
|
|
@ -12,7 +12,6 @@ moltemplate.sh system.lt
|
|||
|
||||
2) Run LAMMPS in this order:
|
||||
|
||||
lmp_mpi -i run.in.min # minimize the energy (to avoid atom overlap) before...
|
||||
lmp_mpi -i run.in.npt # running the simulation at constant pressure
|
||||
lmp_mpi -i run.in.nvt # running the simulation at constant temperature
|
||||
|
||||
|
|
|
|||
|
|
@ -2,9 +2,7 @@
|
|||
|
||||
This software is experimental, and the force-fields and equilbration protocols
|
||||
have not been tested carefully by me. There is no gaurantee that the simulation
|
||||
will reproduce the behavior of real hexadecane molecules,
|
||||
(or even of hexadecane molecules simulated using AMBER, which should
|
||||
be using the same force-field).
|
||||
will reproduce the behavior of real hexadecane molecules
|
||||
|
||||
# -------- REQUEST FOR HELP: --------
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 2.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.7 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 5.1 KiB After Width: | Height: | Size: 7.9 KiB |
|
|
@ -2,7 +2,7 @@
|
|||
|
||||
# First, load the OPLS force field parameters we will need.
|
||||
# These 2 files are located in the "force_fields" subdirectory
|
||||
# of the moltemplate distribution.
|
||||
# distributed with moltemplate.
|
||||
|
||||
import "oplsaa.lt" # <-- defines the standard "OPLSAA" force field
|
||||
import "loplsaa.lt" # <-- custom parameters for long alkane chains taken from
|
||||
|
|
|
|||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue