git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13900 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
@ -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,32 +0,0 @@
|
||||
# -------- REQUIREMENTS: ---------
|
||||
# You must define your MOLTEMPLATE_PATH environment variable
|
||||
# and set it to the "common" subdirectory of your moltemplate distribution.
|
||||
# (See the "Installation" section in the moltemplate manual.)
|
||||
|
||||
# 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
|
||||
0.0 36.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,70 +0,0 @@
|
||||
# -- Init Section --
|
||||
|
||||
include system.in.init
|
||||
|
||||
# -- Atom Definition Section --
|
||||
|
||||
read_data system.data
|
||||
|
||||
# -- Settings Section --
|
||||
|
||||
include system.in.settings
|
||||
|
||||
# -- Run Section --
|
||||
|
||||
|
||||
timestep 1.0
|
||||
dump 1 all custom 500 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 500 # time interval for printing out "thermo" data
|
||||
|
||||
|
||||
velocity groupB zero angular
|
||||
velocity groupB zero linear
|
||||
# (I'm not sure if the two lines above are necessary, but they don't hurt.)
|
||||
|
||||
# Only the groupB atoms are immobile.
|
||||
|
||||
group mobile subtract all groupB
|
||||
|
||||
# ------------------------- NPT ---------------------------
|
||||
|
||||
# Set temp=300K, pressure=500bar, and equilibrate volume only in the z direction
|
||||
|
||||
fix fMoveStuff mobile npt temp 300 300 100 z 500 500 1000.0 dilate mobile
|
||||
|
||||
# ------ QUESTIONABLE (see below): ------
|
||||
|
||||
fix Ffreezestuff groupB rigid/npt single temp 300 300 100 z 500 500 1000.0 force * off off off torque * off off off dilate mobile
|
||||
|
||||
# -- Alternate npt rigid method --
|
||||
# I'm not sure which way is more correct, however
|
||||
# this also seems to behave in a reasonable-looking way:
|
||||
#fix Ffreezestuff groupB rigid single force * off off off torque * off off off dilate mobile
|
||||
#
|
||||
# The use of either "fix rigid" or "fix rigid/npt" to immobilize
|
||||
# an object is somewhat controversial. Feel free to omit it.
|
||||
#(Neither Trung or Steve Plimpton use rigid or rigid/npt 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.)
|
||||
# ----------------------------------------
|
||||
|
||||
# IMPORTANT for NPT: You must use "neigh_modify" to turn off calculation of the
|
||||
# forces between immobilized atoms.
|
||||
neigh_modify exclude group groupB groupB
|
||||
|
||||
# The next two lines recalculate the temperature
|
||||
# using only the mobile degrees of freedom:
|
||||
|
||||
compute tempMobile mobile temp
|
||||
compute pressMobile all pressure tempMobile
|
||||
|
||||
thermo_style custom step c_tempMobile c_pressMobile temp press vol
|
||||
|
||||
fix_modify fMoveStuff temp tempMobile
|
||||
|
||||
run 60000
|
||||
|
||||
write_data system_after_npt.data
|
||||
@ -1,49 +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
|
||||
|
||||
# -- 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).
|
||||
|
||||
@ -1,29 +0,0 @@
|
||||
---------------
|
||||
The average diameter of a mammalian cell nucleus is is 6 micrometers (μm),
|
||||
which occupies about 10% of the total cell volume.
|
||||
|
||||
(See "Molecular Biology of the Cell, Chapter 4, pages 191–234 (4th ed.)", by
|
||||
Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, Peter Walter, 2002)
|
||||
|
||||
... of that, 25% of it is occupied by the nucleolus
|
||||
http://en.wikipedia.org/wiki/Nucleolus
|
||||
("citation needed")
|
||||
---------------
|
||||
|
||||
From the supplemental material for the original HiC paper
|
||||
(Lieberman-Aiden et al., Science 2009)
|
||||
|
||||
Appendix 1.
|
||||
Estimate of the Volume Fraction of Chromatin in Human Cells
|
||||
In the simulations we sought to obtain an ensemble of structures that, in their statistical properties, resemble some of the features of chromatin arrangement in the cell. Below we demonstrate that chromatin occupies a significant fraction of cell volume, a property that we reproduced in simulations. Taking the nuclear diameter of a tissue culture cell to be 5-10um, and assuming close to a spherical shape we obtain the volume in the range 50-500 um^3, with a (geometric) mean of ~160 um^3. If we assume that the chromatin is built of DNA wrapped around nucleosomes, then we have 6x10^9bp/200bp=3x10^7 nucleosomes. Each may be approximated as a cylinder ~10nm in diameter and ~5nm in height, suggesting a volume of about 500nm3 each. The linker DNA after each nucleosome is about 50bps long, suggesting a volume of about 50*0.34nm*3.14*1nm^2=50nm^3. Thus the total volume of chromatin = 550x3x10^7 =16 um^3, or ~10% (3-23%) of the nuclear volume. This strikingly large volume fraction is itself a significant underestimate, since we ignored, among other things, all other DNA-bound proteins. Note that any further packing or localization of chromatin inside the nucleus will increase local density.
|
||||
---- This next section mostly only justifies why they ----
|
||||
---- they did not stop the simulation when the globules ----
|
||||
---- were fully crumpled (ie with uniform density) ----
|
||||
In our simulations, the radius of the final crumpled globule was R≈12.5 and the volume V≈8000 cubic units. The total volume of the 4000 monomers, 1 unit in diameters each, is V≈2000. This implies a volume fraction of about 25%, which is consistent with the volume fraction estimated above.
|
||||
---- ----
|
||||
|
||||
Appendix 2.
|
||||
Monomer length in base pairs
|
||||
Each monomer of the chain corresponds to a fragment of chromatin that equals the Kuhn length of the chromatin fiber, i.e. approximately twice the persistence length of the fiber. Although the persistence length of the chromatin fiber is unknown it can be estimated using the following arguments. DNA is packed into nucleosomes, where 150 bps are wrapped around the histone core and do not contribute to flexibility of the fiber. The linker DNA of about 50 bps that connects consecutive nucleosomes is bendable, and is the source of flexibility in the fiber. Since the persistence length of double-stranded DNA is 150 bps, an equally flexible region of the nucleosomal DNA should contain 3 linkers, i.e. 3 consecutive nucleosomes packing about 600 bps of DNA. The excluded volume of the nucleosomes, nucleosome interactions, and other DNA-bound proteins can make the fiber less flexible or prohibit certain conformation and may tend to increase the persistence length of the fiber. Using this estimated lower bound estimate for the persistence length, we obtain the Kuhn length of the equivalent freely-jointed chain to be 6 nucleosomes, or ~ 1200bp. A simulated chain of 4000 monomers corresponds to 4.8Mb of packed DNA. The size of each monomer was chosen such that its volume is equal to (or slightly above) that of 6 nucleosomes (V=6 x 600 nm^3); thus the radius of the spherical monomer is R=10nm. The diameter of each globule shown in Figure 4 is about 200 nm.
|
||||
|
||||
|
||||
@ -1,7 +0,0 @@
|
||||
# Run lammps using the following 3 commands:
|
||||
# (assuming "lmp_mpi" is the name of your LAMMPS binary)
|
||||
|
||||
lmp_mpi -i run.in.min
|
||||
lmp_mpi -i run.in.stage1
|
||||
lmp_mpi -i run.in.stage2
|
||||
|
||||
@ -1,58 +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
|
||||
|
||||
# First, rescale and interpolate the positions
|
||||
# where the monomers will be located. (This step is
|
||||
# not needed if the coords_orig.raw file already has correct coordinates.)
|
||||
|
||||
./interpolate_coords.py 32768 1.6198059 < coords_orig.raw > coords.raw
|
||||
|
||||
# Then, build the "system.lt" file
|
||||
|
||||
./generate_system_lt.py 32768 51 < coords.raw > system.lt
|
||||
|
||||
# 32768 is the number of monomers in the polymer
|
||||
# (which may be different from the number of coordinates
|
||||
# in the "coords_orig.raw" file) This number will vary
|
||||
# depending on how long you want the polymer to be.
|
||||
# The second argument "51" is the average interval between
|
||||
# condensin anchors (IE the "loop size" in monomers.)
|
||||
|
||||
|
||||
# Run moltemplate
|
||||
|
||||
moltemplate.sh system.lt -a "@bond:stage1 1" \
|
||||
-a "@bond:stage2 2" \
|
||||
-a "@atom:Monomer/A 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).
|
||||
#
|
||||
# We used the "-a" command to set the variable @bond:condensin to "2"
|
||||
# because we will refer to it later in the "run.in" LAMMPS input script.
|
||||
# (Of coarse, LAMMPS knows nothing about moltemplate variables,
|
||||
# so in that file we refer to it as dihedral type "1")
|
||||
|
||||
mv -f system.in* system.data ../
|
||||
|
||||
# We also need the table of bond forces used during "stage 2".
|
||||
# (Like the system.data and the various input scripts, this file is needed by
|
||||
# LAMMPS, so we need to copy it to the directory where we will run the sim.)
|
||||
cp -f table_bonds_stage2.dat ../
|
||||
|
||||
# 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/
|
||||
|
||||
# Optional:
|
||||
# Remove the "system.lt" file created by "generate_system_lt.py"
|
||||
#rm -f system.lt
|
||||
|
||||
cd ../
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user