forked from lijiext/lammps
remove old moltemplate version
This commit is contained in:
Axel Kohlmeyer
parent
b11fe2eddb
commit
57aafba7c3
|
|
@ -1,28 +0,0 @@
|
|||
|
||||
Author: Andrew Jewett, Shea Group, http://www.chem.ucsb.edu/~sheagroup/
|
||||
Copyright (c) 2014, Regents of the University of California
|
||||
All rights reserved.
|
||||
Redistribution and use in source and binary forms, with or without
|
||||
modification, are permitted provided that the following conditions are met:
|
||||
|
||||
* Redistributions of source code must retain the above copyright notice,
|
||||
this list of conditions and the following disclaimer.
|
||||
* Redistributions in binary form must reproduce the above copyright notice,
|
||||
this list of conditions and the following disclaimer in the documentation
|
||||
and/or other materials provided with the distribution.
|
||||
* Neither the name of the University of California, Santa Barbara nor the
|
||||
names of its contributors may be used to endorse or promote products
|
||||
derived from this software without specific prior written permission.
|
||||
|
||||
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
|
||||
THE POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
-- Description: --
|
||||
|
||||
Moltemplate is a cross-platform text-based molecule builder for LAMMPS.
|
||||
|
||||
-- Typical usage: --
|
||||
|
||||
moltemplate.sh [-atomstyle style] [-pdb/-xyz coord_file] [-vmd] system.lt
|
||||
|
||||
-- Web page: --
|
||||
|
||||
Documentation, examples, and supporting code can be downloaded at:
|
||||
|
||||
http://www.moltemplate.org
|
||||
The most up-to-date version of moltemplate can be downloaded here.
|
||||
(After download, you can unpack the archive using:
|
||||
tar xzf moltemplate_2012-3-31.tar.gz
|
||||
The date will vary from version to version.)
|
||||
|
||||
----------------------------------------------------
|
||||
---------- INSTALLATION INSTRUCTIONS: ------------
|
||||
----------------------------------------------------
|
||||
|
||||
This directory should contain two folders:
|
||||
src/ <-- location of all python and bash scripts
|
||||
common/ <-- location of shared force fields and molecules
|
||||
|
||||
The ``moltemplate.sh'' script and the python scripts that it invokes are
|
||||
located in the ``src/'' subdirectory. You should update your PATH environment
|
||||
variable to include this directory.
|
||||
|
||||
If you do not know what a PATH environment variable is, read:
|
||||
http://www.linfo.org/path_env_var.html
|
||||
(I receive this question often.)
|
||||
|
||||
It is also a good idea to set your MOLTEMPLATE_PATH environment variable to
|
||||
point to the ``common/'' subdirectory.
|
||||
(Force fields and commonly used molecules will eventually be located here.)
|
||||
|
||||
-- Installation example ---
|
||||
|
||||
Suppose the directory with this README.TXT file is located at ~/moltemplate.
|
||||
|
||||
If you use the bash shell, typically you would edit your
|
||||
~/.profile, ~/.bash_profile or ~/.bashrc files to contain the following lines:
|
||||
export PATH="$PATH:$HOME/moltemplate/src"
|
||||
export MOLTEMPLATE_PATH="$HOME/moltemplate/common"
|
||||
If you use the tcsh shell, typically you would edit your
|
||||
~/.login, ~/.cshrc, or ~/.tcshrc files to contain the following lines:
|
||||
setenv PATH "$PATH:$HOME/moltemplate/src"
|
||||
setenv MOLTEMPLATE_PATH "$HOME/moltemplate/common"
|
||||
|
||||
-- Requirements: --
|
||||
|
||||
Moltemplate requires the Bourne-shell, and a recent version of python
|
||||
(2.7, 3.0 or higher), and can run on OS X, linux, or windows (if a
|
||||
suitable shell environment has been installed).
|
||||
|
||||
-- License: --
|
||||
|
||||
Moltemplate is available under the terms of the open-source 3-clause BSD
|
||||
license. (See LICENSE.TXT.)
|
||||
|
|
@ -1,24 +0,0 @@
|
|||
These are examples for the "moltemplate" molecule builder for LAMMPS.
|
||||
http://www.moltemplate.org
|
||||
|
||||
Each directory contains one or more examples.
|
||||
|
||||
Each example directory contains:
|
||||
|
||||
images/ This folder has pictures of the molecules in the system
|
||||
moltemplate_files/ This folder contains LT files and other auxiliary files
|
||||
README_setup.sh Instructions for how to use moltemplate (executable)
|
||||
README_visualize.txt Instructions for viewing in DATA/DUMP files in VMD
|
||||
|
||||
...and one or more LAMMPS input scripts with names like
|
||||
|
||||
run.in.min
|
||||
run.in.npt
|
||||
run.in.nvt
|
||||
|
||||
You can run these scripts using
|
||||
lmp_linux -i run.in.npt
|
||||
(The name of your lammps binary, "lmp_linux" in this example, may vary.
|
||||
Sometimes, these scripts must be run in a certain order. For example
|
||||
it may be necessary to run run.in.min to minimize the system before you can use run.in.npt, and later run.in.nvt. The README_run.sh file in each subdirectory
|
||||
specifies indicates the order. These files have not been optimized.)
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
# -------- REQUIREMENTS: ---------
|
||||
# 1) This example requires the "MANYBODY" package.
|
||||
# As of 2012-9, it is included by default, but this may change in the future.
|
||||
# If lammps complains of a missing pair style enter "make yes-MANYBODY"
|
||||
# into the shell before compiling lammps. For details see:
|
||||
# http://lammps.sandia.gov/doc/Section_start.html#start_3
|
||||
This is a relatively complex example containing two different types of
|
||||
molecules, and a hybrid of Lennard-Jones and 3-body SW "pair" styles.
|
||||
|
||||
The cyclododecane molecule uses the
|
||||
TraPPE force field for hydrocarbon chains.
|
||||
The parameters for the TraPPE force field are
|
||||
in a file named "trappe1998.lt" which should be
|
||||
located in the MOLTEMPLATE_PATH.
|
||||
(See moltemplate installation instructions.)
|
||||
|
||||
The water solvent is implemented using the 3-body single-particle
|
||||
coarse-grained "mW" water model:
|
||||
Molinero, V. and Moore, E.B., J. Phys. Chem. B 2009, 113, 4008-4016
|
||||
|
||||
More detailed instructions on how to build LAMMPS input files and
|
||||
run a short simulation are provided in other README files.
|
||||
|
||||
step 1)
|
||||
README_setup.sh
|
||||
|
||||
step 2)
|
||||
README_run.sh
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -- Prerequisites: --
|
||||
# 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, system.in.sw
|
||||
# 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 # 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 are ignored when beginning the simulation at constant volume.
|
||||
# This can be fixed. Read "run.in.nvt" for equilibration instructions.)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# 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)
|
||||
|
|
@ -1,25 +0,0 @@
|
|||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh -a "@atom:/WatMW/mW 1" system.lt
|
||||
|
||||
# Here we just want to make sure that the "mW" atom type is assigned to
|
||||
# number "1". It should be by default, so usually you can leave out
|
||||
# -a "@atom:/WatMW/mW 1".
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
# These files are the input files directly read by LAMMPS. Move them to
|
||||
# the parent directory (or wherever you plan to run the simulation).
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# 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 ../
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
------- 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.
|
Before Width: | Height: | Size: 51 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 38 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 38 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 3.1 KiB |
|
|
@ -1,11 +0,0 @@
|
|||
# Use this command to generate the LAMMPS input files:
|
||||
|
||||
moltemplate.sh -a "@atom:/WatMW/mW 1" system.lt
|
||||
|
||||
# The -a argument insures that the "mW" atom type is assigned to "1".
|
||||
# (This is necessary for the pair_coeff command to work.
|
||||
# See system.lt for details.)
|
||||
|
||||
# Note: To get rid of the annoying "atom_style unspecified warnings,
|
||||
# use the "-atomstyle" command line argument, as in:
|
||||
# moltemplate.sh -atomstyle full -a "@atom:/WatMW/mW 1" system.lt
|
||||
|
|
@ -1,55 +0,0 @@
|
|||
import "trappe1998.lt"
|
||||
|
||||
# The "trappe1998.lt" file is usually located in $MOLTEMPLATE_PATH (and is
|
||||
# distributed with moltemplate. See the "Installation" section in the manual.)
|
||||
# It contains definitions of the atoms "CH2", "CH3", and "CH4", as well
|
||||
# as "saturated" bonds, and the parameters for (bonded/nonbonded)
|
||||
# interactions between these atoms (all enclosed within the "TraPPE" namespace).
|
||||
|
||||
|
||||
Cyclododecane {
|
||||
|
||||
write('Data Atoms') {
|
||||
$atom:C1 $mol:. @atom:TraPPE/CH2 0.0 0.00000 2.94118 0.0
|
||||
$atom:C2 $mol:. @atom:TraPPE/CH2 0.0 0.00000 2.54714 1.47059
|
||||
$atom:C3 $mol:. @atom:TraPPE/CH2 0.0 0.00000 1.47059 2.54714
|
||||
$atom:C4 $mol:. @atom:TraPPE/CH2 0.0 0.00000 0.0 2.94118
|
||||
$atom:C5 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -1.47059 2.54714
|
||||
$atom:C6 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -2.54714 1.47059
|
||||
$atom:C7 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -2.94118 0.0
|
||||
$atom:C8 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -2.54714 -1.47059
|
||||
$atom:C9 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -1.47059 -2.54714
|
||||
$atom:C10 $mol:. @atom:TraPPE/CH2 0.0 0.00000 -0.0 -2.94118
|
||||
$atom:C11 $mol:. @atom:TraPPE/CH2 0.0 0.00000 1.47059 -2.54714
|
||||
$atom:C12 $mol:. @atom:TraPPE/CH2 0.0 0.00000 2.54714 -1.47059
|
||||
}
|
||||
|
||||
# The "." in "$mol:." refers to the current object's molecule ID,
|
||||
# and "@atom:TraPPE/CH2" refers to the "CH2" atom-type defined in TraPPE
|
||||
|
||||
write('Data Bonds') {
|
||||
$bond:bond1 @bond:TraPPE/saturated $atom:C1 $atom:C2
|
||||
$bond:bond2 @bond:TraPPE/saturated $atom:C2 $atom:C3
|
||||
$bond:bond3 @bond:TraPPE/saturated $atom:C3 $atom:C4
|
||||
$bond:bond4 @bond:TraPPE/saturated $atom:C4 $atom:C5
|
||||
$bond:bond5 @bond:TraPPE/saturated $atom:C5 $atom:C6
|
||||
$bond:bond6 @bond:TraPPE/saturated $atom:C6 $atom:C7
|
||||
$bond:bond7 @bond:TraPPE/saturated $atom:C7 $atom:C8
|
||||
$bond:bond8 @bond:TraPPE/saturated $atom:C8 $atom:C9
|
||||
$bond:bond9 @bond:TraPPE/saturated $atom:C9 $atom:C10
|
||||
$bond:bond10 @bond:TraPPE/saturated $atom:C10 $atom:C11
|
||||
$bond:bond11 @bond:TraPPE/saturated $atom:C11 $atom:C12
|
||||
$bond:bond12 @bond:TraPPE/saturated $atom:C12 $atom:C1
|
||||
}
|
||||
|
||||
} # Cyclododecane
|
||||
|
||||
|
||||
# coordinates in the "Data Atoms" section generated by this python code:
|
||||
# from math import *
|
||||
# bond_length=1.54
|
||||
# N=12
|
||||
# R=(N*bond_length)/(2*pi)
|
||||
# for i in range(0,N):
|
||||
# print('$atom:C'+str(i+1)+' $mol:... @atom:TraPPE/CH2 0.0 0.00000 '+
|
||||
# str(round(R*cos(i*2*pi/N),5))+' '+str(round(R*sin(i*2*pi/N),5)))
|
||||
|
|
@ -1,62 +0,0 @@
|
|||
# This is a relatively complex example containing two different types of
|
||||
# molecules, and a hybrid of Lennard-Jones and 3-body SW "pair" styles.
|
||||
|
||||
import "watmw.lt"
|
||||
import "cyclododecane.lt"
|
||||
|
||||
write_once("Data Boundary") {
|
||||
0.000000 48.000 xlo xhi
|
||||
0.000000 48.000 ylo yhi
|
||||
0.000000 48.000 zlo zhi
|
||||
}
|
||||
|
||||
wat = new WatMW [12].move(0, 0, 4.0)
|
||||
[12].move(0, 4.0, 0)
|
||||
[12].move(4.0, 0, 0)
|
||||
|
||||
cyclododecane = new Cyclododecane [4].move(0, 0, 12.0)
|
||||
[4].move(0, 12.0, 0)
|
||||
[4].move(12.0, 0, 0)
|
||||
|
||||
# (Move them by (6.0,6.0,6.0) to avoid overlap with the water.)
|
||||
cyclododecane[*][*][*].move(6.0,6.0,6.0)
|
||||
|
||||
write_once("In Init") {
|
||||
# -- Tell LAMMPS we want to use two different pair styles
|
||||
# -- (This overrides earlier settings.)
|
||||
pair_style hybrid sw lj/charmm/coul/charmm 9.0 11.0 9.0 11.0
|
||||
}
|
||||
|
||||
|
||||
write_once("In Settings") {
|
||||
# -- Now indicate which atom type(s) are simulated using the "sw" pair style
|
||||
# -- In this case only one of the atom types is used (the mW water "atom").
|
||||
|
||||
pair_coeff * * sw system.in.sw mW NULL NULL NULL
|
||||
|
||||
# -- Unfortunately LAMMPS itself does not understand molemlate syntax, so
|
||||
# -- the atoms are identified by order in the list, not by name. (The "mW"
|
||||
# -- refers to to an identifier in the system.in.sw file, not watmw.lt.)
|
||||
# -- This command says that the first atom type corresponds to the "mW"
|
||||
# -- atom in system.in.sw, and to ignore the remaining three atom types
|
||||
# -- (correspond to the CH2, CH3, CH4 atom types defined in trappe1998.lt.
|
||||
# -- We don't want to use the "sw" force field for interactions involving
|
||||
# -- these atom types, so we put "NULL" there.)
|
||||
# -- Note: For this to work, you should probably run moltemplate this way:
|
||||
# -- moltemplate.sh -a "@atom:WatMW/mW 1" system.lt
|
||||
# -- This assigns the atom type named @atom:WatMW/mW to 1 (the first atom)
|
||||
}
|
||||
|
||||
|
||||
|
||||
# -- Somewhere we must eventually define interactions
|
||||
# -- between atoms from different molecule types
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH2 lj/charmm/coul/charmm 0.11914784667210733 3.558
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH3 lj/charmm/coul/charmm 0.17390830404497651 3.458
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH4 lj/charmm/coul/charmm 0.21371654257637612 3.448
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
# This file stores complete LAMMPS data for the TraPPE model of saturated
|
||||
# hydrocarbon chains. In this "united-atom" model, each methyl group is
|
||||
# represented by a single atom. Forces between "atoms" are taken from the
|
||||
# TraPPE force-field. (J Phys Chem B, 1998, volume 102, pp.2569-2577)
|
||||
|
||||
TraPPE {
|
||||
|
||||
write_once("In Init") {
|
||||
# -- Default styles for "TraPPE" --
|
||||
units real
|
||||
atom_style full
|
||||
# (Hybrid force field styles were used for portability.)
|
||||
bond_style hybrid harmonic
|
||||
angle_style hybrid harmonic
|
||||
dihedral_style hybrid opls
|
||||
improper_style none
|
||||
pair_style hybrid lj/charmm/coul/charmm 9.0 11.0 9.0 11.0
|
||||
pair_modify mix arithmetic
|
||||
special_bonds lj 0.0 0.0 0.0
|
||||
}
|
||||
|
||||
write_once("Data Masses") {
|
||||
@atom:CH2 14.1707
|
||||
@atom:CH3 15.2507
|
||||
@atom:CH4 16.3307
|
||||
}
|
||||
|
||||
write_once("Data Angles By Type") {
|
||||
@angle:backbone @atom:CH? @atom:CH? @atom:CH? @bond:saturated @bond:saturated
|
||||
}
|
||||
|
||||
write_once("Data Dihedrals By Type") {
|
||||
@dihedral:backbone @atom:CH? @atom:CH? @atom:CH? @atom:CH? @bond:saturated @bond:saturated @bond:saturated
|
||||
}
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:CH2 @atom:CH2 lj/charmm/coul/charmm 0.091411522 3.95
|
||||
pair_coeff @atom:CH3 @atom:CH3 lj/charmm/coul/charmm 0.194746286 3.75
|
||||
pair_coeff @atom:CH4 @atom:CH4 lj/charmm/coul/charmm 0.294106636 3.73
|
||||
bond_coeff @bond:saturated harmonic 120.0 1.54
|
||||
angle_coeff @angle:backbone harmonic 62.0022 114
|
||||
dihedral_coeff @dihedral:backbone opls 1.411036 -0.271016 3.145034 0.0
|
||||
}
|
||||
|
||||
write_once("In Settings") {
|
||||
group TraPPE type @atom:CH2 @atom:CH3 @atom:CH4
|
||||
}
|
||||
|
||||
} # class TraPPE
|
||||
|
||||
|
|
@ -1,80 +0,0 @@
|
|||
# This is a relatively complex example containing two different types of
|
||||
# molecules, and a hybrid of Lennard-Jones and 3-body SW "pair" styles.
|
||||
|
||||
import "watmw.lt"
|
||||
import "cyclododecane.lt"
|
||||
|
||||
write_once("Data Boundary") {
|
||||
0.000000 48.000 xlo xhi
|
||||
0.000000 48.000 ylo yhi
|
||||
0.000000 48.000 zlo zhi
|
||||
}
|
||||
|
||||
wat = new WatMW [12].move(0, 0, 4.0)
|
||||
[12].move(0, 4.0, 0)
|
||||
[12].move(4.0, 0, 0)
|
||||
|
||||
cyclododecane = new Cyclododecane [4].move(0, 0, 12.0)
|
||||
[4].move(0, 12.0, 0)
|
||||
[4].move(12.0, 0, 0)
|
||||
|
||||
# (Move them by (6.0,6.0,6.0) to avoid overlap with the water.)
|
||||
cyclododecane[*][*][*].move(6.0,6.0,6.0)
|
||||
|
||||
write_once("In Init") {
|
||||
# -- Tell LAMMPS we want to use two different pair styles
|
||||
# -- (This overrides earlier settings.)
|
||||
pair_style hybrid sw lj/charmm/coul/charmm 9.0 11.0 9.0 11.0
|
||||
}
|
||||
|
||||
|
||||
|
||||
write_once("In Settings") {
|
||||
# -- Now indicate which atom type(s) are simulated using the "sw" pair style
|
||||
# -- In this case only one of the atom types is used (the mW water "atom").
|
||||
|
||||
pair_coeff * * sw system.in.sw mW NULL NULL NULL
|
||||
|
||||
# -- Unfortunately LAMMPS itself does not understand molemlate syntax, so
|
||||
# -- the atoms are identified by order in the list, not by name. (The "mW"
|
||||
# -- refers to to an identifier in the system.in.sw file, not watmw.lt.)
|
||||
# -- This command says that the first atom type corresponds to the "mW"
|
||||
# -- atom in system.in.sw, and to ignore the remaining three atom types
|
||||
# -- (correspond to the CH2, CH3, CH4 atom types defined in trappe1998.lt.
|
||||
# -- We don't want to use the "sw" force field for interactions involving
|
||||
# -- these atom types, so we put "NULL" there.)
|
||||
#
|
||||
# For this to work, the first atom type (assigned to "1")
|
||||
# must refer to the "mW" atom type (defined in watmw.lt).
|
||||
# (This is why we included "watmw.lt" first, to insure that the
|
||||
# atom counters in WatMW are assinged first, starting with 1.)
|
||||
# Alternately we can further insure that this happens, it's
|
||||
# a good idea to run moltemplate.sh using the "-a" argument:
|
||||
# moltemplate.sh -a "@atom:/WatMW/mW 1" system.lt
|
||||
# This assigns the atom type named @atom:/WatMW/mW to 1
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
# -- Somewhere we must eventually define interactions
|
||||
# -- between atoms from different molecule types
|
||||
# -- Now define interactions between DIFFERENT molecules
|
||||
# Note: In the SPC/E model, the epsilon,sigma parameters for water is 0.1553
|
||||
# 3.166. As a crude guess, I chose the LJ parameters for the interaction
|
||||
# between water & the CH2,CH3,CH4 atoms using Lorentz-Berthelot mixing rules
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH2 lj/charmm/coul/charmm 0.11914784667210733 3.558
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH3 lj/charmm/coul/charmm 0.17390830404497651 3.458
|
||||
pair_coeff @atom:WatMW/mW @atom:TraPPE/CH4 lj/charmm/coul/charmm 0.21371654257637612 3.448
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,54 +0,0 @@
|
|||
# This file stores LAMMPS data for the "mW" water model.
|
||||
# (Molinero, V. and Moore, E.B., J. Phys. Chem. B 2009, 113, 4008-4016)
|
||||
#
|
||||
# In this model, each water molecule is represented by a single "mW" particle.
|
||||
# These particles interact with their neighbors via 3-body Stillinger-Weber
|
||||
# forces whose parameters are tuned to mimic directional hydrogen-bonding
|
||||
# in liquid water (as well as hexagonal ice, type II ice, and
|
||||
# low-density super-cooled liquid/amorphous water phases).
|
||||
|
||||
WatMW {
|
||||
write("Data Atoms") {
|
||||
$atom:mW $mol:. @atom:mW 0.0 0.0 0.0 0.0
|
||||
}
|
||||
|
||||
write_once("Data Masses") {
|
||||
@atom:mW 18.02
|
||||
}
|
||||
|
||||
write_once("system.in.sw") {
|
||||
mW mW mW 6.189 2.3925 1.8 23.15 1.2 -0.333333333 7.049556277 0.602224558 4 0 0
|
||||
}
|
||||
|
||||
write_once("In Init") {
|
||||
# -- Default styles for "WatMW" --
|
||||
units real
|
||||
pair_style sw
|
||||
}
|
||||
|
||||
write_once("In Settings") {
|
||||
# --Now indicate which atom type(s) are simulated using the "sw" pair style
|
||||
# -- In this case only one of the atom types is used (the mW water "atom").
|
||||
|
||||
pair_coeff * * sw system.in.sw mW NULL NULL NULL
|
||||
|
||||
# -- Unfortunately LAMMPS itself does not understand molemlate syntax, so
|
||||
# -- the atoms are identified by order in the list, not by name. (The "mW"
|
||||
# -- refers to to an identifier in the system.in.sw file, not watmw.lt.)
|
||||
# -- This command says that the first atom type corresponds to the "mW"
|
||||
# -- atom in system.in.sw, and to ignore the remaining three atom types
|
||||
# -- (correspond to the CH2, CH3, CH4 atom types defined in trappe1998.lt.
|
||||
# -- We don't want to use the "sw" force field for interactions involving
|
||||
# -- these atom types, so we put "NULL" there.)
|
||||
# -- Note: For this to work, you should probably run moltemplate this way:
|
||||
# -- moltemplate.sh -a "@atom:WatMW/mW 1" system.lt
|
||||
# -- This assigns the atom type named @atom:WatMW/mW to 1 (the first atom)
|
||||
}
|
||||
|
||||
# -- optional --
|
||||
|
||||
write_once("In Settings") {
|
||||
group WatMW type @atom:mW #(Atoms of this type belong to the "WatMW" group)
|
||||
}
|
||||
|
||||
} # WatMW
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
# run.in.npt
|
||||
#
|
||||
# -- Usage --
|
||||
#
|
||||
# lmp_linux -i run.in.npt
|
||||
# (assuming lmp_linux is the name of your lammps binary)
|
||||
#
|
||||
# -- Prerequisite Input Files: --
|
||||
# systen.data, system.in.init, system.in.settings, system.in.sw
|
||||
#
|
||||
# You can generate these files with this command:
|
||||
# moltemplate.sh -a "@atom:/WatMW/mW 1" system.lt
|
||||
# ---------------------------------
|
||||
|
||||
# ----- Init Section -----
|
||||
|
||||
|
||||
include system.in.init
|
||||
|
||||
|
||||
# ----- Atom Definition Section -----
|
||||
|
||||
|
||||
read_data system.data
|
||||
|
||||
|
||||
# ----- Settings Section -----
|
||||
|
||||
|
||||
include system.in.settings
|
||||
|
||||
|
||||
# ----- 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.
|
||||
|
||||
minimize 1.0e-5 1.0e-7 100000 400000
|
||||
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
|
||||
timestep 2.0 # <- This can be increased to 5.0 or 10.0 for bulk water
|
||||
dump 1 all custom 500 traj_npt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0
|
||||
|
||||
|
||||
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
|
||||
thermo 500 # time interval for printing out "thermo" data
|
||||
|
||||
run 200000
|
||||
|
||||
write_data system_after_npt.data
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also.)
|
||||
|
|
@ -1,81 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
# 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.)
|
||||
|
||||
# run.in.nvt
|
||||
#
|
||||
# -- Usage --
|
||||
#
|
||||
# lmp_g++ -i run.in.nvt
|
||||
# (assuming lmp_g++ is the name of your lammps binary)
|
||||
#
|
||||
# -- Prerequisite Input Files: --
|
||||
# systen.data, system.in.init, system.in.settings, system.in.sw
|
||||
# system_after_npt.data
|
||||
#
|
||||
# You can generate these files using this procedure
|
||||
#
|
||||
# moltemplate.sh -a "@atom:/WatMW/mW 1" system.lt
|
||||
#
|
||||
# lmp_linux -i run.in.npt
|
||||
|
||||
# ---------------------------------
|
||||
|
||||
|
||||
# -- init section --
|
||||
|
||||
|
||||
include system.in.init
|
||||
|
||||
|
||||
|
||||
# -- atom definition section --
|
||||
|
||||
|
||||
# Read the coordinates generated by an earlier NPT simulation
|
||||
|
||||
read_data system_after_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
|
||||
|
||||
|
||||
# -- run section --
|
||||
|
||||
|
||||
timestep 2.0
|
||||
dump 1 all custom 1000 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
dump 2 TraPPE custom 1000 traj_alkane_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
|
||||
|
||||
# The following commands are useful if you want to calculate the distribution
|
||||
# of alkane-chain radius-of-gyration at a given temperature & pressure.
|
||||
#compute cRg TraPPE gyration
|
||||
#variable vRg equal c_cRg
|
||||
#compute cPE all pe
|
||||
#variable vPE equal c_cPE
|
||||
#fix FprintPE all print 1000 "${vPE}" file U.dat
|
||||
#fix FprintRg all print 1000 "${vRg}" file Rg.dat
|
||||
|
||||
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
|
||||
thermo 1000 # time interval for printing out "thermo" data
|
||||
#thermo_modify flush yes
|
||||
|
||||
restart 100000 restart_nvt
|
||||
|
||||
run 1000000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
||||
|
||||
|
|
@ -1,15 +0,0 @@
|
|||
|
||||
This example of the formation of a coarse-grained DPPC lipid-bilayer uses the
|
||||
Martini force-field v2.0 (2013-10), was provided by Saeed Momeni Bashusqeh.
|
||||
In this example, the initial coordinates are generated by PACKMOL.
|
||||
If you prefer, there is also an example of a Martini DPPC bilayer
|
||||
which has been preassembled using moltemplate commands.
|
||||
(That example does not require PACKMOL.)
|
||||
|
||||
step 1)
|
||||
To build the files which LAMMPS needs, follow the instructions in:
|
||||
README_setup.sh
|
||||
|
||||
step 2)
|
||||
To run LAMMPS with these files, follow these instructions:
|
||||
README_run.sh
|
||||
|
|
@ -1,21 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -------- PREREQUISITES: --------
|
||||
# The 2 files "run.in.min", "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.npt # 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)
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
|
||||
# Create the coordinates of the atoms using PACKMOL
|
||||
cd packmol_files
|
||||
|
||||
packmol < mix_lipids+water.inp
|
||||
mv -f system.xyz ../moltemplate_files/
|
||||
|
||||
cd ..
|
||||
|
||||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh -xyz system.xyz system.lt
|
||||
|
||||
# This will generate 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 ../
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
------- 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.0 -0.5}
|
||||
pbc box -shiftcenterrel {0.0 0.0 -0.5} -style tubes -width 0.75
|
||||
|
||||
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.
|
Before Width: | Height: | Size: 13 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 43 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 44 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 44 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 46 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 1.6 KiB |
|
|
@ -1,78 +0,0 @@
|
|||
DPPC {
|
||||
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
bond_style hybrid harmonic
|
||||
angle_style hybrid cosine/squared
|
||||
dihedral_style none
|
||||
improper_style none
|
||||
pair_style hybrid lj/gromacs/coul/gromacs 9 12 0.000001 12
|
||||
special_bonds lj/coul 0 1 1
|
||||
dielectric 15
|
||||
neigh_modify every 10
|
||||
}
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:1 $mol:. @atom:Q0 1.0 9.09 9.83 0.75
|
||||
$atom:2 $mol:. @atom:Qa -1.0 5.68 7.31 0.00
|
||||
$atom:3 $mol:. @atom:Na 0.0 5.50 5.61 3.28
|
||||
$atom:4 $mol:. @atom:Na 0.0 6.65 2.22 3.04
|
||||
$atom:5 $mol:. @atom:C1 0.0 5.15 7.65 7.06
|
||||
$atom:6 $mol:. @atom:C1 0.0 7.91 7.17 10.54
|
||||
$atom:7 $mol:. @atom:C1 0.0 9.24 8.25 14.96
|
||||
$atom:8 $mol:. @atom:C1 0.0 12.19 11.75 16.38
|
||||
$atom:9 $mol:. @atom:C1 0.0 5.52 1.61 7.40
|
||||
$atom:10 $mol:. @atom:C1 0.0 6.53 2.26 12.25
|
||||
$atom:11 $mol:. @atom:C1 0.0 3.51 1.81 16.01
|
||||
$atom:12 $mol:. @atom:C1 0.0 0.00 0.00 18.19
|
||||
}
|
||||
write("Data Bonds") {
|
||||
$bond:b1 @bond:Bo $atom:1 $atom:2
|
||||
$bond:b2 @bond:Bo $atom:2 $atom:3
|
||||
$bond:b3 @bond:Short $atom:3 $atom:4
|
||||
$bond:b4 @bond:Bo $atom:3 $atom:5
|
||||
$bond:b5 @bond:Bo $atom:5 $atom:6
|
||||
$bond:b6 @bond:Bo $atom:6 $atom:7
|
||||
$bond:b7 @bond:Bo $atom:7 $atom:8
|
||||
$bond:b8 @bond:Bo $atom:4 $atom:9
|
||||
$bond:b9 @bond:Bo $atom:9 $atom:10
|
||||
$bond:b10 @bond:Bo $atom:10 $atom:11
|
||||
$bond:b11 @bond:Bo $atom:11 $atom:12
|
||||
}
|
||||
write("Data Angles") {
|
||||
$angle:a1 @angle:An1 $atom:1 $atom:2 $atom:3
|
||||
$angle:a2 @angle:An2 $atom:2 $atom:3 $atom:5
|
||||
$angle:a3 @angle:An2 $atom:2 $atom:3 $atom:4
|
||||
$angle:a4 @angle:An2 $atom:4 $atom:3 $atom:5
|
||||
$angle:a5 @angle:An1 $atom:3 $atom:4 $atom:9
|
||||
$angle:a6 @angle:An1 $atom:3 $atom:5 $atom:6
|
||||
$angle:a7 @angle:An1 $atom:5 $atom:6 $atom:7
|
||||
$angle:a8 @angle:An1 $atom:6 $atom:7 $atom:8
|
||||
$angle:a9 @angle:An1 $atom:4 $atom:9 $atom:10
|
||||
$angle:a10 @angle:An1 $atom:9 $atom:10 $atom:11
|
||||
$angle:a11 @angle:An1 $atom:10 $atom:11 $atom:12
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:Q0 72.0
|
||||
@atom:Qa 72.0
|
||||
@atom:Na 72.0
|
||||
@atom:C1 72.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:Q0 @atom:Q0 lj/gromacs/coul/gromacs 0.8365200764818 4.7
|
||||
pair_coeff @atom:Q0 @atom:Qa lj/gromacs/coul/gromacs 1.0755258126195 4.7
|
||||
pair_coeff @atom:Q0 @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Q0 @atom:C1 lj/gromacs/coul/gromacs 0.4780114722753 6.2
|
||||
pair_coeff @atom:Qa @atom:Qa lj/gromacs/coul/gromacs 1.1950286806883 4.7
|
||||
pair_coeff @atom:Qa @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Qa @atom:C1 lj/gromacs/coul/gromacs 0.4780114722753 6.2
|
||||
pair_coeff @atom:Na @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Na @atom:C1 lj/gromacs/coul/gromacs 0.6453154875717 4.7
|
||||
pair_coeff @atom:C1 @atom:C1 lj/gromacs/coul/gromacs 0.8365200764818 4.7
|
||||
bond_coeff @bond:Bo harmonic 1.4937858508604 4.7
|
||||
bond_coeff @bond:Short harmonic 1.4937858508604 3.7
|
||||
angle_coeff @angle:An1 cosine/squared 2.9875717017208 180
|
||||
angle_coeff @angle:An2 cosine/squared 2.9875717017208 120
|
||||
}
|
||||
} #DPPC
|
||||
|
|
@ -1,24 +0,0 @@
|
|||
import "water.lt"
|
||||
import "lipid.lt"
|
||||
|
||||
# The lipids and water must be listed instantiated in the same order
|
||||
# they appear in the packmol_files/mix_lipids+water.inp file:
|
||||
|
||||
lipids = new DPPC[300]
|
||||
|
||||
waters = new MW[6000]
|
||||
|
||||
|
||||
write_once("Data Boundary") {
|
||||
0 100.0 xlo xhi
|
||||
0 100.0 ylo yhi
|
||||
0 100.0 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Q0 lj/gromacs/coul/gromacs 1.3384321223709 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Qa lj/gromacs/coul/gromacs 1.3384321223709 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/C1 lj/gromacs/coul/gromacs 0.4780114722753 4.7
|
||||
}
|
||||
|
|
@ -1,17 +0,0 @@
|
|||
MW {
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
pair_style hybrid lj/gromacs/coul/gromacs 9 12 0.000001 12
|
||||
}
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:1 $mol:. @atom:P4 0 0 0 0
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:P4 72.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:P4 @atom:P4 lj/gromacs/coul/gromacs 1.1950286806883 4.7
|
||||
}
|
||||
} #MW
|
||||
|
|
@ -1,8 +0,0 @@
|
|||
You can use packmol to create a file containing the atomic coordinates
|
||||
for a system of coarse-grained lipids mixed with water using this command:
|
||||
|
||||
If it takes too long for packmol to run, try lowering the tolerance.
|
||||
(tolerance 2.0 should work)
|
||||
|
||||
packmol < mix_lipids+water.inp
|
||||
|
||||
|
|
@ -1,14 +0,0 @@
|
|||
12
|
||||
DPPC
|
||||
Q0 9.09 9.83 0.75
|
||||
Qa 5.68 7.31 0.00
|
||||
Na 5.50 5.61 3.28
|
||||
Na 6.65 2.22 3.04
|
||||
C1 5.15 7.65 7.06
|
||||
C1 7.91 7.17 10.54
|
||||
C1 9.24 8.25 14.96
|
||||
C1 12.19 11.75 16.38
|
||||
C1 5.52 1.61 7.40
|
||||
C1 6.53 2.26 12.25
|
||||
C1 3.51 1.81 16.01
|
||||
C1 0.00 0.00 18.19
|
||||
|
|
@ -1,34 +0,0 @@
|
|||
#
|
||||
# A mixture of coarse-grained (martini) DPPC (lipid) and water.
|
||||
#
|
||||
|
||||
# All the atoms from diferent molecules will be separated at least 3.0
|
||||
# Anstroms at the solution.
|
||||
|
||||
tolerance 3.0 # minimal distance between atoms in different molecules
|
||||
# (you should also consider changing the "discale"
|
||||
# parameter. I think discale=1.0 by default.)
|
||||
seed 123 # seed for random number generator
|
||||
|
||||
# The file type of input and output files is XYZ
|
||||
|
||||
filetype xyz
|
||||
|
||||
# The name of the output file
|
||||
|
||||
output system.xyz
|
||||
|
||||
# DPPC (lipid) molecules and water molecules will be put in a box
|
||||
# defined by the minimum coordinates x, y and z = 0 0 0. and maximum
|
||||
# coordinates 100 100 100. (Box size: 100x100x100)
|
||||
|
||||
structure lipid.xyz
|
||||
number 300
|
||||
inside box 0.0 0.0 0.0 100.0 100.0 100.0
|
||||
end structure
|
||||
|
||||
structure water.xyz
|
||||
number 6000
|
||||
inside box 0.0 0.0 0.0 100.0 100.0 100.0
|
||||
end structure
|
||||
|
||||
|
|
@ -1,3 +0,0 @@
|
|||
1
|
||||
water
|
||||
W 0 0 0
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data "system.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
thermo 5
|
||||
dump 1 all custom 100 traj_equib0_min.lammpstrj id mol type x y z ix iy iz
|
||||
minimize 1.0e-4 1.0e-6 100000 400000
|
||||
|
||||
write_data system_after_min.data
|
||||
|
|
@ -1,117 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
# 2) You must minimize the coordinates using by running lammps witn
|
||||
# run.in.min
|
||||
#
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
#read_data "system.data"
|
||||
read_data "system_after_min.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
print "---------------------------------------------------------------------------"
|
||||
print "I often use Langevin dynamics initially at high temperatures and small"
|
||||
print "timesteps to relax the system. It seems more stable than Nose-Hoover."
|
||||
print "(This is probably not necessary.)"
|
||||
print "---------------------------------------------------------------------------"
|
||||
print " Note: Turning off 1-3 interactions during equilibration. (Turn them on later)"
|
||||
print "---------------------------------------------------------------------------"
|
||||
special_bonds lj/coul 0 0 1
|
||||
|
||||
fix fxlan all langevin 450.0 450.0 100 123456 # temp: 450 K
|
||||
fix fxnph all nph iso 170.0 170.0 10000.0 # pressure: 170 barr
|
||||
thermo 100
|
||||
dump dmAll all custom 2000 traj_equib1_npt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
timestep 1.0 # (safer to use a small timestep initially)
|
||||
run 1000
|
||||
timestep 2.0
|
||||
run 1000
|
||||
timestep 5.0
|
||||
run 1000
|
||||
timestep 10.0
|
||||
run 2000
|
||||
timestep 20.0
|
||||
run 5000
|
||||
timestep 30.0 # (40.0 should be possible for lipid systems)
|
||||
|
||||
run 50000
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnph
|
||||
undump dmAll
|
||||
|
||||
print "---------------------------------------------------------------------------"
|
||||
print "--- Now continue the simulation using a Nose-Hoover Thermostat/Barostat ---"
|
||||
print "---------------------------------------------------------------------------"
|
||||
|
||||
velocity all zero linear # <- eliminate drift due to non-zero total momentum
|
||||
#fix 1 all momentum 1000 linear 1 1 1 # also works
|
||||
thermo 100
|
||||
#thermo_modify flush yes
|
||||
|
||||
# temperature: 300 K, pressure: 1 barr
|
||||
fix fxnpt all npt temp 300.0 300.0 3000.0 iso 1.0 1.0 30000.0 drag 1.0
|
||||
dump dmAll all custom 5000 traj_equib2_npt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
run 50000
|
||||
unfix fxnpt
|
||||
undump dmAll
|
||||
|
||||
|
||||
# Pressure was equilibrated at 300K, but I initially run the simulation at
|
||||
# higher temperatures temporarily to speed up the proces of bilayer formation.
|
||||
# (Note: The boiling point of Martini-water is
|
||||
# I have to keep the volume constant (NVT) to prevent the water from boiling.
|
||||
dump dmAll all custom 10000 traj_equib3_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 450.0 450.0 1000.0 tchain 1
|
||||
|
||||
run 1000000
|
||||
unfix fxnvt
|
||||
undump dmAll
|
||||
|
||||
|
||||
print "---------------------------------------------------------------------------"
|
||||
print " Note: Turning ON 1-3 interactions again."
|
||||
print "---------------------------------------------------------------------------"
|
||||
special_bonds lj/coul 0 1 1
|
||||
|
||||
dump dmAll all custom 100 traj_equib4_min.lammpstrj id mol type x y z ix iy iz
|
||||
minimize 1.0e-4 1.0e-6 100000 400000
|
||||
undump dmAll
|
||||
|
||||
|
||||
# Hopefully a bilayer has formed at this point.
|
||||
# (If not, run the equilibration simulation for longer.)
|
||||
# Now I lower the temperature back to 300K.
|
||||
# I should probably re-equilibrate the solvent pressure and surface tension
|
||||
# (Simulation under NPT conditions using "anisotropic" boundaries.)
|
||||
# (so that the surface tension in the plane is allowed to relax independently
|
||||
# of the water volume, perpendicular to the plane.) I do that now:
|
||||
# temperature: 300 K, pressure: 1 barr
|
||||
fix fxnpt all npt temp 300.0 300.0 30000.0 couple xy aniso 1.0 1.0 1000.0 drag 1.0
|
||||
dump dmAll all custom 10000 traj_equib5_npt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
timestep 30.0
|
||||
run 100000
|
||||
unfix fxnpt
|
||||
undump dmAll
|
||||
|
||||
|
||||
write_data system_after_npt.data
|
||||
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
# 2) You must minimize the coordinates using by running lammps witn
|
||||
# run.in.min
|
||||
# 3) 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 simulation
|
||||
|
||||
#read_data "system.data"
|
||||
#read_data "system_after_min.data"
|
||||
read_data "system_after_npt.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
velocity all zero linear # <- eliminate drift due to non-zero total momentum
|
||||
#fix 1 all momentum 1000 linear 1 1 1 # also works
|
||||
timestep 30.0 # (40.0 should be possible for lipid systems)
|
||||
thermo 100
|
||||
#thermo_modify flush yes
|
||||
|
||||
|
||||
|
||||
# Continue the simulation at constant volume (NVT) at 300K.
|
||||
dump dmAll all custom 10000 traj_nvt_300K.lammpstrj id mol type x y z ix iy iz
|
||||
fix fxnvt all nvt temp 300.0 300.0 3000.0 tchain 1
|
||||
|
||||
run 10000000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -1,13 +0,0 @@
|
|||
|
||||
This example of the formation of a coarse-grained DPPC lipid-bilayer uses the
|
||||
Martini force-field v2.0 (2013-10), was provided by Saeed Momeni Bashusqeh.
|
||||
It's probably a good idea to run the simulation for a few ns to allow the
|
||||
lipids to reorient themselves.
|
||||
|
||||
step 1)
|
||||
To build the files which LAMMPS needs, follow the instructions in:
|
||||
README_setup.sh
|
||||
|
||||
step 2)
|
||||
To run LAMMPS with these files, follow these instructions:
|
||||
README_run.sh
|
||||
|
|
@ -1,21 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -------- PREREQUISITES: --------
|
||||
# The 2 files "run.in.min", "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.npt # 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)
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
# Use these commands to generate the LAMMPS input script and data file
|
||||
# (and other auxilliary files):
|
||||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
# These files are the input files directly read by LAMMPS. Move them to
|
||||
# the parent directory (or wherever you plan to run the simulation).
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# 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 ../
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
------- 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.0 -0.5}
|
||||
pbc box -shiftcenterrel {0.0 0.0 -0.5} -style tubes -width 0.75
|
||||
|
||||
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.
|
Before Width: | Height: | Size: 13 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 40 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 56 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 1.6 KiB |
|
|
@ -1,78 +0,0 @@
|
|||
DPPC {
|
||||
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
bond_style hybrid harmonic
|
||||
angle_style hybrid cosine/squared
|
||||
dihedral_style none
|
||||
improper_style none
|
||||
pair_style hybrid lj/gromacs/coul/gromacs 9 12 0.000001 12
|
||||
special_bonds lj/coul 0 1 1
|
||||
dielectric 15
|
||||
neigh_modify every 10
|
||||
}
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:1 $mol:. @atom:Q0 1.0 2.67583 4.37417 19.25
|
||||
$atom:2 $mol:. @atom:Qa -1.0 -0.73417 1.85417 20.00
|
||||
$atom:3 $mol:. @atom:Na 0.0 -0.91417 0.15417 16.72
|
||||
$atom:4 $mol:. @atom:Na 0.0 0.23583 -3.23583 16.96
|
||||
$atom:5 $mol:. @atom:C1 0.0 -1.26417 2.19417 12.94
|
||||
$atom:6 $mol:. @atom:C1 0.0 1.49583 1.71417 9.46
|
||||
$atom:7 $mol:. @atom:C1 0.0 2.82583 2.79417 5.04
|
||||
$atom:8 $mol:. @atom:C1 0.0 5.77583 6.29417 3.62
|
||||
$atom:9 $mol:. @atom:C1 0.0 -0.89417 -3.84583 12.6
|
||||
$atom:10 $mol:. @atom:C1 0.0 0.11583 -3.19583 7.75
|
||||
$atom:11 $mol:. @atom:C1 0.0 -2.90417 -3.64583 3.99
|
||||
$atom:12 $mol:. @atom:C1 0.0 -6.41417 -5.45583 1.81
|
||||
}
|
||||
write("Data Bonds") {
|
||||
$bond:b1 @bond:Bo $atom:1 $atom:2
|
||||
$bond:b2 @bond:Bo $atom:2 $atom:3
|
||||
$bond:b3 @bond:Short $atom:3 $atom:4
|
||||
$bond:b4 @bond:Bo $atom:3 $atom:5
|
||||
$bond:b5 @bond:Bo $atom:5 $atom:6
|
||||
$bond:b6 @bond:Bo $atom:6 $atom:7
|
||||
$bond:b7 @bond:Bo $atom:7 $atom:8
|
||||
$bond:b8 @bond:Bo $atom:4 $atom:9
|
||||
$bond:b9 @bond:Bo $atom:9 $atom:10
|
||||
$bond:b10 @bond:Bo $atom:10 $atom:11
|
||||
$bond:b11 @bond:Bo $atom:11 $atom:12
|
||||
}
|
||||
write("Data Angles") {
|
||||
$angle:a1 @angle:An1 $atom:1 $atom:2 $atom:3
|
||||
$angle:a2 @angle:An2 $atom:2 $atom:3 $atom:5
|
||||
$angle:a3 @angle:An2 $atom:2 $atom:3 $atom:4
|
||||
$angle:a4 @angle:An2 $atom:4 $atom:3 $atom:5
|
||||
$angle:a5 @angle:An1 $atom:3 $atom:4 $atom:9
|
||||
$angle:a6 @angle:An1 $atom:3 $atom:5 $atom:6
|
||||
$angle:a7 @angle:An1 $atom:5 $atom:6 $atom:7
|
||||
$angle:a8 @angle:An1 $atom:6 $atom:7 $atom:8
|
||||
$angle:a9 @angle:An1 $atom:4 $atom:9 $atom:10
|
||||
$angle:a10 @angle:An1 $atom:9 $atom:10 $atom:11
|
||||
$angle:a11 @angle:An1 $atom:10 $atom:11 $atom:12
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:Q0 72.0
|
||||
@atom:Qa 72.0
|
||||
@atom:Na 72.0
|
||||
@atom:C1 72.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:Q0 @atom:Q0 lj/gromacs/coul/gromacs 0.8365200764818 4.7
|
||||
pair_coeff @atom:Q0 @atom:Qa lj/gromacs/coul/gromacs 1.0755258126195 4.7
|
||||
pair_coeff @atom:Q0 @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Q0 @atom:C1 lj/gromacs/coul/gromacs 0.4780114722753 6.2
|
||||
pair_coeff @atom:Qa @atom:Qa lj/gromacs/coul/gromacs 1.1950286806883 4.7
|
||||
pair_coeff @atom:Qa @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Qa @atom:C1 lj/gromacs/coul/gromacs 0.4780114722753 6.2
|
||||
pair_coeff @atom:Na @atom:Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:Na @atom:C1 lj/gromacs/coul/gromacs 0.6453154875717 4.7
|
||||
pair_coeff @atom:C1 @atom:C1 lj/gromacs/coul/gromacs 0.8365200764818 4.7
|
||||
bond_coeff @bond:Bo harmonic 1.4937858508604 4.7
|
||||
bond_coeff @bond:Short harmonic 1.4937858508604 3.7
|
||||
angle_coeff @angle:An1 cosine/squared 2.9875717017208 180
|
||||
angle_coeff @angle:An2 cosine/squared 2.9875717017208 120
|
||||
}
|
||||
} #DPPC
|
||||
|
|
@ -1,27 +0,0 @@
|
|||
import "water.lt"
|
||||
import "lipid.lt"
|
||||
|
||||
write_once("Data Boundary") {
|
||||
0.0 100.0 xlo xhi
|
||||
0.0 100.0 ylo yhi
|
||||
-50.0 50.0 zlo zhi
|
||||
}
|
||||
|
||||
lipids = new DPPC [13].move(7.6923, 0, 0)
|
||||
[13].move(0, 7.6923, 0)
|
||||
[2].rot(180, 1, 0, 0)
|
||||
|
||||
waters = new MW [25].move(4.0, 0, 0)
|
||||
[25].move(0, 4.0, 0)
|
||||
[13].move(0, 0, 4.23)
|
||||
|
||||
# Move the waters upwards so that they don't overlap with the lipids.
|
||||
|
||||
waters[*][*][*].move(0, 0, 22.4)
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Q0 lj/gromacs/coul/gromacs 1.3384321223709 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Qa lj/gromacs/coul/gromacs 1.3384321223709 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/Na lj/gromacs/coul/gromacs 0.9560229445507 4.7
|
||||
pair_coeff @atom:MW/P4 @atom:DPPC/C1 lj/gromacs/coul/gromacs 0.4780114722753 4.7
|
||||
}
|
||||
|
|
@ -1,17 +0,0 @@
|
|||
MW {
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
pair_style hybrid lj/gromacs/coul/gromacs 9 12 0.000001 12
|
||||
}
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:1 $mol:. @atom:P4 0 0 0 0
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:P4 72.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:P4 @atom:P4 lj/gromacs/coul/gromacs 1.1950286806883 4.7
|
||||
}
|
||||
} #MW
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data "system.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
thermo 5
|
||||
dump 1 all custom 100 traj_min.lammpstrj id mol type x y z ix iy iz
|
||||
minimize 1.0e-4 1.0e-6 100000 400000
|
||||
|
||||
write_data system_after_min.data
|
||||
|
|
@ -1,66 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
# 2) You must minimize the coordinates using by running lammps witn
|
||||
# run.in.min
|
||||
#
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include "system.in.init"
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
#read_data "system.data"
|
||||
read_data "system_after_min.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
print "---------------------------------------------------------------------------"
|
||||
print "I often use Langevin dynamics initially at high temperatures and small"
|
||||
print "timesteps to relax the system. It seems more stable than Nose-Hoover."
|
||||
print "(This is probably not necessary.)"
|
||||
print "---------------------------------------------------------------------------"
|
||||
|
||||
fix fxlan all langevin 450.0 450.0 100 12345 # temp: 450 K
|
||||
fix fxnph all nph aniso 100.0 100.0 1000.0 couple xy drag 1.0 #pressure:100barr
|
||||
thermo 100
|
||||
dump dmNPTall all custom 5000 traj_npt_step1.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
timestep 1.0 # (safer to use a small timestep initially)
|
||||
run 1000
|
||||
timestep 3.0
|
||||
run 1000
|
||||
timestep 10.0
|
||||
run 1000
|
||||
timestep 30.0 # (40.0 should be possible for lipid systems)
|
||||
run 100000
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnph
|
||||
undump dmNPTall
|
||||
|
||||
print "---------------------------------------------------------------------------"
|
||||
print "--- Now continue the simulation using a Nose-Hoover Thermostat/Barostat ---"
|
||||
print "---------------------------------------------------------------------------"
|
||||
|
||||
velocity all zero linear # <- eliminate drift due to non-zero total momentum
|
||||
#fix 1 all momentum 1000 linear 1 1 1 # also works
|
||||
|
||||
# temperature: 300 K, pressure: 1 barr
|
||||
fix fxnpt all npt temp 300.0 300.0 100.0 aniso 1.0 1.0 1000.0 drag 1.0 couple xy
|
||||
thermo 100
|
||||
#thermo_modify flush yes
|
||||
dump dmNPTall all custom 10000 traj_npt_step2.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
run 100000
|
||||
|
||||
write_data system_after_npt.data
|
||||
|
|
@ -1,47 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# 1) You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details.)
|
||||
# 2) You must minimize the coordinates using by running lammps witn
|
||||
# run.in.min
|
||||
# 3) 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 simulation
|
||||
|
||||
#read_data "system.data"
|
||||
#read_data "system_after_min.data"
|
||||
read_data "system_after_npt.data"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include "system.in.settings"
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
velocity all zero linear # <- eliminate drift due to non-zero total momentum
|
||||
#fix 1 all momentum 1000 linear 1 1 1 # also works
|
||||
|
||||
timestep 30.0 # (40.0 should be possible for lipid systems)
|
||||
dump 1 all custom 20000 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 100
|
||||
#thermo_modify flush yes
|
||||
|
||||
run 10000000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
|
@ -1,17 +0,0 @@
|
|||
|
||||
This directory contains an example of a couarse-grained (vaguely protein-like)
|
||||
heteropolymer consisting of 14 residues, each of which has 2 atoms
|
||||
(one backbone atom, one residue atom.)
|
||||
|
||||
There are two types of residues, H and P.
|
||||
The R-atom for the H residue are attracted to eachother.
|
||||
All other atoms are repulsive.
|
||||
|
||||
Instructions on how to build LAMMPS input files and
|
||||
run a short simulation are provided in other README files.
|
||||
|
||||
step 1)
|
||||
README_setup.sh
|
||||
|
||||
step2)
|
||||
README_run.sh
|
||||
|
|
@ -1,20 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -- Prerequisites: --
|
||||
# The "run.in.nvt" file is a LAMMPS input script containing
|
||||
# references 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.nvt
|
||||
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
# Use these commands to generate the LAMMPS input script and data file
|
||||
# (and other auxilliary files):
|
||||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
# These files are the input files directly read by LAMMPS. Move them to
|
||||
# the parent directory (or wherever you plan to run the simulation).
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# 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 ../
|
||||
|
|
@ -1,86 +0,0 @@
|
|||
|
||||
------- 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 ----
|
||||
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.
|
Before Width: | Height: | Size: 52 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 99 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 29 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 37 KiB |
|
|
@ -1,6 +0,0 @@
|
|||
# Use these commands to generate the LAMMPS input script and data file
|
||||
# (and other auxilliary files):
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
|
||||
|
|
@ -1,139 +0,0 @@
|
|||
# The "2beadFF" is a force-field environment object containing
|
||||
# force-field data, atomic masses, and bond rules.
|
||||
# Later, when we define molecules (such as "H" and "P"), we can inherit
|
||||
# atom types and bond-rules from this force-field. This will automatically
|
||||
# assign bonds and angular interactions according to atom (and bond) type.
|
||||
# (You can also assign charge by atom type. However in this example I assigned
|
||||
# charge to each atom manually (not by type). The OPLSAA examples in the
|
||||
# "all_atoms" directory demonstrate how to assign charge by atom type.)
|
||||
|
||||
|
||||
|
||||
2beadFF {
|
||||
|
||||
# There are 3 atom types:
|
||||
|
||||
write_once("Data Masses") {
|
||||
@atom:CA 13.0
|
||||
@atom:HR 50.0
|
||||
@atom:PR 50.0
|
||||
}
|
||||
|
||||
# 2-body (non-bonded) interactions:
|
||||
#
|
||||
# Uij(r) = 4*eps_ij * ( (sig_ij/r)^12 - (sig_ij/r)^6 )
|
||||
#
|
||||
# Hydrophobic side-chain (R) atoms are attractive (large epsilon parameter).
|
||||
# Polar side-chains and backbone atoms are not attractive (small epsilon).
|
||||
#
|
||||
# i j pairstylename eps sig
|
||||
#
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:CA @atom:CA lj/cut 0.10 2.0
|
||||
pair_coeff @atom:HR @atom:HR lj/cut 2.50 3.6
|
||||
pair_coeff @atom:PR @atom:PR lj/cut 0.10 3.6
|
||||
}
|
||||
|
||||
# (By default, interactions between different AtomTypes use "arithmetic"rules:
|
||||
# eps_ij=sqrt(eps_ii*eps_ij) and sig_ij=0.5*(sig_ii+sig_jj)
|
||||
# Look for the line containing "pair_modify mix arithmetic" below...)
|
||||
|
||||
# Optional: Assign bond types @bond:Backbone or @bond:Sidechain
|
||||
# according to atom type. (This can be overridden.)
|
||||
|
||||
write_once("Data Bonds By Type") {
|
||||
@bond:Backbone @atom:CA @atom:CA
|
||||
@bond:Sidechain @atom:CA @atom:HR
|
||||
@bond:Sidechain @atom:CA @atom:PR
|
||||
}
|
||||
|
||||
# 2-body (bonded) interactions:
|
||||
#
|
||||
# Ubond(r) = (k/2)*(r-0)^2
|
||||
#
|
||||
# The corresponding command is:
|
||||
#
|
||||
# bond_coeff bondType bondstylename k r0
|
||||
#
|
||||
|
||||
write_once("In Settings") {
|
||||
bond_coeff @bond:Sidechain harmonic 30.0 3.4
|
||||
bond_coeff @bond:Backbone harmonic 30.0 3.7
|
||||
}
|
||||
|
||||
|
||||
# Rules for determining 3 and 4-body bonded interactions by type
|
||||
|
||||
# angle-type atomType1 atomType2 atomType3
|
||||
|
||||
write_once("Data Angles By Type") {
|
||||
@angle:Backbone @atom:CA @atom:CA @atom:CA
|
||||
@angle:Sidechain @atom:CA @atom:CA @atom:*R # Note: "*R" <--> "HR" or "PR"
|
||||
}
|
||||
|
||||
# dihedral-type AtomType1 AtomType2 AtomType3 AtomType4
|
||||
|
||||
write_once("Data Dihedrals By Type") {
|
||||
@dihedral:CCCC @atom:CA @atom:CA @atom:CA @atom:CA
|
||||
@dihedral:RCCR @atom:*R @atom:CA @atom:CA @atom:*R #"*R" <--> "HR" or "PR"
|
||||
}
|
||||
|
||||
# 3-body interactions in this example are listed by atomType and bondType
|
||||
# The atomIDs involved are determined automatically. The forumula used is:
|
||||
#
|
||||
# Uangle(theta) = (k/2)*(theta-theta0)^2
|
||||
# (k in kcal/mol/rad^2, theta0 in degrees)
|
||||
#
|
||||
# The corresponding command is:
|
||||
#
|
||||
# angle_coeff angleType anglestylename k theta0
|
||||
|
||||
write_once("In Settings") {
|
||||
angle_coeff @angle:Backbone harmonic 30.00 114
|
||||
angle_coeff @angle:Sidechain harmonic 30.00 123
|
||||
}
|
||||
|
||||
|
||||
# 4-body interactions in this example are listed by atomType and bondType
|
||||
# The atomIDs involved are determined automatically. The forumula used is:
|
||||
#
|
||||
# Udihedral(phi) = K * (1 + cos(n*phi - d))
|
||||
#
|
||||
# The d parameter is in degrees, K is in kcal/mol/rad^2.
|
||||
#
|
||||
# The corresponding command is
|
||||
# dihedral_coeff dihedralType dihedralstylename K n d w (ignore "w")
|
||||
|
||||
write_once("In Settings") {
|
||||
dihedral_coeff @dihedral:CCCC charmm -0.5 1 -180 0.0
|
||||
dihedral_coeff @dihedral:RCCR charmm -1.5 1 -180 0.0
|
||||
} # write_once("In Settings")
|
||||
|
||||
|
||||
# LAMMPS supports a large number of force-field styles. We must select
|
||||
# which ones we need. This information belongs in the "In Init" section.
|
||||
# (Hybrid styles used for portability. These choices can be overridden later.)
|
||||
|
||||
write_once("In Init") {
|
||||
# -- Default styles for "2bead" --
|
||||
# (Hybrid force fields were not necessary but are used for portability.)
|
||||
units real
|
||||
atom_style full
|
||||
bond_style hybrid harmonic
|
||||
angle_style hybrid harmonic
|
||||
dihedral_style hybrid charmm
|
||||
pair_style hybrid lj/cut 11.0
|
||||
|
||||
# If charges are needed, (assuming biopolymers), try one of:
|
||||
#dielectric 80.0
|
||||
#pair_style hybrid lj/cut/coul/debye 0.1 11.0
|
||||
# or (for short distances, below a couple nm)
|
||||
#pair_style hybrid lj/charmm/coul/charmm/implicit 9.0 11.0
|
||||
|
||||
pair_modify mix arithmetic
|
||||
special_bonds lj 0.0 0.0 0.0
|
||||
}
|
||||
|
||||
|
||||
} # 2beadFF
|
||||
|
||||
|
|
@ -1,86 +0,0 @@
|
|||
# In this example, we define two types of molecules: "H" and "P",
|
||||
# both containing two atoms, whose ids (names) are "ca" and "r",
|
||||
# and whose atom-types vary.
|
||||
#
|
||||
# "H" molecules: "P" molecules:
|
||||
#
|
||||
# @HR @PR
|
||||
# | |
|
||||
# @CA @CA
|
||||
#
|
||||
# Eventually, we will connect multiple "H" and "P" molecules
|
||||
# together to form a polymer, as shown below:
|
||||
#
|
||||
# @HR @HR
|
||||
# | |
|
||||
# _@CA_ _@CA_
|
||||
# ... -.@CA-' `-@CA-' ` ...
|
||||
# | |
|
||||
# @PR @PR
|
||||
#
|
||||
# The "H" and "P" molecules both share the same type of
|
||||
# backbone atom ("CA"), but have their own custom "r"
|
||||
# sidechain atoms with different properties:
|
||||
# The "HR" atoms belonging to "H" molecules are attracted to each other.
|
||||
# The "PR" atoms in "P" molecules are not.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
import "forcefield.lt" # defines "2beadFF"
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Define the "H" monomer type ("H" <--> "hydrophobic")
|
||||
|
||||
H inherits 2beadFF {
|
||||
# atom-id(name) mol-id atom-type charge x y z
|
||||
write("Data Atoms") {
|
||||
$atom:ca $mol:... @atom:CA 0.0 0.000 1.0000 0.0000000
|
||||
$atom:r $mol:... @atom:HR 0.0 0.000 4.4000 0.0000000
|
||||
}
|
||||
|
||||
write("Data Bond List") {
|
||||
$bond:cr $atom:ca $atom:r
|
||||
}
|
||||
|
||||
# This will look up the bond-parameters according to atom type.
|
||||
# Use "Data Bonds" instead if you prefer to assign the bond type manually:
|
||||
# write("Data Bonds") {
|
||||
# $bond:cr @bond:Sidechain $atom:ca $atom:r
|
||||
# }
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Define the "P" monomer type ("P" <--> "polar")
|
||||
|
||||
P inherits 2beadFF {
|
||||
|
||||
# atom-id(name) mol-id atom-type charge x y z
|
||||
write("Data Atoms") {
|
||||
$atom:ca $mol:... @atom:CA 0.0 0.000 1.0000 0.0000000
|
||||
$atom:r $mol:... @atom:PR 0.0 0.000 4.4000 0.0000000
|
||||
}
|
||||
|
||||
write("Data Bond List") {
|
||||
$bond:CR $atom:ca $atom:r
|
||||
}
|
||||
|
||||
# This will look up the bond-parameters according to atom type.
|
||||
# Use "Data Bonds" instead if you prefer to assign the bond type manually:
|
||||
# write("Data Bonds") {
|
||||
# $bond:cr @bond:Sidechain $atom:ca $atom:r
|
||||
# }
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Note: The "..." in "$mol:..." tells moltemplate that this molecule may
|
||||
# be a part of a larger molecule, and (if so) to use the larger
|
||||
# molecule's id number as it's own.
|
||||
|
|
@ -1,64 +0,0 @@
|
|||
|
||||
import "monomers.lt" # This defines the monomer types named "H" and "P"
|
||||
|
||||
|
||||
|
||||
Polymer {
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
# (The "Data Atoms" in H and P must use "$mol:..." notation.)
|
||||
|
||||
# This causes mon1,mon2,mon3,...,mon14 to share the same molecule counter
|
||||
# because in the 2bead.lt file, the "..." in "$mol:..." preferentially looks
|
||||
# for a counter of that type in a parent molecule or earlier ancestor.
|
||||
|
||||
|
||||
# A polymer of alternating "H" and "P" monomers:
|
||||
|
||||
mon1 = new P
|
||||
mon2 = new P.rot(180.0, 1,0,0).move(3.2,0,0)
|
||||
mon3 = new H.rot( 0.0, 1,0,0).move(6.4,0,0)
|
||||
mon4 = new H.rot(180.0, 1,0,0).move(9.6,0,0)
|
||||
mon5 = new H.rot( 0.0, 1,0,0).move(12.8,0,0)
|
||||
mon6 = new H.rot(180.0, 1,0,0).move(16.0,0,0)
|
||||
mon7 = new P.rot( 0.0, 1,0,0).move(19.2,0,0)
|
||||
mon8 = new P.rot(180.0, 1,0,0).move(22.4,0,0)
|
||||
mon9 = new P.rot( 0.0, 1,0,0).move(25.6,0,0)
|
||||
mon10 = new H.rot(180.0, 1,0,0).move(28.8,0,0)
|
||||
mon11 = new H.rot( 0.0, 1,0,0).move(32.0,0,0)
|
||||
mon12 = new H.rot(180.0, 1,0,0).move(35.2,0,0)
|
||||
mon13 = new P.rot( 0.0, 1,0,0).move(38.4,0,0)
|
||||
mon14 = new P.rot(180.0, 1,0,0).move(41.6,0,0)
|
||||
|
||||
|
||||
# Now, link the monomers together this way:
|
||||
write("Data Bond List") {
|
||||
$bond:backbone1 $atom:mon1/ca $atom:mon2/ca
|
||||
$bond:backbone2 $atom:mon2/ca $atom:mon3/ca
|
||||
$bond:backbone3 $atom:mon3/ca $atom:mon4/ca
|
||||
$bond:backbone4 $atom:mon4/ca $atom:mon5/ca
|
||||
$bond:backbone5 $atom:mon5/ca $atom:mon6/ca
|
||||
$bond:backbone6 $atom:mon6/ca $atom:mon7/ca
|
||||
$bond:backbone7 $atom:mon7/ca $atom:mon8/ca
|
||||
$bond:backbone8 $atom:mon8/ca $atom:mon9/ca
|
||||
$bond:backbone9 $atom:mon9/ca $atom:mon10/ca
|
||||
$bond:backbone10 $atom:mon10/ca $atom:mon11/ca
|
||||
$bond:backbone11 $atom:mon11/ca $atom:mon12/ca
|
||||
$bond:backbone12 $atom:mon12/ca $atom:mon13/ca
|
||||
$bond:backbone13 $atom:mon13/ca $atom:mon14/ca
|
||||
}
|
||||
|
||||
|
||||
# Use "Data Bonds" instead if you prefer to assign the bond types manually:
|
||||
# write("Data Bonds") {
|
||||
# $bond:backbone1 @bond:Backbone $atom:mon1/ca $atom:mon2/ca
|
||||
# $bond:backbone2 @bond:Backbone $atom:mon2/ca $atom:mon3/ca
|
||||
# : : : :
|
||||
# }
|
||||
|
||||
} # Polymer
|
||||
|
||||
|
||||
|
||||
# Angle, dihedral and improper interactions between monomers will be generated
|
||||
# automatically according to the instructions in the "force_field.lt" file.
|
||||
|
|
@ -1,36 +0,0 @@
|
|||
import "polymer.lt"
|
||||
|
||||
|
||||
# Specify the periodic boundary conditions:
|
||||
write_once("Data Boundary") {
|
||||
0 180.0 xlo xhi
|
||||
0 180.0 ylo yhi
|
||||
0 180.0 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
# Create 27 polymers (=3x3x3) in a rectangular grid
|
||||
|
||||
polymers = new Polymer [3].move(60.0, 0, 0)
|
||||
[3].move(0, 60.0, 0)
|
||||
[3].move(0, 0, 60.0)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# ----- everything below is optional: -----
|
||||
# Shift some of the polymers in the Z direction by a distance of 20.0
|
||||
|
||||
polymers[1][*][*].move(0,0,20)
|
||||
|
||||
# We applied this move command to all the
|
||||
# polymers in the middle slab (with constant X).
|
||||
# More examples of applying the "move" command:
|
||||
|
||||
polymers[*][1][*].move(20,0,0)
|
||||
polymers[*][*][1].move(0,20,0)
|
||||
|
||||
|
|
@ -1,32 +0,0 @@
|
|||
# -- Init Section --
|
||||
|
||||
include system.in.init
|
||||
|
||||
# -- Atom Definition Section --
|
||||
|
||||
read_data system.data
|
||||
|
||||
# -- Settings Section --
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# -- Run Section --
|
||||
|
||||
|
||||
timestep 2.0
|
||||
dump 1 all custom 2500 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
# To use Langevin dynamics in LAMMPS you need both "fix langevin" and "fix nve".
|
||||
# (See http://lammps.sandia.gov/doc/fix_langevin.html for details.)
|
||||
|
||||
fix fxlan all langevin 300.0 300.0 5000.0 48279
|
||||
fix fxnve all nve
|
||||
|
||||
|
||||
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
|
||||
thermo 2500 # time interval for printing out "thermo" data
|
||||
|
||||
run 20000000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
||||
|
|
@ -1,19 +0,0 @@
|
|||
This is an example of how to build a polymer out of randomly-chosen monomers.
|
||||
In this case, monomers will be chosen at random from two types
|
||||
(denoted "2bead" and "3bead", although you can have as many types as you like).
|
||||
You can also constrain the end-caps to be a particular type (eg "3bead").
|
||||
|
||||
The properties of the bonds connecting monomers (ie length, rigidity) will
|
||||
be automatically determined, depending on the type of monomers at that location
|
||||
in the polymer. The same is true for the 3-body angle, and 4-body dihedral
|
||||
interactions.
|
||||
|
||||
|
||||
Instructions on how to build LAMMPS input files and
|
||||
run a short simulation are provided in other README files.
|
||||
|
||||
step 1)
|
||||
README_setup.sh
|
||||
|
||||
step2)
|
||||
README_run.sh
|
||||
|
|
@ -1,14 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -- Prerequisites: --
|
||||
# The "run.in.nvt" file is a LAMMPS input script containing
|
||||
# references 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.nvt
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
# Use these commands to generate the LAMMPS input script and data file
|
||||
# (and other auxilliary files):
|
||||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
# These files are the input files directly read by LAMMPS. Move them to
|
||||
# the parent directory (or wherever you plan to run the simulation).
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# 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 ../
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
------- 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.
|
Before Width: | Height: | Size: 2.6 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 4.3 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 51 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 45 KiB |
|
|
@ -1,107 +0,0 @@
|
|||
# The "ForceField" object contains a list of atom types, masses, force-field
|
||||
# parameters, and rules for generating 3-body & 4-body bonded interactions.
|
||||
# Molecules which use "inherit ForceField" share these rules, and consequently
|
||||
# can usually be written in a much more concise way.
|
||||
|
||||
|
||||
ForceField {
|
||||
|
||||
# LAMMPS supports a large number of force-field styles. We must select
|
||||
# which ones we need. This information belongs in the "In Init" section (and
|
||||
# (you can specify it anywhere in your LT files, multiple times if you like).
|
||||
# If different molecules use different force-field styles, you can use hybrid
|
||||
# styles. (In this example the molecules share the same pair_style.)
|
||||
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
bond_style harmonic
|
||||
angle_style harmonic
|
||||
dihedral_style opls
|
||||
pair_style lj/cut 9.0
|
||||
# If you have charged molecules immersed in a salty implicit
|
||||
# solvent, you might try something like this this instead:
|
||||
# pair_style lj/cut/coul/debye 0.1 9.0
|
||||
pair_modify mix arithmetic
|
||||
dielectric 80.0
|
||||
special_bonds lj 0.0 0.0 0.0
|
||||
}
|
||||
|
||||
# atom-type mass
|
||||
|
||||
write_once("Data Masses") {
|
||||
@atom:CA2 12.0
|
||||
@atom:R2 17.0
|
||||
@atom:CA3 12.0
|
||||
@atom:R3 17.0
|
||||
}
|
||||
|
||||
|
||||
# Connect different monomers together with bonds whose type
|
||||
# (length, rigidity, etc) depend on the type of atom at either end.
|
||||
|
||||
write_once("Data Bonds By Type") {
|
||||
@bond:sidechain @atom:CA* @atom:R*
|
||||
@bond:two_two @atom:CA2 @atom:CA2
|
||||
@bond:two_three @atom:CA2 @atom:CA3
|
||||
@bond:three_three @atom:CA3 @atom:CA3
|
||||
}
|
||||
# Note: The next line is redundant and unnecessary:
|
||||
# @bond:two_three @atom:CA3 @atom:CA2
|
||||
|
||||
|
||||
|
||||
# You can also generate angles and dihedrals and impropers in a similar way:
|
||||
|
||||
# Rules for determining 3 and 4-body bonded interactions by type
|
||||
|
||||
# angle-type atomType1 atomType2 atomType3 bondType1 bondType2
|
||||
|
||||
write_once("Data Angles By Type") {
|
||||
@angle:backbone @atom:CA* @atom:CA* @atom:CA* @bond:* @bond:*
|
||||
@angle:sidechain @atom:CA* @atom:CA* @atom:R* @bond:* @bond:*
|
||||
@angle:RCR @atom:R* @atom:CA* @atom:R* @bond:* @bond:*
|
||||
}
|
||||
# Note: The next line is redundant and unnecessary:
|
||||
# @angle:sidechain @atom:R* @atom:CA* @atom:CA* @bond:* @bond:*
|
||||
|
||||
|
||||
# dihedral-type AtomType1 AtomType2 AtomType3 AtomType4 bondType1 btyp2 btyp3
|
||||
|
||||
write_once("Data Dihedrals By Type") {
|
||||
@dihedral:backbone @atom:CA* @atom:CA* @atom:CA* @atom:CA* * * *
|
||||
@dihedral:two_two @atom:R2 @atom:CA* @atom:CA* @atom:R2 * * *
|
||||
@dihedral:two_three @atom:R2 @atom:CA* @atom:CA* @atom:R3 * * *
|
||||
@dihedral:three_three @atom:R3 @atom:CA* @atom:CA* @atom:R3 * * *
|
||||
}
|
||||
# Note: The next line is redundant and unnecessary:
|
||||
# @dihedral:two_three @atom:R3 @atom:CA* @atom:CA* @atom:R2 * * *
|
||||
|
||||
# Force field parameters:
|
||||
write_once("In Settings") {
|
||||
# atom-type atom-type epsilon sigma
|
||||
pair_coeff @atom:CA2 @atom:CA2 0.05 3.0
|
||||
pair_coeff @atom:R2 @atom:R2 0.60 4.0
|
||||
pair_coeff @atom:CA3 @atom:CA3 0.10 2.0
|
||||
pair_coeff @atom:R3 @atom:R3 0.50 3.0
|
||||
|
||||
# bond-type k r0
|
||||
bond_coeff @bond:sidechain 20.0 3.4
|
||||
bond_coeff @bond:two_two 20.0 3.7
|
||||
bond_coeff @bond:two_three 20.0 3.5
|
||||
bond_coeff @bond:three_three 20.0 3.3
|
||||
|
||||
# angle-type k theta0
|
||||
angle_coeff @angle:backbone 40.00 120
|
||||
angle_coeff @angle:sidechain 40.00 120
|
||||
angle_coeff @angle:RCR 40.00 120
|
||||
|
||||
# dihedral-type K1 K2 K3 K4
|
||||
dihedral_coeff @dihedral:backbone 0.3 0.0 0.0 0.0
|
||||
dihedral_coeff @dihedral:two_two 0.08 0.0 0.0 0.0
|
||||
dihedral_coeff @dihedral:two_three 0.08 0.0 -0.05 0.0
|
||||
dihedral_coeff @dihedral:three_three 0.08 0.0 0.0 0.05
|
||||
}
|
||||
|
||||
} # ForceField
|
||||
|
||||
|
|
@ -1,62 +0,0 @@
|
|||
import "forcefield.lt" #<-- "forcefield.lt" contains atom type definitions,
|
||||
# force-field parameters, and rules for generating
|
||||
# 3-body and 4-body angle, dihedral, and improper
|
||||
# interactions for molecules (polymers) made using
|
||||
# "2bead" and "3bead" objects as building blocks.
|
||||
|
||||
|
||||
|
||||
# ----------------------------------------------------------------------
|
||||
# -- General comment: --
|
||||
# -- The write() and write_once() commands create and append text to --
|
||||
# -- files (replacing variables beginning with @ or $ with counters.) --
|
||||
# -- File names beginning with "In " or "Data " are special. --
|
||||
# -- They will be pasted into the LAMMPS input script and --
|
||||
# -- data files which are generated by moltemplate. The syntax --
|
||||
# -- of these files is exactly the same as the syntax from the --
|
||||
# -- corresponding sections of a LAMMPS input script or data file. --
|
||||
# ----------------------------------------------------------------------
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
3bead inherits ForceField {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:CA $mol:... @atom:CA3 0.0 0.000 1.000 0.000
|
||||
$atom:R1 $mol:... @atom:R3 0.0 0.000 2.700 2.950
|
||||
$atom:R2 $mol:... @atom:R3 0.0 0.000 2.700 -2.950
|
||||
}
|
||||
|
||||
# bond-id atom-id1 atom-id2
|
||||
write("Data Bond List") {
|
||||
$bond:CR1 $atom:CA $atom:R1
|
||||
$bond:CR2 $atom:CA $atom:R2
|
||||
}
|
||||
|
||||
} # 3bead
|
||||
|
||||
|
||||
|
||||
|
||||
2bead inherits ForceField {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:CA $mol:... @atom:CA2 0.0 0.000 1.000 0.0000
|
||||
$atom:R $mol:... @atom:R2 0.0 0.000 4.400 0.0000
|
||||
}
|
||||
|
||||
# bond-id atom-id1 atom-id2
|
||||
write("Data Bond List") {
|
||||
$bond:CR $atom:CA $atom:R
|
||||
}
|
||||
|
||||
} # 2bead
|
||||
|
||||
|
||||
|
|
@ -1,163 +0,0 @@
|
|||
import "forcefield.lt" # <-- Defines force-field parameters, and rules.
|
||||
import "monomers.lt" # <-- Defines the "2bead", "3bead" objects we will use
|
||||
# as building-blocks to build the polymer.
|
||||
|
||||
|
||||
|
||||
|
||||
RandomHeteropolymer inherits ForceField {
|
||||
|
||||
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
# If you want to share molecule-IDs, do this first
|
||||
# (This matches with "$mol:..." in "monomers.lt")
|
||||
|
||||
|
||||
|
||||
monomers[0] = new 3bead # Let the first monomer be of type "3bead"
|
||||
|
||||
# Now, fill the middle of the chain with random monomers (2bead, 3bead):
|
||||
|
||||
monomers[1-48] = new random([2bead,3bead],
|
||||
[30, # <-- 30 "2bead" molecules and
|
||||
18], # <-- 18 "3bead" molecules
|
||||
12345) # <-- optional random seed
|
||||
|
||||
# and place them on the X-axis (separated by 2.95 Angstroms)
|
||||
[48].rot(180,1,0,0).move(2.95, 0, 0)
|
||||
|
||||
# Note: The two numbers (30 and 18) must sum up to the number of
|
||||
# monomers we created ("[48]").
|
||||
|
||||
monomers[49] = new 3bead # Let the last monomer be of type "3bead"
|
||||
|
||||
|
||||
# Afterwards, there should be 20 "3bead" monomers and 30 "2bead" monomers
|
||||
|
||||
|
||||
|
||||
|
||||
# Now, physically move the monomers to make sure the monomers on the end of
|
||||
# the chain (monomers[0] & monomers[49]) don't overlap with monomers[1-48]
|
||||
|
||||
monomers[0].rot(180,1,0,0) #leave monomer[0] where it is, but rotate it
|
||||
monomers[1-48].move(2.95,0,0) #move the remaining monomers to make room for it
|
||||
monomers[49].move(144.55,0,0) #move the last monomer (note:144.55=2.95*50))
|
||||
|
||||
|
||||
|
||||
|
||||
# Now, link the monomers together using "Data Bond List"
|
||||
# (Using "Data Bond List" instead of "Data Bonds", allows you to omit the
|
||||
# bond type. This tells moltemplate to look up the appropriate
|
||||
# bond type according to the type of atom at either end of the bond, which
|
||||
# depends on what type of monomer is on either side: "2bead" or "3bead".
|
||||
# This all happens automatically. The user can control how this is done
|
||||
# by editing the file "forcefield.lt".)
|
||||
|
||||
|
||||
write("Data Bond List") {
|
||||
$bond:bb1 $atom:monomers[0]/CA $atom:monomers[1]/CA
|
||||
$bond:bb2 $atom:monomers[1]/CA $atom:monomers[2]/CA
|
||||
$bond:bb3 $atom:monomers[2]/CA $atom:monomers[3]/CA
|
||||
$bond:bb4 $atom:monomers[3]/CA $atom:monomers[4]/CA
|
||||
$bond:bb5 $atom:monomers[4]/CA $atom:monomers[5]/CA
|
||||
$bond:bb6 $atom:monomers[5]/CA $atom:monomers[6]/CA
|
||||
$bond:bb7 $atom:monomers[6]/CA $atom:monomers[7]/CA
|
||||
$bond:bb8 $atom:monomers[7]/CA $atom:monomers[8]/CA
|
||||
$bond:bb9 $atom:monomers[8]/CA $atom:monomers[9]/CA
|
||||
$bond:bb10 $atom:monomers[9]/CA $atom:monomers[10]/CA
|
||||
$bond:bb11 $atom:monomers[10]/CA $atom:monomers[11]/CA
|
||||
$bond:bb12 $atom:monomers[11]/CA $atom:monomers[12]/CA
|
||||
$bond:bb13 $atom:monomers[12]/CA $atom:monomers[13]/CA
|
||||
$bond:bb14 $atom:monomers[13]/CA $atom:monomers[14]/CA
|
||||
$bond:bb15 $atom:monomers[14]/CA $atom:monomers[15]/CA
|
||||
$bond:bb16 $atom:monomers[15]/CA $atom:monomers[16]/CA
|
||||
$bond:bb17 $atom:monomers[16]/CA $atom:monomers[17]/CA
|
||||
$bond:bb18 $atom:monomers[17]/CA $atom:monomers[18]/CA
|
||||
$bond:bb19 $atom:monomers[18]/CA $atom:monomers[19]/CA
|
||||
$bond:bb20 $atom:monomers[19]/CA $atom:monomers[20]/CA
|
||||
$bond:bb21 $atom:monomers[20]/CA $atom:monomers[21]/CA
|
||||
$bond:bb22 $atom:monomers[21]/CA $atom:monomers[22]/CA
|
||||
$bond:bb23 $atom:monomers[22]/CA $atom:monomers[23]/CA
|
||||
$bond:bb24 $atom:monomers[23]/CA $atom:monomers[24]/CA
|
||||
$bond:bb25 $atom:monomers[24]/CA $atom:monomers[25]/CA
|
||||
$bond:bb26 $atom:monomers[25]/CA $atom:monomers[26]/CA
|
||||
$bond:bb27 $atom:monomers[26]/CA $atom:monomers[27]/CA
|
||||
$bond:bb28 $atom:monomers[27]/CA $atom:monomers[28]/CA
|
||||
$bond:bb29 $atom:monomers[28]/CA $atom:monomers[29]/CA
|
||||
$bond:bb30 $atom:monomers[29]/CA $atom:monomers[30]/CA
|
||||
$bond:bb31 $atom:monomers[30]/CA $atom:monomers[31]/CA
|
||||
$bond:bb32 $atom:monomers[31]/CA $atom:monomers[32]/CA
|
||||
$bond:bb33 $atom:monomers[32]/CA $atom:monomers[33]/CA
|
||||
$bond:bb34 $atom:monomers[33]/CA $atom:monomers[34]/CA
|
||||
$bond:bb35 $atom:monomers[34]/CA $atom:monomers[35]/CA
|
||||
$bond:bb36 $atom:monomers[35]/CA $atom:monomers[36]/CA
|
||||
$bond:bb37 $atom:monomers[36]/CA $atom:monomers[37]/CA
|
||||
$bond:bb38 $atom:monomers[37]/CA $atom:monomers[38]/CA
|
||||
$bond:bb39 $atom:monomers[38]/CA $atom:monomers[39]/CA
|
||||
$bond:bb40 $atom:monomers[39]/CA $atom:monomers[40]/CA
|
||||
$bond:bb41 $atom:monomers[40]/CA $atom:monomers[41]/CA
|
||||
$bond:bb42 $atom:monomers[41]/CA $atom:monomers[42]/CA
|
||||
$bond:bb43 $atom:monomers[42]/CA $atom:monomers[43]/CA
|
||||
$bond:bb44 $atom:monomers[43]/CA $atom:monomers[44]/CA
|
||||
$bond:bb45 $atom:monomers[44]/CA $atom:monomers[45]/CA
|
||||
$bond:bb46 $atom:monomers[45]/CA $atom:monomers[46]/CA
|
||||
$bond:bb47 $atom:monomers[46]/CA $atom:monomers[47]/CA
|
||||
$bond:bb48 $atom:monomers[47]/CA $atom:monomers[48]/CA
|
||||
$bond:bb49 $atom:monomers[48]/CA $atom:monomers[49]/CA
|
||||
}
|
||||
|
||||
} # RandomHeteropolymer
|
||||
|
||||
|
||||
|
||||
|
||||
# COMMENTS:
|
||||
#
|
||||
#
|
||||
# 1)
|
||||
# Angle, dihedral and improper interactions will be generated
|
||||
# automatically according to the instructions in "monomers.lt"
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# 2)
|
||||
# These lines in the "Data Bond List" section can be quickly generated in python
|
||||
# N = 50
|
||||
# for i in range(0,N-1):
|
||||
# print(' $bond:bb'+str(i+1)+' $atom:monomers['
|
||||
# +str(i)+']/CA $atom:monomers['+str(i+1)+']/CA')
|
||||
#
|
||||
#
|
||||
#
|
||||
# 3)
|
||||
# The "[1-50]" in "monomers[1-50] = new random([2bead,3bead], ..."
|
||||
# causes moltemplate to create a list of molecules with names beginning with
|
||||
# "monomers[1]", for example:
|
||||
# "monomers[1], "monomers[2]", "monomers[3], ..., "monomers[50]"
|
||||
# (This causes the indexing to begin at [1] instead of [0].)
|
||||
#
|
||||
#
|
||||
#
|
||||
# 4)
|
||||
# ALTERNATE METHOD: You can also generate a random array this way:
|
||||
#
|
||||
# monomers[1-48] = new random([2bead,3bead],
|
||||
# [0.625, 0.375], # <-- probabilities
|
||||
# 12345) # <-- optional random seed
|
||||
# [48].rot(180,1,0,0).move(2.95, 0, 0)
|
||||
#
|
||||
# The command above also works, but it chooses each molecule (monomer) randomly
|
||||
# (independently of the others). Consequently, this does not gaurantee that
|
||||
# exactly 62.5% and 37.5% of the monomers will be of type 2bead and 3bead.
|
||||
#
|
||||
#
|
||||
#
|
||||
# 5) RandomHeteropolymer uses "2bead" and "3bead" objects as building-blocks,
|
||||
# and these objects inherit the properties from the "ForceField" object
|
||||
# (atom types, bond types, etc, defined in "forcefield.lt" in this example).
|
||||
# So I declared that "RandomHeteropolymer inherits ForceField {" so that
|
||||
# you can easily access those atom types, bond types, etc as well.
|
||||
|
|
@ -1,11 +0,0 @@
|
|||
|
||||
write_once("Data Boundary") {
|
||||
0.0 150.0 xlo xhi
|
||||
0.0 150.0 ylo yhi
|
||||
0.0 150.0 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
import "polymer.lt"
|
||||
|
||||
polymer = new RandomHeteropolymer
|
||||
|
|
@ -1,29 +0,0 @@
|
|||
# -- Init Section --
|
||||
|
||||
include system.in.init
|
||||
|
||||
# -- Atom Definition Section --
|
||||
|
||||
read_data system.data
|
||||
|
||||
# -- Settings Section --
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# -- Run Section --
|
||||
|
||||
timestep 2.0
|
||||
dump 1 all custom 1000 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
|
||||
thermo 1000 # time interval for printing out "thermo" data
|
||||
|
||||
fix fxlan all langevin 300.0 300.0 120 48279
|
||||
fix fxnve all nve
|
||||
|
||||
# Temperature = 500 degrees
|
||||
|
||||
run 500000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
|
||||
This example contains a (crude and somewhat simple) example of
|
||||
the translocation of a (rather short) polymer through a hole in a wall,
|
||||
surrounded by an explicit LJ solvent.
|
||||
|
||||
(I used a short polymer because a longer polymer would require a larger box.
|
||||
But this example looked more impressive visually when I used a smaller box.)
|
||||
----
|
||||
Note: You must compile LAMMPS with the optional "RIGID" package installed. To
|
||||
do this, go to the "src" directory of your lammps installation and type:
|
||||
make yes-RIGID
|
||||
make clean-all
|
||||
make NAME_OF_TARGET #<--("make ubuntu", "make g++", "make linux".)
|
||||
----
|
||||
|
||||
Instructions on how to build LAMMPS input files and
|
||||
run a short simulation are provided in other README files.
|
||||
|
||||
step 1)
|
||||
README_setup.sh
|
||||
|
||||
step2)
|
||||
README_run.sh
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
# --- Running LAMMPS ---
|
||||
# -- Prerequisites: --
|
||||
# The "run.in.nvt" file is a LAMMPS input script containing
|
||||
# references 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.nvt # Run a simulation at constant volume
|
||||
|
||||
#or
|
||||
|
||||
lmp_mpi -i run.in.npt # Run a simulation at constant pressure
|
||||
# (Note: Constant pressure conditions have not been
|
||||
# well tested. The "run.in.npt" script may fail.)
|
||||
|
||||
|
||||
|
||||
# If you have compiled the MPI version of lammps, you can run lammps in parallel
|
||||
#mpirun -np 4 lmp_mpi -i run.in.nvt
|
||||
#or
|
||||
#mpirun -np 4 lmp_mpi -i run.in.npt
|
||||
# (assuming you have 4 processors available)
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
# Use these commands to generate the LAMMPS input script and data file
|
||||
# (and other auxilliary files):
|
||||
|
||||
|
||||
# Create LAMMPS input files this way:
|
||||
cd moltemplate_files
|
||||
|
||||
# run moltemplate
|
||||
|
||||
moltemplate.sh system.lt
|
||||
|
||||
# This will generate various files with names ending in *.in* and *.data.
|
||||
# These files are the input files directly read by LAMMPS. Move them to
|
||||
# the parent directory (or wherever you plan to run the simulation).
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# 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 ../
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
|
||||
------- 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.
|
Before Width: | Height: | Size: 5.4 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 23 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 103 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 13 KiB |
|
|
@ -1,111 +0,0 @@
|
|||
# ----------------------------------------------------------------------
|
||||
# -- General comment: --
|
||||
# -- The write() and write_once() commands create and append text to --
|
||||
# -- files (replacing variables beginning with @ or $ with counters.) --
|
||||
# -- File names beginning with "In " or "Data " are special. --
|
||||
# -- They will be pasted into the LAMMPS input script and --
|
||||
# -- data files which are generated by moltemplate. The syntax --
|
||||
# -- of these files is exactly the same as the syntax from the --
|
||||
# -- corresponding sections of a LAMMPS input script or data file. --
|
||||
# ----------------------------------------------------------------------
|
||||
|
||||
|
||||
Monomer {
|
||||
|
||||
# atom-id mol-id atom-type charge x y z
|
||||
|
||||
write("Data Atoms") {
|
||||
$atom:CA $mol:... @atom:CA 0.0 0.000 0.4000 0.00000
|
||||
$atom:R1 $mol:... @atom:R 0.0 0.000 1.000 1.000
|
||||
$atom:R2 $mol:... @atom:R 0.0 0.000 1.000 -1.000
|
||||
}
|
||||
|
||||
# Note: The "..." in "$mol:..." tells moltemplate that this molecule may
|
||||
# be a part of a larger molecule, and (if so) to use the larger
|
||||
# parent object's molecule id number as it's own
|
||||
|
||||
|
||||
# atom-type mass
|
||||
|
||||
write_once("Data Masses") {
|
||||
@atom:CA 13.0
|
||||
@atom:R 50.0
|
||||
}
|
||||
|
||||
# atom-type atom-type epsilon sigma
|
||||
|
||||
write_once("In Settings") {
|
||||
pair_coeff @atom:CA @atom:CA 0.05 2.0
|
||||
pair_coeff @atom:R @atom:R 0.50 2.0
|
||||
}
|
||||
|
||||
# bond-id bond-type atom-id1 atom-id2
|
||||
|
||||
write("Data Bonds") {
|
||||
$bond:CR1 @bond:sidechain $atom:CA $atom:R1
|
||||
$bond:CR2 @bond:sidechain $atom:CA $atom:R2
|
||||
}
|
||||
|
||||
write_once("In Settings") {
|
||||
# bond-type k r0
|
||||
bond_coeff @bond:sidechain 30.0 1.2
|
||||
bond_coeff @bond:bb 30.0 2.0 # "bb" shorthand for "backbone"
|
||||
}
|
||||
|
||||
# For a compound molecule consisting of smaller building blocks (such as a
|
||||
# polymer built from monomers), it is tedious to explicitly list all of the
|
||||
# angles, dihedrals in the entire molecule. Instead, you can define rules
|
||||
# for automatically generating all the angular interactions between bonded
|
||||
# atoms according to their connectivity and the atom/bond type.
|
||||
# Later, when you connect multiple monomers together to form a polymer,
|
||||
# appropriate bond-angle forces will be applied to these atoms automatically
|
||||
# (as well as dihedral and improper forces, if defined).
|
||||
|
||||
# Rules for determining 3 and 4-body bonded interactions by type
|
||||
|
||||
# angle-type atomType1 atomType2 atomType3 bondType1 bondType2
|
||||
|
||||
write_once("Data Angles By Type") {
|
||||
@angle:backbone @atom:CA @atom:CA @atom:CA @bond:* @bond:*
|
||||
@angle:sidechain @atom:CA @atom:CA @atom:R @bond:* @bond:*
|
||||
@angle:RCR @atom:R @atom:CA @atom:R @bond:* @bond:*
|
||||
}
|
||||
|
||||
# ("@angle:RCR" defines the angle between the R-C-R atoms within a monomer.
|
||||
# The other angular interactions are between atoms in neighboring monomers.)
|
||||
|
||||
|
||||
# dihedral-type AtomType1 AtomType2 AtomType3 AtomType4 bondType1 btyp2 btyp3
|
||||
|
||||
write_once("Data Dihedrals By Type") {
|
||||
@dihedral:backbn @atom:CA @atom:CA @atom:CA @atom:CA @bond:* @bond:* @bond:*
|
||||
}
|
||||
|
||||
# Parameters for these new angular interactions must be defined. (I recommend
|
||||
# putting all force-field parameters (coeffs) in the "In Settings" section.)
|
||||
|
||||
write_once("In Settings") {
|
||||
# angle-type k theta0
|
||||
angle_coeff @angle:backbone 50.00 160
|
||||
angle_coeff @angle:sidechain 50.00 120
|
||||
angle_coeff @angle:RCR 50.00 120
|
||||
# dihedral-type K1 K2 K3 K4
|
||||
dihedral_coeff @dihedral:backbn 1.411036 -0.271016 3.145034 0.0
|
||||
}
|
||||
|
||||
} # Monomer
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# -------------------------------------------------------------------------
|
||||
# Heteropolymers:
|
||||
#
|
||||
# There is a similar example for heteropolymers which is distributed online
|
||||
# bundled with the moltemplate software. It is named "2bead_heteropolymer",
|
||||
# and it demonstrates how to share backbone (CA) atoms, bonds and angles
|
||||
# (so that you don't have to define them seperately for each type of monomer).
|
||||
# -------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -1,32 +0,0 @@
|
|||
import "monomer.lt"
|
||||
|
||||
Polymer {
|
||||
|
||||
create_var {$mol} # optional:force all monomers to share the same molecule-ID
|
||||
# (The "Data Atoms" in Monomer must use "$mol:..." notation.)
|
||||
|
||||
# Make a chain of monomers
|
||||
monomers = new Monomer [12].rot(180, 1,0,0).move(2.0, 0, 0)
|
||||
|
||||
|
||||
# Now, link the monomers together this way:
|
||||
write("Data Bonds") {
|
||||
$bond:bb1 @bond:Monomer/bb $atom:monomers[0]/CA $atom:monomers[1]/CA
|
||||
$bond:bb2 @bond:Monomer/bb $atom:monomers[1]/CA $atom:monomers[2]/CA
|
||||
$bond:bb3 @bond:Monomer/bb $atom:monomers[2]/CA $atom:monomers[3]/CA
|
||||
$bond:bb4 @bond:Monomer/bb $atom:monomers[3]/CA $atom:monomers[4]/CA
|
||||
$bond:bb5 @bond:Monomer/bb $atom:monomers[4]/CA $atom:monomers[5]/CA
|
||||
$bond:bb6 @bond:Monomer/bb $atom:monomers[5]/CA $atom:monomers[6]/CA
|
||||
$bond:bb7 @bond:Monomer/bb $atom:monomers[6]/CA $atom:monomers[7]/CA
|
||||
$bond:bb8 @bond:Monomer/bb $atom:monomers[7]/CA $atom:monomers[8]/CA
|
||||
$bond:bb9 @bond:Monomer/bb $atom:monomers[8]/CA $atom:monomers[9]/CA
|
||||
$bond:bb10 @bond:Monomer/bb $atom:monomers[9]/CA $atom:monomers[10]/CA
|
||||
$bond:bb11 @bond:Monomer/bb $atom:monomers[10]/CA $atom:monomers[11]/CA
|
||||
}
|
||||
|
||||
} # Polymer
|
||||
|
||||
|
||||
|
||||
# Angle, dihedral and improper interactions will be generated
|
||||
# automatically according to the instructions in "monomer.lt"
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
###################### SOLVENT #########################
|
||||
|
||||
import "solvent_single.lt"
|
||||
|
||||
# Fill the simulation box with a solvent.
|
||||
# In this example, the solvent is made of many
|
||||
# copies of "MoleculeA" (which has only one atom).
|
||||
|
||||
solvent = new MoleculeA [12].move(3.0,0,0)
|
||||
[12].move(0,3.0,0)
|
||||
[12].move(0,0,3.0)
|
||||
|
||||
# To start with a reasonable conformation, it's a good idea to delete the
|
||||
# solvent where the walls or the polymer is going to be. Here we do it manually:
|
||||
|
||||
delete solvent[*][*][2] # <-- 1st wall will go here
|
||||
delete solvent[*][*][8] # <-- 2nd wall will go here
|
||||
delete solvent[6-7][0-8][5-6] # <-- polymer will go here
|
||||
|
||||
# Alternate notation:
|
||||
# [a:b] notation also works, however the "b" is a strict upper bound...
|
||||
# ...hence the last line is equivalent to "delete solvent[6:8][0:9][5:7]"
|
||||
# [a*b] notation also works, and is equivalent to [a-b]
|
||||
|
|
@ -1,22 +0,0 @@
|
|||
# The two files "solvent_single.lt" and "wall_single.lt"
|
||||
# define two very simple molecules containing one atom each.
|
||||
# Both atoms have a similar size (the have the same sigma parameter).
|
||||
|
||||
|
||||
MoleculeA {
|
||||
|
||||
# atomID molID atomType charge x y z
|
||||
write("Data Atoms") {
|
||||
$atom:a $mol:. @atom:a 0.0 0.0 0.0 0.0
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:a 10.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
# i j epsilon sigma cutoff
|
||||
pair_coeff @atom:a @atom:a 0.60 3.0 7.5 #<--attractive
|
||||
group groupA type @atom:a #(Atoms of this type belong to the "A" group)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -1,57 +0,0 @@
|
|||
|
||||
# LAMMPS supports a large number of force-field styles. We must select
|
||||
# which ones we need. This information belongs in the "In Init" section (and
|
||||
# (you can specify it anywhere in your LT files, multiple times if you like).
|
||||
# If different molecules use different force-field styles, you can use hybrid
|
||||
# styles. (In this example the molecules share the same pair_style.)
|
||||
|
||||
write_once("In Init") {
|
||||
units real
|
||||
atom_style full
|
||||
bond_style harmonic
|
||||
angle_style harmonic
|
||||
dihedral_style opls
|
||||
pair_style lj/cut 9.0
|
||||
# If you have charged molecules immersed in a salty implicit
|
||||
# solvent, you might try something like this this instead:
|
||||
# pair_style lj/cut/coul/debye 0.1 9.0
|
||||
pair_modify mix arithmetic
|
||||
dielectric 80.0
|
||||
special_bonds lj 0.0 0.0 0.0
|
||||
}
|
||||
|
||||
|
||||
write_once("Data Boundary") {
|
||||
0.0 36.0 xlo xhi
|
||||
0.0 36.0 ylo yhi
|
||||
-18.0 18.0 zlo zhi
|
||||
}
|
||||
|
||||
|
||||
import "solvent.lt"
|
||||
|
||||
import "walls.lt"
|
||||
|
||||
import "polymer.lt"
|
||||
|
||||
polymer = new Polymer
|
||||
polymer.rot(-90.0, 0,0,1) # rotate it -90 degrees around the Y axis
|
||||
polymer.move(19.5,22.5,16.5) # move it near the openning of the hole
|
||||
|
||||
|
||||
|
||||
####################### Notes: #########################
|
||||
#
|
||||
# In this example we deleted solvent and wall molecule objects.
|
||||
# You can also delete a monomer inside the polymer. To do that use:
|
||||
# delete polymer/monomers[6]
|
||||
# You can also delete individual atoms, bonds, angles, dihedrals, & impropers
|
||||
# from existing molecules. For example to delete an atom in the middle
|
||||
# of the polymer try this. (Bonds and other interactions will also be removed.)
|
||||
# delete polymer/monomers[6]/CA
|
||||
# To delete a bond, try this
|
||||
# delete polymer/bb6
|
||||
# Note: This will not delete the angular interactions if they were explicitly
|
||||
# defined (ie, using "Data Angles" instead of "Data Angles By Type").
|
||||
# Delete explicit angle, dihedral, and improper interactions manually.
|
||||
# Note: In both cases the two molecule fragments will keep the same mol counter.
|
||||
|
|
@ -1,21 +0,0 @@
|
|||
# The two files "solvent_single.lt" and "wall_single.lt"
|
||||
# define two very simple molecules containing one atom each.
|
||||
# Both atoms have a similar size (the have the same sigma parameter).
|
||||
|
||||
MoleculeB {
|
||||
|
||||
# atomID molID atomType charge x y z
|
||||
write("Data Atoms") {
|
||||
$atom:b $mol:. @atom:b 0.0 0.0 0.0 0.0
|
||||
}
|
||||
write_once("Data Masses") {
|
||||
@atom:b 10.0
|
||||
}
|
||||
write_once("In Settings") {
|
||||
# i j epsilon sigma cutoff
|
||||
pair_coeff @atom:b @atom:b 0.05 3.0 7.5 #<--repulsive (approximately)
|
||||
group groupB type @atom:b #(Atoms of this type belong to the "B" group)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -1,23 +0,0 @@
|
|||
####################### WALLS ##########################
|
||||
|
||||
import "wall_single.lt"
|
||||
|
||||
# Create a wall at position z=6.0 (6.0 = 2*3.0)
|
||||
|
||||
wall1 = new MoleculeB [12].move(3.0, 0, 0)
|
||||
[12].move(0, 3.0, 0)
|
||||
|
||||
wall1[*][*].move(0,0,6.0)
|
||||
|
||||
# Create a second wall at position z=24.0 (24.0 = 8*3.0)
|
||||
|
||||
wall2 = new MoleculeB [12].move(3.0, 0, 0)
|
||||
[12].move(0, 3.0, 0)
|
||||
|
||||
wall2[*][*].move(0,0,24.0)
|
||||
|
||||
# Now delete some of the molecules in "wall2" to create a hole.
|
||||
|
||||
delete wall2[6-7][6-9]
|
||||
delete wall2[5-8][7-8]
|
||||
|
||||
|
|
@ -1,120 +0,0 @@
|
|||
# THIS EXAMPLE HAS NOT BEEN RIGOROUSLY TESTED.
|
||||
# (This simulation may fail.
|
||||
# However the "run.in.nvt" example in this directory should work.)
|
||||
#
|
||||
# Requirements:
|
||||
# To run this system at constant pressure, it might help to compile LAMMPS with
|
||||
# the optional RIGID package, and use "fix rigid" on the carbon. (Optional.)
|
||||
# The use of fix rigid is controversial. This method is demonstrated below.
|
||||
|
||||
# ------------------------------- Initialization Section --------------------
|
||||
|
||||
include system.in.init
|
||||
|
||||
# ------------------------------- Atom Definition Section -------------------
|
||||
|
||||
read_data system.data
|
||||
|
||||
# ------------------------------- Settings Section --------------------------
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# ------------------------------- Run Section -------------------------------
|
||||
|
||||
|
||||
# Only the groupB atoms are immobile.
|
||||
group mobile subtract all groupB
|
||||
|
||||
# Unfortunately you can not use the LAMMPS "minimize" command on this system
|
||||
# because there is no way to immobilize the wall atoms during minimization.
|
||||
# Instead, we can use langevin dynamics with a fast
|
||||
# damping parameter and a small timestep.
|
||||
|
||||
print "--------- beginning minimization (using fix langevin) ---------"
|
||||
|
||||
timestep 0.1
|
||||
fix fxlan mobile langevin 1.0 1.0 100.0 48279
|
||||
fix fxnve mobile nve # <-- needed by fix langevin (see lammps documentation)
|
||||
thermo 100
|
||||
run 2500
|
||||
|
||||
unfix fxlan
|
||||
unfix fxnve
|
||||
|
||||
# -- simulation protocol --
|
||||
|
||||
print "--------- beginning simulation (using fix nvt) ---------"
|
||||
|
||||
dump 1 all custom 1000 traj_npt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
|
||||
thermo 200 # time interval for printing out "thermo" data
|
||||
|
||||
|
||||
# ------------------------- NPT ---------------------------
|
||||
|
||||
|
||||
# ------ QUESTIONABLE (see below): ------
|
||||
|
||||
fix Ffreezestuff groupB rigid single force * off off off torque * off off off
|
||||
|
||||
# Comment:
|
||||
# The use of "fix rigid" to immobilize an object is somewhat controversial.
|
||||
# Feel free to omit it.
|
||||
# (Neither Trung or Steve Plimpton use fix rigid for immobilizing
|
||||
# molecules, but I noticed that at NPT, it does a better job of maintaining
|
||||
# the correct volume. However "fix rigid" has changed since then (2011),
|
||||
# so this may no longer be true. Please use this example with caution.)
|
||||
|
||||
|
||||
|
||||
# Thermostat+Barostat
|
||||
# Set temp=300K, pressure=200bar, and equilibrate volume only in the z direction
|
||||
|
||||
fix fxMoveStuff mobile npt temp 300 300 100 z 200 200 1000.0 dilate mobile drag 2.0
|
||||
|
||||
# ----------------------------------------
|
||||
|
||||
# The next two lines recalculate the temperature using
|
||||
# only the mobile degrees of freedom (ie. water atom velocities):
|
||||
|
||||
compute tempMobile mobile temp
|
||||
compute pressMobile all pressure tempMobile
|
||||
|
||||
thermo_style custom step c_tempMobile c_pressMobile temp press vol
|
||||
|
||||
fix_modify fxMoveStuff temp tempMobile
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
timestep 0.5
|
||||
|
||||
run 100000
|
||||
|
||||
timestep 1.0
|
||||
|
||||
run 100000
|
||||
|
||||
|
||||
write_data system_after_npt.data
|
||||
|
||||
# (The "write_restart" and "read_restart" commands were buggy in 2012,
|
||||
# but they should work also.)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# ----- Comment: Avoid using fix rigid/npt on large single rigid objects -----
|
||||
#
|
||||
# Use of the following is not recommended:
|
||||
#
|
||||
# fix Ffreezestuff groupB rigid/npt single temp 300 300 100 z 200 200 1000.0 force * off off off torque * off off off dilate mobile
|
||||
# (temp=300K, pressure=200bar, and equilibrate volume only in the z direction)
|
||||
#
|
||||
# In my experience, the system becomes unstable when applying "fix rigid/npt"
|
||||
# to the immobile atoms, while also applying "fix npt" on the solvent atoms.
|
||||
# (It is probably a bad idea to use two barostats simultaneously.)
|
||||
# ----------------------------------------------------------------------------
|
||||
|
|
@ -1,53 +0,0 @@
|
|||
# PREREQUISITES:
|
||||
#
|
||||
# You must use moltemplate.sh to create 3 files:
|
||||
# system.data system.in.init system.in.settings
|
||||
# (See README_setup.sh for details)
|
||||
|
||||
# -- Init Section --
|
||||
|
||||
include system.in.init
|
||||
|
||||
# -- Atom Definition Section --
|
||||
|
||||
read_data system.data
|
||||
|
||||
# -- Settings Section --
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# Use "neigh_modify" to turn off calculation of interactions between immobilized
|
||||
# atoms. (Note: The "groupB" group was defined in the file "system.insettings")
|
||||
neigh_modify exclude group groupB groupB
|
||||
|
||||
# -- Run Section --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 500 traj_nvt.lammpstrj id mol type x y z ix iy iz
|
||||
|
||||
thermo_style custom step temp pe etotal vol epair ebond eangle edihed
|
||||
thermo 500 # time interval for printing out "thermo" data
|
||||
|
||||
|
||||
# Optional: Improve efficiency by omitting the calcuation of interactions
|
||||
# between immobile atoms. (Note: This is not optional under NPT conditions.)
|
||||
neigh_modify exclude group groupB groupB
|
||||
|
||||
# Only the groupB atoms are immobile.
|
||||
group mobile subtract all groupB
|
||||
|
||||
# The next two lines recalculate the temperature
|
||||
# using only the mobile degrees of freedom:
|
||||
|
||||
compute tempMobile mobile temp
|
||||
|
||||
# Integrate the equations of motion:
|
||||
fix fMoveStuff mobile nvt temp 300.0 300.0 100.0
|
||||
fix_modify fMoveStuff temp tempMobile
|
||||
|
||||
|
||||
run 100000
|
||||
|
||||
write_data system_after_nvt.data
|
||||
|
||||
|
|
@ -1,48 +0,0 @@
|
|||
This is an implementation of the "two-stage" model used by Maxim Imakaev
|
||||
in the Naumova et Al 2013 Science paper on metaphase chromatin.
|
||||
(Download the supplemental materials section and scroll down to the section:
|
||||
"Two-stage process: linear compaction - axial compression")
|
||||
|
||||
---- SMALL MODIFICATION ----
|
||||
|
||||
Unlike that study, I did not use "softened" Lennard-Jones potentials
|
||||
(which allow the chains to pass through each other).
|
||||
|
||||
--- Why use moltemplate? ---
|
||||
|
||||
Honestly, you don't need to use moltemplate to build this polymer.
|
||||
It is almost counter-productive to use moltemplate to build this kind of
|
||||
polymer because it is so simple. (The polymer has only 1 bead per atom.
|
||||
It just makes it more complicated to introduce all these extra
|
||||
files including monomer.lt, condensin.lt and system.lt, especially considering
|
||||
that system.lt is a complex file which is generated by a separate script.)
|
||||
|
||||
However building the sytem using moltemplate may pay off if you
|
||||
replace each point-like monomer with a multi-atom molecule later on.
|
||||
(Right now, using moltemplate to build this system is sort of overkill.
|
||||
I'll post an example of building more complex models of chromatin eventually.)
|
||||
|
||||
Anyway, the two-stage model at the end of Naumova et al Science 2013 uses the "30nm-fiber" model, whose details are (somewhat vaguely) described in the supplemental materials section.
|
||||
|
||||
---- 10-nm fiber model: ----
|
||||
|
||||
For the 10nm model,
|
||||
n=128000,
|
||||
L=200,
|
||||
U(alpha)=5*(1 - cos(alpha))
|
||||
bond_length=1.0 (=10nm)
|
||||
sigma=1.0 (particle radius = 10nm)
|
||||
|
||||
---- 30-nm fiber model: ----
|
||||
|
||||
"The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level."
|
||||
|
||||
I interpret this to mean that, for the 30nm model,
|
||||
n=128000/4.25~=30117 (however I rounded up to 32768=2^15)
|
||||
L=200/4.25~=47 (however I rounded up to 51)
|
||||
U(alpha)=1.17647*(1 - cos(alpha)) (5/4.25=1.17647)
|
||||
|
||||
To increase the volume by a factor o 4.25, I increase both the diameter of each
|
||||
bead (the "sigma" parameter), and the bond-lengths connecting them from
|
||||
1.0 (corresponding to 10nm) to 4.25^(1/3)~=1.6198 (corresponding to 16.198nm).
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue