jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9808 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-05-07 21:14:00 +00:00
parent 03aeb17625
commit b7da4d691c
335 changed files with 0 additions and 76240 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
tools/moltemplate/LICENSE.TXT
View File

@ -1,28 +0,0 @@
Author: Andrew Jewett, Shea Group, http://www.chem.ucsb.edu/~sheagroup/
Copyright (c) 2012, Regents of the University of California
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the name of the University of California, Santa Barbara nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

61
tools/moltemplate/README.TXT
View File

@ -1,61 +0,0 @@
-- Description: --
Moltemplate is a cross-platform text-based molecule builder for LAMMPS.
-- Typical usage: --
moltemplate.sh [-atomstyle style] [-pdb/-xyz coord_file] file.lt
-- Web page: --
Documentation, examples, and supporting code can be downloaded at:
http://www.moltemplate.org
The most up-to-date version of moltemplate can be downloaded here.
(After download, you can unpack the archive using:
tar xzf moltemplate_2012-3-31.tar.gz
The date will vary from version to version.)
----------------------------------------------------
---------- INSTALLATION INSTRUCTIONS: ------------
----------------------------------------------------
This directory should contain two folders:
src/ <-- location of all python and bash scripts
common/ <-- location of shared force fields and molecules
The ``moltemplate.sh'' script and the python scripts that it invokes are
located in the ``src/'' subdirectory. You should update your PATH environment
variable to include this directory.
If you do not know what a PATH environment variable is, read:
http://www.linfo.org/path_env_var.html
(I receive this question often.)
It is also a good idea to set your MOLTEMPLATE_PATH environment variable to
point to the ``common/'' subdirectory.
(Force fields and commonly used molecules will eventually be located here.)
-- Installation example ---
Suppose the directory with this README.TXT file is located at ~/moltemplate.
If you use the bash shell, typically you would edit your
~/.profile, ~/.bash_profile or ~/.bashrc files to contain the following lines:
export PATH="$PATH:$HOME/moltemplate/src"
export MOLTEMPLATE_PATH="$HOME/moltemplate/common"
If you use the tcsh shell, typically you would edit your
~/.login, ~/.cshrc, or ~/.tcshrc files to contain the following lines:
setenv PATH "$PATH:$HOME/moltemplate/src"
setenv MOLTEMPLATE_PATH "$HOME/moltemplate/common"
-- Requirements: --
Moltemplate requires the Bourne-shell, and a recent version of python
(2.7, 3.0 or higher), and can run on OS X, linux, or windows (if a
suitable shell environment has been installed).
-- License: --
Moltemplate is available under the terms of the open-source 3-clause BSD
license. (See LICENSE.TXT.)

70
tools/moltemplate/common/graphene.lt
View File

@ -1,70 +0,0 @@
# You can use the definitions in this file to
# create graphene, graphite, or nanotubes.
Graphene {
2AtomCellAlignX {
# The 2AtomCellAlignX "molecule" is a minimal basis cell
# for any hexagonal arrangement of atoms in 2-dimensions.
# The distance between these two atoms is equal to "l", where "l" is
# the length of each side of a hexegon, which I set to 1.420 Angstroms.
#
# l = length of each hexagonal side = 1.42 Angstroms
# L = length of each hexagon = 2*l = 2.84 Angstroms
# W = width of each hexagon = 2*l*sqrt(3)/2 = 2.4595121467478056 Angstroms
#
# The Lattice-cell vectors for graphene (graphite) are
# (2.4595121467478, 0, 0) (aligned with X axis)
# (1.2297560733739, 2.13, 0) (2.13 = 1.5*l)
# ( 0 1.00, 3.35) (3.35 = distance between layers of
# graphene sheets in graphite.)
# atomID molID atomType charge x y z
write("Data Atoms") {
$atom:C1 $mol:... @atom:../C 0.0 -0.61487803668695 -0.355 0.0
$atom:C2 $mol:... @atom:../C 0.0 0.61487803668695 0.355 0.0
}
} # 2AtomCellAlignX
# Notice that the two atoms in the unit-cell above lie in the XY plane.
# (Their z-coord is zero). It's also handy to have a version of this
# object which lies in the XZ plan, so we make this below.
2AtomCellAlignXZ = 2AtomCellAlignX.rot(90,1,0,0)
# Define properties of the Carbon-Hydrogen graphene atom
write_once("Data Masses") {
@atom:C 12.0
}
write_once("In Init") {
pair_style hybrid lj/charmm/coul/charmm 9.0 10.0
}
write_once("In Settings") {
# Define a group consisting of only carbon atoms in graphene molecules
group Cgraphene type @atom:C
# i j epsilon sigma
pair_coeff @atom:C @atom:C lj/charmm/coul/charmm 0.068443 3.407
# These Lennard-Jones parameters come from
# R. Saito, R. Matsuo, T. Kimura, G. Dresselhaus, M.S. Dresselhaus,
# Chem Phys Lett, 348:187 (2001)
}
} # Graphene

52
tools/moltemplate/common/spce.lt
View File

@ -1,52 +0,0 @@
# file "spce.lt"
#
# H1 H2
# \ /
# O
SPCE {
write_once("In Init") {
# -- Default styles (for solo "SPCE" water) --
units real
atom_style full
# (Hybrid force fields were not necessary but are used for portability.)
pair_style hybrid lj/charmm/coul/long 9.0 10.0 10.0
bond_style hybrid harmonic
angle_style hybrid harmonic
kspace_style pppm 0.0001
pair_modify mix arithmetic
}
write("Data Atoms") {
$atom:O $mol:. @atom:O -0.8476 0.0000000 0.00000 0.000000
$atom:H1 $mol:. @atom:H 0.4238 0.8164904 0.00000 0.5773590
$atom:H2 $mol:. @atom:H 0.4238 -0.8164904 0.00000 0.5773590
}
write_once("Data Masses") {
@atom:O 15.9994
@atom:H 1.008
}
write("Data Bonds") {
$bond:OH1 @bond:OH $atom:O $atom:H1
$bond:OH2 @bond:OH $atom:O $atom:H2
}
write("Data Angles") {
$angle:HOH @angle:HOH $atom:H1 $atom:O $atom:H2
}
write_once("In Settings") {
bond_coeff @bond:OH harmonic 200.0 1.0
angle_coeff @angle:HOH harmonic 200.0 109.47
pair_coeff @atom:O @atom:O lj/charmm/coul/long 0.1553 3.166
pair_coeff @atom:H @atom:H lj/charmm/coul/long 0.0 2.058
group spce type @atom:O @atom:H
fix fSHAKE spce shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fSHAKE during minimization.)
}
} # end of definition of "SPCE" water molecule type

13
tools/moltemplate/common/tip3p_attempts/README.TXT
View File

@ -1,13 +0,0 @@
This directory contains two LT files corresponding to
different versions of TIP3P:
tip3pcharmm.lt # The implementation of TIP3P used by CHARMM (I think).
tip3p2004.lt # The newer Price & Brooks, J. Chem Phys 2004 model
# which uses long-range coulombics
I have not tested these files so I moved them here.
(If you have tested these files, and they work, or if you have other comments
or suggestions, feel free to email me at jewett.aij at gmail dot com.)
Andrew
2012-10-20

88
tools/moltemplate/common/tip3p_attempts/tip3p2004.lt
View File

@ -1,88 +0,0 @@
# file "tip3p2004.lt"
#
# H1 H2
# \ /
# O
#
# I think this is the TIP3P water described in the paper by
# Daniel J. Price and Charles L. Brooks III
# J. Chem. Phys., 121(20): 10096 (2004)
# Specifically I think it refers to the "Model B" version of long-range TIP3P
# described in the 3rd-to-last column of "Table I", on p.10099.
TIP3P2004 {
write_once("In Init") {
# -- Default styles (for solo "TIP3P2004" water) --
units real
atom_style full
pair_style hybrid lj/charmm/coul/long 10.0 10.5 10.5
bond_style hybrid harmonic
angle_style hybrid harmonic
kspace_style pppm 0.0001
pair_modify mix arithmetic
}
write("Data Atoms") {
$atom:O $mol:. @atom:O -0.830 0.0000000 0.00000 0.000000
$atom:H1 $mol:. @atom:H 0.415 0.756950327 0.00000 0.5858822766
$atom:H2 $mol:. @atom:H 0.415 -0.756950327 0.00000 0.5858822766
}
write_once("Data Masses") {
@atom:O 15.9994
@atom:H 1.008
}
write("Data Bonds") {
$bond:OH1 @bond:OH $atom:O $atom:H1
$bond:OH2 @bond:OH $atom:O $atom:H2
}
write("Data Angles") {
$angle:HOH @angle:HOH $atom:H1 $atom:O $atom:H2
}
write_once("In Settings") {
bond_coeff @bond:OH harmonic 450.0 0.9572
angle_coeff @angle:HOH harmonic 55.0 104.52
#########################################################################
#### There are two choices for for the O-O interactions
#########################################################################
#### O-O nonbonded interactions
# For the 1983 Jorgensen version of TIP3P use:
# pair_coeff @atom:O @atom:O lj/charmm/coul/charmm 0.1521 3.1507
# For the 2004 Price & Brooks version of TIP3P use:
pair_coeff @atom:O @atom:O lj/charmm/coul/charmm 0.102 3.188
#########################################################################
#### There are three choices for for the O-H and H-H interactions
#########################################################################
#### 1) CHARMM uses an arithmetic mixing-rule for the O-H sigma parameter
pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.0460 0.4000
pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.0836 1.7753 #arithmetic
#########################################################################
#### 2) OPLS-AA uses geometric a mixing-fule for the O-H sigma parameter,
#### If you want to use this, uncomment the following two lines:
# pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.0460 0.4000
# pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.0836 1.1226 #geometric
#########################################################################
#### 3) The original Jorgensen 1983 parameterization has no OH or HH
# lennard-jones interactions. For this behavior, uncomment these lines:
# pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.00 0.4000
# pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.00 1.7753
#########################################################################
# Define a group for the tip3p water molecules:
group tip3p type @atom:O @atom:H
# Optional: Constrain the angles and distances.
# (Most implementations use this, but it is optional.)
fix fSHAKE tip3p shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fSHAKE during minimization.)
}
} # "TIP3P2004" water molecule type

91
tools/moltemplate/common/tip3p_attempts/tip3pcharmm.lt
View File

@ -1,91 +0,0 @@
# file "tip3p_charmm.lt"
#
# H1 H2
# \ /
# O
#
# I think this is the TIP3P water model used by CHARMM (and probably AMBER)
# It is (mostly) based on this paper:
# Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem Phys, 79, 926 (1983)
TIP3Pcharmm {
write_once("In Init") {
# -- Default styles (for solo "TIP3Pcharmm" water) --
units real
atom_style full
# I'm not sure exactly which cutoffs distances are traditionally used in
# the TIP3P water model used by CHARMM.
# (See the Price JCP 2004 paper for a review.)
# pair_style hybrid lj/charmm/coul/charmm 7.5 8.0 10.0 10.5
# Try this instead:
pair_style hybrid lj/charmm/coul/charmm 10.0 10.5 10.0 10.5
bond_style hybrid harmonic
angle_style hybrid harmonic
pair_modify mix arithmetic
}
write("Data Atoms") {
$atom:O $mol:. @atom:O -0.834 0.0000000 0.00000 0.000000
$atom:H1 $mol:. @atom:H 0.417 0.756950327 0.00000 0.5858822766
$atom:H2 $mol:. @atom:H 0.417 -0.756950327 0.00000 0.5858822766
}
write_once("Data Masses") {
@atom:O 15.9994
@atom:H 1.008
}
write("Data Bonds") {
$bond:OH1 @bond:OH $atom:O $atom:H1
$bond:OH2 @bond:OH $atom:O $atom:H2
}
write("Data Angles") {
$angle:HOH @angle:HOH $atom:H1 $atom:O $atom:H2
}
write_once("In Settings") {
bond_coeff @bond:OH harmonic 450.0 0.9572
angle_coeff @angle:HOH harmonic 55.0 104.52
#########################################################################
#### There are two choices for for the O-O interactions
#########################################################################
#### O-O nonbonded interactions
# For the 1983 Jorgensen version of TIP3P use:
pair_coeff @atom:O @atom:O lj/charmm/coul/charmm 0.1521 3.1507
# For the 2004 Price & Brooks version of TIP3P use:
# pair_coeff @atom:O @atom:O lj/charmm/coul/charmm 0.102 3.188
#########################################################################
#### There are three choices for for the O-H and H-H interactions
#########################################################################
#### 1) CHARMM uses an arithmetic mixing-rule for the O-H sigma parameter
pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.0460 0.4000
pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.0836 1.7753 #arithmetic
#########################################################################
#### 2) OPLS-AA uses geometric a mixing-fule for the O-H sigma parameter,
#### If you want to use this, uncomment the following two lines:
# pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.0460 0.4000
# pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.0836 1.1226 #geometric
#########################################################################
#### 3) The original Jorgensen 1983 parameterization has no OH or HH
# lennard-jones interactions. For this behavior, uncomment these lines:
# pair_coeff @atom:H @atom:H lj/charmm/coul/charmm 0.00 0.4000
# pair_coeff @atom:O @atom:H lj/charmm/coul/charmm 0.00 1.7753
#########################################################################
# Define a group for the tip3p water molecules:
group tip3p type @atom:O @atom:H
# Optional: Constrain the angles and distances.
# (Most implementations use this, but it is optional.)
fix fSHAKE tip3p shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fSHAKE during minimization.)
}
} # "TIP3Pcharmm" water molecule type

50
tools/moltemplate/common/trappe1998.lt
View File

@ -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

33
tools/moltemplate/common/watmw.lt
View File

@ -1,33 +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") {
group WatMW type @atom:mW #(Atoms of this type belong to the "WatMW" group)
}
} # WatMW

56
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/README.TXT
View File

@ -1,56 +0,0 @@
This example shows how to put a protein (inclusion) in a
lipid bilayer mixture composed of two different lipids (DPPC and DLPC).
The DPPC lipid model is described here:
G. Brannigan, P.F. Philips, and F.L.H. Brown,
Physical Review E, Vol 72, 011915 (2005)
(The new DLPC model is a truncated version of DPPC.)
The protein model is described here:
G. Bellesia, AI Jewett, and J-E Shea,
Protein Science, Vol19 141-154 (2010)
Note:
This example may require additional features to be added to LAMMPS.
If LAMMPS complains about an "Invalid pair_style", then copy the code
in the "additional_lammps_code" directory into your LAMMPS "src" directory
and recompile LAMMPS.
----- Details --------
This example contains a coarse-grained model of a 4-helix bundle protein
inserted into a lipid bilayer (made from a mixture of DPPC and DLPC).
-- Protein Model: --
The coarse-grained protein is described in:
G. Bellesia, AI Jewett, and J-E Shea, Protein Science, Vol19 141-154 (2010)
Here we use the "AUF2" model described in that paper.
(The hydrophobic beads face outwards.)
-- Memebrane Model: --
The DPPC lipid bilayer described in:
G. Brannigan, P.F. Philips, and F.L.H. Brown,
Physical Review E, Vol 72, 011915 (2005)
and:
M.C. Watson, E.S. Penev, P.M. Welch, and F.L.H. Brown
J. Chem. Phys. 135, 244701 (2011)
As in Watson(JCP 2011), rigid bond-length constraints
have been replaced by harmonic bonds.
A truncated version of this lipid (named "DLPC") has also been added.
Unlike the original "DPPC" molecule model, "DLPC" has not been carefully
parameterized to reproduce the correct behavior in a lipid bilayer mixture.
-------------
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

33
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/README_run.sh
View File

@ -1,33 +0,0 @@
# You would probably run lammps this way:
#
# lmp_ubuntu -i run.in.nvt
# The files "run.in.min", and "run.in.nvt" are LAMMPS input scripts which refer
# to the input scripts & data files you created earlier when you ran moltemplate
# system.in.init, system.in.settings, system.data
# -----------------------------------
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
"$LAMMPS_COMMAND" -i run.in.npt # pressure equilibration
"$LAMMPS_COMMAND" -i run.in.nvt # production run
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.npt # pressure equilibration
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt # production run

24
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/README_setup.sh
View File

@ -1,24 +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 ../
cp -f table_int.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/
cd ../

87
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/README_visualise.txt
View File

@ -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.)

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/4HelixOrig+Lipid2005_t=1290ps_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 83 KiB

2611
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/for_visualization/psf_file_created_by_topotools/system.psf
View File

File diff suppressed because it is too large Load Diff

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/membrane+protein_t=0ps_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 69 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/membrane+protein_t=0ps_no_pbc_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 50 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/membrane+protein_t=0ps_top_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 66 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/images/membrane+protein_t=5000ps_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 69 KiB

237
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/1beadProtSci2010.lt
View File

@ -1,237 +0,0 @@
# This file defines a 4-helix bundle coarse-grained protein model (AUF2) used in
# G. Bellesia, AI Jewett, and J-E Shea,
# Protein Science, Vol19 141-154 (2010)
#
# Strategy:
#
#1) First I'll define some building blocks (A16, B16, T3)
# which are helices, sheets and turns of a predetermined length)
#
#2) Then I'll copy and paste them together to build
# a 4-helix bundle (or a 4-strand beta-barrel).
# This approach is optional. If your protein has helices which are not
# identical, you should probably just include all 4 helices in a single
# "Data Atoms" section and don't try to subdivide the protein into pieces.)
1beadProtSci2010 { # <-- enclose definitions in a namespace for portability
# A16 is a coarse-grained alpha-helix containing 16 residues (one "atom" each)
A16 {
# AtomID MoleculeID AtomType Charge X Y Z
write('Data Atoms') {
$atom:a1 $mol:... @atom:../sL 0.0 -2.4 -2.4 0.0
$atom:a2 $mol:... @atom:../sL 0.0 2.4 -2.4 3.6
$atom:a3 $mol:... @atom:../sH 0.0 2.4 2.4 7.2
$atom:a4 $mol:... @atom:../sH 0.0 -2.4 2.4 10.8
$atom:a5 $mol:... @atom:../sL 0.0 -2.4 -2.4 14.4
$atom:a6 $mol:... @atom:../sL 0.0 2.4 -2.4 18.0
$atom:a7 $mol:... @atom:../sH 0.0 2.4 2.4 21.6
$atom:a8 $mol:... @atom:../sH 0.0 -2.4 2.4 25.2
$atom:a9 $mol:... @atom:../sL 0.0 -2.4 -2.4 28.8
$atom:a10 $mol:... @atom:../sL 0.0 2.4 -2.4 32.4
$atom:a11 $mol:... @atom:../sH 0.0 2.4 2.4 36.0
$atom:a12 $mol:... @atom:../sH 0.0 -2.4 2.4 39.6
$atom:a13 $mol:... @atom:../sL 0.0 -2.4 -2.4 43.2
$atom:a14 $mol:... @atom:../sL 0.0 2.4 -2.4 46.8
$atom:a15 $mol:... @atom:../sH 0.0 2.4 2.4 50.4
$atom:a16 $mol:... @atom:../sH 0.0 -2.4 2.4 54.0
}
write('Data Bonds') {
$bond:b1 @bond:../backbone $atom:a1 $atom:a2
$bond:b2 @bond:../backbone $atom:a2 $atom:a3
$bond:b3 @bond:../backbone $atom:a3 $atom:a4
$bond:b4 @bond:../backbone $atom:a4 $atom:a5
$bond:b5 @bond:../backbone $atom:a5 $atom:a6
$bond:b6 @bond:../backbone $atom:a6 $atom:a7
$bond:b7 @bond:../backbone $atom:a7 $atom:a8
$bond:b8 @bond:../backbone $atom:a8 $atom:a9
$bond:b9 @bond:../backbone $atom:a9 $atom:a10
$bond:b10 @bond:../backbone $atom:a10 $atom:a11
$bond:b11 @bond:../backbone $atom:a11 $atom:a12
$bond:b12 @bond:../backbone $atom:a12 $atom:a13
$bond:b13 @bond:../backbone $atom:a13 $atom:a14
$bond:b14 @bond:../backbone $atom:a14 $atom:a15
$bond:b15 @bond:../backbone $atom:a15 $atom:a16
}
} # A16
T3 { # T3 is a "turn" region consisting of 3 beads
# AtomID MoleculeID AtomType Charge X Y Z
write('Data Atoms') {
$atom:a1 $mol:... @atom:../tN 0.0 -4.8 0.0 0.0
$atom:a2 $mol:... @atom:../tN 0.0 0.0 3.3 -1.44
$atom:a3 $mol:... @atom:../tN 0.0 4.8 0.0 0.0
}
write('Data Bonds') {
$bond:b1 @bond:../backbone $atom:a1 $atom:a2
$bond:b2 @bond:../backbone $atom:a2 $atom:a3
}
} # T3
# ----- Now build a larger molecule using A16 and T3 -------
# Create a 4-Helix bundle.
# In this version, the hydrophobic beads are poing outward.
# I oriented them this way because I want to place this protein in a membrane.
# (There is another file in this directory containing alternate version
# of this same molecule with the hydrophobic beads pointing inward.)
4HelixInsideOut {
helix1 = new A16.rot(-225, 0,0,1).move(-5.70,-5.70,-32.4)
helix2 = new A16.rot(-135, 0,0,1).move( 5.70,-5.70,-28.8)
helix3 = new A16.rot( -45, 0,0,1).move( 5.70, 5.70,-25.2)
helix4 = new A16.rot( 45, 0,0,1).move(-5.70, 5.70,-21.6)
turn1 = new T3.rot(180,1,0,0).rot(-20,0,1,0).rot( 10,0,0,1).move(0.78,-4.2, 27.9)
turn2 = new T3.rot(-10,1,0,0).rot( 20,0,1,0).rot(-70,0,0,1).move(4.55, 2.4,-33.0)
turn3 = new T3.rot(180,1,0,0).rot(-20,0,1,0).rot(190,0,0,1).move(-0.78,4.2, 34.2)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:helix1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:helix2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:helix3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:helix2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:helix3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:helix4/a16
}
create_var { $mol } # molecule ID number shared by all atoms in this protein
} # 4HelixInsideOut
# -------------- Force-Field Parameters ------------
# Units and force-field styles for this protein model
# (These can be overridden later.)
write_once("In Init") {
units real
atom_style full
bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid class2
pair_style hybrid lj/charmm/coul/charmm/inter es4k4l maxmax 21.0 24.0
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 #(turn on "1-4" interactions)
}
# --- Distance Units ---
# In this version of the model, sigma (the bond-length
# and particle diameter) is rounded to 6.0 Angstroms.
#
# --- Energy & Temperature Units ---
# In this protein model, "epsilon" represents the free energy
# bonus for bringing two hydrophobic amino acids together.
# Here I choose to set epsilon to 1.806551818181818 kCal/mole.
# This value was chosen so that a temperature of 300 Kelvin lies at
# 0.33 epsilon, which is the unfolding temperature of the marginally stable
# "ASF1" protein model from the Bellesia et al 2010 paper.
# This choice insures that both the "ASF1" model from that paper,
# as well as the much more stable "AUF2" protein we use here (which
# unfolds at 0.42*eps) should definitely remain stable at 300 degrees Kelvin,
# in the bulk at least. (However it's not clear that these energy
# parameters will work well for a protein in membrane. Perhaps I'll
# run some tests and fine tune these parameters for this scenario.)
# 2-body (non-bonded) interactions:
#
# Uij(r) = 4*eps_ij * (K*(sig_ij/r)^12 + L*(sig_ij/r)^6)
#
# i j pairstylename eps sig K L
#
write_once("In Settings") {
# Override any pair_coeff and bond_coeff
# statements from the 1beadProteSci2010 model
# (Because the energies &distances in that model
# were in the wrong units.)
pair_coeff @atom:sH @atom:sH lj/charmm/coul/charmm/inter 1.8065518 6.0 1 -1
pair_coeff @atom:sL @atom:sL lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:tN @atom:tN lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
}
# The exact value of the bond_coeff does not matter too much as long as
# it is "stiff enough". Here I use a softer bond spring than the one
# used in the paper so that I can increase the time step.
# I also use a relatively soft spring to constrain the bond angles.
# bond_coeff bondType bondstylename k r0
write_once("In Settings") {
bond_coeff @bond:1beadProtSci2010/backbone harmonic 10.0 6.0
}
# angleType atomtypes1 2 3 bondtypes1 2
write_once("Data Angles By Type") {
@angle:backbone @atom:* @atom:* @atom:* @bond:* @bond:*
}
# angle_coeff angleType anglestylename k theta0
write_once("In Settings") {
angle_coeff @angle:backbone harmonic 100.0 105.0
}
# dihedralType atomtypes1 2 3 4 bondtypes1 2 3
write_once("Data Dihedrals By Type") {
# For a chain of sH and sL atoms, use the @dihedral:delta65_0
# parameters. (This corresponds to the "AUF2" model from the
# Bellesia et. al 2010 paper. Note: If you want the "AF" model,
# then change this to "delta60_0".)
@dihedral:delta65_0 @atom:s* @atom:s* @atom:s* @atom:s* * * *
# If "tN" (turn) atoms are present, use the @dihedral:turn parameters
@dihedral:turn @atom:tN @atom:* @atom:* @atom:* * * *
}
write_once("In Settings") {
dihedral_coeff @dihedral:delta60_0 class2 2.167862 120.0 0 0 2.167862 180.0
dihedral_coeff @dihedral:delta65_0 class2 2.167862 115.0 0 0 2.167862 180.0
dihedral_coeff @dihedral:turn class2 0 0 0 0 0.361310 180.0
# Note: 2.167862=1.2*epsilon and 0.361310=0.2*epsilon.
# All of the cross-terms (for the class2 force-field) are zero (see docs):
dihedral_coeff @dihedral:delta60_0 class2 mbt 0 0 0 0
dihedral_coeff @dihedral:delta60_0 class2 ebt 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:delta60_0 class2 at 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:delta60_0 class2 aat 0 0 0
dihedral_coeff @dihedral:delta60_0 class2 bb13 0 0 0
dihedral_coeff @dihedral:delta65_0 class2 mbt 0 0 0 0
dihedral_coeff @dihedral:delta65_0 class2 ebt 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:delta65_0 class2 at 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:delta65_0 class2 aat 0 0 0
dihedral_coeff @dihedral:delta65_0 class2 bb13 0 0 0
dihedral_coeff @dihedral:turn class2 mbt 0 0 0 0
dihedral_coeff @dihedral:turn class2 ebt 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:turn class2 at 0 0 0 0 0 0 0 0
dihedral_coeff @dihedral:turn class2 aat 0 0 0
dihedral_coeff @dihedral:turn class2 bb13 0 0 0
}
# --- Mass Units ---
# Typical amino acids weigh approximately 110.0 grams/mole. (Rounding down):
write_once("Data Masses") {
@atom:1beadProtSci2010/sH 100.0
@atom:1beadProtSci2010/sL 100.0
@atom:1beadProtSci2010/tN 100.0
}
} # 1beadProtSci2010 (namespace)

208
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/1beadProtSci2010_variations.lt
View File

@ -1,208 +0,0 @@
### THIS FILE IS OPTIONAL AND IS NOT NECESSARY. IN THIS FILE, I DEFINED SOME ##
### ADDITIONAL PROTEIN TYPES FROM THE PAPER THAT I DID NOT USE IN THIS EXAMPLE##
#
# This file defines a family of coarse-grained protein models used in:
# G. Bellesia, AI Jewett, and J-E Shea,
# Protein Science, Vol19 141-154 (2010)
#
# Strategy:
#
#1) First I'll define some building blocks (A16, B16, T3)
# which are helices, sheets and turns of a predetermined length)
import "1beadProtSci2010.lt"
# We defined A16 and T3 earlier in "1beadPRotSci2010.lt" Will define B16 below
#
#2) Then I'll copy and paste them together to build
# a 4-helix bundle or a 4-strand beta-barrel.
1beadProtSci2010 { #<-- Add new molecules to existing namespace defined earlier
# This way we don't have to start from scratch. We can
# use all the atom types and angle settings defined earlier
# B16 is a coarse-grained beta-strand containing 16 residues (one "atom" each)
B16 {
# AtomID MoleculeID AtomType Charge X Y Z
write('Data Atoms') {
$atom:a1 $mol:... @atom:../sL 0.0 -1.8 0.0 0.0
$atom:a2 $mol:... @atom:../sH 0.0 1.8 0.0 4.8
$atom:a3 $mol:... @atom:../sL 0.0 -1.8 0.0 9.6
$atom:a4 $mol:... @atom:../sH 0.0 1.8 0.0 14.4
$atom:a5 $mol:... @atom:../sL 0.0 -1.8 0.0 19.2
$atom:a6 $mol:... @atom:../sH 0.0 1.8 0.0 24.0
$atom:a7 $mol:... @atom:../sL 0.0 -1.8 0.0 28.8
$atom:a8 $mol:... @atom:../sH 0.0 1.8 0.0 33.6
$atom:a9 $mol:... @atom:../sL 0.0 -1.8 0.0 38.4
$atom:a10 $mol:... @atom:../sH 0.0 1.8 0.0 43.2
$atom:a11 $mol:... @atom:../sL 0.0 -1.8 0.0 48.0
$atom:a12 $mol:... @atom:../sH 0.0 1.8 0.0 52.8
$atom:a13 $mol:... @atom:../sL 0.0 -1.8 0.0 57.6
$atom:a14 $mol:... @atom:../sH 0.0 1.8 0.0 62.4
$atom:a15 $mol:... @atom:../sL 0.0 -1.8 0.0 67.2
$atom:a16 $mol:... @atom:../sH 0.0 1.8 0.0 72.0
}
write('Data Bonds') {
$bond:b1 @bond:../backbone $atom:a1 $atom:a2
$bond:b2 @bond:../backbone $atom:a2 $atom:a3
$bond:b3 @bond:../backbone $atom:a3 $atom:a4
$bond:b4 @bond:../backbone $atom:a4 $atom:a5
$bond:b5 @bond:../backbone $atom:a5 $atom:a6
$bond:b6 @bond:../backbone $atom:a6 $atom:a7
$bond:b7 @bond:../backbone $atom:a7 $atom:a8
$bond:b8 @bond:../backbone $atom:a8 $atom:a9
$bond:b9 @bond:../backbone $atom:a9 $atom:a10
$bond:b10 @bond:../backbone $atom:a10 $atom:a11
$bond:b11 @bond:../backbone $atom:a11 $atom:a12
$bond:b12 @bond:../backbone $atom:a12 $atom:a13
$bond:b13 @bond:../backbone $atom:a13 $atom:a14
$bond:b14 @bond:../backbone $atom:a14 $atom:a15
$bond:b15 @bond:../backbone $atom:a15 $atom:a16
}
} # B16
# ----- Now build larger molecules using B16 and T3 -------
4SheetBarrel {
sheet1 = new B16.rot( 45, 0,0,1).move(-4.762203156,-4.762203156, -36.0)
sheet2 = new B16.rot( 135, 0,0,1).move( 4.762203156,-4.762203156, -36.0)
sheet3 = new B16.rot( 225, 0,0,1).move( 4.762203156, 4.762203156, -36.0)
sheet4 = new B16.rot( 315, 0,0,1).move(-4.762203156, 4.762203156, -36.0)
turn1 = new T3.rot(180,1,0,0).rot( 0, 0,0,1).move( 0, -7.8, 39.6)
turn2 = new T3.rot( 0,1,0,0).rot(-90,0,0,1).move(4.2, 0.0,-41.4)
turn3 = new T3.rot(180,1,0,0).rot(-180,0,0,1).move( 0, 7.8, 39.6)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:sheet1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:sheet2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:sheet3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:sheet2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:sheet3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:sheet4/a16
}
create_var { $mol } # molecule ID number shared by all atoms in this protein
}
# Below I define several alternate conformations of the"4HelixBundleInsideOut"
# molecule I defined earlier in "1beadProtSci2010.lt". Same molecule however.
4HelixBundle {
helix1 = new A16.rot( -45, 0,0,1).move(-5.70,-5.70,-32.4)
helix2 = new A16.rot( 45, 0,0,1).move( 5.70,-5.70,-28.8)
helix3 = new A16.rot( 135, 0,0,1).move( 5.70, 5.70,-25.2)
helix4 = new A16.rot( 225, 0,0,1).move(-5.70, 5.70,-21.6)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:helix1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:helix2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:helix3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:helix2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:helix3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:helix4/a16
}
turn1 = new T3.rot(150,1,0,0).rot(-23,0,1,0).rot( 8,0,0,1).move(-3.6,-4.8,28.2)
turn2 = new T3.rot(-5,1,0,0).rot( 21,0,1,0).rot(-100,0,0,1).move(4.2,-0.66,-30.9)
turn3 = new T3.rot(150,1,0,0).rot(-23,0,1,0).rot(188,0,0,1).move(3.6,4.8,35.4)
create_var { $mol } # molecule ID number shared by all atoms in this protein
} # 4HelixBundle
# --- alternate conformations (same molecule) ----
# In the following version, the helices are oriented in a similar way,
# but they are separated a little further away from eachother.
4HelixBundleLoose {
helix1 = new A16.rot( -45, 0,0,1).move(-6.7347723,-6.7347723, -27.0)
helix2 = new A16.rot( 45, 0,0,1).move( 6.7347723,-6.7347723, -27.0)
helix3 = new A16.rot( 135, 0,0,1).move( 6.7347723, 6.7347723, -27.0)
helix4 = new A16.rot( 225, 0,0,1).move(-6.7347723, 6.7347723, -27.0)
turn1 = new T3.rot(180,1,0,0).rot(-17,0,0,1).move(-1.2,-4.2,32.4)
turn2 = new T3.rot( 0,1,0,0).rot(-100,0,0,1).move(4.2,-0.9,-28.8)
turn3 = new T3.rot(180,1,0,0).rot(163,0,0,1).move(1.2,4.2,32.4)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:helix1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:helix2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:helix3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:helix2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:helix3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:helix4/a16
}
create_var { $mol } # molecule ID number shared by all atoms in this protein
}
# In following version, the helices are oriented in a similar way,
# but they are separated a little further away from eachother.
4HelixInsideOutLoose {
helix1 = new A16.rot(-225, 0,0,1).move(-6.7347723,-6.7347723, -27.0)
helix2 = new A16.rot(-135, 0,0,1).move( 6.7347723,-6.7347723, -27.0)
helix3 = new A16.rot( -45, 0,0,1).move( 6.7347723, 6.7347723, -27.0)
helix4 = new A16.rot( 45, 0,0,1).move(-6.7347723, 6.7347723, -27.0)
turn1 = new T3.rot(180,1,0,0).rot( 10,0,0,1).move( 0.78,-4.2,28.8)
turn2 = new T3.rot( 70,1,0,0).rot(-70,0,0,1).move( 10.8,2.4,-28.2)
turn3 = new T3.rot(180,1,0,0).rot(190,0,0,1).move(-0.78,4.2,28.8)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:helix1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:helix2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:helix3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:helix2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:helix3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:helix4/a16
}
create_var { $mol } # molecule ID number shared by all atoms in this protein
} # 4HelixInsideOutLoose
# In the following version, the 4 helices are arranged next to each other,
# side-by-side, in a planar conformation (instead of a compact bundle).
4HelixPlanar {
helix1 = new A16.rot(-00, 0,0,1).move(0, 0, -27.0)
helix2 = new A16.rot( 00, 0,0,1).move(14.4, 0, -27.0)
helix3 = new A16.rot(-00, 0,0,1).move(28.8, 0, -27.0)
helix4 = new A16.rot( 00, 0,0,1).move(43.2, 0, -27.0)
turn1 = new T3.rot(180,1,0,0).rot( 0,0,0,1).move( 4.8, 0, 31.8)
turn2 = new T3.rot( 0,1,0,0).rot(180,0,0,1).move(19.2, 0,-31.8)
turn3 = new T3.rot(180,1,0,0).rot( 0,0,0,1).move(34.6, 0, 31.8)
write('Data Bonds') {
$bond:turn1a @bond:../backbone $atom:turn1/a1 $atom:helix1/a16
$bond:turn1b @bond:../backbone $atom:turn1/a3 $atom:helix2/a16
$bond:turn2a @bond:../backbone $atom:turn2/a1 $atom:helix3/a1
$bond:turn2b @bond:../backbone $atom:turn2/a3 $atom:helix2/a1
$bond:turn3a @bond:../backbone $atom:turn3/a1 $atom:helix3/a16
$bond:turn3b @bond:../backbone $atom:turn3/a3 $atom:helix4/a16
}
create_var { $mol } # molecule ID number shared by all atoms in this protein
} # 4HelixPlanar
} # 1beadProtSci2010 (namespace)

186
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/CGLipidBr2005.lt
View File

@ -1,186 +0,0 @@
# Note:
#
# This example may require additional features to be added to LAMMPS.
# If LAMMPS complains about an "Invalid pair_style", then copy the code
# in the "additional_lammps_code" directory into your LAMMPS "src" directory
# and recompile LAMMPS.
#
# -------- Description --------
#
# This example contains an implementation of the DPPC lipid bilayer described in
# G. Brannigan, P.F. Philips, and F.L.H. Brown,
# Physical Review E, Vol 72, 011915 (2005)
# and:
# M.C. Watson, E.S. Penev, P.M. Welch, and F.L.H. Brown
# J. Chem. Phys. 135, 244701 (2011)
#
# As in Watson(JCP 2011), rigid bond-length constraints have been replaced
# by harmonic bonds.
#
# A truncated version of this lipid (named "DLPC") has also been added.
# Unlike the original "DPPC" molecule model, "DLPC" has not been carefully
# parameterized to reproduce the correct behavior in a lipid bilayer mixture.
CGLipidBr2005 {
write_once("In Init") {
# -- Default styles for "CGLipidBr2005" --
units real
atom_style full
# (Hybrid force field styles were used for portability.)
bond_style hybrid harmonic
angle_style hybrid cosine/delta
dihedral_style none
improper_style none
pair_style hybrid table linear 1001 lj/charmm/coul/charmm/inter es4k4l 14.5 15
# Unfortunately, using both lj/charmm and "table" pair styles in the same
# simulation seems to be very inneficient. (The simulation is twice as slow
# as using only the "lj/charmm" pair styles for every pairwise interaction,
# ...and about 25% slower than using "table" for every pairwise interaction.
# However the lennard-jones pair styles support mixing, so we use them to
# make it easier to run these molecules with other molecules which don't use
# pair_table. I felt that portability was worth the extra 25% slow down.)
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 # turn off pairs if "less than 3 bonds"
# Optional: The "multi" style could be useful to efficiently simulate
# particles different sizes together. The tabular potential has range
# 22.5 Angstoms, but the Lennard-Jones forces have range of only
# 15.0 Angstroms.
neighbor 2.0 multi
communicate multi
}
# Particles shared by all lipid types
write_once("Data Masses") {
@atom:int 200.0
@atom:tail 200.0
@atom:head 200.0 #<- Default "head" type. We will override it later
}
# All force-field parameters appear in the "In Settings" section,
# which is later included in the LAMMPS input script.
# These force-field parameters are shared between all lipids by default.
# These can be customized and overridden for specific lipid types later.
# The "epsilon" parameter in their model is approximately 2.75 kJ/mole
# ( = 0.657265774378585 kCal/mole, using 1kCal=4.184kJ)
# The "sigma" parameter corresponds to 7.5 angstroms.
write_once("In Settings") {
# -- Default settings/parameters for "CGLipidBr2005" --
# (Hybrid bond & angle styles were used for portability.)
# As in Watson(JCP 2011), rigid bond-length constraints
# have been replaced by harmonic bonds.
# The k_theta parameter should lie in between 5*epsilon and 10*epsilon.
bond_coeff @bond:backbone harmonic 116.847 7.5 #<--2*5000*eps/sig^2
}
write_once("In Settings") {
#angle_coeff @angle:backbone cosine/delta 4.60086042 180 #<-- 7*eps
angle_coeff @angle:backbone cosine/delta 6.57265774 180 #<-- 10*eps
}
write_once("In Settings") {
# The interaction of "atom:int" with other "atom:int" atoms is given by
# epsilon*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2), shifted and cutoff at
# r=3*sigma. This was implemented using pair_style table.
# Unfortunately, mixing lj/charmm and "table" pair styles in the same
# simulation is very inneficient.
pair_coeff @atom:int @atom:int table table_int.dat INT
# The interaction of tail beads with eachother is given by the formula below
# and with other atoms ...using Lorenz-Berthelot and "repulsive wins" rules:
# epsilon*(0.4*(sigma/r)^12 - 1.0*(sigma/r)^6),
pair_coeff @atom:tail @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
pair_coeff @atom:int @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
# The interaction of head beads which all other beads is given by:
# epsilon*(0.4*(sigma/r)^12 - 0.0*(sigma/r)^6),
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
} # write_once("In Settings")
# Note: I divided epsilon by 4 to get "0.1643" because we are using the
# "es4k4l" coeffstyle, corresponding to U(r)=eps(4*K*(s/r)^12 + 4*L*(s/r)^6)
# The es4k4l coeffstyle is popular. Using this convention makes it easier
# to mix this coarse-grained lipid model with other molecular models.
DPPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 33.75 # custom "head" atom
$atom:i $mol:. @atom:../int 0.0 -1.00 0.00 26.25
$atom:t1 $mol:. @atom:../tail 0.0 1.00 0.00 18.75
$atom:t2 $mol:. @atom:../tail 0.0 -1.00 0.00 11.25
$atom:t3 $mol:. @atom:../tail 0.0 1.00 0.00 3.75
}
write("Data Bonds") {
$bond:b1 @bond:../backbone $atom:h $atom:i
$bond:b2 @bond:../backbone $atom:i $atom:t1
$bond:b3 @bond:../backbone $atom:t1 $atom:t2
$bond:b4 @bond:../backbone $atom:t2 $atom:t3
}
write("Data Angles") {
$angle:a1 @angle:../backbone $atom:h $atom:i $atom:t1
$angle:a2 @angle:../backbone $atom:i $atom:t1 $atom:t2
$angle:a3 @angle:../backbone $atom:t1 $atom:t2 $atom:t3
}
# Define properties of the local (lipid-specific) atom:head type atom:
write_once("Data Masses") {
@atom:head 200.0
}
write_once("In Settings") {
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
}
} #DPPC
DLPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 30.00 # custom "head" atom
$atom:i $mol:. @atom:../int 0.0 -1.00 0.00 22.50
$atom:t1 $mol:. @atom:../tail 0.0 1.00 0.00 15.00
$atom:t2 $mol:. @atom:../tail 0.0 -1.00 0.00 7.50
}
write("Data Bonds") {
$bond:b1 @bond:../backbone $atom:h $atom:i
$bond:b2 @bond:../backbone $atom:i $atom:t1
$bond:b3 @bond:../backbone $atom:t1 $atom:t2
}
write("Data Angles") {
$angle:a1 @angle:../backbone $atom:h $atom:i $atom:t1
$angle:a2 @angle:../backbone $atom:i $atom:t1 $atom:t2
}
# Define properties of the local (lipid-specific) atom:head type atom:
write_once("Data Masses") {
@atom:head 200.0
}
write_once("In Settings") {
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
}
} #DLPC
} # CGLipidBr2005

34
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/calc_table/calc_CGLipidTableINTvsINT.py
View File

@ -1,34 +0,0 @@
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between "INT" atoms
# in the lipid membrane model described in
# Brannigan et al, Phys Rev E, 72, 011915 (2005)
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
# However it is truncated at rc2 = 22.5 (shifted upwards to maintain continuity)
# The previous version included the repulsive core term
def U(r, eps, sigma):
return eps* (0.4*pow((sigma/r),12) - 3.0*sigma*sigma/(r*r))
def F(r, eps, sigma):
return eps*(12*0.4*pow((sigma/r),13)/sigma - 2*3.0*sigma*sigma/(r*r*r))
# We don't want to do that. Instead compute the core repulsion using a
# different pair_style and add the attractive term on top of it using the table.
# This way it the core repulsion acts as a default interaction with other atom
# types (using the new repulsive mixing rules).
epsilon = 2.75/4.184 # kCal/mole
sigma = 7.5
Rmin = 2.6
Rmax = 22.6
rcut = 22.5
N = 1001
for i in range(0,N):
r = Rmin + i*(Rmax-Rmin)/(N-1)
U_r = U(r, epsilon, sigma) - U(rcut, epsilon, sigma)
F_r = F(r, epsilon, sigma)
if r > rcut:
U_r = 0.0
F_r = 0.0
print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))

70
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/calc_table/version_charmm_cutoff/calc_table.py
View File

@ -1,70 +0,0 @@
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between "INT" atoms
# in the lipid membrane model described in
# Brannigan et al, Phys Rev E, 72, 011915 (2005)
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
# I realized later this is not what we want because although energy is conserved
# all enrgies are shifted with respect to energies used in the Brannigan paper
# (by 0.27 kCal/mole) and the later Watson JCP 2011 paper (by 0.224 kCal/mole).
# (So don't use this.)
# Calculate and print a
def S(r, rc1, rc2, derivative=False):
"""
Calculate the switching function S(r) which decays continuously
between 1 and 0 in the range from rc1 to rc2 (rc2>rc1):
S(r) = (rc2^2 - r^2)^2 * (rc2^2 + 2*r^2 - 3*rc1^2) / (rc2^2-rc1^2)^3
I'm using the same smoothing/switching cutoff function used by the CHARMM
force-fields. (I'm even using the same code to implement it, taken
from lammps charmm/coul/charmm pair style, rewritten in python.)
"""
assert(rc2>rc1)
rsq = r*r
rc1sq = rc1*rc1
rc2sq = rc2*rc2
denom_lj_inv = (1.0 / ((rc2sq-rc1sq)*
(rc2sq-rc1sq)*
(rc2sq-rc1sq)))
if rsq > rc2sq:
return 0.0
elif rsq < rc1sq:
if derivative:
return 0.0
else:
return 1.0
else:
rc2sq_minus_rsq = (rc2sq - rsq)
rc2sq_minus_rsq_sq = rc2sq_minus_rsq * rc2sq_minus_rsq
if derivative:
return (12.0 * rsq * rc2sq_minus_rsq * (rsq-rc1sq) * denom_lj_inv)
else:
return (rc2sq_minus_rsq_sq *
(rc2sq + 2.0*rsq - 3.0*rc1sq) * denom_lj_inv)
def U(r, eps, sigma):
return eps* (0.4*pow((sigma/r),12) - 3.0*sigma*sigma/(r*r))
def F(r, eps, sigma):
return eps*(12*0.4*pow((sigma/r),13)/sigma - 2*3.0*sigma*sigma/(r*r*r))
epsilon = 2.75/4.184 # kCal/mole
sigma = 7.5
Rmin = 2.6
Rmax = 22.6
Rc1 = 22.0
Rc2 = 22.5
N = 1001
for i in range(0,N):
r = Rmin + i*(Rmax-Rmin)/(N-1)
U_r = U(r, epsilon, sigma)
F_r = F(r, epsilon, sigma)
# Multiply U(r) & F(r) by the smoothing/switch function
U_r = U_r * S(r, Rc1, Rc2)
F_r = U_r * S(r, Rc1, Rc2, True) + F_r * S(r, Rc1, Rc2, False)
print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))

178
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/system.lt
View File

@ -1,178 +0,0 @@
# Description:
# This example shows how to put a protein (inclusion) in a
# lipid bilayer mixture composed of two different lipids (DPPC and DLPC).
# The DPPC lipid model is described here:
# G. Brannigan, P.F. Philips, and F.L.H. Brown,
# Physical Review E, Vol 72, 011915 (2005)
# The protein model is described here:
# G. Bellesia, AI Jewett, and J-E Shea,
# Protein Science, Vol19 141-154 (2010)
# The new DLPC model is a truncated version of DPPC,
# (Its behaviour has not been rigorously tested.)
# Note that 50%/50% mixtures of DPPC & DLPC are commonly used to
# build liposomes http://www.ncbi.nlm.nih.gov/pubmed/10620293
# Note:
# This example may require additional features to be added to LAMMPS.
# If LAMMPS complains about an "Invalid pair_style", then copy the code
# in the "additional_lammps_code" directory into your LAMMPS "src" directory
# and recompile LAMMPS.
import "CGLipidBr2005.lt"
using namespace CGLipidBr2005
# The "= new random" syntax chooses one of several molecules at random
lipids = new random([DPPC, DLPC], [0.5,0.5], 1234) #"1234"=random_seed
[13].move(7.5, 0, 0)
[15].move(3.75, 6.49519, 0) # <-- hexagonal lattice
[2].rot(180, 1, 0, 0) # <-- 2 monolayers
# Move all the lipds up to the center of the box
lipids[*][*][*].move(0,0,75.0)
# Although this patch of lipids is not square or rectangular, (it looks
# like a parallelogram), this is no longer the case after rectangular
# periodic boundary conditions are applied. (Check by visualising in VMD.)
write_once("Data Boundary") {
0 97.5 xlo xhi
0 97.42785792 ylo yhi
0 150.0 zlo zhi
}
# A note on geometry:
# We want to create a bilayer arranged in a hexagonal lattice consisting of
# 15 rows (each row is aligned with the x-axis)
# 13 columns (aligned at a 60 degree angle from the x axis)
# The lattice spacing is 7.5 Angstroms.
# When wrapped onto a rectangular box, the dimensions of the system are:
# 13 * 7.5 Angstroms in the X direction
# 15 * 7.5*sqrt(3)/2 Angstroms in the Y direction
# ------------------- protein inclusion ---------------------
import "1beadProtSci2010.lt"
using namespace 1beadProtSci2010
protein = new 4HelixInsideOut
protein.move(45.0, 25.98076211, 75.0)
# Delete a hole in the membrane to create space for the protein.
# (In the future moltemplate will be able to avoid occlusion automatically.)
delete lipids[4][2][*]
delete lipids[6][2][*]
delete lipids[3-6][3][*]
delete lipids[3-5][4][*]
delete lipids[2-5][5][*]
delete lipids[2][6][*]
delete lipids[4][6][*]
# -------- interactions between protein and lipids ----------
# Note: All atom types must include the full path (the name of
# the namespace which defined them as well as the atom type name).
# (This is because we are no longer inside that namespace.)
write_once("In Settings") {
# Interactions between the protein and lipid atoms are usually
# determined by mixing rules. (However this is not possible some
# for atoms, such as the "int" atoms in the lipid model which
# interact using -1/r^2 attraction.) Mixing rules do not make
# sense for these atoms so we must explicitly define their
# interaction with all other atoms.
# i j pairstylename eps sig K L
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 1.8065518 6.0 1 -1
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
# We want the interactions between hydrophobic residues and atoms in
# the interior of the lipd to be energetically similar to the attractive
# interactions between hydrophobic residues. (See 1beadProtSci2010.)
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 1.8065518 6.0 1 -1
# We also add an artificial attractive interaction between the
# turn residues of the protein and the lipid head groups in
# order to keep the protein upright. Probably not necessary.
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 6.0 1 -1
pair_coeff @atom:CGLipidBr2005/DLPC/head @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 6.0 1 -1
# All other interactions between proteins and lipids are steric.
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/DLPC/head @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
pair_coeff @atom:CGLipidBr2005/DLPC/head @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 1.8065518 6.0 1 0
# Add a weak attractive interaction between hydrophilic "sL" beads
# (Whose strength mimics the strength of interaction between tail beads
# in the lipid. This was absent from the original protein model.
# However without some kind of weak attraction between residues,
# the negative pressure in the interior of the bilayer membrane
# allways pulls the protein apart. Recall that in the membrane,
# the hydrophobic beads in the protein will face outwards towards the lipid
# tails leaving the hydrophilic amino acids of the protein in the interior.
# In reality, these polar groups form hydrogen bonds with each other.)
pair_coeff @atom:1beadProtSci2010/sL @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.3286 7.5 0.4 -1
# However these hydrophilic amino acids are not attrected to
# the bilayer interior.
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
}
# Finally, we must combine the two force-field styles which were used for
# the coarse-grained lipid and protein. To do that, we write one last time
# to the "In Init" section. When reading the "Init" section LAMMPS will
# read these commands last and this will override any earlier settings.
write_once("In Init") {
# -- These styles override earlier settings --
units real
atom_style full
# (Hybrid force field styles were used for portability.)
bond_style hybrid harmonic
angle_style hybrid cosine/delta harmonic
dihedral_style hybrid class2
improper_style none
pair_style hybrid table linear 1001 lj/charmm/coul/charmm/inter es4k4l 14.5 15
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 # turn off pairs if "less than 3 bonds"
}

1010
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/moltemplate_files/table_int.dat
View File

File diff suppressed because it is too large Load Diff

19
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/run.in.min
View File

@ -1,19 +0,0 @@
# -- Init section --
include system.in.init
# -- Atom definition section --
read_data system.data
# -- Settings Section --
include system.in.settings
# -- Run section --
dump 1 all custom 50 traj_min.lammpstrj id mol type x y z ix iy iz
minimize 1.0e-5 1.0e-7 500 2000
write_restart system_after_min.rst

60
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/run.in.npt
View File

@ -1,60 +0,0 @@
# -------- REQUIREMENTS: ---------
# 1) This example requires the CLASS2 package.
# http://lammps.sandia.gov/doc/Section_start.html#start_3
# 2) It also may require additional features and bug fixes for LAMMPS.
# So, after typing "make yes-class2" in to the shell, ...
# be sure to download and copy the "additional_lammps_code" from
# http://moltemplate.org (upper-left corner menu)
# 3) Unpack it
# 4) copy the .cpp and .h files to the src folding of your lammps installation.
# 5) Compile LAMMPS.
#
# If LAMMPS complains about an "Invalid pair_style", or "Invalid dihedral_style"
# then you made a mistake in the instructions above.
#
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
read_data system.data
# -- Settings Section --
include system.in.settings
# -- Run Section --
timestep 10.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 10000 traj_npt.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp pe etotal vol epair ebond eangle
thermo 1000 # time interval for printing out "thermo" data
fix fxlan all langevin 300.0 300.0 120 48279
fix fxnph all nph x 0 0 1000 y 0 0 1000
# Note: The temperature 300.0 K corresponds to 0.907033536873*epsilon
# (for the "epsilon" used by the coarse-grained lipid), and
# to 0.33*epsilon (for the "epsilon" used in the coarse-grained protein)
# Note: The langevin damping parameter "120" corresponds to
# the 0.12ps damping time used in Watson et. al JCP 2011.
# Note: We maintain the system system at constant (zero) tention
# using a barostat damping parameter Pdamp=1000 ("0 0 1000")
# optional (not sure if this helps):
# balance x uniform y uniform
#restart 1000000
run 2000000
write_restart system_after_npt.rst

40
tools/moltemplate/examples/CG_membrane_examples/membrane+protein/run.in.nvt
View File

@ -1,40 +0,0 @@
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
# read_data system.data <-- commenting out
# Use the pressure-equilibrated restart file instead:
read_restart system_after_npt.rst
# -- Settings Section --
include system.in.settings
# -- Run Section --
timestep 10.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 20000 traj_nvt.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp pe etotal vol epair ebond eangle
thermo 1000 # time interval for printing out "thermo" data
fix fxlan all langevin 300.0 300.0 120 48279
fix fxnve all nve
# Note: The energy scale "epsilon" = 2.75kJ/mole = 330.7485200981 Kelvin*kB.
# So a temperature of 300.0 Kelvin corresponds to 0.907033536873*epsilon.
# Note: The langevin damping parameter "120" corresponds to
# the 0.12ps damping time used in Watson et. al JCP 2011.
#restart 1000000
run 20000000
write_restart system_after_nvt.rst

33
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/README.TXT
View File

@ -1,33 +0,0 @@
Note:
This example may require additional features to be added to LAMMPS.
If LAMMPS complains about an "Invalid pair_style", then copy the code
in the "additional_lammps_code" directory into your LAMMPS "src" directory
and recompile LAMMPS.
----- Description --------
This example contains an implementation of the DPPC lipid bilayer described in:
G. Brannigan, P.F. Philips, and F.L.H. Brown,
Physical Review E, Vol 72, 011915 (2005)
and:
M.C. Watson, E.S. Penev, P.M. Welch, and F.L.H. Brown
J. Chem. Phys. 135, 244701 (2011)
As in Watson(JCP 2011), rigid bond-length constraints
have been replaced by harmonic bonds.
A truncated version of this lipid (named "DLPC") has also been added.
Unlike the original "DPPC" molecule model, "DLPC" has not been carefully
parameterized to reproduce the correct behavior in a lipid bilayer mixture.
-------------
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

33
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/README_run.sh
View File

@ -1,33 +0,0 @@
# You would probably run lammps this way:
#
# lmp_ubuntu -i run.in.nvt
# The files "run.in.min", and "run.in.nvt" are LAMMPS input scripts which refer
# to the input scripts & data files you created earlier when you ran moltemplate
# system.in.init, system.in.settings, system.data
# -----------------------------------
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
"$LAMMPS_COMMAND" -i run.in.npt # pressure equilibration
"$LAMMPS_COMMAND" -i run.in.nvt # production run
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.npt # pressure equilibration
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt # production run

24
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/README_setup.sh
View File

@ -1,24 +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 ../
cp -f table_int.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/
cd ../

87
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/README_visualise.txt
View File

@ -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.)

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/images/DLPC.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 4.8 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/images/DPPC+DLPC_bilayer32x37_t=0ps_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 57 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/images/DPPC+DLPC_bilayer32x37_t=0ps_no_pbc_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 36 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/images/DPPC+DLPC_bilayer32x37_t=500ps_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 61 KiB

BIN
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/images/DPPC.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 5.7 KiB

186
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/CGLipidBr2005.lt
View File

@ -1,186 +0,0 @@
# Note:
#
# This example may require additional features to be added to LAMMPS.
# If LAMMPS complains about an "Invalid pair_style", then copy the code
# in the "additional_lammps_code" directory into your LAMMPS "src" directory
# and recompile LAMMPS.
#
# -------- Description --------
#
# This example contains an implementation of the DPPC lipid bilayer described in
# G. Brannigan, P.F. Philips, and F.L.H. Brown,
# Physical Review E, Vol 72, 011915 (2005)
# and:
# M.C. Watson, E.S. Penev, P.M. Welch, and F.L.H. Brown
# J. Chem. Phys. 135, 244701 (2011)
#
# As in Watson(JCP 2011), rigid bond-length constraints have been replaced
# by harmonic bonds.
#
# A truncated version of this lipid (named "DLPC") has also been added.
# Unlike the original "DPPC" molecule model, "DLPC" has not been carefully
# parameterized to reproduce the correct behavior in a lipid bilayer mixture.
CGLipidBr2005 {
write_once("In Init") {
# -- Default styles for "CGLipidBr2005" --
units real
atom_style full
# (Hybrid force field styles were used for portability.)
bond_style hybrid harmonic
angle_style hybrid cosine/delta
dihedral_style none
improper_style none
pair_style hybrid table linear 1001 lj/charmm/coul/charmm/inter es4k4l 14.5 15
# Unfortunately, using both lj/charmm and "table" pair styles in the same
# simulation seems to be very inneficient. (The simulation is twice as slow
# as using only the "lj/charmm" pair styles for every pairwise interaction,
# ...and about 25% slower than using "table" for every pairwise interaction.
# However the lennard-jones pair styles support mixing, so we use them to
# make it easier to run these molecules with other molecules which don't use
# pair_table. I felt that portability was worth the extra 25% slow down.)
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 # turn off pairs if "less than 3 bonds"
# Optional: The "multi" style could be useful to efficiently simulate
# particles different sizes together. The tabular potential has range
# 22.5 Angstoms, but the Lennard-Jones forces have range of only
# 15.0 Angstroms.
neighbor 2.0 multi
communicate multi
}
# Particles shared by all lipid types
write_once("Data Masses") {
@atom:int 200.0
@atom:tail 200.0
@atom:head 200.0 #<- Default "head" type. We will override it later
}
# All force-field parameters appear in the "In Settings" section,
# which is later included in the LAMMPS input script.
# These force-field parameters are shared between all lipids by default.
# These can be customized and overridden for specific lipid types later.
# The "epsilon" parameter in their model is approximately 2.75 kJ/mole
# ( = 0.657265774378585 kCal/mole, using 1kCal=4.184kJ)
# The "sigma" parameter corresponds to 7.5 angstroms.
write_once("In Settings") {
# -- Default settings/parameters for "CGLipidBr2005" --
# (Hybrid bond & angle styles were used for portability.)
# As in Watson(JCP 2011), rigid bond-length constraints
# have been replaced by harmonic bonds.
# The k_theta parameter should lie in between 5*epsilon and 10*epsilon.
bond_coeff @bond:backbone harmonic 116.847 7.5 #<--2*5000*eps/sig^2
}
write_once("In Settings") {
#angle_coeff @angle:backbone cosine/delta 4.60086042 180 #<-- 7*eps
angle_coeff @angle:backbone cosine/delta 6.57265774 180 #<-- 10*eps
}
write_once("In Settings") {
# The interaction of "atom:int" with other "atom:int" atoms is given by
# epsilon*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2), shifted and cutoff at
# r=3*sigma. This was implemented using pair_style table.
# Unfortunately, mixing lj/charmm and "table" pair styles in the same
# simulation is very inneficient.
pair_coeff @atom:int @atom:int table table_int.dat INT
# The interaction of tail beads with eachother is given by the formula below
# and with other atoms ...using Lorenz-Berthelot and "repulsive wins" rules:
# epsilon*(0.4*(sigma/r)^12 - 1.0*(sigma/r)^6),
pair_coeff @atom:tail @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
pair_coeff @atom:int @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
# The interaction of head beads which all other beads is given by:
# epsilon*(0.4*(sigma/r)^12 - 0.0*(sigma/r)^6),
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
} # write_once("In Settings")
# Note: I divided epsilon by 4 to get "0.1643" because we are using the
# "es4k4l" coeffstyle, corresponding to U(r)=eps(4*K*(s/r)^12 + 4*L*(s/r)^6)
# The es4k4l coeffstyle is popular. Using this convention makes it easier
# to mix this coarse-grained lipid model with other molecular models.
DPPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 33.75 # custom "head" atom
$atom:i $mol:. @atom:../int 0.0 -1.00 0.00 26.25
$atom:t1 $mol:. @atom:../tail 0.0 1.00 0.00 18.75
$atom:t2 $mol:. @atom:../tail 0.0 -1.00 0.00 11.25
$atom:t3 $mol:. @atom:../tail 0.0 1.00 0.00 3.75
}
write("Data Bonds") {
$bond:b1 @bond:../backbone $atom:h $atom:i
$bond:b2 @bond:../backbone $atom:i $atom:t1
$bond:b3 @bond:../backbone $atom:t1 $atom:t2
$bond:b4 @bond:../backbone $atom:t2 $atom:t3
}
write("Data Angles") {
$angle:a1 @angle:../backbone $atom:h $atom:i $atom:t1
$angle:a2 @angle:../backbone $atom:i $atom:t1 $atom:t2
$angle:a3 @angle:../backbone $atom:t1 $atom:t2 $atom:t3
}
# Define properties of the local (lipid-specific) atom:head type atom:
write_once("Data Masses") {
@atom:head 200.0
}
write_once("In Settings") {
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
}
} #DPPC
DLPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 30.00 # custom "head" atom
$atom:i $mol:. @atom:../int 0.0 -1.00 0.00 22.50
$atom:t1 $mol:. @atom:../tail 0.0 1.00 0.00 15.00
$atom:t2 $mol:. @atom:../tail 0.0 -1.00 0.00 7.50
}
write("Data Bonds") {
$bond:b1 @bond:../backbone $atom:h $atom:i
$bond:b2 @bond:../backbone $atom:i $atom:t1
$bond:b3 @bond:../backbone $atom:t1 $atom:t2
}
write("Data Angles") {
$angle:a1 @angle:../backbone $atom:h $atom:i $atom:t1
$angle:a2 @angle:../backbone $atom:i $atom:t1 $atom:t2
}
# Define properties of the local (lipid-specific) atom:head type atom:
write_once("Data Masses") {
@atom:head 200.0
}
write_once("In Settings") {
pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
}
} #DLPC
} # CGLipidBr2005

34
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/generate_tables/calc_table.py
View File

@ -1,34 +0,0 @@
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between "INT" atoms
# in the lipid membrane model described in
# Brannigan et al, Phys Rev E, 72, 011915 (2005)
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
# However it is truncated at rc2 = 22.5 (shifted upwards to maintain continuity)
# The previous version included the repulsive core term
def U(r, eps, sigma):
return eps* (0.4*pow((sigma/r),12) - 3.0*sigma*sigma/(r*r))
def F(r, eps, sigma):
return eps*(12*0.4*pow((sigma/r),13)/sigma - 2*3.0*sigma*sigma/(r*r*r))
# We don't want to do that. Instead compute the core repulsion using a
# different pair_style and add the attractive term on top of it using the table.
# This way it the core repulsion acts as a default interaction with other atom
# types (using the new repulsive mixing rules).
epsilon = 2.75/4.184 # kCal/mole
sigma = 7.5
Rmin = 2.6
Rmax = 22.6
rcut = 22.5
N = 1001
for i in range(0,N):
r = Rmin + i*(Rmax-Rmin)/(N-1)
U_r = U(r, epsilon, sigma) - U(rcut, epsilon, sigma)
F_r = F(r, epsilon, sigma)
if r > rcut:
U_r = 0.0
F_r = 0.0
print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))

90
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/system.lt
View File

@ -1,90 +0,0 @@
# Description:
# This example shows how to build a lipid bilayer composed of a
# 50%-50% mixture of two different lipids (DPPC and DLPC).
# The DPPC model is described here:
# G. Brannigan, P.F. Philips, and F.L.H. Brown,
# Physical Review E, Vol 72, 011915 (2005)
# The new DLPC model is a truncated version of DPPC,
# (Its behaviour has not been rigorously tested.)
# Note: 50%/50% mixtures of DPPC & DLPC are commonly used to build liposomes
# http://www.ncbi.nlm.nih.gov/pubmed/10620293
import "CGLipidBr2005.lt"
using namespace CGLipidBr2005
lipids = new random([DPPC,DLPC], [0.5,0.5], 12345)
[32].move(7.5, 0, 0)
[37].move(3.75, 6.49519, 0)
[2].rot(180, 1, 0, 0)
# Move the lipds up to the center of the box
lipids[*][*][*].move(0,0,75.0)
# Although this patch of lipids is not square or rectangular, (it looks
# like a parallelogram), this is no longer the case after rectangular
# periodic boundary conditions are applied.
write_once("Data Boundary") {
0 240 xlo xhi
0 240.322 ylo yhi
0 150.0 zlo zhi
}
# File ends here. Only comments below.
# ------------------------------------
# ------------- COMMENTS: ------------
# ------------------------------------
#
# A note on geometry:
# We want to create a bilayer arranged in a hexagonal lattice consisting of
# 32 rows (each row is aligned with the x-axis)
# 37 columns (aligned at a 60 degree angle from the x axis)
# The lattice spacing is 8.0 Angstroms ( ~= 0.95*sigma*2^(1/6), where sigma=7.5)
# When wrapped onto a rectangular box, the dimensions of the system are:
# 32 * 7.5 Angstroms in the X direction
# 37 * 7.5*sqrt(3)/2 Angstroms in the Y direction
# ------------------------------------
#
# Below I show simple ways to create a lipid bilayer:
#
# 1) If you just want to make lipid bilayer out of DPPC,
# without specifying the location of each lipid, you could use this syntax:
# lipids = new DPPC [32][37][2] # 3-D array
# Later you can load in the coordinates of the lipds from a PDB file.
# Alternately you could also use a 1-dimensional array:
# lipids = new DPPC [2368] # 1-D array. Note: 2368 = 32 x 37 x 2
# It does not matter as long as the number of lipids is correct.
# Multidimensional arrays are only useful if you plan to apply independent
# coordinate transformations to each row and column and monolayer. See below:
#
# 2) Instead of loading a PDB file later, we can directly specify the location
# of each DPPC lipid in the LT file itself. For lipid bilayers, this is
# easy, because the bilayer structure resembles 2 planar lattices.
# We can use "move" commands to place each lipid, and the "rot" command
# to turn the lipids in one of the monolayers upside down.
#
# lipids = new DPPC [32].move(7.5, 0, 0)
# [37].move(3.75, 6.49519, 0)
# [2].rot(180, 1, 0, 0)
#
# 3) The example we use here is a lipid mixture of DPPC and DLPC, so we must
# replace "DPPC" in the command above with random([DPPC,DLPC],[0.5,0.5],12345)
# Here "0.5,0.5" are the probabilities for each molecule type, and "12345"
# is an optional random seed.
# lipids = new random([DPPC,DLPC], [0.5,0.5], 12345)
# [32].move(7.5, 0, 0)
# [37].move(3.75, 6.49519, 0)
# [2].rot(180, 1, 0, 0)
#

1010
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/table_int.dat
View File

File diff suppressed because it is too large Load Diff

70
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/version_charmm_cutoff/calc_table.py
View File

@ -1,70 +0,0 @@
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between "INT" atoms
# in the lipid membrane model described in
# Brannigan et al, Phys Rev E, 72, 011915 (2005)
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
# I realized later this is not what we want because although energy is conserved
# all enrgies are shifted with respect to energies used in the Brannigan paper
# (by 0.27 kCal/mole) and the later Watson JCP 2011 paper (by 0.224 kCal/mole).
# (So don't use this.)
# Calculate and print a
def S(r, rc1, rc2, derivative=False):
"""
Calculate the switching function S(r) which decays continuously
between 1 and 0 in the range from rc1 to rc2 (rc2>rc1):
S(r) = (rc2^2 - r^2)^2 * (rc2^2 + 2*r^2 - 3*rc1^2) / (rc2^2-rc1^2)^3
I'm using the same smoothing/switching cutoff function used by the CHARMM
force-fields. (I'm even using the same code to implement it, taken
from lammps charmm/coul/charmm pair style, rewritten in python.)
"""
assert(rc2>rc1)
rsq = r*r
rc1sq = rc1*rc1
rc2sq = rc2*rc2
denom_lj_inv = (1.0 / ((rc2sq-rc1sq)*
(rc2sq-rc1sq)*
(rc2sq-rc1sq)))
if rsq > rc2sq:
return 0.0
elif rsq < rc1sq:
if derivative:
return 0.0
else:
return 1.0
else:
rc2sq_minus_rsq = (rc2sq - rsq)
rc2sq_minus_rsq_sq = rc2sq_minus_rsq * rc2sq_minus_rsq
if derivative:
return (12.0 * rsq * rc2sq_minus_rsq * (rsq-rc1sq) * denom_lj_inv)
else:
return (rc2sq_minus_rsq_sq *
(rc2sq + 2.0*rsq - 3.0*rc1sq) * denom_lj_inv)
def U(r, eps, sigma):
return eps* (0.4*pow((sigma/r),12) - 3.0*sigma*sigma/(r*r))
def F(r, eps, sigma):
return eps*(12*0.4*pow((sigma/r),13)/sigma - 2*3.0*sigma*sigma/(r*r*r))
epsilon = 2.75/4.184 # kCal/mole
sigma = 7.5
Rmin = 2.6
Rmax = 22.6
Rc1 = 22.0
Rc2 = 22.5
N = 1001
for i in range(0,N):
r = Rmin + i*(Rmax-Rmin)/(N-1)
U_r = U(r, epsilon, sigma)
F_r = F(r, epsilon, sigma)
# Multiply U(r) & F(r) by the smoothing/switch function
U_r = U_r * S(r, Rc1, Rc2)
F_r = U_r * S(r, Rc1, Rc2, True) + F_r * S(r, Rc1, Rc2, False)
print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))

1020
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/moltemplate_files/version_charmm_cutoff/table_int.dat
View File

File diff suppressed because it is too large Load Diff

28
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/run.in.min
View File

@ -1,28 +0,0 @@
# -- Init section --
include system.in.init
# -- Atom definition section --
read_data system.data
# -- Settings Section --
include system.in.settings
# Optional: Make sure the pairwise energies look reasonable:
#pair_write 2 2 1001 r 2.6 16.0 test_tail-tail.dat t-t 0 0
#pair_write 2 3 1001 r 2.6 16.0 test_tail-head.dat t-h 0 0
#pair_write 1 2 1001 r 2.6 16.0 test_int-tail.dat i-t 0 0
#pair_write 1 1 2573 r 2.6 16.0 test_int-int.dat i-i 0 0
#pair_write 1 3 1001 r 2.6 16.0 test_int-head.dat i-h 0 0
#pair_write 3 3 1001 r 2.6 16.0 test_head-head.dat h-h 0 0
# -- Run section --
dump 1 all custom 50 traj_min.lammpstrj id mol type x y z ix iy iz
minimize 1.0e-5 1.0e-7 500 2000
write_restart system_after_min.rst

56
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/run.in.npt
View File

@ -1,56 +0,0 @@
# -------- REQUIREMENTS: ---------
# 1) This example may require additional features and bug fixes for LAMMPS.
# Be sure to download and copy the "additional_lammps_code" from
# http://moltemplate.org (upper-left corner menu)
# 2) Unpack it
# 3) copy the .cpp and .h files to the src folding of your lammps installation.
# 4) Compile LAMMPS.
#
# (If LAMMPS complains about an "Invalid pair_style"
# then you made a mistake in the instructions above.)
#
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
read_data system.data
# -- Settings Section --
include system.in.settings
# -- Run Section --
timestep 10.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 10000 traj_npt.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp pe etotal vol epair ebond eangle
thermo 1000 # time interval for printing out "thermo" data
fix fxlan all langevin 300.0 300.0 120 48279
fix fxnph all nph x 0 0 1000 y 0 0 1000
# Note: The temperature 300.0 K corresponds to 0.907033536873*epsilon
# for the "epsilon" used by the coarse-grained lipid.
# Note: The langevin damping parameter "120" corresponds to
# the 0.12ps damping time used in Watson et. al JCP 2011.
# Note: We maintain the system system at constant (zero) tention
# using a barostat damping parameter Pdamp=1000 ("0 0 1000")
# optional (not sure if this helps):
# balance x uniform y uniform
#restart 1000000
run 2000000
write_restart system_after_npt.rst

40
tools/moltemplate/examples/CG_membrane_examples/membrane_BranniganPRE2005/run.in.nvt
View File

@ -1,40 +0,0 @@
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
# read_data system.data #<-- commenting out
# Use the pressure-equilibrated restart file instead:
read_restart system_after_npt.rst
# -- Settings Section --
include system.in.settings
# -- Run Section --
timestep 6.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 5000 traj_nvt.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp pe etotal vol epair ebond eangle
thermo 1000 # time interval for printing out "thermo" data
fix fxlan all langevin 300.0 300.0 120 48279
fix fxnve all nve
# Note: The energy scale "epsilon" = 2.75kJ/mole = 330.7485200981 Kelvin*kB.
# So a temperature of 300.0 Kelvin corresponds to 0.907033536873*epsilon.
# Note: The langevin damping parameter "120" corresponds to
# the 0.12ps damping time used in Watson et. al JCP 2011.
#restart 1000000
run 1000000
write_restart system_after_nvt.rst

17
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/README.TXT
View File

@ -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

22
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/README_run.sh
View File

@ -1,22 +0,0 @@
# This is just an example.
#
# Note: The files "run.in.min", and "run.in.nvt" are LAMMPS input scripts which
# refer to the input scripts and data files you created earlier:
# system.in.init, system.in.settings, system.data
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
"$LAMMPS_COMMAND" -i run.in.nvt # production run
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt # production run

23
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/README_setup.sh
View File

@ -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 ../

86
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/README_visualise.txt
View File

@ -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.)

BIN
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/images/2bead_heteropolymer_t=0.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 74 KiB

BIN
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/images/2bead_heteropolymer_t=1000ps.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 53 KiB

1728
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/images/for_visualization/psf_file_created_by_topotools/system.psf
View File

File diff suppressed because it is too large Load Diff

173
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/2bead.lt
View File

@ -1,173 +0,0 @@
# In this example, we define two types of molecules: "H" and "P",
# both containing two atoms, named "CA" and "R".
#
# @R
# |
# @CA
#
# Eventually, we will connect multiple "H" and "P" molecules
# together to form a polymer, as shown below:
#
# @R @R
# | |
# _@CA_ _@CA_
# ... -.@CA-' `-@CA-' ` ...
# | |
# @R @R
#
# Suppose that the "H" and "P" molecules both use the same
# type of backbone atom ("CA"), but have their own custom "R"
# sidechain atoms with different properties:
# The "R" atoms belonging to "H" molecules are attracted to each other.
# The "R" atoms in "P" molecules are not.
#
# (Note: There is no reason the "H" and "P" molecules in this example need
# to contain the same number of atoms, or the same atom names.
# The point of this example is to illustrate how to share atom types
# and also the difference between local and global atom types.)
#
# By default, all counter variables are local. This means that whenever
# an atom type (or other counter variable) appears inside a molecule
# definition using the normal syntax "@atom:R", that atom type is a local
# property of the the "H" or "P" molecule in which it appears.
# Any properties assigned to the "R" atoms in either molecule are unique
# to that molecule.
#
# However in order to share "CA" atom types, we will override this behavior
# using the "@atom:../CA" syntax instead of "@atom:CA". This will define
# the "CA" atom in the outer (global) environment, and any properties
# (mass, radius, etc...) assigned to the "@atom:../CA" atom apply to
# "CA" atoms in both the the "H" and "P" molecules (and everywhere else).
#
# This might not be a desirable. "CA" is a popular name for carbon
# atoms in different types of molecules. We don't want to prevent other
# molecules from using this atom name. So we enclose the "CA" atom,
# (along with the definitions of the "H" and "P" molecules) within a
# namespace/environment object ("2bead"). This makes the definition
# of "H" and "P" more portable. Later on we can combine "H" and "P"
# molecules with other molecules without worrying whether they contain
# "CA" atoms with different properties.
#
# Note: In this example "2bead" is NOT a molecule. (Because it contains no
# "write("Data Atoms")" section of its own.) "2bead" is simply the
# name of an environment in which other molecules (H,P) are defined.
2bead {
# 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" --
units real
atom_style full
# (Hybrid force fields were not necessary but are used for portability.)
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
}
# Define H (the "hydrophobic" residue)
H {
# atom-id 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:R 0.0 0.000 4.4000 0.0000000
}
write("Data Bonds") {
$bond:CR @bond:../sidechain $atom:CA $atom:R
}
}
# Define P (the "polar" residue)
P {
# atom-id 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:R 0.0 0.000 4.4000 0.0000000
}
write("Data Bonds") {
$bond:CR @bond:../sidechain $atom:CA $atom:R
}
}
# (Note: Again, atom types, bond-types, (dihedral-types, any variable, etc)
# can be shared. The ".." in "@atom:../CA" tells moltemplate that
# atom type CA is defined in the parent's environment. (We are
# sharing the CA atom type between both the H and P residues.
# The same is true of the ".." in "@bond:../sidechain".
#
#
# 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.
# There are 3 atom types: the R sidechain (belonging to the H and P residues)
# and the CA sidechain atom (shared by both residues)
write_once("Data Masses") {
@atom:CA 13.0
@atom:H/R 50.0
@atom:P/R 50.0
}
# --------------------------------------------------------------------
# -- In this example, all force field parameters are stored in the --
# -- file named "In Settings". They can also go in sections like --
# -- "Data Pair Coeffs", "Data Bond Coeffs", "Data Angle Coeffs"... --
# --------------------------------------------------------------------
# 2-body (non-bonded) interactions:
#
# Uij(r) = 4*eps_ij * ( (sig_ij/r)^12 - (sig_ij/r)^6 )
#
# Hydrophobic side-chains 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:H/R @atom:H/R lj/cut 2.50 3.6
pair_coeff @atom:P/R @atom:P/R lj/cut 0.10 3.6
}
#
# (Interactions between different atom types use "arithmetic" mixing rules.)
# 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
}
} # 2bead

61
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/2bead_angles.lt
View File

@ -1,61 +0,0 @@
# Although there's no need to define angular interactions (because the "H"
# and "P" molecules only contains two atoms), we define the settings for angles
# or dihedrals which might be present later when we connect multiple "H" and "P"
# molecules together to build a polymer.
2bead {
# OPTIONAL: Enclose these angle settings within the 2bead {...} environment.
# We do this for the same reason we enclosed "H" and "P" in "2bead".
# This does not overwrite the definition of 2bead. Here we are only
# augmenting "2bead" to include definitions of the following angles
# 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:*
}
# dihedral-type AtomType1 AtomType2 AtomType3 AtomType4 bondType1 btyp2 btyp3
write_once("Data Dihedrals By Type") {
@dihedral:CCCC @atom:CA @atom:CA @atom:CA @atom:CA @bond:* @bond:* @bond:*
@dihedral:RCCR @atom:*/R @atom:CA @atom:CA @atom:*/R @bond:* @bond:* @bond:*
}
# 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
# "w" is the weight for 1-4 pair interactions, which we set to 0.0
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")
}

20
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/README.sh
View File

@ -1,20 +0,0 @@
# The version of "2bead.lt" in this directory is defines two types of
# molecules (named "2bead/H" and "2bead/P").
#
# However, there is another version of this file which is easier to understand.
# I recommend reading that file first.
# It is located at "simplified_version_one_residue/2bead.lt".
# It defines only one type of molecule (named "2bead")
# It is much simpler.
# ------
# Use these commands to generate the LAMMPS input script and data file
# (and other auxilliary files):
# run moltemplate
moltemplate.sh system.lt
# This will generate various files with names ending in *.in* and *.data.

51
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/peptide.lt
View File

@ -1,51 +0,0 @@
import "2bead.lt"
import "2bead_angles.lt"
Peptide {
# A polymer of alternating "H" and "P" beads:
res1 = new 2bead/P
res2 = new 2bead/P.rot(180.0, 1,0,0).move(3.2,0,0)
res3 = new 2bead/H.rot( 0.0, 1,0,0).move(6.4,0,0)
res4 = new 2bead/H.rot(180.0, 1,0,0).move(9.6,0,0)
res5 = new 2bead/H.rot( 0.0, 1,0,0).move(12.8,0,0)
res6 = new 2bead/H.rot(180.0, 1,0,0).move(16.0,0,0)
res7 = new 2bead/P.rot( 0.0, 1,0,0).move(19.2,0,0)
res8 = new 2bead/P.rot(180.0, 1,0,0).move(22.4,0,0)
res9 = new 2bead/P.rot( 0.0, 1,0,0).move(25.6,0,0)
res10 = new 2bead/H.rot(180.0, 1,0,0).move(28.8,0,0)
res11 = new 2bead/H.rot( 0.0, 1,0,0).move(32.0,0,0)
res12 = new 2bead/H.rot(180.0, 1,0,0).move(35.2,0,0)
res13 = new 2bead/P.rot( 0.0, 1,0,0).move(38.4,0,0)
res14 = new 2bead/P.rot(180.0, 1,0,0).move(41.6,0,0)
# Now, link the residues together this way:
write("Data Bonds") {
$bond:backbone1 @bond:2bead/backbone $atom:res1/CA $atom:res2/CA
$bond:backbone2 @bond:2bead/backbone $atom:res2/CA $atom:res3/CA
$bond:backbone3 @bond:2bead/backbone $atom:res3/CA $atom:res4/CA
$bond:backbone4 @bond:2bead/backbone $atom:res4/CA $atom:res5/CA
$bond:backbone5 @bond:2bead/backbone $atom:res5/CA $atom:res6/CA
$bond:backbone6 @bond:2bead/backbone $atom:res6/CA $atom:res7/CA
$bond:backbone7 @bond:2bead/backbone $atom:res7/CA $atom:res8/CA
$bond:backbone8 @bond:2bead/backbone $atom:res8/CA $atom:res9/CA
$bond:backbone9 @bond:2bead/backbone $atom:res9/CA $atom:res10/CA
$bond:backbone10 @bond:2bead/backbone $atom:res10/CA $atom:res11/CA
$bond:backbone11 @bond:2bead/backbone $atom:res11/CA $atom:res12/CA
$bond:backbone12 @bond:2bead/backbone $atom:res12/CA $atom:res13/CA
$bond:backbone13 @bond:2bead/backbone $atom:res13/CA $atom:res14/CA
}
create_var { $mol } # <--create a molecule ID number for this peptide
# This causes res1,res2,res3,...,res14 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.
} # Peptide
# Angle, dihedral and improper interactions will be generated
# according to the instructions in "2bead_angles.lt"

122
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/simplified_version_one_residue/2bead_simple.lt
View File

@ -1,122 +0,0 @@
# This file contains the definition of "2-bead", a simple molecule
# containing two atoms:
#
# @R
# |
# @CA
#
#
# Later on, we may connect them together to form a polymer:
#
#
# @R @R
# | |
# _@CA_ _@CA_
# ... -.@CA-' `-@CA-' ` ...
# | |
# @R @R
#
# THERE ARE TWO VERSIONS OF THIS FILE:
# In the simple version, shown here, there is only one type of side-chain "R".
# More complex examples have multiple types of molecules
# each with different side-chain properties.
2bead {
# atom-id mol-id atom-type charge x y z
write("Data Atoms") {
$atom:CA $mol:... @atom:CA 0.0 0.000 0.0000 0.0000000
$atom:R $mol:... @atom:R 0.0 0.000 3.4000 0.0000000
}
write_once("Data Masses") {
@atom:C 14.0
@atom:R 50.0
}
write("Data Bonds") {
$bond:CR @bond:sidechain $atom:CA $atom:R
}
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:*
}
write_once("Data Dihedrals By Type") {
@dihedral:CCCC @atom:CA @atom:CA @atom:CA @atom:CA @bond:* @bond:* @bond:*
@dihedral:RCCR @atom:R @atom:CA @atom:CA @atom:R @bond:* @bond:* @bond:*
}
write_once("In Settings") {
# 2-body (bonded) interaction parameters are listed by bondType:
#
# Ubond(r) = (k/2)*(r-0)^2
#
# The corresponding command is:
#
# bond_coeff bondType harmonic k r0
#
bond_coeff @bond:sidechain harmonic 30.0 3.4
bond_coeff @bond:backbone harmonic 30.0 3.7
# 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
#
# The corresponding command is:
#
# angle_coeff angleType harmonic k theta0
# (The theta0 parameter is in degrees, k is in kcal/mol/rad^2)
angle_coeff @angle:backbone harmonic 30.0 104
angle_coeff @angle:sidechain harmonic 30.0 127
# 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 corresponding command is
# dihedral_coeff dihedralType charmm K d w
#
# The d parameter is in degrees, K is in kcal/mol/rad^2.
# "w" is the weight for 1-4 pair interactions, which we set to 0.0.
# This should turn off pairwise interactions between 1-4 bonded atoms
# which are also involved in any of the dihedral interactions below:
#
dihedral_coeff @dihedral:CCCC charmm -0.5 1 180 0.0
dihedral_coeff @dihedral:RCCR charmm -1.5 1 -90 0.0
# 2-body (non-bonded) interactions are listed by atom type:
# interaction between atoms of the same type:
# epsilon signma
pair_coeff @atom:CA @atom:CA lj/charmm/coul/charmm 0.010 2.0
pair_coeff @atom:R @atom:R lj/charmm/coul/charmm 2.50 3.6
#
# interactions between different atom types use arithmetic mixing rules
} # write_once("In Settings")
write_once("In Init") {
# -- Default styles for "2bead" --
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
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 0.0
}
} # 2bead

25
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/moltemplate_files/system.lt
View File

@ -1,25 +0,0 @@
import "peptide.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 peptides in a rectangular grid
peptides = new Peptide [3].move(0, 0, 60.0)
[3].move(0, 60.0, 0)
[3].move(60.0, 0, 0)
# Now (for fun) shift some of the peptides
# in the x direction by a distance of 25.0
# Suppose we want to move the middle slice
# (which has constant Z). We do that this way:
peptides[1][*][*].move(25,0,0)
# more examples:
peptides[*][1][*].move(0,0,25)
peptides[*][*][1].move(0,25,0)

19
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/run.in.min
View File

@ -1,19 +0,0 @@
# -- Init section --
include system.in.init
# -- Atom definition section --
read_data system.data
# -- Settings Section --
include system.in.settings
# -- Run section --
dump 1 all custom 2500 traj_min.lammpstrj id mol type x y z ix iy iz
minimize 1.0e-5 1.0e-7 500 2000
write_restart system_after_min.rst

38
tools/moltemplate/examples/CG_misc/2bead_heteropolymer/run.in.nvt
View File

@ -1,38 +0,0 @@
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
# I you want to be careful, you can minimize the system first. (Try using
# "run.in.min" and uncomment the read_restart command in this file below.)
# read_restart system_after_min.rst
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
restart 1000000 restart_nvt
run 1000000
write_restart system_after_nvt.rst

36
tools/moltemplate/examples/CG_misc/random_heteropolymer/README.TXT
View File

@ -1,36 +0,0 @@
This is a very crude attempt to create a polymer melt
consisting of 150 copies of the same polymer.
THE FORCE FIELD PARAMETERS AND CONFORMATION FOR THIS SYSTEM ARE COMPLETELY WRONG
The purpose of this example is to demonstrate one way to create a
large number of randomly generated polymers, and to use
an NPT simulation to pack them all into a small box.
(Smaller than their initial outstretched length.)
Each polymer is a random heteropolymer of length 200 monomers.
This polymer is a polymoer of PVDF and PCTFE monomers
selected randomly in a 1:3 ratio.
Again, this is not a realistic simulation of PVDV or PCTFE polymers,
Furthermore, even after the simulation is done, the arrangement of the
polymers in the box is not characteristic of a truly random polymer melt.
---- A note on size ---
This is a large system with nearly 200000 atoms.
It takes several minutes to compile this example
and (currently requiring at least 4.0 Gb of memory).
(I'm working on reducing that requirement. -Andrew 2012-9-12)
-----------------------
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

19
tools/moltemplate/examples/CG_misc/random_heteropolymer/README_run.sh
View File

@ -1,19 +0,0 @@
# This is just an example.
#
# Note: The "run.in.nvt" file is a LAMMPS input script which attempts to read
# the input scripts and data files you created with moltemplate:
# system.in.init, system.in.settings, system.data
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.nvt
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt

23
tools/moltemplate/examples/CG_misc/random_heteropolymer/README_setup.sh
View File

@ -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 ../

87
tools/moltemplate/examples/CG_misc/random_heteropolymer/README_visualise.txt
View File

@ -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.)

BIN
tools/moltemplate/examples/CG_misc/random_heteropolymer/images/2bead.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 2.6 KiB

BIN
tools/moltemplate/examples/CG_misc/random_heteropolymer/images/3bead.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 4.3 KiB

277
tools/moltemplate/examples/CG_misc/random_heteropolymer/images/for_visualization/psf_file_created_by_topotools/system.psf
View File

@ -1,277 +0,0 @@
PSF
1 !NTITLE
REMARKS VMD generated structure x-plor psf file
116 !NATOM
1 1 3 3 0.000000 12.0000 0
2 1 4 4 0.000000 17.0000 0
3 1 4 4 0.000000 17.0000 0
4 1 3 3 0.000000 12.0000 0
5 1 4 4 0.000000 17.0000 0
6 1 4 4 0.000000 17.0000 0
7 1 1 1 0.000000 12.0000 0
8 1 2 2 0.000000 17.0000 0
9 1 1 1 0.000000 12.0000 0
10 1 2 2 0.000000 17.0000 0
11 1 1 1 0.000000 12.0000 0
12 1 2 2 0.000000 17.0000 0
13 1 3 3 0.000000 12.0000 0
14 1 4 4 0.000000 17.0000 0
15 1 4 4 0.000000 17.0000 0
16 1 1 1 0.000000 12.0000 0
17 1 2 2 0.000000 17.0000 0
18 1 1 1 0.000000 12.0000 0
19 1 2 2 0.000000 17.0000 0
20 1 1 1 0.000000 12.0000 0
21 1 2 2 0.000000 17.0000 0
22 1 1 1 0.000000 12.0000 0
23 1 2 2 0.000000 17.0000 0
24 1 1 1 0.000000 12.0000 0
25 1 2 2 0.000000 17.0000 0
26 1 1 1 0.000000 12.0000 0
27 1 2 2 0.000000 17.0000 0
28 1 1 1 0.000000 12.0000 0
29 1 2 2 0.000000 17.0000 0
30 1 3 3 0.000000 12.0000 0
31 1 4 4 0.000000 17.0000 0
32 1 4 4 0.000000 17.0000 0
33 1 1 1 0.000000 12.0000 0
34 1 2 2 0.000000 17.0000 0
35 1 3 3 0.000000 12.0000 0
36 1 4 4 0.000000 17.0000 0
37 1 4 4 0.000000 17.0000 0
38 1 1 1 0.000000 12.0000 0
39 1 2 2 0.000000 17.0000 0
40 1 3 3 0.000000 12.0000 0
41 1 4 4 0.000000 17.0000 0
42 1 4 4 0.000000 17.0000 0
43 1 3 3 0.000000 12.0000 0
44 1 4 4 0.000000 17.0000 0
45 1 4 4 0.000000 17.0000 0
46 1 1 1 0.000000 12.0000 0
47 1 2 2 0.000000 17.0000 0
48 1 1 1 0.000000 12.0000 0
49 1 2 2 0.000000 17.0000 0
50 1 3 3 0.000000 12.0000 0
51 1 4 4 0.000000 17.0000 0
52 1 4 4 0.000000 17.0000 0
53 1 1 1 0.000000 12.0000 0
54 1 2 2 0.000000 17.0000 0
55 1 1 1 0.000000 12.0000 0
56 1 2 2 0.000000 17.0000 0
57 1 3 3 0.000000 12.0000 0
58 1 4 4 0.000000 17.0000 0
59 1 4 4 0.000000 17.0000 0
60 1 1 1 0.000000 12.0000 0
61 1 2 2 0.000000 17.0000 0
62 1 3 3 0.000000 12.0000 0
63 1 4 4 0.000000 17.0000 0
64 1 4 4 0.000000 17.0000 0
65 1 1 1 0.000000 12.0000 0
66 1 2 2 0.000000 17.0000 0
67 1 1 1 0.000000 12.0000 0
68 1 2 2 0.000000 17.0000 0
69 1 1 1 0.000000 12.0000 0
70 1 2 2 0.000000 17.0000 0
71 1 3 3 0.000000 12.0000 0
72 1 4 4 0.000000 17.0000 0
73 1 4 4 0.000000 17.0000 0
74 1 3 3 0.000000 12.0000 0
75 1 4 4 0.000000 17.0000 0
76 1 4 4 0.000000 17.0000 0
77 1 1 1 0.000000 12.0000 0
78 1 2 2 0.000000 17.0000 0
79 1 1 1 0.000000 12.0000 0
80 1 2 2 0.000000 17.0000 0
81 1 3 3 0.000000 12.0000 0
82 1 4 4 0.000000 17.0000 0
83 1 4 4 0.000000 17.0000 0
84 1 1 1 0.000000 12.0000 0
85 1 2 2 0.000000 17.0000 0
86 1 1 1 0.000000 12.0000 0
87 1 2 2 0.000000 17.0000 0
88 1 1 1 0.000000 12.0000 0
89 1 2 2 0.000000 17.0000 0
90 1 1 1 0.000000 12.0000 0
91 1 2 2 0.000000 17.0000 0
92 1 3 3 0.000000 12.0000 0
93 1 4 4 0.000000 17.0000 0
94 1 4 4 0.000000 17.0000 0
95 1 1 1 0.000000 12.0000 0
96 1 2 2 0.000000 17.0000 0
97 1 1 1 0.000000 12.0000 0
98 1 2 2 0.000000 17.0000 0
99 1 1 1 0.000000 12.0000 0
100 1 2 2 0.000000 17.0000 0
101 1 3 3 0.000000 12.0000 0
102 1 4 4 0.000000 17.0000 0
103 1 4 4 0.000000 17.0000 0
104 1 1 1 0.000000 12.0000 0
105 1 2 2 0.000000 17.0000 0
106 1 1 1 0.000000 12.0000 0
107 1 2 2 0.000000 17.0000 0
108 1 3 3 0.000000 12.0000 0
109 1 4 4 0.000000 17.0000 0
110 1 4 4 0.000000 17.0000 0
111 1 1 1 0.000000 12.0000 0
112 1 2 2 0.000000 17.0000 0
113 1 1 1 0.000000 12.0000 0
114 1 2 2 0.000000 17.0000 0
115 1 1 1 0.000000 12.0000 0
116 1 2 2 0.000000 17.0000 0
115 !NBOND: bonds
1 2 1 3 1 4 4 5
4 6 4 7 7 8 7 9
9 10 9 11 11 12 11 13
13 14 13 15 13 16 16 17
16 18 18 19 18 20 20 21
20 22 22 23 22 24 24 25
24 26 26 27 26 28 28 29
28 30 30 31 30 32 30 33
33 34 33 35 35 36 35 37
35 38 38 39 38 40 40 41
40 42 40 43 43 44 43 45
43 46 46 47 46 48 48 49
48 50 50 51 50 52 50 53
53 54 53 55 55 56 55 57
57 58 57 59 57 60 60 61
60 62 62 63 62 64 62 65
65 66 65 67 67 68 67 69
69 70 69 71 71 72 71 73
71 74 74 75 74 76 74 77
77 78 77 79 79 80 79 81
81 82 81 83 81 84 84 85
84 86 86 87 86 88 88 89
88 90 90 91 90 92 92 93
92 94 92 95 95 96 95 97
97 98 97 99 99 100 99 101
101 102 101 103 101 104 104 105
104 106 106 107 106 108 108 109
108 110 108 111 111 112 111 113
113 114 113 115 115 116
193 !NTHETA: angles
1 4 7 38 40 43 40 43 46
69 71 74 71 74 77 30 33 35
35 38 40 57 60 62 4 7 9
9 11 13 13 16 18 26 28 30
43 46 48 46 48 50 50 53 55
53 55 57 62 65 67 67 69 71
74 77 79 77 79 81 81 84 86
88 90 92 92 95 97 97 99 101
101 104 106 104 106 108 108 111 113
7 9 11 16 18 20 18 20 22
20 22 24 22 24 26 24 26 28
65 67 69 84 86 88 86 88 90
95 97 99 111 113 115 11 13 16
28 30 33 33 35 38 48 50 53
55 57 60 60 62 65 79 81 84
90 92 95 99 101 104 106 108 111
5 4 7 6 4 7 11 13 14
11 13 15 14 13 16 15 13 16
28 30 31 28 30 32 31 30 33
32 30 33 33 35 36 33 35 37
36 35 38 37 35 38 38 40 41
38 40 42 44 43 46 45 43 46
48 50 51 48 50 52 51 50 53
52 50 53 55 57 58 55 57 59
58 57 60 59 57 60 60 62 63
60 62 64 63 62 65 64 62 65
69 71 72 69 71 73 75 74 77
76 74 77 79 81 82 79 81 83
82 81 84 83 81 84 90 92 93
90 92 94 93 92 95 94 92 95
99 101 102 99 101 103 102 101 104
103 101 104 106 108 109 106 108 110
109 108 111 110 108 111 7 9 10
8 7 9 9 11 12 10 9 11
16 18 19 17 16 18 18 20 21
19 18 20 20 22 23 21 20 22
22 24 25 23 22 24 24 26 27
25 24 26 26 28 29 27 26 28
46 48 49 47 46 48 53 55 56
54 53 55 65 67 68 66 65 67
67 69 70 68 67 69 77 79 80
78 77 79 84 86 87 85 84 86
86 88 89 87 86 88 88 90 91
89 88 90 95 97 98 96 95 97
97 99 100 98 97 99 104 106 107
105 104 106 111 113 114 112 111 113
113 115 116 114 113 115 1 4 5
1 4 6 2 1 4 3 1 4
40 43 44 40 43 45 41 40 43
42 40 43 71 74 75 71 74 76
72 71 74 73 71 74 4 7 8
12 11 13 13 16 17 29 28 30
30 33 34 34 33 35 35 38 39
39 38 40 43 46 47 49 48 50
50 53 54 56 55 57 57 60 61
61 60 62 62 65 66 70 69 71
74 77 78 80 79 81 81 84 85
91 90 92 92 95 96 100 99 101
101 104 105 107 106 108 108 111 112
2 1 3 5 4 6 14 13 15
31 30 32 36 35 37 41 40 42
44 43 45 51 50 52 58 57 59
63 62 64 72 71 73 75 74 76
82 81 83 93 92 94 102 101 103
109 108 110
47 !NPHI: dihedrals
38 40 43 46 69 71 74 77
1 4 7 9 40 43 46 48
67 69 71 74 71 74 77 79
35 38 40 43 30 33 35 38
28 30 33 35 33 35 38 40
57 60 62 65 55 57 60 62
4 7 9 11 7 9 11 13
13 16 18 20 24 26 28 30
62 65 67 69 65 67 69 71
81 84 86 88 86 88 90 92
92 95 97 99 95 97 99 101
108 111 113 115 16 18 20 22
18 20 22 24 20 22 24 26
22 24 26 28 84 86 88 90
9 11 13 16 11 13 16 18
26 28 30 33 46 48 50 53
53 55 57 60 48 50 53 55
60 62 65 67 77 79 81 84
79 81 84 86 88 90 92 95
90 92 95 97 97 99 101 104
104 106 108 111 99 101 104 106
106 108 111 113 43 46 48 50
50 53 55 57 74 77 79 81
101 104 106 108
0 !NIMPHI: impropers
0 !NDON: donors
0 !NACC: acceptors
0 !NNB
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0
1 0 !NGRP
0 0 0

BIN
tools/moltemplate/examples/CG_misc/random_heteropolymer/images/random_heteropolymer_t=0ps.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 30 KiB

BIN
tools/moltemplate/examples/CG_misc/random_heteropolymer/images/random_heteropolymer_t=1ns.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 29 KiB

123
tools/moltemplate/examples/CG_misc/random_heteropolymer/moltemplate_files/monomers.lt
View File

@ -1,123 +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. --
# ----------------------------------------------------------------------
MonomerTypes {
2bead {
# atom-id mol-id atom-type charge x y z
write("Data Atoms") {
$atom:CA $mol:... @atom:CA 0.0 0.000 1.000 0.0000
$atom:R $mol:... @atom:R 0.0 0.000 4.400 0.0000
}
# bond-id bond-type atom-id1 atom-id2
write("Data Bonds") {
$bond:CR @bond:../sidechain $atom:CA $atom:R
}
# atom-type mass
write_once("Data Masses") {
@atom:CA 12.0
@atom:R 17.0
}
# atom-type atom-type epsilon sigma
write_once("In Settings") {
pair_coeff @atom:CA @atom:CA 0.10 2.0
pair_coeff @atom:R @atom:R 0.50 3.0
}
} # 2bead
3bead {
# atom-id mol-id atom-type charge x y z
write("Data Atoms") {
$atom:CA $mol:... @atom:CA 0.0 0.000 1.000 0.000
$atom:R1 $mol:... @atom:R 0.0 0.000 2.700 2.950
$atom:R2 $mol:... @atom:R 0.0 0.000 2.700 -2.950
}
# 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
}
# atom-type mass
write_once("Data Masses") {
@atom:CA 12.0
@atom:R 17.0
}
# atom-type atom-type epsilon sigma
write_once("In Settings") {
pair_coeff @atom:CA @atom:CA 0.10 2.0
pair_coeff @atom:R @atom:R 0.50 3.0
}
} # 3bead
# bond-type k r0
write_once("Data Bond Coeffs") {
@bond:sidechain 20.0 3.4
@bond:bb 20.0 3.4 # "bb" shorthand for "backbone"
}
# Although there's no need to define angular interactions (because this
# "molecule" only contains two atoms), we define the settings for angles
# or dihedrals which might be present later when we build a polymer.
# angle-type k theta0
write_once("Data Angle Coeffs") {
@angle:backbone 10.00 160
@angle:sidechain 10.00 120
@angle:RCR 10.00 120
}
# dihedral-type K1 K2 K3 K4
write_once("Data Dihedral Coeffs") {
@dihedral:backbn 0.10 -0.271016 3.145034 0.0
}
# 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:*
}
# 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 * * *
}
} # MonomerTypes

82
tools/moltemplate/examples/CG_misc/random_heteropolymer/moltemplate_files/polymer.lt
View File

@ -1,82 +0,0 @@
import "monomers.lt"
using namespace MonomerTypes
RandomHeteropolymer {
# Make a chain of monomers with random composition:
monomers = new random([2bead,3bead],
[0.6,0.4],
123456) # <-- optional random seed
[50].rot(180,1,0,0).move(2.95, 0, 0)
# Now, link the monomers together this way:
write("Data Bonds") {
$bond:bb1 @bond:MonomerTypes/bb $atom:monomers[0]/CA $atom:monomers[1]/CA
$bond:bb2 @bond:MonomerTypes/bb $atom:monomers[1]/CA $atom:monomers[2]/CA
$bond:bb3 @bond:MonomerTypes/bb $atom:monomers[2]/CA $atom:monomers[3]/CA
$bond:bb4 @bond:MonomerTypes/bb $atom:monomers[3]/CA $atom:monomers[4]/CA
$bond:bb5 @bond:MonomerTypes/bb $atom:monomers[4]/CA $atom:monomers[5]/CA
$bond:bb6 @bond:MonomerTypes/bb $atom:monomers[5]/CA $atom:monomers[6]/CA
$bond:bb7 @bond:MonomerTypes/bb $atom:monomers[6]/CA $atom:monomers[7]/CA
$bond:bb8 @bond:MonomerTypes/bb $atom:monomers[7]/CA $atom:monomers[8]/CA
$bond:bb9 @bond:MonomerTypes/bb $atom:monomers[8]/CA $atom:monomers[9]/CA
$bond:bb10 @bond:MonomerTypes/bb $atom:monomers[9]/CA $atom:monomers[10]/CA
$bond:bb11 @bond:MonomerTypes/bb $atom:monomers[10]/CA $atom:monomers[11]/CA
$bond:bb12 @bond:MonomerTypes/bb $atom:monomers[11]/CA $atom:monomers[12]/CA
$bond:bb13 @bond:MonomerTypes/bb $atom:monomers[12]/CA $atom:monomers[13]/CA
$bond:bb14 @bond:MonomerTypes/bb $atom:monomers[13]/CA $atom:monomers[14]/CA
$bond:bb15 @bond:MonomerTypes/bb $atom:monomers[14]/CA $atom:monomers[15]/CA
$bond:bb16 @bond:MonomerTypes/bb $atom:monomers[15]/CA $atom:monomers[16]/CA
$bond:bb17 @bond:MonomerTypes/bb $atom:monomers[16]/CA $atom:monomers[17]/CA
$bond:bb18 @bond:MonomerTypes/bb $atom:monomers[17]/CA $atom:monomers[18]/CA
$bond:bb19 @bond:MonomerTypes/bb $atom:monomers[18]/CA $atom:monomers[19]/CA
$bond:bb20 @bond:MonomerTypes/bb $atom:monomers[19]/CA $atom:monomers[20]/CA
$bond:bb21 @bond:MonomerTypes/bb $atom:monomers[20]/CA $atom:monomers[21]/CA
$bond:bb22 @bond:MonomerTypes/bb $atom:monomers[21]/CA $atom:monomers[22]/CA
$bond:bb23 @bond:MonomerTypes/bb $atom:monomers[22]/CA $atom:monomers[23]/CA
$bond:bb24 @bond:MonomerTypes/bb $atom:monomers[23]/CA $atom:monomers[24]/CA
$bond:bb25 @bond:MonomerTypes/bb $atom:monomers[24]/CA $atom:monomers[25]/CA
$bond:bb26 @bond:MonomerTypes/bb $atom:monomers[25]/CA $atom:monomers[26]/CA
$bond:bb27 @bond:MonomerTypes/bb $atom:monomers[26]/CA $atom:monomers[27]/CA
$bond:bb28 @bond:MonomerTypes/bb $atom:monomers[27]/CA $atom:monomers[28]/CA
$bond:bb29 @bond:MonomerTypes/bb $atom:monomers[28]/CA $atom:monomers[29]/CA
$bond:bb30 @bond:MonomerTypes/bb $atom:monomers[29]/CA $atom:monomers[30]/CA
$bond:bb31 @bond:MonomerTypes/bb $atom:monomers[30]/CA $atom:monomers[31]/CA
$bond:bb32 @bond:MonomerTypes/bb $atom:monomers[31]/CA $atom:monomers[32]/CA
$bond:bb33 @bond:MonomerTypes/bb $atom:monomers[32]/CA $atom:monomers[33]/CA
$bond:bb34 @bond:MonomerTypes/bb $atom:monomers[33]/CA $atom:monomers[34]/CA
$bond:bb35 @bond:MonomerTypes/bb $atom:monomers[34]/CA $atom:monomers[35]/CA
$bond:bb36 @bond:MonomerTypes/bb $atom:monomers[35]/CA $atom:monomers[36]/CA
$bond:bb37 @bond:MonomerTypes/bb $atom:monomers[36]/CA $atom:monomers[37]/CA
$bond:bb38 @bond:MonomerTypes/bb $atom:monomers[37]/CA $atom:monomers[38]/CA
$bond:bb39 @bond:MonomerTypes/bb $atom:monomers[38]/CA $atom:monomers[39]/CA
$bond:bb40 @bond:MonomerTypes/bb $atom:monomers[39]/CA $atom:monomers[40]/CA
$bond:bb41 @bond:MonomerTypes/bb $atom:monomers[40]/CA $atom:monomers[41]/CA
$bond:bb42 @bond:MonomerTypes/bb $atom:monomers[41]/CA $atom:monomers[42]/CA
$bond:bb43 @bond:MonomerTypes/bb $atom:monomers[42]/CA $atom:monomers[43]/CA
$bond:bb44 @bond:MonomerTypes/bb $atom:monomers[43]/CA $atom:monomers[44]/CA
$bond:bb45 @bond:MonomerTypes/bb $atom:monomers[44]/CA $atom:monomers[45]/CA
$bond:bb46 @bond:MonomerTypes/bb $atom:monomers[45]/CA $atom:monomers[46]/CA
$bond:bb47 @bond:MonomerTypes/bb $atom:monomers[46]/CA $atom:monomers[47]/CA
$bond:bb48 @bond:MonomerTypes/bb $atom:monomers[47]/CA $atom:monomers[48]/CA
$bond:bb49 @bond:MonomerTypes/bb $atom:monomers[48]/CA $atom:monomers[49]/CA
}
# These lines of moltemplate script above were generated in python:
# for i in range(0,49):
# print(' $bond:bb'+str(i+1)+' @bond:MonomerTypes/bb $atom:monomers['
# +str(i)+']/CA $atom:monomers['+str(i+1)+']/CA')
# Finally, create a molecule ID number for this large polymer object
create_var { $mol }
} # RandomHeteropolymer
# Angle, dihedral and improper interactions will be generated
# automatically according to the instructions in "monomers.lt"

33
tools/moltemplate/examples/CG_misc/random_heteropolymer/moltemplate_files/system.lt
View File

@ -1,33 +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 150.0 xlo xhi
0.0 150.0 ylo yhi
0.0 150.0 zlo zhi
}
import "polymer.lt"
polymer = new RandomHeteropolymer

29
tools/moltemplate/examples/CG_misc/random_heteropolymer/run.in.nvt
View File

@ -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_restart system_after_nvt.rst

23
tools/moltemplate/examples/CG_misc/translocation/README.TXT
View File

@ -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

19
tools/moltemplate/examples/CG_misc/translocation/README_run.sh
View File

@ -1,19 +0,0 @@
# This is just an example.
#
# Note: The "run.in.nvt" file is a LAMMPS input script which attempts to read
# the input scripts and data files you created with moltemplate:
# system.in.init, system.in.settings, system.data
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.nvt
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt

23
tools/moltemplate/examples/CG_misc/translocation/README_setup.sh
View File

@ -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 ../

87
tools/moltemplate/examples/CG_misc/translocation/README_visualise.txt
View File

@ -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.)

1995
tools/moltemplate/examples/CG_misc/translocation/images/for_visualization/psf_file_created_by_topotools/system.psf
View File

File diff suppressed because it is too large Load Diff

BIN
tools/moltemplate/examples/CG_misc/translocation/images/walls+solvent+polymer_t=0.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 103 KiB

105
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/monomer.lt
View File

@ -1,105 +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
}
# bond-type k r0
write_once("Data Bond Coeffs") {
@bond:sidechain 30.0 1.2
@bond:bb 30.0 2.0 # "bb" shorthand for "backbone"
}
# Although there's no need to define angular interactions (because this
# "molecule" only contains two atoms), we define the settings for angles
# or dihedrals which might be present later when we build a polymer.
# angle-type k theta0
write_once("Data Angle Coeffs") {
@angle:backbone 50.00 160
@angle:sidechain 50.00 120
@angle:RCR 50.00 120
}
# dihedral-type K1 K2 K3 K4
write_once("Data Dihedral Coeffs") {
@dihedral:backbn 1.411036 -0.271016 3.145034 0.0
}
# 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:*
}
# 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:*
}
} # 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).
# -------------------------------------------------------------------------

35
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/polymer.lt
View File

@ -1,35 +0,0 @@
import "monomer.lt"
Polymer {
# 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
}
create_var { $mol } # Create a molecule ID number for this polymer
# This causes monomer[0], monomer[1], ... 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.
} # Polymer
# Angle, dihedral and improper interactions will be generated
# automatically according to the instructions in "monomer.lt"

23
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/solvent.lt
View File

@ -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]

22
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/solvent_single.lt
View File

@ -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)
}
}

53
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/system.lt
View File

@ -1,53 +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 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 bonded interactions. Delete them too.
# Note: In both cases the two molecule fragments will keep the same mol counter.

21
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/wall_single.lt
View File

@ -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)
}
}

23
tools/moltemplate/examples/CG_misc/translocation/moltemplate_files/walls.lt
View File

@ -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]

61
tools/moltemplate/examples/CG_misc/translocation/run.in.npt
View File

@ -1,61 +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.)
# This next line is somewhat controversial. Feel free to delete this next line
fix Ffreezestuff groupB rigid single force * off off off torque * off off off
#(Neither Trung or Steve Plimpton use fix rigid for immobilizing objects, but
# I noticed that at NPT, it does a better job of maintaining the correct volume)
# IMPORTANT for NPT: You must use "neigh_modify" to turn off calculation of the
# forces between immobilized atoms.
neigh_modify exclude group groupB groupB
# ------------------- NPT -----------------------
# 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
compute pressMobile all pressure tempMobile
thermo_style custom step c_tempMobile c_pressMobile temp press vol
# 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
fix_modify fMoveStuff temp tempMobile
run 60000
write_restart system_after_npt.rst

44
tools/moltemplate/examples/CG_misc/translocation/run.in.nvt
View File

@ -1,44 +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 100 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 100 # 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
compute pressMobile all pressure tempMobile
# Integrate the equations of motion:
fix fMoveStuff mobile nvt temp 300.0 300.0 100.0
fix_modify fMoveStuff temp tempMobile
run 200000
write_restart system_after_nvt.rst

33
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/README.TXT
View File

@ -1,33 +0,0 @@
# This directory contains examples of how to run a short simulation of a
# coarse-grained protein-like polymer, folding in the presence and absence of
# a chaperone (modeled as an attractive or repulsie spherical shell).
#
# The protein models and the chaperone models are described and used here:
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
# (http://www.pnas.org/content/101/36/13192)
# ...and also here:
# AI Jewett and J-E Shea, J. Mol. Biol, Vol 363(5), (2006)
#
# (In the "frustrated+minichaperone" directory, the protein is
# placed outside the chaperone sphere, as opposed to inside.)
#
# -------- REQUIREMENTS: ---------
# 1) These examples require the "USER-MISC" package. (Use "make yes-USER-MISC")
# http://lammps.sandia.gov/doc/Section_start.html#start_3
# 2) They also may require additional features and bug fixes for LAMMPS.
# be sure to download and copy the "additional_lammps_code" from
# http://moltemplate.org (upper-left corner menu)
# 3) Unpack it
# 4) copy the .cpp and .h files to the src folding of your lammps installation.
# 5) Compile LAMMPS.
#
-------------
Instructions on how to build LAMMPS input files and
run a short simulation are provided in other README files in each directory.
step 1)
README_setup.sh
step2)
README_run.sh

32
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/README.TXT
View File

@ -1,32 +0,0 @@
# This directory demonstrates how to run a long simulation of
# the "frustrated" coarse-grained protein confined in a frustrated
# coarse-grained chaperonin (R=6, h=0.475) as described in:
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
# (http://www.pnas.org/content/101/36/13192)
#
# Note: If you want to use a "hydrophilic" chaperone (with h=0.0
# instead of h=0.475), then replace the word "CHAP_INTERIOR_H0.475"
# (at the end of "system.lt") with "CHAP_INTERIOR_H0"
#
# Because this process takes a long time (even with the help of the chaperone)
# I save the data relatively infrequently.
#
# -------- REQUIREMENTS: ---------
# 1) This example requires the "USER-MISC" package. (Use "make yes-USER-MISC")
# http://lammps.sandia.gov/doc/Section_start.html#start_3
# 2) It also may require additional features and bug fixes for LAMMPS.
# be sure to download and copy the "additional_lammps_code" from
# http://moltemplate.org (upper-left corner menu)
# 3) Unpack it
# 4) copy the .cpp and .h files to the src folding of your lammps installation.
# 5) Compile LAMMPS.
-------------
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

31
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/README_run.sh
View File

@ -1,31 +0,0 @@
# You would probably run lammps this way:
#
# lmp_ubuntu -i run.in.nvt
# The files "run.in.min", and "run.in.nvt" are LAMMPS input scripts which refer
# to the input scripts & data files you created earlier when you ran moltemplate
# system.in.init, system.in.settings, system.data
# -----------------------------------
LAMMPS_COMMAND="lmp_ubuntu"
# Here "$LAMMPS_BINARY" is the name of the command you use to invoke lammps
# (such as lmp_linux, lmp_g++, lmp_mac, lmp_cygwin, etc...) Change if necessary.
# Run lammps using the following 3 commands:
"$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
"$LAMMPS_COMMAND" -i run.in.nvt # production run
# Alternately, if you have MPI installed, try something like this:
#NUMPROCS=4
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.min # minimize (OPTIONAL)
#mpirun -np $NUMPROCS "$LAMMPS_COMMAND" -i run.in.nvt # production run

24
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/README_setup.sh
View File

@ -1,24 +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 -overlay-dihdedrals 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 ../
cp -r table*.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/
cd ../

87
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/README_visualise.txt
View File

@ -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.)

86
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/images/for_visualization/psf_file_created_by_topotools/system.psf
View File

@ -1,86 +0,0 @@
PSF
1 !NTITLE
REMARKS VMD generated structure x-plor psf file
28 !NATOM
1 1 2 2 0.000000 1.0000 0
2 1 1 1 0.000000 1.0000 0
3 1 2 2 0.000000 1.0000 0
4 1 1 1 0.000000 1.0000 0
5 1 2 2 0.000000 1.0000 0
6 1 1 1 0.000000 1.0000 0
7 1 3 3 0.000000 1.0000 0
8 1 3 3 0.000000 1.0000 0
9 1 1 1 0.000000 1.0000 0
10 1 2 2 0.000000 1.0000 0
11 1 1 1 0.000000 1.0000 0
12 1 2 2 0.000000 1.0000 0
13 1 1 1 0.000000 1.0000 0
14 1 2 2 0.000000 1.0000 0
15 1 3 3 0.000000 1.0000 0
16 1 3 3 0.000000 1.0000 0
17 1 3 3 0.000000 1.0000 0
18 1 1 1 0.000000 1.0000 0
19 1 1 1 0.000000 1.0000 0
20 1 2 2 0.000000 1.0000 0
21 1 2 2 0.000000 1.0000 0
22 1 1 1 0.000000 1.0000 0
23 1 1 1 0.000000 1.0000 0
24 1 2 2 0.000000 1.0000 0
25 1 2 2 0.000000 1.0000 0
26 1 1 1 0.000000 1.0000 0
27 1 2 2 0.000000 1.0000 0
28 2 4 4 0.000000 100.0000 0
26 !NBOND: bonds
1 2 2 3 3 4 4 5
5 6 6 7 7 8 8 9
9 10 10 11 11 12 12 13
13 14 14 15 15 16 16 17
17 18 18 19 19 20 20 21
21 22 22 23 23 24 24 25
25 26 26 27
25 !NTHETA: angles
13 14 15 7 8 9 6 7 8
16 17 18 15 16 17 2 3 4
4 5 6 9 10 11 11 12 13
14 15 16 1 2 3 3 4 5
10 11 12 12 13 14 25 26 27
5 6 7 8 9 10 17 18 19
18 19 20 22 23 24 21 22 23
19 20 21 20 21 22 23 24 25
24 25 26
19 !NPHI: dihedrals
1 2 3 4 2 3 4 5
3 4 5 6 4 5 6 7
8 9 10 11 9 10 11 12
10 11 12 13 11 12 13 14
12 13 14 15 15 16 17 18
16 17 18 19 17 18 19 20
18 19 20 21 19 20 21 22
20 21 22 23 21 22 23 24
22 23 24 25 23 24 25 26
24 25 26 27
0 !NIMPHI: impropers
0 !NDON: donors
0 !NACC: acceptors
0 !NNB
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0
1 0 !NGRP
0 0 0

BIN
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/images/misfolded+chaperonin_t=0tau_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 15 KiB

BIN
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/images/unfolded+chaperonin_t=508750tau_LR.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 17 KiB

216
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/moltemplate_files/1beadFrustrated.lt
View File

@ -1,216 +0,0 @@
# This file defines the "frustrated" coarse-grained protein model used in:
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
# (http://www.pnas.org/content/101/36/13192)
1beadFrustrated {
# There are 3 atom types (referred to above as B, L, and N)
# Define their masses:
write_once("Data Masses") {
@atom:B 1.0
@atom:L 1.0
@atom:N 1.0
}
# AtomID MoleculeID AtomType Charge X Y Z
write('Data Atoms') {
$atom:a1 $mol @atom:L 0.0 -0.92636654 -1.8409904 -2.1482679
$atom:a2 $mol @atom:B 0.0 -0.57313354 -1.0670787 -1.6182341
$atom:a3 $mol @atom:L 0.0 -0.85707399 -1.2358703 -0.69350966
$atom:a4 $mol @atom:B 0.0 -0.44231274 -0.4584993 -0.23418709
$atom:a5 $mol @atom:L 0.0 -0.75081182 -0.62868078 0.69786737
$atom:a6 $mol @atom:B 0.0 -0.36201977 0.11619615 1.2249098
$atom:a7 $mol @atom:N 0.0 -0.63708237 -0.15973084 2.1723919
$atom:a8 $mol @atom:N 0.0 0.20516047 0.10417157 2.624901
$atom:a9 $mol @atom:B 0.0 0.57223743 0.44728103 1.7695617
$atom:a10 $mol @atom:L 0.0 0.77646279 -0.40630393 1.3168043
$atom:a11 $mol @atom:B 0.0 0.45475664 -0.2077937 0.40045721
$atom:a12 $mol @atom:L 0.0 0.72712495 -1.0397637 -0.087614951
$atom:a13 $mol @atom:B 0.0 0.36971183 -0.85840501 -0.9933019
$atom:a14 $mol @atom:L 0.0 0.74784336 -1.5700415 -1.5859217
$atom:a15 $mol @atom:N 0.0 0.43423905 -1.2758917 -2.4853429
$atom:a16 $mol @atom:N 0.0 0.70583191 -0.30726921 -2.4987711
$atom:a17 $mol @atom:N 0.0 -0.091688915 0.23323014 -2.2051358
$atom:a18 $mol @atom:B 0.0 -0.34243283 -0.035822049 -1.2644719
$atom:a19 $mol @atom:B 0.0 0.41961247 0.18475451 -0.65971014
$atom:a20 $mol @atom:L 0.0 0.51968465 1.1546791 -0.77877053
$atom:a21 $mol @atom:L 0.0 -0.40827985 1.2765273 -0.52550748
$atom:a22 $mol @atom:B 0.0 -0.368141 0.58090904 0.19152224
$atom:a23 $mol @atom:B 0.0 0.40327249 0.86101769 0.7336255
$atom:a24 $mol @atom:L 0.0 0.22707289 1.8326235 0.89673346
$atom:a25 $mol @atom:L 0.0 -0.66500182 1.7285809 1.2783166
$atom:a26 $mol @atom:B 0.0 -0.39205603 1.0475436 1.9328097
$atom:a27 $mol @atom:L 0.0 0.25339027 1.5246265 2.5388463
}
# bond-ID bond-Type atom-ID atom-ID
write('Data Bonds') {
$bond:b1 @bond:backbone $atom:a1 $atom:a2
$bond:b2 @bond:backbone $atom:a2 $atom:a3
$bond:b3 @bond:backbone $atom:a3 $atom:a4
$bond:b4 @bond:backbone $atom:a4 $atom:a5
$bond:b5 @bond:backbone $atom:a5 $atom:a6
$bond:b6 @bond:backbone $atom:a6 $atom:a7
$bond:b7 @bond:backbone $atom:a7 $atom:a8
$bond:b8 @bond:backbone $atom:a8 $atom:a9
$bond:b9 @bond:backbone $atom:a9 $atom:a10
$bond:b10 @bond:backbone $atom:a10 $atom:a11
$bond:b11 @bond:backbone $atom:a11 $atom:a12
$bond:b12 @bond:backbone $atom:a12 $atom:a13
$bond:b13 @bond:backbone $atom:a13 $atom:a14
$bond:b14 @bond:backbone $atom:a14 $atom:a15
$bond:b15 @bond:backbone $atom:a15 $atom:a16
$bond:b16 @bond:backbone $atom:a16 $atom:a17
$bond:b17 @bond:backbone $atom:a17 $atom:a18
$bond:b18 @bond:backbone $atom:a18 $atom:a19
$bond:b19 @bond:backbone $atom:a19 $atom:a20
$bond:b20 @bond:backbone $atom:a20 $atom:a21
$bond:b21 @bond:backbone $atom:a21 $atom:a22
$bond:b22 @bond:backbone $atom:a22 $atom:a23
$bond:b23 @bond:backbone $atom:a23 $atom:a24
$bond:b24 @bond:backbone $atom:a24 $atom:a25
$bond:b25 @bond:backbone $atom:a25 $atom:a26
$bond:b26 @bond:backbone $atom:a26 $atom:a27
}
# (3-body) Angles are specified below
# (4-body) Dihedrals must be defined explicitly for every quartet of atoms.
# (These interactions are not determined by atom type.)
# dihedral-ID dihedral-Type atom-ID atom-ID atom-ID atom-ID
write('Data Dihedrals') {
$dihedral:d1 @dihedral:beta $atom:a1 $atom:a2 $atom:a3 $atom:a4
$dihedral:d2 @dihedral:beta $atom:a2 $atom:a3 $atom:a4 $atom:a5
$dihedral:d3 @dihedral:beta $atom:a3 $atom:a4 $atom:a5 $atom:a6
$dihedral:d4 @dihedral:beta $atom:a4 $atom:a5 $atom:a6 $atom:a7
# Dihedral angle forces in the turn regions were switched off
# (in this model) so just I comment them out (and \ the variable names).
# \$dihedral:d5 \@dihedral:turn $atom:a5 $atom:a6 $atom:a7 $atom:a8
# \$dihedral:d6 \@dihedral:turn $atom:a6 $atom:a7 $atom:a8 $atom:a9
# \$dihedral:d7 \@dihedral:turn $atom:a7 $atom:a8 $atom:a9 $atom:a10
$dihedral:d8 @dihedral:beta $atom:a8 $atom:a9 $atom:a10 $atom:a11
$dihedral:d9 @dihedral:beta $atom:a9 $atom:a10 $atom:a11 $atom:a12
$dihedral:d10 @dihedral:beta $atom:a10 $atom:a11 $atom:a12 $atom:a13
$dihedral:d11 @dihedral:beta $atom:a11 $atom:a12 $atom:a13 $atom:a14
$dihedral:d12 @dihedral:beta $atom:a12 $atom:a13 $atom:a14 $atom:a15
# Dihedral angle forces in the turn regions were switched off
# (in this model) so just I comment them out (and \ the variable names).
# \$dihedral:d13 \@dihedral:turn $atom:a13 $atom:a14 $atom:a15 $atom:a16
# \$dihedral:d14 \@dihedral:turn $atom:a14 $atom:a15 $atom:a16 $atom:a17
$dihedral:d15 @dihedral:alpha $atom:a15 $atom:a16 $atom:a17 $atom:a18
$dihedral:d16 @dihedral:alpha $atom:a16 $atom:a17 $atom:a18 $atom:a19
$dihedral:d17 @dihedral:alpha $atom:a17 $atom:a18 $atom:a19 $atom:a20
$dihedral:d18 @dihedral:alpha $atom:a18 $atom:a19 $atom:a20 $atom:a21
$dihedral:d19 @dihedral:alpha $atom:a19 $atom:a20 $atom:a21 $atom:a22
$dihedral:d20 @dihedral:alpha $atom:a20 $atom:a21 $atom:a22 $atom:a23
$dihedral:d21 @dihedral:alpha $atom:a21 $atom:a22 $atom:a23 $atom:a24
$dihedral:d22 @dihedral:alpha $atom:a22 $atom:a23 $atom:a24 $atom:a25
$dihedral:d23 @dihedral:alpha $atom:a23 $atom:a24 $atom:a25 $atom:a26
$dihedral:d24 @dihedral:alpha $atom:a24 $atom:a25 $atom:a26 $atom:a27
}
# All consecutively bonded triplets of atoms same 3-body bond-angle
# interaction parameters. Of coarse, we could specify them all explicitly
# (as we did for the dihedrals above), but I wanted to show how to specify
# angles by atom type instead. (You can do this for dihedrals & impropers
# also.)
# angle-Type atom-Type atom-Type atom-Type bond-Type bond-Type
write_once('Data Angles By Type') {
@angle:backbone @atom:* @atom:* @atom:* @bond:* @bond:*
}
# (The "*" is a wildcard character. I use "*" to denote any atom-type or
# bond-type which is defined within the current namespace: 1beadFrustrated)
# 2-body (non-bonded) interactions:
#
# Uij(r) = 4*eps_ij * (K*(sig_ij/r)^12 + L*(sig_ij/r)^6)
#
# i j pairstylename eps sig K L
#
write_once("In Settings") {
pair_coeff @atom:B @atom:B lj/charmm/coul/charmm/inter 1.0 1.0 1 -1
pair_coeff @atom:B @atom:L lj/charmm/coul/charmm/inter 0.5833333333 1.0 1 0
pair_coeff @atom:B @atom:N lj/charmm/coul/charmm/inter 0.6666666667 1.0 1 0
pair_coeff @atom:L @atom:L lj/charmm/coul/charmm/inter 0.1666666667 1.0 1 1
pair_coeff @atom:L @atom:N lj/charmm/coul/charmm/inter 0.25 1.0 1 0
pair_coeff @atom:N @atom:N lj/charmm/coul/charmm/inter 0.3333333333 1.0 1 0
}
# 2-body (bonded) interactions:
#
# Ubond(r) = (k/2)*(r-0)^2
#
# The corresponding command is:
#
# bond-Type bondstylename k r0
write_once("In Settings") {
bond_coeff @bond:backbone harmonic 100.0 1.0
}
# 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)
#
# angle-Type anglestylename k theta0
write_once("In Settings") {
angle_coeff @angle:backbone harmonic 13.3333333333 105.0
}
# We use tabular dihedral potentials to implement the dihedral forces.
# (Actually there is a way to use Fourier series, using multiple charmm
# style dihedral interactions, but it's slower and messier.)
write_once("In Settings") {
# style file keyword
dihedral_coeff @dihedral:alpha table table_dihedral_frustrated.dat FRUSTRATED_ALPHA
dihedral_coeff @dihedral:beta table table_dihedral_frustrated.dat FRUSTRATED_BETA
# No need to specify dihedral interactions in the turn regions. (none exist)
}
write_once("In Settings") {
# Optional: define the atoms in the "proteins" group
group proteins type @atom:B
group proteins type @atom:L
group proteins type @atom:N
}
# LAMMPS has many available force field styles (and atom styles).
# Here, we pick the ones which work well for this molecular model:
write_once("In Init") {
# --- Default options for the "1BeadFrustrated" protein model ---
# --- (These can be overridden later.) ---
units lj
atom_style full
bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid table spline 360
pair_style hybrid lj/charmm/coul/charmm/inter es4k4l maxmax 3.5 4.0
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 #(turn on "1-4" interactions)
}
} # 1beadFrustrated

85
tools/moltemplate/examples/CG_protein_folding_examples/1bead+chaperone/frustrated+chaperonin/moltemplate_files/1beadFrustrated_variants.lt
View File

@ -1,85 +0,0 @@
import "1beadFrustrated.lt"
# Alternate starting conformation (same molecule):
1beadMisfolded inherits 1beadFrustrated {
# This molecule "inherits" all of its features from "1beadFrustrated".
# Here we override the atomic positions with new coordinates:
# AtomID MoleculeID AtomType Charge X Y Z
write("Data Atoms") {
$atom:a1 $mol @atom:L 0.0 -0.69801399 -0.22114168 -1.9464876
$atom:a2 $mol @atom:B 0.0 -0.40921658 -0.027063664 -1.0033251
$atom:a3 $mol @atom:L 0.0 0.10259348 0.80836418 -1.0737085
$atom:a4 $mol @atom:B 0.0 0.25857916 1.0054984 -0.11621451
$atom:a5 $mol @atom:L 0.0 0.8258629 1.8325549 -0.18529135
$atom:a6 $mol @atom:B 0.0 0.91366257 2.1950317 0.74175977
$atom:a7 $mol @atom:N 0.0 1.4399539 1.554238 1.2994409
$atom:a8 $mol @atom:N 0.0 0.73372573 1.0161012 1.7397275
$atom:a9 $mol @atom:B 0.0 0.26608782 0.65302497 0.94353938
$atom:a10 $mol @atom:L 0.0 0.97442305 0.13574211 0.50586398
$atom:a11 $mol @atom:B 0.0 0.35889617 -0.18247555 -0.1764186
$atom:a12 $mol @atom:L 0.0 0.87151735 -0.77260824 -0.75240916
$atom:a13 $mol @atom:B 0.0 0.047726486 -1.0530682 -1.1902704
$atom:a14 $mol @atom:L 0.0 0.34530697 -1.7476773 -1.8393331
$atom:a15 $mol @atom:N 0.0 0.65865186 -2.45948 -1.2167056
$atom:a16 $mol @atom:N 0.0 -0.16534524 -2.6219442 -0.67112167
$atom:a17 $mol @atom:N 0.0 -0.010590421 -2.2445242 0.24748633
$atom:a18 $mol @atom:B 0.0 0.18135771 -1.2564919 0.1767523
$atom:a19 $mol @atom:B 0.0 -0.57472665 -0.82852797 -0.27027791
$atom:a20 $mol @atom:L 0.0 -1.3967448 -1.0516787 0.24247346
$atom:a21 $mol @atom:L 0.0 -1.003428 -0.85642681 1.1107555
$atom:a22 $mol @atom:B 0.0 -0.25156735 -0.3182346 0.74262946
$atom:a23 $mol @atom:B 0.0 -0.61751956 0.30115562 0.070426493
$atom:a24 $mol @atom:L 0.0 -1.3347934 0.83310182 0.52625934
$atom:a25 $mol @atom:L 0.0 -0.83315257 1.270904 1.2564086
$atom:a26 $mol @atom:B 0.0 -0.10469759 1.6988523 0.72597181
$atom:a27 $mol @atom:L 0.0 -0.57854905 2.3367737 0.11206868
}
} # 1beadMisfolded
1beadUnfolded inherits 1beadFrustrated {
# This molecule "inherits" all of its features from "1beadFrustrated"
# Here we override the atomic positions with new coordinates:
# AtomID MoleculeID AtomType Charge X Y Z
write('Data Atoms') {
$atom:a1 $mol @atom:L 0.0 -2.4 1.7 -0.0
$atom:a2 $mol @atom:B 0.0 -1.8 1.7 0.8
$atom:a3 $mol @atom:L 0.0 -1.2 2.5 0.8
$atom:a4 $mol @atom:B 0.0 -0.6 2.5 -0.0
$atom:a5 $mol @atom:L 0.0 0.0 1.7 -0.0
$atom:a6 $mol @atom:B 0.0 0.6 1.7 0.8
$atom:a7 $mol @atom:N 0.0 1.2 2.5 0.8
$atom:a8 $mol @atom:N 0.0 1.8 2.5 -0.0
$atom:a9 $mol @atom:B 0.0 2.4 1.7 -0.0
$atom:a10 $mol @atom:L 0.0 3.0 1.7 -0.8
$atom:a11 $mol @atom:B 0.0 3.0 0.7 -0.8
$atom:a12 $mol @atom:L 0.0 3.0 0.1 -0.0
$atom:a13 $mol @atom:B 0.0 3.8 -0.5 -0.0
$atom:a14 $mol @atom:L 0.0 3.8 -1.1 -0.8
$atom:a15 $mol @atom:N 0.0 3.0 -1.7 -0.8
$atom:a16 $mol @atom:N 0.0 3.0 -1.7 0.2
$atom:a17 $mol @atom:N 0.0 2.4 -2.5 0.2
$atom:a18 $mol @atom:B 0.0 1.8 -2.5 -0.6
$atom:a19 $mol @atom:B 0.0 1.2 -1.7 -0.6
$atom:a20 $mol @atom:L 0.0 0.6 -1.7 0.2
$atom:a21 $mol @atom:L 0.0 -0.0 -2.5 0.2
$atom:a22 $mol @atom:B 0.0 -0.6 -2.5 -0.6
$atom:a23 $mol @atom:B 0.0 -1.2 -1.7 -0.6
$atom:a24 $mol @atom:L 0.0 -1.8 -1.7 0.2
$atom:a25 $mol @atom:L 0.0 -2.4 -2.5 0.2
$atom:a26 $mol @atom:B 0.0 -3.0 -2.5 -0.6
$atom:a27 $mol @atom:L 0.0 -3.6 -1.7 -0.6
}
} # 1beadUnfolded

Some files were not shown because too many files have changed in this diff Show More