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

This commit is contained in:
sjplimp 2014-05-02 19:17:28 +00:00
parent 9a0d1c765b
commit 6883b645fa
496 changed files with 136255 additions and 0 deletions

View File

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

View File

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

View File

@ -0,0 +1,74 @@
This directory contains scripts used for converting AMBER parameter files
into moltemplate (.LT) format. When a newer version of the AMBER parameters
is eventually published, you can use these scripts to convert the new files
again. (Some tinkering may be necessary.)
The main bash script is a wrapper which simply splits up the parameter (".dat")
file into fragments which (it thinks) correspond to the mass, pair, bond,
angle, dihedral, and improper section of the original .dat file.
(However sometimes it gets this wrong and you have to split it up manually!)
Then this bash script invokes the relevant python script to convert
each section into .LT format:
amberparm_to_mass.py
amberparm_to_pair.py
amberparm_to_bond.py
amberparm_to_angle.py
amberparm_to_dihedral.py
amberparm_to_improper.py
In case this goes wrong, you may have to run these scripts manaully.
Find out how to run this bash script by invoking it without any arguments:
./amberparm2lt.sh
------------ IMPORTANT ------------
BEFORE YOU RUN THIS SCRIPT, BE SURE TO CHANGE THE ORDER OF THE IMPROPER DIHEDRAL
PARAMETERS SO THAT THE "SPECIFIC" IMPROPER DIHEDRALS APPEAR LAST, AND THE
"GENERIC" IMPROPER DIHEDRALS APPEAR FIRST.
For example replace these two lines:
X -o -c -o 1.1 180. 2. JCC,7,(1986),230
X -X -c -o 10.5 180. 2. JCC,7,(1986),230
with these two lines:
X -X -c -o 10.5 180. 2. JCC,7,(1986),230
X -o -c -o 1.1 180. 2. JCC,7,(1986),230
Why:
This is the order that moltemplate expects: generic first. specific last.
So far only the improper dihedral parameters in the gaff.dat file seem
to violate this order. The bonds, angles and dihedrals seem to obey this,
but check to make sure.
There is a discussion of these parameters here:
http://structbio.vanderbilt.edu/archives/amber-archive/2005/3444.php
excerpt:
> > In the parm99 file (for example), sometimes the wild-card is used, as it
> > is done in the following example:
> >
> > X -X -C -O 10.5 180. 2. JCC,7,(1986),230
> >
> > The first example is the specific case while the second one is the generic
> > case. In page # 257 of the AMBER Manual, it is talking about Dihedral
> > Angle, and how these dihedral parameters are used to calculate the
> > energies. I am wondering what the difference between generic and specific
> > case is for improper torsions.
>
> "specific" torsions are search for first, and used if a match is found. If
> no match is found, then a search is made to see if a "generic" (aka wild-card)
> torsion with match.
> ...good luck...dac
Good luck
-Andrew
2014-4-19

View File

@ -0,0 +1,203 @@
#!/bin/sh
SYNTAX_MSG=$(cat <<EOF
Typical Usage:
amberparm2lt.sh gaff.dat GAFF > gaff.lt
You can also try:
amberparm2lt.sh parm94.dat "AMBERFF94 inherits GAFF" > amberff94.lt
(However, this later usage may not work.
You may need to manually split the .dat file and run these scripts instead:
amberparm_pair_to_lt.py, amberparm_bond_to_lt.py, amberparm_angle_to_lt.py...)
Be sure that all of these .py files are in your PATH as well.)
EOF
)
if [ "$#" != "2" ]; then
echo "${SYNTAX_MSG}" >&2
echo "" >&2
echo "Error: This script requires two arguments," >&2
echo " 1) the name of the amber parm file to be converted (eg \"gaff.dat\")" >&2
echo " 2) the name of the moltemplate object to be created (eg \"GAFF\")" >&2
echo " (This may include the \"inherits\" keyword and parent classes.)" >&2
exit 1
fi
MOLTEMPLATE_USAGE_MSG=$(cat <<EOF
# Background information and usage explanation:
# This file contanis a list of atom types and rules for generating bonded
# interactions between these atoms (hopefully) according to AMBER conventions.
# By using the atom types shown below in your own molecules, bonds and angular
# interactions will be automatically generated.
# AMBER (GAFF) force-field parameters will also be assigned to each angle
# interaction (according to these atom types).
# One way to apply the GAFF force field to a particular type of molecule, is
# to use the "inherits" keyword when you define that molecule. For example:
# import("gaff.lt")
# MoleculeType inherits GAFF {
# write_once("Data Atoms") {
# \$atom:C1 \$mol:... @atom:cx 0.0 4.183 3.194 13.285
# \$atom:C2 \$mol:... @atom:cx 0.0 4.291 4.618 13.382
# : : :
# }
# }
#(See "Inheritance" and "short names vs. full names" in the moltemplate manual.)
EOF
)
# (Note that the full name of the atom type in this example is "@atom:/GAFF/cx"
# You can always refer to atom types this way as well. Using "inherits GAFF"
# allows you to use more conventient "@atom:cx" shorthand notation instead.)
echo "####################################################################"
echo "# To use this, LAMMPS currently must be compiled with the USER-MISC package."
echo "# (Type \"make yes-user-misc\" into the shell before compiling LAMMPS.)"
echo "####################################################################"
echo "# This moltemplate (LT) file was generated automatically using"
echo "# amberparm2lt.sh $1 $2"
echo "####################################################################"
echo "$MOLTEMPLATE_USAGE_MSG"
echo "####################################################################"
echo "# Moltemplate can not assign atom charge. You must assign atomic"
echo "# charges yourself. (Moltemplate is only a simple text manipulation tool.)"
echo "####################################################################"
echo ""
echo ""
if ! which ./amberparm_mass_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_mass_to_lt.py\" not found.\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
if ! which ./amberparm_pair_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_pair_to_lt.py\" not found.\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
if ! which ./amberparm_bond_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_bond_to_lt.py\" not found.\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
if ! which ./amberparm_angle_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_angle_to_lt.py\" not found.\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
if ! which ./amberparm_dihedral_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_dihedral_to_lt.py\" not found.\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
if ! which ./amberparm_improper_to_lt.py > /dev/null; then
echo "\nError: \"amberparm_improper_to_lt.py\" not found. (Update your PATH?)\n" >&2
echo " (Try running this script from the directory containing amberparm2lt.sh)" >&2
exit 2
fi
#PARM_FILE='gaff.dat'
PARM_FILE=$1
# sections are separated by blank lines
# some sections have comment lines at the beginning
# The 1st section is the mass (note: skip the first line)
tail -n +2 < "$PARM_FILE" | \
awk -v n=1 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
> "${PARM_FILE}.mass"
# The 2nd section has the list of 2-body bond force-field params
awk -v n=2 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
| tail -n +2 \
> "${PARM_FILE}.bond"
# The 3rd section has the list of 3-body angle force-field params
awk -v n=3 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
> "${PARM_FILE}.angle"
# The 4th section has the list of 4-body dihedral force-field params
awk -v n=4 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
> "${PARM_FILE}.dihedral"
# The 5th section has the list of 4-body improper force-field params
awk -v n=5 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
> "${PARM_FILE}.improper"
# The 6th section has the hbond-parameters (no-longer used. ignore)
awk -v n=6 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
> "${PARM_FILE}.hbond"
# The 7th "section" is just a blank line. (skip that)
# The 8th section has the list of non-bonded ("pair") force-field parameters
awk -v n=8 '{if (NF==0) nblanks++; else {if (nblanks+1==n) print $0}}' \
< "$PARM_FILE" \
| tail -n +2 \
> "${PARM_FILE}.pair"
./amberparm_mass_to_lt.py < "${PARM_FILE}.mass" > "${PARM_FILE}.mass.lt"
./amberparm_pair_to_lt.py < "${PARM_FILE}.pair" > "${PARM_FILE}.pair.lt"
./amberparm_bond_to_lt.py < "${PARM_FILE}.bond" > "${PARM_FILE}.bond.lt"
./amberparm_angle_to_lt.py < "${PARM_FILE}.angle" > "${PARM_FILE}.angle.lt"
./amberparm_dihedral_to_lt.py \
< "${PARM_FILE}.dihedral" > "${PARM_FILE}.dihedral.lt"
./amberparm_improper_to_lt.py \
< "${PARM_FILE}.improper" > "${PARM_FILE}.improper.lt"
echo "$2 {"
echo ""
echo " # ----------------------------------------------------------------------"
#echo " # This file was automatically generated by \"common/amber/amberparm2lt.sh\""
echo " # The basic atom nomenclature and conventions are explained here:"
echo " # http://ambermd.org/antechamber/gaff.pdf"
echo " # For reference, the original gaff.dat file and format documenation are here:"
echo " # http://ambermd.org/AmberTools-get.html"
echo " # http://ambermd.org/formats.html#parm.dat"
echo " # ----------------------------------------------------------------------"
echo ""
cat "$PARM_FILE.mass.lt" \
"$PARM_FILE.pair.lt" \
"$PARM_FILE.bond.lt" \
"$PARM_FILE.angle.lt" \
"$PARM_FILE.dihedral.lt" \
"$PARM_FILE.improper.lt"
AMBER_STYLES_INIT=$(cat <<EOF
write_once("In Init") {
# Default styles and settings for AMBER based force-fields:
units real
atom_style full
bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
improper_style hybrid cvff
pair_style hybrid lj/charmm/coul/long 9.0 10.0 10.0
kspace_style pppm 0.0001
# NOTE: If you do not want to use long-range coulombic forces,
# comment out the two lines above and uncomment this line:
# pair_style hybrid lj/charmm/coul/charmm 9.0 10.0
pair_modify mix arithmetic
special_bonds amber
}
EOF
)
echo "$AMBER_STYLES_INIT"
echo ""
echo "}"
echo ""
echo ""

View File

@ -0,0 +1,50 @@
#!/usr/bin/env python
import sys
lines_gaff = sys.stdin.readlines()
angle_style_name = 'harmonic'
sys.stdout.write(' write_once("In Settings") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:8].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
angletype = '@angle:'+atype1+'-'+atype2+'-'+atype3
tokens= line[8:].split()
keq = tokens[0]
req = tokens[1]
comments=' '.join(tokens[2:])
sys.stdout.write(' angle_coeff '+angletype+' '+angle_style_name+' '+keq+' '+req+' # '+comments+'\n')
sys.stdout.write(' } # (end of angle_coeffs)\n')
sys.stdout.write('\n')
sys.stdout.write(' write_once("Data Angles By Type") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:8].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
angletype = '@angle:'+atype1+'-'+atype2+'-'+atype3
#tokens= line[8:].split()
#keq = tokens[0]
#req = tokens[1]
#comments=' '.join(tokens[2:])
sys.stdout.write(' '+angletype+' @atom:'+at1+' @atom:'+at2+' @atom:'+at3+'\n')
sys.stdout.write(' } # (end of Angles By Type)\n')
sys.stdout.write('\n')

View File

@ -0,0 +1,47 @@
#!/usr/bin/env python
import sys
lines_gaff = sys.stdin.readlines()
bond_style_name = 'harmonic'
sys.stdout.write(' write_once("In Settings") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
tokens= line.split()
atypes = line[:6].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
bondtype = '@bond:'+atype1+'-'+atype2
tokens= line[5:].split()
keq = tokens[0]
req = tokens[1]
comments=' '.join(tokens[2:])
sys.stdout.write(' bond_coeff '+bondtype+' '+bond_style_name+' '+keq+' '+req+' # '+comments+'\n')
sys.stdout.write(' } # (end of bond_coeffs)\n')
sys.stdout.write('\n')
sys.stdout.write(' write_once("Data Bonds By Type") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:6].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
bondtype = '@bond:'+atype1+'-'+atype2
#tokens= line[5:].split()
#keq = tokens[0]
#req = tokens[1]
#comments=' '.join(tokens[2:])
sys.stdout.write(' '+bondtype+' @atom:'+at1+' @atom:'+at2+'\n')
sys.stdout.write(' } # (end of Bonds By Type)\n')
sys.stdout.write('\n')

View File

@ -0,0 +1,157 @@
#!/usr/bin/env python
# SOME UGLY CODE HERE
import sys
lines_gaff = sys.stdin.readlines()
dihedral_style_name = 'fourier'
in_dihedral_coeffs = []
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:11].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
atype4 = atypes[3].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
at4 = atype4.replace('X','*')
dihedraltype = '@dihedral:'+atype1+'-'+atype2+'-'+atype3+'-'+atype4
tokens= line[11:].split()
npth = float(tokens[0])
Kn = float(tokens[1])
Kn /= npth # The coeff for each fourier term is Kn/npth
# ...I THINK (?). (Very confusing. See documentation below...)
dn = float(tokens[2])
n = int(float(tokens[3]))
comments=' # '+(' '.join(tokens[4:]))
in_dihedral_coeffs.append([dihedraltype, Kn, n, dn, comments])
#print(Kn, n, dn)
#for entry in in_dihedral_coeffs:
# print(entry)
#exit()
# ---- processing dihedral fourier series ----
# ---- (negative "n" values means the
# ---- Fourier series is not yet complete.
i = 0
while i < len(in_dihedral_coeffs):
type_str = in_dihedral_coeffs[i][0]
Kn = in_dihedral_coeffs[i][1]
n = in_dihedral_coeffs[i][2]
dn = in_dihedral_coeffs[i][3]
#if (i>0):
# sys.stderr.write('prev_n='+str(in_dihedral_coeffs[i-1][-3])+'\n')
#sys.stderr.write('n='+str(n)+'\n')
if ((i>0) and (in_dihedral_coeffs[i-1][-3] < 0)):
#sys.stdout.write('interation_before_append: '+str(in_dihedral_coeffs[i-1])+'\n')
assert(in_dihedral_coeffs[i-1][0] == in_dihedral_coeffs[i][0])
in_dihedral_coeffs[i-1][-3] = -in_dihedral_coeffs[i-1][-3]
comments = in_dihedral_coeffs[i-1][-1]
in_dihedral_coeffs[i-1][-1] = Kn
in_dihedral_coeffs[i-1].append(n)
in_dihedral_coeffs[i-1].append(dn)
in_dihedral_coeffs[i-1].append(comments)
#sys.stdout.write('interation_after_append: '+str(in_dihedral_coeffs[i-1])+'\n')
del in_dihedral_coeffs[i]
#elif len(in_dihedral_coeffs) < 3:
# del in_dihedral_coeffs[i]
else:
i += 1
for i in range(0, len(in_dihedral_coeffs)):
type_str = in_dihedral_coeffs[i][0]
params = in_dihedral_coeffs[i][1:]
params = map(str, params)
num_fourier_terms = (len(params)-1)/3
dihedral_coeff_str = 'dihedral_coeff '+type_str+' '+\
dihedral_style_name+' '+\
str(num_fourier_terms)+' '+ \
' '.join(params)
in_dihedral_coeffs[i] = dihedral_coeff_str
# ---- finished processing dihedral fourier series ----
sys.stdout.write(' write_once(\"In Settings\") {\n ')
sys.stdout.write('\n '.join(in_dihedral_coeffs)+'\n')
sys.stdout.write(' } # (end of dihedral_coeffs)\n')
sys.stdout.write('\n')
sys.stdout.write(' write_once("Data Dihedrals By Type") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:11].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
atype4 = atypes[3].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
at4 = atype4.replace('X','*')
dihedraltype = '@dihedral:'+atype1+'-'+atype2+'-'+atype3+'-'+atype4
sys.stdout.write(' '+dihedraltype+' @atom:'+at1+' @atom:'+at2+' @atom:'+at3+' @atom:'+at4+'\n')
sys.stdout.write(' } # (end of Dihedrals By Type)\n')
sys.stdout.write('\n')
"""
- 6 - ***** INPUT FOR DIHEDRAL PARAMETERS *****
IPT , JPT , KPT , LPT , IDIVF , PK , PHASE , PN
FORMAT(A2,1X,A2,1X,A2,1X,A2,I4,3F15.2)
IPT, ... The atom symbols for the atoms forming a dihedral
angle. If IPT .eq. 'X ' .and. LPT .eq. 'X ' then
any dihedrals in the system involving the atoms "JPT" and
and "KPT" are assigned the same parameters. This is
called the general dihedral type and is of the form
"X "-"JPT"-"KPT"-"X ".
IDIVF The factor by which the torsional barrier is divided.
Consult Weiner, et al., JACS 106:765 (1984) p. 769 for
details. Basically, the actual torsional potential is
(PK/IDIVF) * (1 + cos(PN*phi - PHASE))
PK The barrier height divided by a factor of 2.
PHASE The phase shift angle in the torsional function.
The unit is degrees.
PN The periodicity of the torsional barrier.
NOTE: If PN .lt. 0.0 then the torsional potential
is assumed to have more than one term, and the
values of the rest of the terms are read from the
next cards until a positive PN is encountered. The
negative value of pn is used only for identifying
the existence of the next term and only the
absolute value of PN is kept.
The input is terminated by a blank card.
"""

View File

@ -0,0 +1,90 @@
#!/usr/bin/env python
import sys
lines_gaff = sys.stdin.readlines()
improper_style_name = 'cvff'
sys.stdout.write(' write_once("In Settings") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:11].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
atype4 = atypes[3].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
at4 = atype4.replace('X','*')
#impropertype = '@improper:'+atype1+'-'+atype2+'-'+atype3+'-'+atype4
#sys.stdout.write(' '+impropertype+' @atom:'+at1+' @atom:'+at2+' @atom:'+at3+' @atom:'+at4+'\n')
# Oops. This is incorrect.
# In moltemplate, the central atom is the first atom,
# In "gaff.dat", the central atom is the third atom
# http://archive.ambermd.org/201307/0519.html
impropertype = '@improper:'+atype3+'-'+atype1+'-'+atype2+'-'+atype4
tokens= line[11:].split()
Kn = float(tokens[0])
dn = float(tokens[1])
n = int(float(tokens[2]))
comments=' '.join(tokens[3:])
if (dn < 0.001):
sys.stdout.write(' improper_coeff '+impropertype+' '+improper_style_name+' '+str(Kn)+' 1 '+str(n)+' # '+comments+'\n')
elif (179.999 < abs(dn) < 180.001):
sys.stdout.write(' improper_coeff '+impropertype+' '+improper_style_name+' '+str(Kn)+' -1 '+str(n)+' # '+comments+'\n')
else:
sys.stderr.write('Error: Illegal bondImproper parameters:\n'
' As of 2013-8-03, LAMMPS doens hot have an improper style\n'
' which can handle impropers with gamma != 0 or 180\n')
exit(-1)
sys.stdout.write(' } # (end of improper_coeffs)\n')
sys.stdout.write('\n')
sys.stdout.write(' write_once("Data Impropers By Type") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
atypes = line[:11].split('-')
atype1 = atypes[0].strip()
atype2 = atypes[1].strip()
atype3 = atypes[2].strip()
atype4 = atypes[3].strip()
at1 = atype1.replace('X','*')
at2 = atype2.replace('X','*')
at3 = atype3.replace('X','*')
at4 = atype4.replace('X','*')
#impropertype = '@improper:'+atype1+'-'+atype2+'-'+atype3+'-'+atype4
#sys.stdout.write(' '+impropertype+' @atom:'+at1+' @atom:'+at2+' @atom:'+at3+' @atom:'+at4+'\n')
# Oops. This is incorrect.
# In moltemplate, the central atom is the first atom,
# In "gaff.dat", the central atom is the third atom
# http://archive.ambermd.org/201307/0519.html
impropertype = '@improper:'+atype3+'-'+atype1+'-'+atype2+'-'+atype4
sys.stdout.write(' '+impropertype+' @atom:'+at3+' @atom:'+at1+' @atom:'+at2+' @atom:'+at4+'\n')
sys.stdout.write(' } # (end of Impropers By Type)\n')
sys.stdout.write('\n')
# NOTE: AMBER documentation is not clear how the improper angle is defined.
# It's not clear if we should be using the dihedral angle between
# planes I-J-K and J-K-L. As of 2014-4, improper_style cvff does this.
# Even if we create improper interactions with the angle defined between
# the wrong planes, at least the minima should be the same
# (0 degrees or 180 degrees).
# So I'm not too worried we are getting this detail wrong long as
# we generate new impropers realizing that the 3rd atom (K) is the
# central atom (according to AMBER conventions).
#
# http://structbio.vanderbilt.edu/archives/amber-archive/2007/0408.php
#
# Currently, we only apply improper torsional angles for atoms
# in a planar conformations. Is it clear?
# Junmei

View File

@ -0,0 +1,19 @@
#!/usr/bin/env python
import sys
lines_gaff = sys.stdin.readlines()
sys.stdout.write(' write_once(\"Data Masses\") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
tokens= line.split()
atype = tokens[0]
mass=tokens[1]
# what is the next number? (the one in tokens[2]?)
comments=' '.join(tokens[3:])
sys.stdout.write(' @atom:'+atype+' '+mass+' # '+comments+'\n')
sys.stdout.write(' } # (end of masses)\n')
sys.stdout.write('\n')

View File

@ -0,0 +1,32 @@
#!/usr/bin/env python
import sys
lines_gaff = sys.stdin.readlines()
#pair_style = 'lj/charmm/coul/long'
# NOTE: Long-range coulombic forces were disabled intentionally. (See below)
# If you want to use long-range electrostatics, uncomment these lines:
# Instead I use hybrid lj/charmm/coul/charmm by default, because
# LAMMPS complains if you attempt to use lj/charmm/coul/long on a
# system if it does not contain any charged particles.
# Currently, moltemplate does not assign atomic charge,
# so this problem occurs frequently.
#pair_style = 'lj/charmm/coul/charmm'
pair_style = 'lj/charmm/coul/long'
sys.stdout.write(' write_once(\"In Settings\") {\n')
for i in range(0, len(lines_gaff)):
line = lines_gaff[i]
tokens= line.split()
atype = tokens[0]
sig=tokens[1]
eps=tokens[2]
comments=' '.join(tokens[3:])
sys.stdout.write(' pair_coeff @atom:'+atype+' @atom:'+atype+' '+pair_style+' '+eps+' '+sig+' # '+comments+'\n')
sys.stdout.write(' } # (end of pair_coeffs)\n')
sys.stdout.write('\n')

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,106 @@
# This file contains a unit cell for building graphene and nanotubes
#
#
# The 2AtomCellAlignX "molecule" defined below is a minimal unit cell for any
# hexagonal tesselation in 2-dimensions. (See "graphene_unit_cell.jpg")
# Surfaces constructed with this unit cell can be flat or curved into tubes.
# The distance between nearest-neighbor carbon atoms (ie the length of a
# carbon-carbon bond) is equal to "d" which I set to 1.420 Angstroms.
#
# d = length of each hexagon's side = 1.42 Angstroms
# L = length of each hexagon = 2*d = 2.84 Angstroms
# W = width of each hexagon = 2*d*sqrt(3)/2 = 2.4595121467478056 Angstroms
#
# Consequently, the Lattice-cell vectors for singe-layer graphene are:
# (2.4595121467478, 0, 0) (aligned with X axis)
# (1.2297560733739, 2.13, 0) (2.13 = 1.5*d)
# So, to build a sheet of graphite, you could use:
# sheet = new Graphene/2AtomCellAlignX [10].move(2.4595121467478,0,0)
# [10].move(1.2297560733739,2.13,0)
Graphene {
2AtomCellAlignX
{
# 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
}
}
# Now define properties of the Carbon graphene atom
write_once("In Init") {
pair_style hybrid lj/charmm/coul/charmm 9.0 10.0
}
write_once("Data Masses") {
@atom:C 12.0
}
write_once("In Settings") {
# 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)
# Define a group consisting of only carbon atoms in graphene molecules
group Cgraphene type @atom:C
}
# Notice that the two atoms in the unit-cell above lie in the XY plane.
# (Their z-coordinate is zero). It's also useful to have a version of
# this object which lies in the XZ plan. So we define this below:
2AtomCellAlignXZ = 2AtomCellAlignX.rot(90,1,0,0)
} # Graphene
# ------------ Graphite -----------
#
# Note: For graphite: sheets stacked in the Z direction are separated by a
# distance of 3.35 Angstroms, and shifted in an alternating +/-Y direction
# by a distance of d (1.42 Angstroms). To add additional graphene layers
# you could use:
# sheet2 = new Graphene/2AtomCellAlignX [10].move(2.4595121467478,0,0)
# [10].move(1.2297560733739,2.13,0)
# sheet2[*][*].move(0, 1.42, 3.35)
# sheet3 = new Graphene/2AtomCellAlignX [10].move(2.4595121467478,0,0)
# [10].move(1.2297560733739,2.13,0)
# sheet3[*][*].move(0, -1.42, 6.70)
# etc...
# However, to build a thick sheet of graphite, it would
# be more efficient to use a 4-atom unit cell:
#
#Graphene {
# GraphiteCell {
# # 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
# $atom:C3 $mol:... @atom:../C 0.0 -0.61487803668695 1.065 3.35
# $atom:C4 $mol:... @atom:../C 0.0 0.61487803668695 1.775 3.35
# }
# } # GraphiteCell
#}
#
# Then you could create a thick sheet of graphite this way:
#
# graphite = new Graphene/GraphiteCell [10].move(2.4595121467478,0,0)
# [10].move(1.2297560733739,2.13,0)
# [5].move(0,0,6.70)

Binary file not shown.

After

Width:  |  Height:  |  Size: 13 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 21 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 5.5 KiB

View File

@ -0,0 +1,3 @@
OPLSAA force-field conversion tools provided by Jason Lambert.

View File

@ -0,0 +1,118 @@
Unfortunately, moltemplate does not come with a file containing OPLSAA
paramters which is ready to use. You must build it yourself.
This directory has tools and instructions to explain how to do this.
-----------------------------
When you want to run a new simulation, you must download the full
"oplsaa.prm" force-field file (from the TINKER web site) and manually
delete all the atom types which you do not need. (See below.)
Then you must use the "oplsaa_moltemplate.py" script to create
"oplsaa.lt" file which moltemplate.sh needs. Then you must run moltemplate.sh.
You must do this for each new simulation you plan to run which uses OPLSAA.
---- Details: ----
Download the original "oplsaa.prm" file here:
http://dasher.wustl.edu/tinker/distribution/params/oplsaa.prm
or here:
http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm
and save this file as "oplsaa_subset.prm". Then you must EDIT THIS FILE
so that it only contains atom types you plan to have in your simulation
(see below). Finally, you must run the opls_moltemplate.py script this way:
python oplsaa_moltemplate.py oplsaa_subset.prm
This will create a file named "oplsaa.lt"
Look over the newly created "oplsaa.lt" file.
Then, if necessary, move this file to wherever you plan to run moltemplate.
There is a directory containing an example how to use this code located here:
examples/all_atom_examples/OPLSAA_force_field_examples/ethylene/moltemplate_files/oplsaa_lt_generator/README.TXT
----- DETAILS: Editing the "oplsaa_subset.prm" file -------
Again, before you run "oplsaa_moltemplate.py", you must edit the "oplsaa.prm"
file (or "oplsaa_subset.prm file) and eliminate atom types which do not
correspond to any of the atoms in your simulation. This means you must
look for lines near the beginning of this file which begin with the word "atom"
and refer to atom types which appear in the simulation you plan to run. All
other lines (beginning with the word "atom") must be deleted or commented out.
(Leave the rest of the file alone.)
For example:
If you were working with ethylene and benzene you would delete every line
beginning with the word "atom", except for these four lines:
# for ethylene:
atom 88 47 CM "Alkene H2-C=" 6 12.011 3
atom 89 46 HC "Alkene H-C=" 1 1.008 1
# for benzene:
atom 90 48 CA "Aromatic C" 6 12.011 3
atom 91 49 HA "Aromatic H-C" 1 1.008 1
Then you are ready to run oplsaa_moltemplate.py on this file.
(You can try to delete more irrelevant information, but be careful.
It is not necessary, and it is easy to make mistakes.)
----- Using the "oplsaa.lt" file -----
Once you have created the "oplsaa.lt" file, you can create files (like
ethylene.lt) which define molecules that refer to these atom types.
Here is an excerpt from "ethyelene.lt":
Ethylene inherits OPLSAA {
write('Data Atoms') {
list of atoms goes here ...
}
write('Data Bond List') {
list of bonds goes here ...
}
}
And then run moltemplate.
----------- CHARGE: -----------
By default, the OPLSAA force-field assigns atom charge according to atom type.
When you run moltemplate, it will create a file named "system.in.charges",
containing commands like:
set type 2 charge -0.42
set type 3 charge 0.21
(This assumes your main moltemplate file is named "system.lt". If it was
named something else, eg "polymer.lt", then the file created by moltemplate
will be named "polymer.in.charges".)
Include these commands somewhere in your LAMMPS input script
(or use the LAMMPS "include" command to load the commands in system.in.charges)
Note that the atom numbers (eg "2", "3") in this file will not match the
OPLS atom numbers. (Check the output_ttree/ttree_assignments.txt file,
created by moltemplate, to see a table of "@atom" type numbers translated
from OPLSAA into LAMMPS.)
----------- CREDIT -----------
If you use these tools and you publish a paper using OPLSAA, please also cite
the TINKER program. (Because the original "oplsaa.prm" file which we depend on
is distributed with TINKER.) I think these are the relevant citations:
1) Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
2) Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
-------------------------------
Jason Lambert and Andrew Jewett
April, 2014
Please email bugs to jewett.aij@gmail.com and jlamber9@gmail.com

View File

@ -0,0 +1,52 @@
MOST USERS CAN IGNORE THIS DIRECTORY.
This is the original version of oplsaa.txt released with
moltemplate version 1.20. Later on I decided to direct
users to the (hopefully current) version of the "oplsaa.prm"
file distributed with TINKER's force-field-explorer program.
http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm
http://dasher.wustl.edu/tinker/distribution/params/oplsaa.prm
The original oplsaa.txt from moltemplate version 1.20 file is located here.
Feel free to pick whichever one works well for your application, but
if you publish a paper using these parameters please cite TINKER,
since that appears to be the source of all of these files.
I think these are the relevant citations:
Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
Here is a discussion of a difference between these two files:
>> On Fri, Apr 11, 2014 at 6:44 PM, Andrew Jewett <jewett.aij@gmail.com> wrote:
>> Jason
>> Lambert graciously wrote the "oplsaa_moltemplate.py" script and
>> provided the "oplsaa.txt" file which it reads as an input file (after
>> some minor manual editing on the part of the user). I did some
>> googling, and I noticed that this file is very similar to the
>> "oplsaa.prm" file located at this web site:
>>
>> http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm
>>
>> I suspect it uses the same parameter conventions.
>>
>> However I just realized that the two files are not identical. (I am
>> emailing Jason now to ask about the difference.)
>>
>> Where did you get the data for the "oplsaa.txt" file? More
>> specifically, where do the atom numbers come from?
>On Tue, Apr 15, 2014 at 6:06 PM, jason Lambert <jlamber9@gmail.com> wrote:
>To be honest I do not know the original source of the OPLSA force
> field file. The file was provided by a graduate student from a group
> I worked with and I did not question it other than spot checking with
> values that were published. I am surprised that the files are different.
> Anyways, I am not aware of any logic to the observed switches in
> the parameters. I did a diff -y on the two files and the difference is
> primarily atom switching. The code works with either one though
> so it is not a problem. We can provide a link to the prm file as well.
> I know I added some dihedrals explicitly but the ones I added were
> already encompassed by the dihedrals containing wild cards.

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,386 @@
#! /usr/bin/env python
#The purpose of this script is to create a moltemplate lt file for the opls-aa forcefield.
#This will assist researchers in building complex simulations using this OPLS-UA and the OPLS-AA forcefields.
__author__="Jason Lambert"
__version__="0.15"
import sys
import os
from operator import itemgetter
print("""
Warning:
Run this program on a SUBSET of the OPLS atoms which are relevant
to your problem. If not, this program (and moltemplate) may crash
your computer and/or generate enormous files that you do not need.
""")
# To do that, first make a copy of the \"oplsaa.prm\" file
# (which can be downloaded from the TINKER web site).
# The lines in this file beginning with the word \"atoms\" should
# define the atoms which you plan to put in your simulation. All other
# lines beginning with the word \"atoms\" should be deleted.
# (Leave the other sections of this file alone.)
#""")
#input data from file containing opls aa force field parameters.
try:
f=open(sys.argv[1],"r")
except:
print("need to specify file name as an input argument:")
print("python oplsaa_moltemplate.py <forcefield file name>")
print("or file name is specified incorrectly")
sys.exit()
#output lt file
g=open("oplsaa.lt","w")
lines = f.readlines()
# Ignore/Comment out lines before the "## Atom Type Definitions ##" section.
for i in range(0, len(lines)):
if (lines[i].find("## Atom Type Definitions ##") != -1):
break
else:
lines[i] = '# ' + lines[i]
# As of 2014-4-19, there appear to be 906 atom types, but we don't assume this.
# First try to infer out how many atom types there were in the original
# oplsaa.prm file, or at least find an upper bound on the atom-type numbers.
# (Keep track of the maximum value of the first column in the "atom" section.)
max_atomType = 0
for line in lines:
# skip over text after a # comment character
ic = line.find('#')
if ic != -1:
line = (line[:ic]).strip()
else:
line = line.strip()
# now look for lines beginning with the word "atom"
tokens = line.split()
if ((len(tokens)>2) and
(tokens[0] == "atom") and
(int(tokens[1]) > max_atomType)):
max_atomType = int(tokens[1])
#temporary storage file
h=open("temp.txt","w+")
atom_lookup={} #this dictionary contains all the atom ffid's as a key and the number of atoms with that key
#atom=[[10000,10000] for i in range(906)] <- don't assume there are 906 atoms
atom=[[-10000,-10000] for i in range(0,max_atomType+1)]
#charge_by_type={} # lookup charge by atom type
#vdw_by_type={} # lookup epsilon & sigma paramters by atom type
charge_by_type=[0.0 for i in range(0,max_atomType+1)] # lookup charge by atom
vdw_by_type=[(0.0,0.0) for i in range(0,max_atomType+1)] # lookup epsilon & sigma
#atom is declared this way so for sorting purposes.
#atom contains the following data upon allocation
#atom[][0]=atom_id( Important for partial charges and non_bonded interactions)
#atom[][1]=atom_ffid( Important for stretches, bending, torsions and impropers)
#atom[][2]=atom_mass
#atom[][3]=partial charge
#atom[][4]=non_bonding sigma
#atom[][5]=non_bonding epsilon
#atom[][6]=atom comment
bond=[]
#bond contains the following data
#bond[0]=atom 1 ffid
#bond[1]=atom 2 ffid
#bond[2]=bond spring constant(OPLS-aa compatible)
#bond[3]=equilibrium bond distance(Angstrom)
angle=[]
#angle contains the following data
#angle[0]=atom 1 ffid
#angle[1]=atom 2 ffid
#angle[2]=atom 3 ffid
#angle[3]=spring constant
#angle[4]=equilibrium angle (degrees)
dihedral=[]
#dihedral contains the following data
#dihedral[0]=atom 1 ffid
#dihedral[1]=atom 2 ffid
#dihedral[2]=atom 3 ffid
#dihedral[3]=atom 4 ffid
#dihedral[4]=v1
#dihedral[5]=v2
#dihedral[6]=v3
#dihedral[7]=v4
improper=[]
#improper[0]=atom 1 ffid
#improper[1]=atom 2 ffid(central atom)
#improper[2]=atom 3 ffid
#improper[3]=atom 4 ffid
#improper[4]=spring coefficient
#improper[5]=equilibrium angle
#This section gets all the parameters from the force field file
for line in lines:
# skip over text after a # comment character
ic = line.find('#')
if ic != -1:
line = (line[:ic]).strip()
else:
line = line.strip()
if line.find("atom") == 0:
line=line.split()
atom[int(line[1])-1]=[int(line[1]),int(line[2]),float(line[-2]),
0.0,0.0,0.0," ".join(line[3:-2])]
elif line.find("vdw") == 0:
line=line.split()
#vdw_temp.append([float(line[1]),float(line[2]),float(line[3])])
if (int(line[1]) <= max_atomType):
vdw_by_type[int(line[1])] = (float(line[2]),float(line[3]))
elif line.find("bond") == 0:
line=line.split()
bond.append([int(line[1]),int(line[2]),float(line[3]),float(line[4])])
elif line.find("angle") == 0:
line=line.split()
angle.append([int(line[1]),int(line[2]),int(line[3]),
float(line[4]),float(line[5])])
elif line.find("torsion") == 0:
line=line.split()
dihedral.append([int(line[1]),int(line[2]),int(line[3]),int(line[4]),
float(line[5]),float(line[8]), float(line[11]), 0.0])
elif line.find("charge") == 0:
line=line.split()
#charge_temp.append([int(line[1]),float(line[2])])
if (int(line[1]) <= max_atomType):
charge_by_type[int(line[1])] = float(line[2])
elif line.find("imptors") == 0:
line=line.split()
improper.append([int(line[1]), int(line[2]),
int(line[3]), int(line[4]), float(line[5]), float(line[6])])
if len(atom) > 600:
sys.stderr.write("WARNING: The number of atom types in your file exceeds 600\n"
" (You were supposed to edit out the atoms you don't need.\n"
" Not doing this may crash your computer.)\n"
"\n"
" Proceed? (Y/N): ")
reply = sys.stdin.readline()
if find(reply.strip().lower(), 'y') != 0:
exit(0)
#adding the charge and Lennard Jones parameters to
#to each atom type.
#----------------------------------------------#
for i in range(0,len(atom)):
atom_type_num = atom[i][0]
#q = charge_by_type.get(atomTypeNum)
#if q:
# atom[i][3] = q
if atom_type_num != -10000:
q = charge_by_type[atom_type_num]
atom[i][3] = q
for i in range(0,len(atom)):
atom_type_num = atom[i][0]
#vdw_params = vdw_by_type.get(atomTypeNum)
#if vdw_params:
# atom[i][4] = vdw_params[0]
# atom[i][5] = vdw_params[1]
if atom_type_num != -10000:
vdw_params = vdw_by_type[atom_type_num]
atom[i][4] = vdw_params[0]
atom[i][5] = vdw_params[1]
del(charge_by_type)
del(vdw_by_type)
#----------------------------------------------------------#
#begin writing content to lt file
g.write("OPLSAA {\n\n" )
#write out the atom masses
#----------------------------------------------------------#
g.write(" write_once(\"Data Masses\"){\n")#checked with gaff
for i,x in enumerate(atom):
if x[0] != -10000:
g.write(" @atom:{} {} #{} partial charge={}\n".format(
x[0],x[2],x[6],x[3]))
g.write(" } #(end of atom masses)\n\n")
#----------------------------------------------------------#
#write out the pair coefficients
#----------------------------------------------------------#
g.write(" write_once(\"In Settings\"){\n")#checked with gaff
for i,x in enumerate(atom):
if x[0] != -10000:
g.write(" pair_coeff @atom:{0} @atom:{0} lj/cut/coul/long {1} {2}\n".format(x[0],x[5],x[4]))
g.write(" } #(end of pair coeffs)\n\n")
g.write(" write_once(\"In Charges\"){\n")#checked with gaff
for i,x in enumerate(atom):
if x[0] != -10000:
g.write(" set type @atom:{0} charge {1}\n".format(x[0],x[3]))
g.write(" } #(end of atom charges)\n\n")
#-----------------------------------------------------------#
# This part of the code creates a lookup dictionary
# that allows you to find every type of atom by its
# force field id. force field id is the id number
# relevant to bonds, angles, dihedrals, and impropers.
# This greatly increases the speed of angle, bond, dihedral
# and improper assignment.
#------------------------------------------------------------#
atom=sorted(atom,key=itemgetter(1))
atom_ffid=0
for x in atom:
if x[1]==atom_ffid:
atom_lookup[x[1]].append(x[0])
elif x[1]>atom_ffid:
atom_lookup[x[1]]=[x[0]]
atom_ffid=x[1]
atom_lookup[0]=["*"]
#-------------------------------------------------------------#
#writing out the bond coefficients and bond parameters#
#-------------------------------------------------------------#
g.write(" write_once(\"In Settings\") {\n")
index1=0
for x in bond:
for y in atom_lookup.get(x[0],[]):
for z in atom_lookup.get(x[1],[]):
g.write(" bond_coeff @bond:{}-{} harmonic {} {}\n".format(y,z,x[2]/2,x[3]))
h.write(" @bond:{0}-{1} @atom:{0} @atom:{1}\n".format(y,z))
g.write(" } #(end of bond_coeffs)\n\n")
h.seek(0,0)
g.write(" write_once(\"Data Bonds By Type\") {\n")
for line in h.readlines():
g.write(line)
g.write(" } #(end of bonds by type)\n\n")
del(bond)
h.close()
#-----------------------------------------------------------#
h=open("temp.txt","w+")
#writing out angle coefficients and angles by type.---------#
#-----------------------------------------------------------#
g.write(" write_once(\"Data Angles By Type\"){\n")
for x in angle:
for y in atom_lookup.get(x[0],[]):
for z in atom_lookup.get(x[1],[]):
for u in atom_lookup.get(x[2],[]):
#print(y,z,u,x)
h.write(" angle_coeff @angle:{}-{}-{} harmonic {} {}\n".format(y,z,u,
x[3]/2.0,x[4]))
g.write(" @angle:{0}-{1}-{2} @atom:{0} @atom:{1} @atom:{2}\n".format(
y,z,u))
g.write(" } #(end of angles by type)\n\n")
h.seek(0,0)
g.write(" write_once(\"In Settings\" ){\n")
for line in h.readlines():
g.write(line)
g.write(" } #(end of angle_coeffs)\n\n")
del(angle)
h.close()
#----------------------------------------------------------#
#writing dihedrals by type and dihedral coefficients-------#
h=h=open("temp.txt","w+")
g.write(" write_once(\"Data Dihedrals By Type\") {\n")
#print(atom_lookup)
for x in dihedral:
for y in atom_lookup.get(x[0],[]):
for z in atom_lookup.get(x[1],[]):
for u in atom_lookup.get(x[2],[]):
for v in atom_lookup.get(x[3],[]):
if x[0]!=0 and x[3]!=0:
g.write(" @dihedral:{0}-{1}-{2}-{3} @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
y,z,u,v))
h.write(" dihedral_coeff @dihedral:{}-{}-{}-{} opls {} {} {} {}\n".format(
y,z,u,v,x[4],x[5],x[6],x[7]))
elif x[0]==0 and x[3]!=0:
g.write(" @dihedral:0-{1}-{2}-{3} @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
y,z,u,v))
h.write(" dihedral_coeff @dihedral:0-{}-{}-{} opls {} {} {} {}\n".format(
z,u,v,x[4],x[5],x[6],x[7]))
elif x[0]==0 and x[3]==0:
g.write(" @dihedral:0-{1}-{2}-0 @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
y,z,u,v))
#h.write(" dihedral_coeff @dihedral:0-{}-{}-0 harmonic {} {} {} {}\n".format(
h.write(" dihedral_coeff @dihedral:0-{}-{}-0 opls {} {} {} {}\n".format(
z,u,x[4],x[5],x[6],x[7]))
del(dihedral)
g.write(" } #(end of Dihedrals by type)\n\n")
h.seek(0,0)
g.write(" write_once(\"In Settings\") {\n")
for line in h.readlines():
g.write(line)
g.write(" } #(end of dihedral_coeffs)\n\n")
h.close()
#-----------------------------------------------------------------------#
#----writing out improper coefficients and impropers by type------------#
h=open("temp.txt","w+")
g.write(" write_once(\"Data Impropers By Type\") {\n")
for x in improper:
for y in atom_lookup.get(x[0],[]):
for z in atom_lookup.get(x[1],[]):
for u in atom_lookup.get(x[2],[]):
for v in atom_lookup.get(x[3],[]):
if x[0]==0 and x[1]==0 and x[3]==0:
g.write(" @improper:{2}-0-0-0 @atom:{2} @atom:{0} @atom:{1} @atom:{3}\n".format(
y,z,u,v))
h.write(" improper_coeff @improper:{2}-0-0-0 harmonic {4} {5} \n".format(
y,z,u,v,x[4]/2,0))
else:
g.write(" @improper:{2}-0-0-{3} @atom:{2} @atom:{0} @atom:{1} @atom:{3}\n".format(
y,z,u,v))
h.write(" improper_coeff @improper:{2}-0-0-{3} harmonic {4} {5} \n".format(
y,z,u,v,x[4]/2,0))
g.write(" } #(end of impropers by type)\n\n")
h.seek(0,0)
g.write(" write_once(\"In Settings\") {\n")
for line in h.readlines():
g.write(line)
g.write(" } #(end of improp_coeffs)\n\n")
#-----------------------------------------------------------------------#
#This section writes out the input parameters required for an opls-aa simulation
# lammps.
g.write(" write_once(\"In Init\") {\n")
g.write(" units real\n")
g.write(" atom_style full\n")
g.write(" bond_style hybrid harmonic\n")
g.write(" angle_style hybrid harmonic\n")
g.write(" dihedral_style hybrid opls\n")
g.write(" improper_style hybrid harmonic\n")
#g.write(" pair_style hybrid lj/cut/coul/cut 10.0 10.0\n")
g.write(" pair_style hybrid lj/cut/coul/long 10.0 10.0\n")
g.write(" pair_modify mix arithmetic\n")
g.write(" special_bonds lj 0.0 0.0 0.5\n")
g.write(" kspace_style pppm 0.0001\n")
g.write(" } #end of init parameters\n\n")
g.write("} # OPLSAA\n")
f.close()
g.close()
h.close()
os.remove("temp.txt")

View File

@ -0,0 +1,52 @@
# 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 1000.0 1.0
angle_coeff @angle:HOH harmonic 1000.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 fShakeSPCE spce shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fShakeSPCE during minimization.)
}
} # end of definition of "SPCE" water molecule type

View File

@ -0,0 +1,57 @@
# This ice (1h) unit cell is rectangular and contains 8 water molecules.
# (Coordinates and cell dimensions converted were from a PDB file.)
# The dimensions of the unit cell (in Angstroms) are: 4.521 7.832 7.362
import "spce.lt" # <-- define the "SPCE" molecule
SpceIceRect8 {
# Create a 3-dimensional array of 8 water molecules
wat = new SPCE[2][2][2]
# Array indices will be correlated with position [xindex][yindex][zindex]
# You can overwrite coordinates of atoms after they were created this way:
# (Order is not important)
# atom-ID molecule-ID atomType charge newX newY newZ
write("Data Atoms") {
$atom:wat[1][0][0]/O $mol:wat[1][0][0] @atom:SPCE/O -0.8476 3.391 1.305 1.381
$atom:wat[1][0][0]/H1 $mol:wat[1][0][0] @atom:SPCE/H 0.4238 3.391 0.370 1.710
$atom:wat[1][0][0]/H2 $mol:wat[1][0][0] @atom:SPCE/H 0.4238 2.582 1.772 1.710
$atom:wat[1][0][1]/O $mol:wat[1][0][1] @atom:SPCE/O -0.8476 3.391 1.305 5.981
$atom:wat[1][0][1]/H1 $mol:wat[1][0][1] @atom:SPCE/H 0.4238 3.391 1.305 6.970
$atom:wat[1][0][1]/H2 $mol:wat[1][0][1] @atom:SPCE/H 0.4238 4.200 1.772 5.652
$atom:wat[0][0][0]/O $mol:wat[0][0][0] @atom:SPCE/O -0.8476 1.131 2.611 2.300
$atom:wat[0][0][0]/H1 $mol:wat[0][0][0] @atom:SPCE/H 0.4238 1.131 2.611 3.289
$atom:wat[0][0][0]/H2 $mol:wat[0][0][0] @atom:SPCE/H 0.4238 0.320 2.143 1.971
$atom:wat[0][0][1]/O $mol:wat[0][0][1] @atom:SPCE/O -0.8476 1.131 2.611 5.061
$atom:wat[0][0][1]/H1 $mol:wat[0][0][1] @atom:SPCE/H 0.4238 1.940 2.143 5.391
$atom:wat[0][0][1]/H2 $mol:wat[0][0][1] @atom:SPCE/H 0.4238 1.131 3.546 5.391
$atom:wat[0][1][0]/O $mol:wat[0][1][0] @atom:SPCE/O -0.8476 1.131 5.221 1.381
$atom:wat[0][1][0]/H1 $mol:wat[0][1][0] @atom:SPCE/H 0.4238 1.131 4.286 1.710
$atom:wat[0][1][0]/H2 $mol:wat[0][1][0] @atom:SPCE/H 0.4238 0.320 5.688 1.710
$atom:wat[0][1][1]/O $mol:wat[0][1][1] @atom:SPCE/O -0.8476 1.131 5.221 5.981
$atom:wat[0][1][1]/H1 $mol:wat[0][1][1] @atom:SPCE/H 0.4238 1.131 5.221 6.970
$atom:wat[0][1][1]/H2 $mol:wat[0][1][1] @atom:SPCE/H 0.4238 1.940 5.688 5.652
$atom:wat[1][1][0]/O $mol:wat[1][1][0] @atom:SPCE/O -0.8476 3.391 6.526 2.300
$atom:wat[1][1][0]/H1 $mol:wat[1][1][0] @atom:SPCE/H 0.4238 3.391 6.526 3.289
$atom:wat[1][1][0]/H2 $mol:wat[1][1][0] @atom:SPCE/H 0.4238 2.582 6.058 1.971
$atom:wat[1][1][1]/O $mol:wat[1][1][1] @atom:SPCE/O -0.8476 3.391 6.526 5.061
$atom:wat[1][1][1]/H1 $mol:wat[1][1][1] @atom:SPCE/H 0.4238 4.200 6.058 5.391
$atom:wat[1][1][1]/H2 $mol:wat[1][1][1] @atom:SPCE/H 0.4238 3.391 7.462 5.391
}
} # IceRect8
# Credit goes to Martin Chaplin.
# These coordinates were orignally downloaded from Martin Chaplin's
# website: http://www.btinternet.com/~martin.chaplin/ice1h.html
# ... and then they were stretched independently in the xy and z
# directions in order to match the lattice parameters measured by
# Rottger et al.,
# "Lattice constants and thermal expansion of H2O and D2O ice Ih"
# between 10 and 265K", Acta Crystallogr. B, 50 (1994) 644-648
# I am using the lattice constants measured at temperature 265K
# (and pressure=100Torr).

View File

@ -0,0 +1,123 @@
#############################################################
# WARNING: THIS FILE HAS NOT BEEN TESTED!
# (If you use this file in a simulation, please email me to let me know
# if it worked. -Andrew 2012-10, (jewett dot aij at gmail dot com))
#########################################################
# There are two different versions of TIP3P:
#
# tip3p_1983.lt # The implementation of TIP3P used by CHARMM (I think).
# tip3p_2004.lt # The newer Price & Brooks, J. Chem Phys 2004 model
# # which uses long-range coulombics
#########################################################
# file "tip3p_1983.lt"
#
# H1 H2
# \ /
# O
#
# I think this is the TIP3P water model used by CHARMM (and AMBER?). See:
# Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem Phys, 79, 926 (1983)
TIP3P_1983 {
write_once("In Init") {
# -- Default styles (for solo "TIP3P_1983" water) --
units real
atom_style full
# I'm not sure exactly which cutoffs distances are traditionally used in
# in the 1983 "TIP3P" model by Jorgensen model, (used by CHARMM).
# (See the Price JCP 2004 paper for a review.)
# My first guess was this:
# pair_style hybrid lj/charmm/coul/charmm 7.5 8.0 10.0 10.5
# However, in the LAMMPS "peptide" example, they use these parameters:
pair_style hybrid lj/charmm/coul/long 8.0 10.0 10.0
# ...alternately, if this does not work, try this:
# pair_style hybrid lj/charmm/coul/long 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/long 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. (?? I think ??)
# 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 fShakeTIP3P tip3p shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fShakeTIP3P during minimization.)
}
} # "TIP3P_1983" water molecule type
# (note to self:)
# In the LAMMPS "peptide" example, these (nearly identical) parameters were used
# and they left the O-H parameters to be determined by the default mixing rules
#pair_style lj/charmm/coul/long 8.0 10.0 10.0
#pair_coeff @atom:H @atom:H 0.046 0.400014 0.046 0.400014
#pair_coeff @atom:O @atom:O 0.1521 3.15057 0.1521 3.15057
#angle_style charmm
#angle_coeff @angle:HOH 55.0 104.52 0.0 0.0

View File

@ -0,0 +1,97 @@
#########################################################
# Warning: This file has not been carefully tested.
#########################################################
# There are two different versions of TIP3P:
#
# tip3p_1983.lt # The implementation of TIP3P used by CHARMM (I think).
# tip3p_2004.lt # The newer Price & Brooks, J. Chem Phys 2004 model
# # which uses long-range coulombics
#########################################################
# file "tip3p_2004.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.
TIP3P_2004 {
write_once("In Init") {
# -- Default styles (for solo "TIP3P_2004" 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:ow -0.830 0.0000000 0.00000 0.000000
$atom:H1 $mol:. @atom:hw 0.415 0.756950327 0.00000 0.5858822766
$atom:H2 $mol:. @atom:hw 0.415 -0.756950327 0.00000 0.5858822766
}
write_once("Data Masses") {
@atom:ow 15.9994
@atom:hw 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:ow @atom:ow lj/charmm/coul/charmm 0.1521 3.1507
# For the 2004 Price & Brooks version of TIP3P use:
pair_coeff @atom:ow @atom:ow lj/charmm/coul/long 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:hw @atom:hw lj/charmm/coul/long 0.0460 0.4000
pair_coeff @atom:ow @atom:hw lj/charmm/coul/long 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:hw @atom:hw lj/charmm/coul/charmm 0.0460 0.4000
# pair_coeff @atom:ow @atom:hw 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:hw @atom:hw lj/charmm/coul/charmm 0.00 0.4000
# pair_coeff @atom:ow @atom:hw lj/charmm/coul/charmm 0.00 1.7753
#########################################################################
# Define a group for the tip3p water molecules:
group tip3p type @atom:ow @atom:hw
# Optional: Constrain the angles and distances.
# (Most implementations use this, but it is optional.)
fix fShakeTIP3P tip3p shake 0.0001 10 100 b @bond:OH a @angle:HOH
# (Remember to "unfix" fShakeTIP3P during minimization.)
}
} # "TIP3P_2004" water molecule type

View File

@ -0,0 +1,53 @@
# 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("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
}
# Optional: Create a group corresponding to atoms used by the TraPPE force-
# field. (This is useful if you mix force-fields together.)
write_once("In Settings") {
group TraPPE type @atom:CH2 @atom:CH3 @atom:CH4
}
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
}
} # class TraPPE

View File

@ -0,0 +1,54 @@
# This file stores LAMMPS data for the "mW" water model.
# (Molinero, V. and Moore, E.B., J. Phys. Chem. B 2009, 113, 4008-4016)
#
# In this model, each water molecule is represented by a single "mW" particle.
# These particles interact with their neighbors via 3-body Stillinger-Weber
# forces whose parameters are tuned to mimic directional hydrogen-bonding
# in liquid water (as well as hexagonal ice, type II ice, and
# low-density super-cooled liquid/amorphous water phases).
WatMW {
write("Data Atoms") {
$atom:mW $mol:. @atom:mW 0.0 0.0 0.0 0.0
}
write_once("Data Masses") {
@atom:mW 18.02
}
write_once("system.in.sw") {
mW mW mW 6.189 2.3925 1.8 23.15 1.2 -0.333333333 7.049556277 0.602224558 4 0 0
}
write_once("In Init") {
# -- Default styles for "WatMW" --
units real
pair_style sw
}
write_once("In Settings") {
# --Now indicate which atom type(s) are simulated using the "sw" pair style
# -- In this case only one of the atom types is used (the mW water "atom").
pair_coeff * * sw system.in.sw mW NULL NULL NULL
# -- Unfortunately LAMMPS itself does not understand molemlate syntax, so
# -- the atoms are identified by order in the list, not by name. (The "mW"
# -- refers to to an identifier in the system.in.sw file, not watmw.lt.)
# -- This command says that the first atom type corresponds to the "mW"
# -- atom in system.in.sw, and to ignore the remaining three atom types
# -- (correspond to the CH2, CH3, CH4 atom types defined in trappe1998.lt.
# -- We don't want to use the "sw" force field for interactions involving
# -- these atom types, so we put "NULL" there.)
# -- Note: For this to work, you should probably run moltemplate this way:
# -- moltemplate.sh -a "@atom:WatMW/mW 1" system.lt
# -- This assigns the atom type named @atom:WatMW/mW to 1 (the first atom)
}
# -- optional --
write_once("In Settings") {
group WatMW type @atom:mW #(Atoms of this type belong to the "WatMW" group)
}
} # WatMW

View File

@ -0,0 +1,70 @@
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 DLPC model is a truncated version of DPPC. Modifications discussed below.)
The protein model is described here:
G. Bellesia, AI Jewett, and J-E Shea,
Protein Science, Vol19 141-154 (2010)
--- PREREQUISITES: ---
1) This example requires the "dihedral_style fourier", which is currently
in the USER-MISC package. Build LAMMPS with this package enabled using
make yes-user-misc
before compiling LAMMPS.
(See http://lammps.sandia.gov/doc/Section_start.html#start_3 for details.)
2) This example may require additional features to be added to LAMMPS.
If LAMMPS complains about an "Invalid pair_style", then
a) download the "additional_lammps_code" from
http://moltemplate.org (upper-left corner menu)
b) unpack it
c) copy the .cpp and .h files to the src folding of your lammps installation.
d) (re)compile 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.
The bending stiffness of each lipid has been increased to compensate
for the additional disorder resulting from mixing two different types
of lipids together. (Otherwise pores appear.)
Unlike the original "DPPC" molecule model, the new "DPPC" and "DLPC" models
have 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

View File

@ -0,0 +1,33 @@
# --- Running LAMMPS ---
# -- Prerequisites: --
# The "run.in.nvt" file is a LAMMPS input script containing
# references to the input scripts and data files
# you hopefully have created earlier with moltemplate.sh:
# system.in.init, system.in.settings, system.data, and table_int.dat
# If not, carry out the instructions in "README_setup.sh".
#
# -- Instructions: --
# If "lmp_linux" is the name of the command you use to invoke lammps,
# then you would run lammps on these files this way:
lmp_linux -i run.in.npt # Run a simulation at constant pressure (tension)
#or
lmp_linux -i run.in.nvt # Run a simulation at constant volume
#(Note: The constant volume simulation lacks pressure equilibration. These are
# completely separate simulations. The results of the constant pressure
# simulation are ignored when beginning the simulation at constant volume.
# This can be fixed. Read "run.in.nvt" for equilibration instructions.)
# If you have compiled the MPI version of lammps, you can run lammps in parallel
#mpirun -np 4 lmp_linux -i run.in.npt
#or
#mpirun -np 4 lmp_linux -i run.in.nvt
# (assuming you have 4 processors available)

View File

@ -0,0 +1,28 @@
# 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 ../
# The "table_int.dat" file contains tabular data for the lipid INT-INT atom
# 1/r^2 interaction. We need it too. (This slows down the simulation by x2,
# so I might look for a way to get rid of it later.)
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 ../

View File

@ -0,0 +1,87 @@
------- To view a lammps trajectory in VMD --------
1) Build a PSF file for use in viewing with VMD.
This step works with VMD 1.9 and topotools 1.2.
(Older versions, like VMD 1.8.6, don't support this.)
a) Start VMD
b) Menu Extensions->Tk Console
c) Enter:
(I assume that the the DATA file is called "system.data")
topo readlammpsdata system.data full
animate write psf system.psf
2)
Later, to Load a trajectory in VMD:
Start VMD
Select menu: File->New Molecule
-Browse to select the PSF file you created above, and load it.
(Don't close the window yet.)
-Browse to select the trajectory file.
If necessary, for "file type" select: "LAMMPS Trajectory"
Load it.
---- A note on trajectory format: -----
If the trajectory is a DUMP file, then make sure the it contains the
information you need for pbctools (see below. I've been using this
command in my LAMMPS scripts to create the trajectories:
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
It's a good idea to use an atom_style which supports molecule-ID numbers
so that you can assign a molecule-ID number to each atom. (I think this
is needed to wrap atom coordinates without breaking molecules in half.)
Of course, you don't have to save your trajectories in DUMP format,
(other formats like DCD work fine) I just mention dump files
because these are the files I'm familiar with.
3) ----- Wrap the coordinates to the unit cell
(without cutting the molecules in half)
a) Start VMD
b) Load the trajectory in VMD (see above)
c) Menu Extensions->Tk Console
d) Try entering these commands:
pbc wrap -compound res -all
pbc box
----- Optional ----
Sometimes the solvent or membrane obscures the view of the solute.
It can help to shift the location of the periodic boundary box
To shift the box in the y direction (for example) do this:
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
pbc box -shiftcenterrel {0.0 0.15 0.0}
Distances are measured in units of box-length fractions, not Angstroms.
Alternately if you have a solute whose atoms are all of type 1,
then you can also try this to center the box around it:
pbc wrap -sel type=1 -all -centersel type=2 -center com
4)
You should check if your periodic boundary conditions are too small.
To do that:
select Graphics->Representations menu option
click on the "Periodic" tab, and
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
5) Optional: If you like, change the atom types in the PSF file so
that VMD recognizes the atom types, use something like:
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
(If you do this, it might effect step 2 above.)

Binary file not shown.

After

Width:  |  Height:  |  Size: 83 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 69 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 50 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 66 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 69 KiB

View File

@ -0,0 +1,254 @@
# 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
# -------- Minor coordinates adjustment: -----------
# Those coordinates in the commands above are a little too large.
# To make it easier to type them in, I was using sigma=6.0 Angstroms.
# Instead, here I'll try using sigma=5.5 Angstroms. 5.5/6 = 0.916667)
4HelixInsideOut.scale(0.9166666666666666)
# Note: "scale()" only effects the initial coordinates of
# the molecule, not the force field parameters.
# (If you plan to minimize the molecule, you don't need to
# be so careful about the initial coordinates. In that case,
# you don't have worry about "scale()". Feel free to remove.)
# -------------- 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 fourier
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 5.5 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") {
pair_coeff @atom:sH @atom:sH lj/charmm/coul/charmm/inter 1.8065518 5.5 1 -1
pair_coeff @atom:sL @atom:sL lj/charmm/coul/charmm/inter 1.8065518 5.5 1 0
pair_coeff @atom:tN @atom:tN lj/charmm/coul/charmm/inter 1.8065518 5.5 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 5.5
}
# We use the same bond-angle forces whenever
# there are 3 consecutively-bonded atoms:
#
# 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
}
# Dihedral interactions are also determined by atom-type name in this case.
# I chose atoms whose type-names begin with "t" to be "turn" atoms.
# Atoms whose type-names begin with "s" could be either helices or sheets.
# (In this case, helices. In this example, we use the @dihedral:delta65_0
# parameters. This corresponds to the "AUF2" model from the
# Bellesia et. al 2010 paper.)
# dihedralType atomtypes1 2 3 4 bondtypes1 2 3
write_once("Data Dihedrals By Type") {
@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:* * * *
}
# From the Bellesia et al 2010 paper:
# for helices: U_{dih}(\phi) = 1.2*(cos(3\phi) + cos(\phi+\delta))
# for turns: U_{dih}(\phi) = 0.2*cos(3\phi)
#
# General formula used for "dihedral_style fourier":
# U_{dih}(\phi) = \Sum_{i=1}^m K_i [ 1.0 + cos(n_i \phi - d_i) ]
#
# Syntax:
# dihedralType dihedralstyle m K_1 n_1 d_1 K2 n_2 d_2 ...
write_once("In Settings") {
dihedral_coeff @dihedral:delta60_0 fourier 2 2.167862 3 0 2.167862 1 -60.0
dihedral_coeff @dihedral:delta62_5 fourier 2 2.167862 3 0 2.167862 1 -62.5
dihedral_coeff @dihedral:delta65_0 fourier 2 2.167862 3 0 2.167862 1 -65.0
dihedral_coeff @dihedral:turn fourier 1 0.361310 3 0
# Note: 2.167862=1.2*epsilon and 0.361310=0.2*epsilon.
}
# --- 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)

View File

@ -0,0 +1,225 @@
### 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
# -------- Minor coordinates adjustment: -----------
# Those coordinates in the commands above are a little too large.
# To make it easier to type them in, I was using sigma=6.0 Angstroms.
# Instead, here I'll try using sigma=5.5 Angstroms. 5.5/6 = 0.916667)
4SheetBarrel.scale(0.9166666666666666)
4HelixBundle.scale(0.9166666666666666)
4HelixBundleLoose.scale(0.9166666666666666)
4HelixInsideOutLoose.scale(0.9166666666666666)
4HelixPlanar.scale(0.9166666666666666)
# Note: "scale()" only effects the initial coordinates of
# the molecule, not the force field parameters.
# (If you plan to minimize the molecule, you don't need to
# be so careful about the initial coordinates. In that case,
# you don't have worry about "scale()". Feel free to remove.)
} # 1beadProtSci2010 (namespace)

View File

@ -0,0 +1,198 @@
# Note:
#
# This example may require additional features to be added to LAMMPS. If
# LAMMPS complains about an "Invalid pair_style", then download copy the
# "additional_lammps_code" from moltemplate.org, unpack it 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.
#
# Units:
#
# 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.
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 # <- used in the original article
angle_style hybrid harmonic # <- prevents unphysical acute angle turns
# Explanation:
# angle_style cosine/delta: U(theta) = k*(1-cos(theta-theta0))
# angle_style harmonic: U(theta) = k*(theta-theta0)^2
dihedral_style none
improper_style none
pair_style hybrid table linear 1130 &
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"
neighbor 2.0 multi # <- perhaps unnecessary
communicate multi # <- perhaps unnecessary
}
DPPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 33.75 # DPPC 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 # DLPC 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
# Particles and properties 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 may override it later.
}
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") {
# cosine/delta: U(theta) = k*(1-cos(theta-theta0))
#angle_coeff @angle:backbone cosine/delta 4.60086042 180 #<-- 7*eps
# harmonic: U(theta) = k*(theta-theta0)^2 not (k/2)*(theta-theta0)^2
angle_coeff @angle:backbone harmonic 9.85898661 180 #<-->30*eps
}
# I use a stiffer bond-angle than the original Brannigan & Brown 2005 paper
# to attempt to compensate for the fact that here we are using a lipid
# mixture of DPPC and DLPC. (The mixture of lipids introduces a great deal
# of disorder into the bilayer which would not be present in a DPPC bilayer.
# This causes pores to form. Increasing the angle stiffness prevents this.)
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 the default.) Using this convention makes it
# easier to mix this coarse-grained lipid model with other molecular models.
} # CGLipidBr2005
# Note: This example has not been optimized for speed.
#
# 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.)

View File

@ -0,0 +1,29 @@
#!/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)
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 = 0.02
Rmax = 22.6
rcut = 22.5
N = 1130
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))

View File

@ -0,0 +1,70 @@
#!/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 = 0.02
Rmax = 22.6
Rc1 = 22.0
Rc2 = 22.5
N = 1130
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))

View File

@ -0,0 +1,207 @@
# 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
# 13 rows (each row is aligned with the x-axis)
# 15 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.
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][*]
# 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 protein and lipids ----------
# -----------------------------------------------------------
# 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.
# We want the hydrophobic interactions between hydrophobic residues in
# the protein and beads the interior of the lipid to be energetically
# similar to the attractive interactions between the lipid tails.
#
# Note: I made the width of the outward-facing protein beads slightly larger
# ("12.5") whenever they interact with the "tail" beads in each lipid
# (in order to make the protein wider there).
# This hopefully relieves some of the internal negative pressure in the center
# of the bilayer which can otherwise rip apart the protein or suck it into
# the bilaer. (This is a hack, and I'm not sure if it is necessary.
# For different protein or lipid models, you probably don't need this.)
#
# i j pairstylename eps sig K L
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 0.1643 12.5 0.4 -1
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
# To help keep the protein from tilting 90 degrees and burying itself
# within the lipid bilayer, we make the turn regions at either
# end of the protein (strongly) attracted to the head groups
# of the lipid. (In reality, they would probably be attracted
# to the water as well.)
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 5.5 1 -1
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 1.8065518 5.5 1 -1
# All other interactions between proteins and lipids are steric.
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/tail @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/int @atom:1beadProtSci2010/tN lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/DPPC/head @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/DLPC/head @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
pair_coeff @atom:CGLipidBr2005/DLPC/head @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
# -----------------------------------------------------------
# -------- Modifications to the protein model: --------------
# -----------------------------------------------------------
#
# Turn off attraction between the hydrophobic "@atom:sH" beads:
# (These beads are located in the outside of a trans-membrane protein.)
pair_coeff @atom:1beadProtSci2010/sH @atom:1beadProtSci2010/sH lj/charmm/coul/charmm/inter 1.8065518 5.5 1 0
# (Why: These beads are only attracted to
# each other in an aqueous environment)
# ... and
# Turn ON attraction between the hydrophilic "@atom:sL" beads.
# (These beads are located in the interior of a trans-membrane protein.)
pair_coeff @atom:1beadProtSci2010/sL @atom:1beadProtSci2010/sL lj/charmm/coul/charmm/inter 1.8065518 5.5 1 -1
# Why?
# In reality, polar groups in the interior of trans-membrane
# proteins do form hydrogen bonds with each other. This was
# absent from the original protein model because, in an aqueous
# environment, these groups preferentially interact with the water.
#
# Why is this necessary?
# Shouldn't attraction between lipid tails and the protein create
# an effective force which brings the hydrophilic beads together?
# (similar to the hydrophobic effect, but in reverse?).
# Answer:
# Unlike an aqueous environment (~zero pressure, or +1atm), there is
# a large negative pressure in the interior of some bilayer membrane
# models (such as this one). Without some kind of direct attraction
# between interior residues, the protein will get pulled apart.
# (Perhaps the attractive force I am using is too strong?)
}
# 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 fourier
improper_style none
pair_style hybrid table linear 1130 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"
}

View File

@ -0,0 +1,35 @@
# -------- 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.
# So, after typing "make yes-user-misc" 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 --
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

View File

@ -0,0 +1,65 @@
# PREREQUISITES:
#
# You must use moltemplate.sh to create 3 files:
# system.data system.in.init system.in.settings
# (Follow the instructions in README_setup.sh, or run it using ./README_sh.)
#
# -------- 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.
# So, after typing "make yes-user-misc" 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 --
#minimize 1.0e-5 1.0e-7 500 2000
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 couple xy
# 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
run 100000
write_data system_after_npt.data

View File

@ -0,0 +1,71 @@
# PREREQUISITES:
#
# 1) You must use moltemplate.sh to create 3 files:
# system.data system.in.init system.in.settings
# (Follow the instructions in README_setup.sh, or run it using ./README_sh.)
# 2) You must equilibrate the system beforehand using "run.in.npt".
# This will create the file "system_after_npt.data" which this file reads.
# (Note: I have not verified that this equilibration protocol works well.)
#
# -------- LAMMPS 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.
# So, after typing "make yes-user-misc" 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.
#
# ------------------------------- Initialization Section --------------------
include system.in.init
# ------------------------------- Atom Definition Section -------------------
# Read the coordinates generated by an earlier NPT simulation
read_data system_after_npt.data
# (The "write_restart" and "read_restart" commands were buggy in 2012,
# but they should work also. I prefer "write_data" and "read_data".)
# ------------------------------- Settings Section --------------------------
include system.in.settings
# ------------------------------- Run Section -------------------------------
timestep 10.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 10000 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 500000
run 10000000
write_data system_after_nvt.data
# (The "write_restart" and "read_restart" commands were buggy in 2012,
# but they should work also.)

View File

@ -0,0 +1,33 @@
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.
-------------
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

View File

@ -0,0 +1,33 @@
# --- Running LAMMPS ---
# -- Prerequisites: --
# The "run.in.nvt" file is a LAMMPS input script containing
# references to the input scripts and data files
# you hopefully have created earlier with moltemplate.sh:
# system.in.init, system.in.settings, system.data, and table_int.dat
# If not, carry out the instructions in "README_setup.sh".
#
# -- Instructions: --
# If "lmp_linux" is the name of the command you use to invoke lammps,
# then you would run lammps on these files this way:
lmp_linux -i run.in.npt # Run a simulation at constant pressure (tension)
#or
lmp_linux -i run.in.nvt # Run a simulation at constant volume
#(Note: The constant volume simulation lacks pressure equilibration. These are
# completely separate simulations. The results of the constant pressure
# simulation are ignored when beginning the simulation at constant volume.
# This can be fixed. Read "run.in.nvt" for equilibration instructions.)
# If you have compiled the MPI version of lammps, you can run lammps in parallel
#mpirun -np 4 lmp_linux -i run.in.npt
#or
#mpirun -np 4 lmp_linux -i run.in.nvt
# (assuming you have 4 processors available)

View File

@ -0,0 +1,28 @@
# 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 ../
# The "table_int.dat" file contains tabular data for the lipid INT-INT atom
# 1/r^2 interaction. We need it too. (This slows down the simulation by x2,
# so I might look for a way to get rid of it later.)
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 ../

View File

@ -0,0 +1,87 @@
------- To view a lammps trajectory in VMD --------
1) Build a PSF file for use in viewing with VMD.
This step works with VMD 1.9 and topotools 1.2.
(Older versions, like VMD 1.8.6, don't support this.)
a) Start VMD
b) Menu Extensions->Tk Console
c) Enter:
(I assume that the the DATA file is called "system.data")
topo readlammpsdata system.data full
animate write psf system.psf
2)
Later, to Load a trajectory in VMD:
Start VMD
Select menu: File->New Molecule
-Browse to select the PSF file you created above, and load it.
(Don't close the window yet.)
-Browse to select the trajectory file.
If necessary, for "file type" select: "LAMMPS Trajectory"
Load it.
---- A note on trajectory format: -----
If the trajectory is a DUMP file, then make sure the it contains the
information you need for pbctools (see below. I've been using this
command in my LAMMPS scripts to create the trajectories:
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
It's a good idea to use an atom_style which supports molecule-ID numbers
so that you can assign a molecule-ID number to each atom. (I think this
is needed to wrap atom coordinates without breaking molecules in half.)
Of course, you don't have to save your trajectories in DUMP format,
(other formats like DCD work fine) I just mention dump files
because these are the files I'm familiar with.
3) ----- Wrap the coordinates to the unit cell
(without cutting the molecules in half)
a) Start VMD
b) Load the trajectory in VMD (see above)
c) Menu Extensions->Tk Console
d) Try entering these commands:
pbc wrap -compound res -all
pbc box
----- Optional ----
Sometimes the solvent or membrane obscures the view of the solute.
It can help to shift the location of the periodic boundary box
To shift the box in the y direction (for example) do this:
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
pbc box -shiftcenterrel {0.0 0.15 0.0}
Distances are measured in units of box-length fractions, not Angstroms.
Alternately if you have a solute whose atoms are all of type 1,
then you can also try this to center the box around it:
pbc wrap -sel type=1 -all -centersel type=2 -center com
4)
You should check if your periodic boundary conditions are too small.
To do that:
select Graphics->Representations menu option
click on the "Periodic" tab, and
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
5) Optional: If you like, change the atom types in the PSF file so
that VMD recognizes the atom types, use something like:
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
(If you do this, it might effect step 2 above.)

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.8 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 57 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 5.7 KiB

View File

@ -0,0 +1,191 @@
# Note:
#
# This example may require additional features to be added to LAMMPS. If
# LAMMPS complains about an "Invalid pair_style", then download copy the
# "additional_lammps_code" from moltemplate.org, unpack it 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.
#
# --- DLPC lipids ---
# A truncated version of the DPPC 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.
# (You may need to stiffen the bond-angle forces to make it behave correctly,
# but I did not do this here.)
#
# Units:
#
# 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.
#
#
# The new DLPC model is a truncated version of DPPC,
# (Its behaviour has not been rigorously tested.)
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 1130 &
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"
neighbor 2.0 multi # <- perhaps unnecessary
communicate multi # <- perhaps unnecessary
}
DPPC {
write("Data Atoms") {
$atom:h $mol:. @atom:head 0.0 0.00 0.00 33.75 # DPPC 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 # DLPC 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
# 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 may override it later.
}
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 the default.) Using this convention makes it
# easier to mix this coarse-grained lipid model with other molecular models.
} # CGLipidBr2005
# Note: This example has not been optimized for speed.
#
# 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.)

View File

@ -0,0 +1,29 @@
#!/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)
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 = 0.02
Rmax = 22.6
rcut = 22.5
N = 1130
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))

View File

@ -0,0 +1,70 @@
#!/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 = 0.02
Rmax = 22.6
Rc1 = 22.0
Rc2 = 22.5
N = 1130
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))

View File

@ -0,0 +1,94 @@
# Description:
# This constructs a bilayer constructed from coarse-grained DPPC lipids
# (implicit solvent). The DPPC lipid model is described here:
# G. Brannigan, P.F. Philips, and F.L.H. Brown,
# Physical Review E, Vol 72, 011915 (2005)
#
# NOTE: There is an example of a 50%/50% DPPC & DLPC mixture
# in the "membrane+protein" and "vesicle" examples.
import "CGLipidBr2005Orig.lt"
using namespace CGLipidBr2005
lipids = new DPPC [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. We apply them below:
# width: 240 = 32*7.5
# height: 240.322 = 37*6.49519
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) If you want to create a bilayer from a mixture of DPPC and DLPC, you 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)
#

View File

@ -0,0 +1,28 @@
# -- 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

View File

@ -0,0 +1,57 @@
# -------- 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 6.0 # The time-step in Watson et. al 2011 was 0.002*3ps = 6fs
dump 1 all custom 5000 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 couple xy
# 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

View File

@ -0,0 +1,53 @@
# ------------------------------- Initialization Section --------------------
include system.in.init
# ------------------------------- Atom Definition Section -------------------
# Normally, I would minimize the system and equilibrate the system at constant
# pressure and temperature beforehand. If you run lammps with "run.in.npt",
# it will generate a restart file "system_after_npt.rst" with reasonable
# coordinates at that temperature and pressure. Then we could load it now:
#
#read_restart system_after_npt.rst
#
# Unfortunately the LAMMPS "read_restart" command has been undependable over
# the past year (2012), and I feel it is safer to remove it from the examples.
# Instead, for this example, I just read the raw coordinates generated by
# moltemplate (and the default volume). (I get fewer questions this way.)
# However you should never run any liquid simulations at constant volume without
# pressure equilibration first. Hopefully in the future "read_restart" will
# work. Until then, try "read_dump", "dump2data.py", or "restart2data".
read_data system.data
# ------------------------------- 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

View File

@ -0,0 +1,47 @@
This is an implementation of the "two-stage" model used by Maxim Imakaev
in the Naumova et Al 2013 Science paper on metaphase chromatin.
(Download the supplemental materials section and scroll down to the section:
"Two-stage process: linear compaction - axial compression")
---- SMALL MODIFICATION ----
Unlike that study, I did not use "softened" Lennard-Jones potentials
(which allow the chains to pass through each other).
--- Why use moltemplate? ---
Honestly, you don't need to use moltemplate to build this polymer.
It is almost counter-productive to use moltemplate to build this kind of
polymer because it is so simple. (The polymer has only 1 bead per atom.
It just makes it more complicated to introduce all these extra
files including monomer.lt, condensin.lt and system.lt, especially considering
that system.lt is a complex file which is generated by a separate script.)
However building the sytem using moltemplate may pay off if you
replace each point-like monomer with a multi-atom molecule later on.
(Right now, using moltemplate to build this system is sort of overkill.
I'll post an example of building more complex models of chromatin eventually.)
---- 30-nm fiber model: ----
Anyway, the two-stage model at the end of Naumova et al Science 2013 uses the "30nm-fiber" model, whose details are (somewhat vaguely) described in the supplemental materials section.
For the 10nm model,
n=128000,
L=200,
U(alpha)=5*(1 - cos(alpha))
bond_length=1.0 (=10nm)
sigma=1.0 (particle radius = 10nm)
30nm-fiber model details:
"The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level."
I interpret this to mean that, for the 30nm model,
n=128000/4.25~=30117 (however I rounded up to 32768=2^15)
L=200/4.25~=47 (however I rounded up to 51)
U(alpha)=1.17647*(1 - cos(alpha)) (5/4.25=1.17647)
To increase the volume by a factor o 4.25, I increase both the diameter of each
bead (the "sigma" parameter), and the bond-lengths connecting them from
1.0 (corresponding to 10nm) to 4.25^(1/3)~=1.6198 (corresponding to 16.198nm).

View File

@ -0,0 +1,29 @@
---------------
The average diameter of a mammalian cell nucleus is is 6 micrometers (μm),
which occupies about 10% of the total cell volume.
(See "Molecular Biology of the Cell, Chapter 4, pages 191234 (4th ed.)", by
Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, Peter Walter, 2002)
... of that, 25% of it is occupied by the nucleolus
http://en.wikipedia.org/wiki/Nucleolus
("citation needed")
---------------
From the supplemental material for the original HiC paper
(Lieberman-Aiden et al., Science 2009)
Appendix 1.
Estimate of the Volume Fraction of Chromatin in Human Cells
In the simulations we sought to obtain an ensemble of structures that, in their statistical properties, resemble some of the features of chromatin arrangement in the cell. Below we demonstrate that chromatin occupies a significant fraction of cell volume, a property that we reproduced in simulations. Taking the nuclear diameter of a tissue culture cell to be 5-10um, and assuming close to a spherical shape we obtain the volume in the range 50-500 um^3, with a (geometric) mean of ~160 um^3. If we assume that the chromatin is built of DNA wrapped around nucleosomes, then we have 6x10^9bp/200bp=3x10^7 nucleosomes. Each may be approximated as a cylinder ~10nm in diameter and ~5nm in height, suggesting a volume of about 500nm3 each. The linker DNA after each nucleosome is about 50bps long, suggesting a volume of about 50*0.34nm*3.14*1nm^2=50nm^3. Thus the total volume of chromatin = 550x3x10^7 =16 um^3, or ~10% (3-23%) of the nuclear volume. This strikingly large volume fraction is itself a significant underestimate, since we ignored, among other things, all other DNA-bound proteins. Note that any further packing or localization of chromatin inside the nucleus will increase local density.
---- This next section mostly only justifies why they ----
---- they did not stop the simulation when the globules ----
---- were fully crumpled (ie with uniform density) ----
In our simulations, the radius of the final crumpled globule was R≈12.5 and the volume V≈8000 cubic units. The total volume of the 4000 monomers, 1 unit in diameters each, is V≈2000. This implies a volume fraction of about 25%, which is consistent with the volume fraction estimated above.
---- ----
Appendix 2.
Monomer length in base pairs
Each monomer of the chain corresponds to a fragment of chromatin that equals the Kuhn length of the chromatin fiber, i.e. approximately twice the persistence length of the fiber. Although the persistence length of the chromatin fiber is unknown it can be estimated using the following arguments. DNA is packed into nucleosomes, where 150 bps are wrapped around the histone core and do not contribute to flexibility of the fiber. The linker DNA of about 50 bps that connects consecutive nucleosomes is bendable, and is the source of flexibility in the fiber. Since the persistence length of double-stranded DNA is 150 bps, an equally flexible region of the nucleosomal DNA should contain 3 linkers, i.e. 3 consecutive nucleosomes packing about 600 bps of DNA. The excluded volume of the nucleosomes, nucleosome interactions, and other DNA-bound proteins can make the fiber less flexible or prohibit certain conformation and may tend to increase the persistence length of the fiber. Using this estimated lower bound estimate for the persistence length, we obtain the Kuhn length of the equivalent freely-jointed chain to be 6 nucleosomes, or ~ 1200bp. A simulated chain of 4000 monomers corresponds to 4.8Mb of packed DNA. The size of each monomer was chosen such that its volume is equal to (or slightly above) that of 6 nucleosomes (V=6 x 600 nm^3); thus the radius of the spherical monomer is R=10nm. The diameter of each globule shown in Figure 4 is about 200 nm.

View File

@ -0,0 +1,7 @@
# Run lammps using the following 3 commands:
# (assuming "lmp_linux" is the name of your LAMMPS binary)
lmp_linux -i run.in.min
lmp_linux -i run.in.stage1
lmp_linux -i run.in.stage2

View File

@ -0,0 +1,58 @@
# Use these commands to generate the LAMMPS input script and data file
# (and other auxilliary files):
# Create LAMMPS input files this way:
cd moltemplate_files
# First, rescale and interpolate the positions
# where the monomers will be located. (This step is
# not needed if the coords_orig.raw file already has correct coordinates.)
./interpolate_coords.py 32768 1.6198059 < coords_orig.raw > coords.raw
# Then, build the "system.lt" file
./generate_system_lt.py 32768 51 < coords.raw > system.lt
# 32768 is the number of monomers in the polymer
# (which may be different from the number of coordinates
# in the "coords_orig.raw" file) This number will vary
# depending on how long you want the polymer to be.
# The second argument "51" is the average interval between
# condensin anchors (IE the "loop size" in monomers.)
# Run moltemplate
moltemplate.sh system.lt -a "@bond:stage1 1" \
-a "@bond:stage2 2" \
-a "@atom:Monomer/A 1"
# This will generate various files with names ending in *.in* and *.data.
# These files are the input files directly read by LAMMPS. Move them to
# the parent directory (or wherever you plan to run the simulation).
#
# We used the "-a" command to set the variable @bond:condensin to "2"
# because we will refer to it later in the "run.in" LAMMPS input script.
# (Of coarse, LAMMPS knows nothing about moltemplate variables,
# so in that file we refer to it as dihedral type "1")
mv -f system.in* system.data ../
# We also need the table of bond forces used during "stage 2".
# (Like the system.data and the various input scripts, this file is needed by
# LAMMPS, so we need to copy it to the directory where we will run the sim.)
cp -f table_bonds_stage2.dat ../
# Optional:
# The "./output_ttree/" directory is full of temporary files generated by
# moltemplate. They can be useful for debugging, but are usually thrown away.
rm -rf output_ttree/
# Optional:
# Remove the "system.lt" file created by "generate_system_lt.py"
#rm -f system.lt
cd ../

View File

@ -0,0 +1,131 @@
NOTE: VMD DOES NOT ALLOW YOU TO VISUALIZE SYSTEMS WITH MANY BONDS ATTACHED
TO EACH ATOM. (IF IT DID, THE RESULTS WOULD BE UGLY ANWAY.)
HOWEVER THIS MODEL ATTACHES APPROXIMATELY 60 BONDS TO EACH CONDENSIN ATOM.
IN ORDER TO PULL THE CONDENSIN MONOMERS TOGETHER. YOU MUST DELETE THOSE
BONDS (of type "1" or "2") FROM THE "system.data" FILE BEFORE YOU CARRY
OUT THE COMMANDS BELOW. (...And backup your "system.data" file. You'll need
all the bonds when you run the simulations.)
-------------- COLORS ---------------
In order to show how the polymer is distributed along the length of the
cylinder, I recommend to select the
Graphics->Graphical Representations
menu option, and select "Index" from the "Coloring Method" pull-down menu.
After doing this, you can switch from a red-white-blue scheme, to a
rainbow ("jet") scheme, by selecting the Extensions->Tk Console menu option
and loading the "vmd_colorscale_jet.tcl" file located in the "images" directory.
-------------------------------------------
First, if you have not done so, download and install VMD:
http://www.ks.uiuc.edu/Research/vmd/
http://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD
------- To view a lammps trajectory in VMD --------
The system coordinates are initialy stored in a LAMMPS' ".data" file.
(If that file was built with moltemplate, it will be named "system.data".)
The first step is to view that file.
Then you should create a ".psf" file
(The .psf file is necessary after you run the simulation
for viewing trajectories.)
1) Build a PSF file for use in viewing with VMD
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
You will see a snapshot of the system on the screen.
(presumably the initial conformation at t=0)
2)
Later once you have run a simulation,
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
(It usually has names like "traj.lammpstrj". It depends on how you saved it.)
If necessary, for "file type" select: "LAMMPS Trajectory".
(However VMD should recognize the file type by the file extension.)
Load it.
##################### PERIODIC BOUNDARY CONDITIONS #####################
If you are only simulating a single molecule and you are not
using periodic boundary conditions, then ignore everything below.
########################################################################
---- A note on trajectory format: -----
If the trajectory is the standard LAMMPS format, (aka a "DUMP" file with
a ".lammpstrj" extension), then it's a good idea when you run the simulation
to tell LAMMPS you want to store the information needed for showing periodic
boundary conditions. (Even if you are not using periodic boundaries.
It never hurts to include a tiny bit of extra information.) To do that,
I've been using this command in my LAMMPS scripts to create the trajectories:
dump 1 all custom 5000 traj.lammpstrj id mol type x y z ix iy iz
(Also: 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 visually without breaking molecules in half. Again
you don't need to worry about this if you are not using periodic boundaries.)
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.5 -0.5 -0.5}
pbc box -shiftcenterrel {-0.5 -0.5 -0.5}
Distances are measured in units of box-length fractions, not Angstroms.
Alternately if you have a solute whose atoms are all of type 1,
then you can also try this to center the box around it:
pbc wrap -sel type=1 -all -centersel type=2 -center com
4)
You should check if your periodic boundary conditions are too small.
To do that:
select Graphics->Representations menu option
click on the "Periodic" tab, and
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
5) Optional: If you like, change the atom types in the PSF file so
that VMD recognizes the atom types, use something like:
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
(If you do this, it might effect step 2 above.)

Binary file not shown.

After

Width:  |  Height:  |  Size: 40 KiB

View File

@ -0,0 +1,87 @@
proc lerpcolor { col1 col2 alpha } {
set dc [vecsub $col2 $col1]
set nc [vecadd $col1 [vecscale $dc $alpha]]
return $nc
}
proc coltogs { col } {
foreach {r g b} $col {}
set gray [expr ($r + $g + $b) / 3.0]
return [list $gray $gray $gray]
}
proc jet_tricolor_scale {} {
display update off
set mincolorid [expr [colorinfo num] - 1]
set maxcolorid [expr [colorinfo max] - 1]
set colrange [expr $maxcolorid - $mincolorid]
set colhalf [expr $colrange / 2]
for {set i $mincolorid} {$i < $maxcolorid} {incr i} {
set colpcnt [expr ($i - $mincolorid) / double($colrange)]
# The following color definitions for "jet" sort-of came from:
#http://stackoverflow.com/questions/7706339/grayscale-to-red-green-blue-matlab-jet-color-scale
# but it was missing "green", so I inserted a green somewhere in the middle.
# darkblue/violet 0.0
set color0 { 0.08 0.0 0.77 }
# blue 0.19
set color1 { 0.0 0.0 1.0 }
# cyan 0.34
set color2 { 0.0 1.0 1.0 }
# turquoise 0.4001
set color3 { 0.0 1.0 0.78 }
# green 0.445
set color4 { 0.0 1.0 0.0 }
# chartreuse 0.535
set color5 { 0.875 1.0 0.0 }
# yellow 0.69
set color6 { 1.0 1.0 0.0 }
# orange 0.73
set color7 { 1.0 0.25 0.0 }
# red 0.755
set color8 { 1.0 0.0 0.0 }
# darkred 1.0
set color9 { 0.93 0.0 0.0 }
if { $colpcnt < 0.19 } {
set nc [lerpcolor $color0 $color1 [expr $colpcnt/(0.19-0.0)]]
} elseif { $colpcnt < 0.34 } {
set nc [lerpcolor $color1 $color2 [expr ($colpcnt-0.19)/(0.34-0.19)]]
} elseif { $colpcnt < 0.4001 } {
set nc [lerpcolor $color2 $color3 [expr ($colpcnt-0.34)/(0.4001-0.34)]]
} elseif { $colpcnt < 0.445 } {
set nc [lerpcolor $color2 $color3 [expr ($colpcnt-0.4001)/(0.445-0.4001)]]
} elseif { $colpcnt < 0.535 } {
set nc [lerpcolor $color3 $color4 [expr ($colpcnt-0.445)/(0.535-0.445)]]
} elseif { $colpcnt < 0.69 } {
set nc [lerpcolor $color4 $color5 [expr ($colpcnt-0.535)/(0.69-0.535)]]
} elseif { $colpcnt < 0.73} {
set nc [lerpcolor $color5 $color6 [expr ($colpcnt-0.69)/(0.73-0.69)]]
} elseif { $colpcnt < 0.755} {
set nc [lerpcolor $color6 $color7 [expr ($colpcnt-0.73)/(0.755-0.73)]]
} else {
set nc [lerpcolor $color7 $color8 [expr ($colpcnt-0.755)/(1.0-0.755)]]
}
# set nc [coltogs $nc]
foreach {r g b} $nc {}
puts "index: $i $r $g $b -- $colpcnt"
display update ui
color change rgb $i $r $g $b
}
display update on
}
jet_tricolor_scale

View File

@ -0,0 +1,19 @@
# First, rescale and interpolate the positions
# where the monomers will be located. (This step is not needed
# if the coords_orig.raw file already has correct coordinates.)
# The first argument, 32768 is the number of atoms in the desired file.
# The second argument, 1.6198059 = 4.25^(1/3), tells interpolate_coords.py
# to multiply all the coordinates (scale them up) by 1.6198059.
./interpolate_coords.py 32768 1.6198059 < coords_orig.raw > coords.raw
# Then, build the "system.lt" file
./generate_system_lt.py 32768 51 < coords.raw > system.lt
# 32768 is the number of monomers in the polymer
# (which may be different from the number of coordinates
# in the "coords_orig.raw" file) This number will vary
# depending on how long you want the polymer to be.
# The second argument "51" is the average interval between
# condensin anchors (IE the "loop size" in monomers.)

View File

@ -0,0 +1,57 @@
---- Andrew's comments ----
The two-stage model at the end of Naumova et al Science 2013 uses the "30nm-fiber" model, whose details are (somewhat vaguely) described in the supplemental materials section.
For the 10nm model,
n=128000,
L=200,
U(alpha)=5*(1 - cos(alpha))
bond_length=1.0 (=10nm)
sigma=1.0 (particle radius = 10nm)
30nm-fiber model details:
"The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level."
I interpret this to mean that, for the 30nm model,
n=128000/4.25~=30117
L=200/4.25~=47
U(alpha)=1.17647*(1 - cos(alpha)) (5/4.25=1.17647)
To increase the volume by a factor o 4.25, originally I thought I should increase the "sigma" parameter from 1.0 to 4.25^(1/3)~=1.6198. But I suspect that the bond-lengths between monomers should be fixed at 1.0. If that is the case, then, perhaps I should increase "sigma" from 1.0 to 4.25^(1/2)~=2.061552, and keep the bond-length fixed at 1.0 (which in the units used by thsi paper, corresponds to 10nm). (That would increase the volume of a cylinder of radius "sigma" and length="bond-length" by a factor of 4.25)
bond_length=1.0 (10nm again)
sigma=2.061552 (Yes, this is less than 3.0<-->30nm. See below.)
--- Excerpts from the Supplemental section of Naumova et al Science 2013 ---
From p. 18 of the supplemental materials section of Naumova et al Science 2013.
(This section was probably written by Maxim Imakaev.)
In vivo, the structure of the chromatin fiber can be complicated and many details remain unknown, particularly in metaphase. Given this uncertainty, we simulated chromatin as a homogeneous “beads-on-a-string” polymer fiber. We consider a 10nm fiber, as the pervasiveness of the 30nm fiber in vivo has become increasingly contested. In our simulations, 77Mb is represented by a densely-packed 10nm fiber of 128,000 monomers. Each monomer represents a 10nm-sized DNA-histone complex containing 3 nucleosomes (around 600bp). The fiber has a persistence length of 4 monomers (~2.4Kb), which is based on earlier estimates of 5-10 nucleosomes for interphase (14). Those estimates arise from the assumption that 5-10 linker DNA fragments, each of 20-40bp, can collectively provide flexibility equal to that of the 150bp persistence length of DNA. Binding of proteins to the linker DNA (e.g. histone H1) and interactions between neighboring nucleosomes can further constrain dynamics, requiring more linkers to provide the persistence length. Due to the tight packing of nucleosomes in metaphase, we use the upper limit of this range, i.e. 12 nucleosomes.
For the consecutive loops on a scaffold model (the final folded state model with the best agreement with Hi-C data), we also performed simulations with a more flexible 10nm fiber, or with a 30nm fiber, and found similar results. The more flexible 10nm fiber was modeled by decreasing persistence length to 1.8 monomers. The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level. We note that a classic model of a 30nm fiber is much less dense than a compact metaphase chromosome. A textbook model of a 30nm fiber assumes packing of about 6 nucleosomes per 10nm of fiber length. This model predicts that only 28% of the volume of the fiber (a 30nm-diameter cylinder) is occupied by nucleosomes, assuming a nucleosome shell volume of 328 nm^3. This is much less than the estimated 30-50% density of nucleosomes in a metaphase chromosome, assuming a diameter of 600nm, a packing density of 50-70 Mb/um, and the same nucleosome volume. (See also (15), which gives an estimate of 0.14-0.18 pg/μm for DNA only, and would give about twice the density if DNA is counted with nucleosomes). As follows, these fibers would have to interdigitate, and fill in gaps within each other. We account for this overlap by assuming the effective diameter of the fiber to be less than 30nm. The effective diameter was chosen to make the volume of the fiber equal the total volume of all the nucleosomes.
We accounted for topoisomerase II activity by allowing chromatin fibers to pass through each other, while still having excluded volume interactions. This was achieved by using a soft-core Lennard-Jones potential with 1kT energy cost for monomer overlap (see below). This allows for changes in the topological state of a chromosome that are known to occur during compaction in vivo.
Our simulations of a two-step folding process show that Hi-C data for mitotic chromosomes is consistent with a linearly compressed array of consecutive chromatin loops. Whereas mechanisms for formation of consecutive chromatin loops have been proposed, the process of axial compression is less understood. Chromatid compression cannot be accomplished by increased chromatin-chromatin affinity alone, as this would lead to condensation into a globular geometry (14, 16, 17). However, mechanisms which locally compress the fiber of loop bases naturally allow for anisotropic compression into a shorter and thicker fiber, with the same width regardless of chromosome length (18). Differences in the duration or efficiency of the first and second stages of chromosomal condensation provide a natural mechanism for condensation-related proteins to separately affect mitotic chromosome length and width (19). We also note that the axis of loop-bases in our two-stage model does not necessarily form a continuous and rigid scaffold (Figure S26). As follows, we remain agnostic about the molecular details of the chromosomal scaffold, which might for example be formed by a network consisting of protein-protein and/or protein-DNA interactions (20).
1. Polymer simulations
To perform Langevin dynamics polymer simulations we used OpenMM, a high-performance GPU-assisted molecular dynamics API (21, 22). To represent chromatin fibers as polymers, we used a sequence of spherical monomers of 1 unit of length in diameter. Here and below all distances are measured in monomer sizes, set to be 10nm unless specified otherwise. Neighboring monomers are connected by harmonic bonds, with a potential U = 100*(r - 1)^2 (here and below in units of kT). Polymer stiffness is modeled with a three point interaction term, with the potential U = 5*(1 - cos(alpha)), where alpha is the angle between neighboring bonds. All monomers interact via either a shifted Lennard-Jones (LJ) repulsive potential, or an attractive Lennard-Jones potential. At high densities in a confined volume, the details of the inter-monomer interactions become negligible due to screening (23), and we therefore used the computationally efficient shifted LJ potential. The shifted LJ potential allows for a short-range repulsion by truncating the LJ potential at its minimum and shifting the minimum to zero: U = 4 * (1/r^12 - 1/r^6) + 1, for r<2^(1/6); U=0 for r > 2^(1/6). The shifted LJ potential is one of the most computationally efficient repulsive potentials due to a very short cutoff radius.
To allow chain passing, which represents activity of topoisomerase II, we softened the shifted LJ potential by truncating the interaction energy at Ecutoff = 1 kT. At energies more than 0.5 Ecutoff, the LJ potential was softened via: Usoftened = 0.5 * Ecutoff * (1 + tanh(2*U/Ecutoff - 1)). To avoid numerical 19instabilities in the calculation of U at r ~ 0, the interaction radius r was truncated at r=0.3 via: rtruncated = (r^10 + (0.3)^10)^0.1, which introduced negligible shift in a final softened potential. For an attractive LJ potential, we used: U = 4 * e * (1/r^12 - 1/r^6), with e = 0.46 kT, slightly below the theta-temperature. The attractive potential was similarly softened at 2 kT and cut off at r=2.5. Unless noted, we used a softened shifted LJ repulsive potential.
Polymer models were visualized using Pymol and Rasmol. For images with loop bases highlighted, a base of each loop and 3 monomers surrounding it in each direction were labeled in red.
SECTIONS 2-5 SKIPPED
6. Two-stage process: linear compaction - axial compression
To simulate the two-stage process of metaphase chromosome folding, we used the 30nm fiber representation described above for its computational efficiency. Simulations were initialized from 30000 monomer fractal globule conformations; fractal globule is a model for interphase chromatin organization. First, random consecutive loops with L=100 monomers (see above) were introduced, and anchors of neighboring loops were brought together using harmonic springs with a potential U = k * (r r0)2; r0=0.5. To avoid abrupt motion of the loop anchors, the force was gradually turned on over the first
300000 timesteps, with k linearly increasing in time from 0 to 10 kT. We used softened shifted repulsive LJ potential for inter-monomer interaction.
Upon completion of linear compaction, axial compression was initiated. This involves following changes: the repulsive LJ force is replaced with an attractive LJ force for all monomers, and the chromosomal core of loop anchors is homogeneously compressed. To achieve the latter, all anchor pairs separated by less than 30 anchors were attracted via a potential U = step(d-3) * abs(d-3) * 10 kT, which implements a constant attractive force between two anchors if they are separated by a distance larger than 3. The interactions between neighboring loop anchors were kept throughout this process.
To obtain the contact map from this simulation, 50 independent runs of 1.5e7 timesteps were performed, and 250 conformations were collected from the second half of each run. The contact map was calculated from all conformations of all runs at a 30-monomer resolution, and was further averaged over three 10000-monomer blocks along the diagonal of the heatmap. The latter was done to show contact map at a relevant length scale (0 to 25 Mb), and to achieve a better averaging of the contact map.

View File

@ -0,0 +1,8 @@
for ((i=0; i<60; i++)); do echo "$((i+1)) " `echo "$i*0.05" | bc` 0 0; done
echo 61 3.0 0 -5
for ((i=61; i<=4000; i++)); do echo "$((i+1)) " `echo "$i*0.05" | bc` `echo "($i-60)*0.5"|bc` -10; done

View File

@ -0,0 +1,47 @@
# This file contains the definition of a molecule named "CondensinMonomer".
# (This particular molecule contain only one atom, but that is up to you.)
# Later, multiple CondensinMonomers can be connected together to build a molecule.
CondensinMonomer {
# atom-id mol-id(ignore) atom-type q x y z
write("Data Atoms") {
$atom:a $mol @atom:A 0.000 0.00000 0.00000 0.00000
}
# (The x y z positions will be changed later with move commands
# You can spedify charge and other properties by changing the atom_style.)
# atom-type mass
write_once("Data Masses") {
@atom:A 1.0
}
# pairwise interactions (between non-bonded atoms):
#
# U(r) = 4*eps*((r/sig)^12 - (r/sig)^6)
#
# Note: when sigma=0.8908987181403393=2^(1/6), the minimia is at r=1.0
#
# atom-type atom-type pair_style epsilon sigma
write_once("In Settings") {
# I usually use sigma = 2^(-(1/6)), with a cutoff of 1
#pair_coeff @atom:A @atom:A lj/cut 1.0 0.8908987181403393 1.0
# In the 2013 Science (metaphase) paper, Imakaev used sigma=1.0
# with a cutoff of 2^(1/6). Here we are trying to reproduce his results.
# 10nm fiber
#pair_coeff @atom:A @atom:A lj/cut 1.0 1.0 1.122462048309373
# 30nm fiber (4.25^(1/2)=2.0615528128088303)
#pair_coeff @atom:A @atom:A lj/cut 1.0 2.0615528128088303 2.314014792963349
# 30nm fiber (4.25^(1/3)=1.6198059006387417)
pair_coeff @atom:A @atom:A lj/cut 1.0 1.6198059006387417 1.8181706490945708
}
} # CondensinMonomer

View File

@ -0,0 +1,301 @@
#!/usr/bin/env python
err_msg = """
Usage:
generate_system_lt.py n < monomer_coords.raw > system.lt
Example:
generate_system_lt.py 30118 47 < coords.raw > system.lt
Explanation:
ARGUMENTS:
n = total length of the polymer (in monomers)
L = the (average) length of each condensin interval (Poisson-
distributed) This is also 1/probability that each monomer
is a "condensin monomer".
(Note: 30117 ~= 128000/4.25, but using 30118 makes interpolation cleaner,
and 47 = 200/4.25. Note that 128000 and 200 are for the 10nm model.
See the supplemental section of Naumova et al Science 2013, p 18.)
"""
import sys
import random
from math import *
# Parse the argument list:
if len(sys.argv) <= 2:
sys.stderr.write("Error:\n\nTypical Usage:\n\n"+err_msg+"\n")
exit(1)
N=int(sys.argv[1])
L=float(sys.argv[2])
if len(sys.argv) > 3:
delta_x = float(sys.argv[3])
else:
delta_x = 2.0
if len(sys.argv) > 4:
x_offset = float(sys.argv[4])
else:
x_offset = -((N-1.0)/2) * delta_x
coords = [[0.0, 0.0, 0.0] for i in range(0,N)]
lines = sys.stdin.readlines()
if len(lines) != N:
sys.stderr.write("Error: Number of lines in input file ("+str(len(lines))+")\n"
" does not match first argument ("+str(N)+")\n")
exit(1)
for i in range(0, N):
coords[i] = map(float, lines[i].split())
# Now calculate the box_boundaries:
box_bounds_min = [0.0, 0.0, 0.0]
box_bounds_max = [0.0, 0.0, 0.0]
for i in range(0, N):
for d in range(0, 3):
if i == 0:
box_bounds_min[d] = coords[i][d]
box_bounds_max[d] = coords[i][d]
else:
if coords[i][d] > box_bounds_max[d]:
box_bounds_max[d] = coords[i][d]
if coords[i][d] < box_bounds_min[d]:
box_bounds_min[d] = coords[i][d]
# Now scale the box boundaries outward by 50%
box_scale = 1.5
for d in range(0,3):
box_bounds_cen = 0.5*(box_bounds_max[d] + box_bounds_min[d])
box_bounds_width = box_bounds_max[d] - box_bounds_min[d]
box_bounds_min[d] = box_bounds_cen - 0.5*box_scale*box_bounds_width
box_bounds_max[d] = box_bounds_cen + 0.5*box_scale*box_bounds_width
# Now calculate the direction each molecule should be pointing at:
direction_vects = [[0.0, 0.0, 0.0] for i in range(0,N)]
for d in range(0, 3):
direction_vects[0][d] = coords[1][d] - coords[0][d]
direction_vects[N-1][d] = coords[N-1][d] - coords[N-2][d]
for i in range(1, N-1):
for d in range(0, 3):
direction_vects[i][d] = coords[i+1][d] - coords[i-1][d]
# Optional: normalize the direction vectors
for i in range(1, N-1):
direction_len = 0.0
for d in range(0, 3):
direction_len += (direction_vects[i][d])**2
direction_len = sqrt(direction_len)
for d in range(0, 3):
direction_vects[i][d] /= direction_len
# Now, begin writing the text for the system.lt file:
sys.stdout.write(
"""
import "monomer.lt" # <-- defines "Monomer"
import "condensin.lt" # <-- defines "CondensinMonomer"
"""
)
# Figure out which monomers are "Monomers" and which monomers are
# "CondensinMonomers"
ic = 0 # count the number of condensins added so far
condensin_is_here = [False for i in range(0, N)]
for i in range(0, N):
#add_condensin_here = random.random() < (1.0 / L)
add_condensin_here = random.random() < (1.0 / (L-2.0))
# We do not allow condensin at successive sites separated by less than 2
# subunits (the "L-2.0" above is to compensate for this)
if (((i > 0) and condensin_is_here[i-1]) or
((i > 1) and condensin_is_here[i-2])):
add_condensin_here = False
if add_condensin_here:
condensin_is_here[i] = True
ic += 1
Nc = ic
ic = 0
for i in range(0, N):
if condensin_is_here[i]:
sys.stdout.write("condensins["+str(ic)+"] = new CondensinMonomer.scale(0.5,0.8,0.8).rotvv(1,0,0,")
ic+=1
else:
sys.stdout.write("monomers["+str(i)+"] = new Monomer.scale(0.5,0.8,0.8).rotvv(1,0,0,")
sys.stdout.write(str(direction_vects[i][0])+","
+str(direction_vects[i][1])+","
+str(direction_vects[i][2])+
").move("
+str(coords[i][0])+","
+str(coords[i][1])+","
+str(coords[i][2])+")\n")
#if condensin_is_here[i]:
# if i < N-1:
# sys.stdout.write("\n"
# "#(override the dihedral angle for this monomer)\n"
# "write(\"Data Dihedrals\") {\n"
# " $dihedral:twistor"+str(i+1)+" @dihedral:CondensinMonomer/TWISTOR $atom:monomers["+str(i)+"]/t $atom:monomers["+str(i)+"]/c $atom:monomers["+str(i+1)+"]/c $atom:monomers["+str(i+1)+"]/t\n"
# "}\n"
# "\n")
sys.stdout.write(
"""
# ---------------- simulation box -----------------
# Now define a box big enough to hold a polymer with this (initial) shape
"""
)
sys.stdout.write("write_once(\"Data Boundary\") {\n"
+str(box_bounds_min[0])+" "+str(box_bounds_max[0])+" xlo xhi\n"
+str(box_bounds_min[1])+" "+str(box_bounds_max[1])+" ylo yhi\n"
+str(box_bounds_min[2])+" "+str(box_bounds_max[2])+" zlo zhi\n"
"}\n\n\n")
sys.stdout.write(
"""
# What kind of boundary conditions are we using?
write_once("In Init") {
boundary s s s # <-- boundary conditions in x y z directions
#boundary p p p # <-- boundary conditions in x y z directions
}
# "p" stands for "periodic"
# "s" stands for "shrink-wrapped" (non-periodic)
# ---- Bonds ----
write_once("In Settings") {
# 10nm model:
#bond_coeff @bond:backbone harmonic 100.0 1.0
# 30nm fiber (4.25^(1/3)=1.6198059006387417)
bond_coeff @bond:backbone harmonic 100.0 1.6198059006387417
}
"""
)
sys.stdout.write("write(\"Data Bonds\") {\n")
# Old bond-loop was simple:
#for i in range(0, N-1):
# sys.stdout.write(" $bond:b"+str(i+1)+" @bond:backbone $atom:monomers["+str(i)+"]/a $atom:monomers["+str(i+1)+"]/a\n")
ic = 0
for i in range(0, N-1):
#sys.stderr.write("i="+str(i)+", ic="+str(ic)+", Nc="+str(Nc)+"\n")
# Figure out if the first atom in the bond pair
# belongs to a regular Monomer or a CondensinMonomer
if condensin_is_here[i]:
sys.stdout.write(" $bond:b"+str(i+1)+" @bond:backbone $atom:condensins["+str(ic)+"]/a")
ic+=1
else:
sys.stdout.write(" $bond:b"+str(i+1)+" @bond:backbone $atom:monomers["+str(i)+"]/a")
# Do the same thing for the second atom in the bond pair
if condensin_is_here[i+1]:
assert(ic<Nc)
sys.stdout.write(" $atom:condensins["+str(ic)+"]/a\n")
else:
sys.stdout.write(" $atom:monomers["+str(i+1)+"]/a\n")
sys.stdout.write("}\n\n\n")
sys.stdout.write("""
write_once("Data Angles By Type") {
@angle:backbone @atom:* @atom:* @atom:* @bond:backbone @bond:backbone
}
write_once("In Settings") {
# Most parameters here were taken from the supplemental material of
# Naumova et al. Science 2013 (simulations by Maxim Imakaev, see Supp Mat)
#angle_coeff @angle:backbone cosine 5.0 #<-10nm fiber
angle_coeff @angle:backbone cosine 1.1764705882352942 #<-30nm fiber
}
""")
sys.stdout.write(
"""
# ---- Condensins randomly located on the polymer ----
# Stage 1:
# Add bonds between consecutive condensin anchors.
# Imakaev calls this "stage 1: linear compaction":
"""
)
sys.stdout.write("write(\"Data Bonds\") {\n")
ic = 0
for i in range(0, N):
if condensin_is_here[i]:
if (0 < ic):
sys.stdout.write(" $bond:bstage1_"+str(ic-1)+" @bond:stage1 $atom:condensins["+str(ic-1)+"]/a $atom:condensins["+str(ic)+"]/a\n")
ic += 1
sys.stdout.write("}\n\n")
sys.stdout.write("""
# Stage 2:
# Add additional bonds between all pairs of condensin anchors
# in a window of |ic-jc| <= 30 anchors (along the chain).
# In the paper, they call this stage 2 axial compression.
write("Data Bonds") {
""")
jcwindowsize = 30
for ic in range(0, Nc):
#jcmin = max(ic-jcwindowsize, 0)
#jcmax = min(ic+jcwindowsize, Nc-1)
jcmax = min(ic+jcwindowsize, Nc-1)
for jc in range(ic+2, jcmax+1):
sys.stdout.write(" $bond:bstage2_"+str(ic)+"_"+str(jc)+" @bond:stage2 $atom:condensins["+str(ic)+"]/a $atom:condensins["+str(jc)+"]/a\n")
sys.stdout.write("}\n")
sys.stdout.write("""
write_once("In Settings") {
# stage 1 bonds are initially off
bond_coeff @bond:stage1 harmonic 0.0 0.5 # <--(we can override this later)"
# stage 2 bonds are initially off
bond_coeff @bond:stage2 harmonic 0.0 0.0 # <--(we can override this later)"
}
""")
sys.stdout.write("\n\n# "+str(Nc)+" condensin molecules added\n\n")
sys.stderr.write("\n\n# "+str(Nc)+" condensin molecules added\n\n")

View File

@ -0,0 +1,74 @@
#!/usr/bin/env python
err_msg = """
Usage:
interpolate_coords.py Ndesired [scale] < coords_orig.raw > coords.raw
Example:
interpolate_coords.py 30118 3.0 < coords_orig.raw > coords.raw
# (Note: 30117 ~= 128000/4.25, but using 30118 makes interpolation cleaner.
# See the supplemental section of Naumova et al Science 2013, p 18.)
"""
import sys
from math import *
# Parse the argument list:
if len(sys.argv) <= 1:
sys.stderr.write("Error:\n\nTypical Usage:\n\n"+err_msg+"\n")
exit(1)
n_new = int(sys.argv[1])
if len(sys.argv) > 2:
scale = float(sys.argv[2])
else:
scale = 1.0
coords_orig = []
lines = sys.stdin.readlines()
for line in lines:
tokens = line.split()
if (len(tokens) > 0):
coords_orig.append(map(float, tokens))
g_dim = len(tokens)
n_orig = len(coords_orig)
if n_orig < 2:
sys.stderr.write("Error:\n\nInput file contains less than two lines of coordinates\n")
exit(1)
if n_new < 2:
sys.stderr.write("Error:\n\nOutput file will contain less than two lines of coordinates\n")
exit(1)
coords_new = [[0.0 for d in range(0, g_dim)] for i in range(0, n_new)]
for i_new in range(0, n_new):
I_orig = (i_new) * (float(n_orig-1) / float(n_new-1))
i_orig = int(floor(I_orig))
i_remainder = I_orig - i_orig
if (i_new < n_new-1):
for d in range(0, g_dim):
coords_new[i_new][d] = scale*(coords_orig[i_orig][d]
+
i_remainder*(coords_orig[i_orig+1][d]-
coords_orig[i_orig][d]))
else:
for d in range(0, g_dim):
coords_new[i_new][d] = scale*coords_orig[n_orig-1][d]
# print the coordates
for d in range(0, g_dim-1):
sys.stdout.write(str(coords_new[i_new][d]) + ' ')
sys.stdout.write(str(coords_new[i_new][g_dim-1]) + "\n")

View File

@ -0,0 +1,84 @@
# This file contains the definition of a molecule named "Monomer".
# (This particular molecule contain only one atom, but that is up to you.)
# Later, multiple monomers can be connected together to build a molecule.
Monomer {
# atom-id mol-id(ignore) atom-type q x y z
write("Data Atoms") {
$atom:a $mol @atom:A 0.000 0.00000 0.00000 0.00000
}
# (The x y z positions will be changed later with move commands
# You can spedify charge and other properties by changing the atom_style.)
# atom-type mass
write_once("Data Masses") {
@atom:A 1.0
}
# pairwise interactions (between non-bonded atoms):
#
# U(r) = 4*eps*((r/sig)^12 - (r/sig)^6)
#
# Note: when sigma=0.8908987181403393=2^(1/6), the minimia is at r=1.0
#
# atom-type atom-type pair_style epsilon sigma rcutoff
write_once("In Settings") {
# I usually use sigma = 2^(-(1/6)), with a cutoff of 1
#pair_coeff @atom:A @atom:A lj/cut 1.0 0.8908987181403393 1.0
# In the 2013 Science (metaphase) paper, Imakaev used sigma=1.0
# with a cutoff of 2^(1/6). Here we are trying to reproduce his results.
# 10nm fiber
#pair_coeff @atom:A @atom:A lj/cut 1.0 1.0 1.122462048309373
# 30nm fiber (4.25^(1/2)=2.0615528128088303)
#pair_coeff @atom:A @atom:A lj/cut 1.0 2.0615528128088303 2.314014792963349
# 30nm fiber (4.25^(1/3)=1.6198059006387417)
pair_coeff @atom:A @atom:A lj/cut 1.0 1.6198059006387417 1.8181706490945708
}
} # Monomer
# --------------------------------------------------------------------
#
# At some point we need to specify which force-field styles we want.
# LAMMPS also allows you to customize the kinds of properties you want
# each atom to have (the "atom_style"), such as charge, molecule-id, dipole etc.
# I typically specify this here. Doing it this way means that all systems built
# from "Monomers" (ie which import "monomer.lt") share these atom-styles
# and force-field styles by default. You can override these settings later.
write_once("In Init") {
# Default styles for molecules built out of "Monomers"
units lj
atom_style full
bond_style hybrid harmonic table linear 4001
angle_style hybrid cosine
dihedral_style none
# If you need angles, dihedrals and impropers, uncomment or replace:
# angle_style hybrid harmonic
# dihedral_style hybrid fourier
pair_style hybrid lj/cut 2.5
# If you are using gpu acceleration uncomment these lines:
# package gpu force/neigh 0 0 1.0
# pair_style hybrid lj/cut/gpu 4.0
pair_modify mix arithmetic
special_bonds lj/coul 1 1 1
}

View File

@ -0,0 +1,41 @@
#########################################################
# Example how to run this file:
#
# 1) Choose a ransom seed (in this case 141203)
# (or use `bash -c 'echo $RANDOM'`)
#
# 2) Then, from the shell, invoke LAMMPS to collapse the polymer:
#
# lmp_ubuntu_parallel -i run.in -var seed 141203
#
#########################################################
# eg:
# time mpirun -np 2 lmp_ubuntu_parallel -i run.in.min
#########################################################
# -- 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 10000 traj_min.lammpstrj id mol type x y z ix iy iz
thermo_style custom step pe etotal vol epair ebond eangle edihed
thermo_modify norm no #(report total energy not energy / num_atoms)
thermo 100 #(time interval for printing out "thermo" data)
# Now minimize the system:
min_style quickmin
min_modify dmax 0.05
minimize 1.0e-7 1.0e-8 30000 100000000
write_data system_after_min.data

View File

@ -0,0 +1,110 @@
# PREREQUISITES: You must run LAMMPS using "run.in.min" beforehand.
# (This will create the "system_after_min0.data" file needed below.)
#########################################################
# Run using:
#
# lmp_ubuntu_parallel -i run.in.stage1
#
#########################################################
# GPUs:
# To enable gpu acceleration, make sure settings.in.init includes this line:
# package gpu force/neigh 0 0 1.0 (make sure it is not commented out.)
# ...and replace "lj/cut" in the "settings.in.init" and "settings.in.settings"
# files with "lj/cut/gpu"
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
#read_data system.data
read_data system_after_min.data
# -- Settings Section --
include system.in.settings
# -- Run Section --
# DON'T MINIMIZE FIRST UNLESS YOU CHOOSE THE CORRECT INITIAL KbondC FORMULA
#thermo_style custom step pe etotal vol epair ebond eangle edihed
#thermo_modify norm no #(report total energy not energy / num_atoms)
#thermo 20 #(time interval for printing out "thermo" data)
#min_style sd
#min_modify dmax 0.05
#minimize 1.0e-7 1.0e-8 20000 1000000
#write_data system_after_min_t=0.data
mass * 1.0
timestep 0.005 # "dt"
dump 1 all custom 25000 traj_stage1.lammpstrj id mol type x y z ix iy iz
reset_timestep 0
# --- run the simulation ---
# set the velocity to zero
velocity all create 0.0 123456
# 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.)
# Tstart Tstop tdamp randomseed
fix fxlan all langevin 1.0 1.0 10.0 123456
# pstart pstop pdamp(time-units, 2000 iters usually)
#fix fxnph all nph x -0.000 -0.000 1.0
fix fxnve all nve
# (See http://lammps.sandia.gov/doc/fix_langevin.html)
fix fxcenter all recenter 0.0 0.0 0.0
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
thermo_modify norm no #(report total energy not energy / num_atoms)
thermo 1000 #(time interval for printing out "thermo" data)
#balance dynamic x 20 1.0 -out tmp.balance
#balance x uniform
variable nloop1 loop 300
label loop1
print "############### LOOP ${nloop1} ###############"
# Now, change the bond-strength between condensin monomers.
# From the Naumova et al Science 2013 paper (supp materials)
# "Two-stage model: linear compaction + axial compression"
# "First, random consecutive loops with L=100 monomers (see above) were
# introduced, and anchors of neighboring loops were brought together
# using harmonic springs with a potential U = k * (r r0)2; r0=0.5.
# To avoid abrupt motion of the loop anchors, the force was gradually
# turned on over the first 300000 timesteps, with k linearly increasing
# in time from 0 to 10 kT."
# Do this by changing the parameters in the force-field for these
# bonds.
#
# Formula used for "bond_style harmonic":
# Ubond(r) = k*(r-r0)^2
# bondType style
#bond_coeff 1 harmonic 0.1 0.5
variable time equal step
variable KbondC equal $((v_time+1)*(10.0/300000.0))
print "timestep = ${time}, KbondC = ${KbondC}" file KbondC_vs_time.dat
#bond_coeff 1 harmonic ${KbondC} 0.5
bond_coeff 1 harmonic ${KbondC} 0.5
run 1000
next nloop1
jump SELF loop1
write_data system_after_stage1.data

View File

@ -0,0 +1,90 @@
# PREREQUISITES: You must run LAMMPS using "run.in.stage1" beforehand.
# (This will create the "system_after_stage1.data" file.)
#########################################################
# Run using:
#
# lmp_ubuntu_parallel -i run.in.stage2
#
#########################################################
# eg:
# time mpirun -np 2 lmp_ubuntu_parallel -i run.in -var seed 1
#########################################################
# GPUs:
# To enable gpu acceleration, make sure settings.in.init includes this line:
# package gpu force/neigh 0 0 1.0 (make sure it is not commented out.)
# ...and replace "lj/cut" in the "settings.in.init" and "settings.in.settings"
# files with "lj/cut/gpu"
# -- Init Section --
include system.in.init
# -- Atom Definition Section --
#read_data system.data
read_data system_after_stage1.data
# -- Settings Section --
include system.in.settings
# -- Run Section --
mass * 1.0
timestep 0.005 # "dt"
dump 1 all custom 50000 traj_stage2.lammpstrj id mol type x y z ix iy iz
reset_timestep 300000
# --- run the simulation ---
# set the velocity to zero
velocity all create 0.0 123456
# 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.)
# Tstart Tstop tdamp randomseed
fix fxlan all langevin 1.0 1.0 10.0 123456
# pstart pstop pdamp(time-units, 2000 iters usually)
fix fxnve all nve
# (See http://lammps.sandia.gov/doc/fix_langevin.html)
fix fxcenter all recenter 0.0 0.0 0.0
thermo_style custom step temp pe etotal press vol epair ebond eangle edihed
thermo_modify norm no #(report total energy not energy / num_atoms)
thermo 1000 #(time interval for printing out "thermo" data)
#balance dynamic x 20 1.0 -out tmp.balance
#balance x uniform
# atomTypes pairStyle epsilon sigma rcutoff
# 10nm-fiber
#pair_coeff 1 1 lj/cut 1.0 1.0 2.5
#pair_coeff 2 2 lj/cut 1.0 1.0 2.5
# 30nm fiber (4.25^(1/3)=1.6198059006387417)
pair_coeff 1 1 lj/cut 1.0 1.6198059006387417 4.049514751596854
pair_coeff 2 2 lj/cut 1.0 1.6198059006387417 4.049514751596854
# During stage 2, add attractive forces between all pairs of non-consecutive
# condensin anchors. These forces are stored in the table file below:
# bondType bondStyle filename label
bond_coeff 2 table table_bonds_stage2.dat STAGE2
# During stage 2, I assume the stage-1 bonds remain in place
# (They have length 0.5.
# After 300000 timesteps during stage 1, the "k" value should be 10.0.)
# bondType bondStyle k r0
bond_coeff 1 harmonic 10.0 0.5
timestep 0.005
run 1700000
write_data system_after_stage2.data

View File

@ -0,0 +1,33 @@
# 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

View File

@ -0,0 +1,32 @@
# 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

View File

@ -0,0 +1,31 @@
# 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

View File

@ -0,0 +1,24 @@
# 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 ../

View File

@ -0,0 +1,87 @@
------- To view a lammps trajectory in VMD --------
1) Build a PSF file for use in viewing with VMD.
This step works with VMD 1.9 and topotools 1.2.
(Older versions, like VMD 1.8.6, don't support this.)
a) Start VMD
b) Menu Extensions->Tk Console
c) Enter:
(I assume that the the DATA file is called "system.data")
topo readlammpsdata system.data full
animate write psf system.psf
2)
Later, to Load a trajectory in VMD:
Start VMD
Select menu: File->New Molecule
-Browse to select the PSF file you created above, and load it.
(Don't close the window yet.)
-Browse to select the trajectory file.
If necessary, for "file type" select: "LAMMPS Trajectory"
Load it.
---- A note on trajectory format: -----
If the trajectory is a DUMP file, then make sure the it contains the
information you need for pbctools (see below. I've been using this
command in my LAMMPS scripts to create the trajectories:
dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
It's a good idea to use an atom_style which supports molecule-ID numbers
so that you can assign a molecule-ID number to each atom. (I think this
is needed to wrap atom coordinates without breaking molecules in half.)
Of course, you don't have to save your trajectories in DUMP format,
(other formats like DCD work fine) I just mention dump files
because these are the files I'm familiar with.
3) ----- Wrap the coordinates to the unit cell
(without cutting the molecules in half)
a) Start VMD
b) Load the trajectory in VMD (see above)
c) Menu Extensions->Tk Console
d) Try entering these commands:
pbc wrap -compound res -all
pbc box
----- Optional ----
Sometimes the solvent or membrane obscures the view of the solute.
It can help to shift the location of the periodic boundary box
To shift the box in the y direction (for example) do this:
pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
pbc box -shiftcenterrel {0.0 0.15 0.0}
Distances are measured in units of box-length fractions, not Angstroms.
Alternately if you have a solute whose atoms are all of type 1,
then you can also try this to center the box around it:
pbc wrap -sel type=1 -all -centersel type=2 -center com
4)
You should check if your periodic boundary conditions are too small.
To do that:
select Graphics->Representations menu option
click on the "Periodic" tab, and
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
5) Optional: If you like, change the atom types in the PSF file so
that VMD recognizes the atom types, use something like:
sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf
(If you do this, it might effect step 2 above.)

View File

@ -0,0 +1,86 @@
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

View File

@ -0,0 +1,216 @@
# 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

View File

@ -0,0 +1,85 @@
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

View File

@ -0,0 +1,41 @@
# Here we define a trivial molecule containing only one particle.
Chaperonin {
# atomID molID atomType charge x y z
write("Data Atoms") {
$atom:C $mol @atom:C 0.0 0.0 0.0 0.0
}
write_once("Data Masses") {
@atom:C 100.0
}
write_once("In Settings") {
# If for some reason there are multiple chaperones present,
# I assume that they interact repulsively (hence, L=0)
# i j epsilon sigma K L
pair_coeff @atom:C @atom:C lj/charmm/coul/charmm/inter 1.0 6.0 1 0
# Optional: define the atoms in the "chaperonins" group:
# (Defining a group for the chaperone makes it easy to immobilize it later.)
group chaperonins type @atom:C
}
# Specify which pair_styles, and atom styles work well with
# this model. (Again this can be overridden later.)
write_once("In Init") {
units lj
atom_style full
pair_style hybrid lj/charmm/coul/charmm/inter es4k4l maxmax 11.0 12.0
}
} # Chaperonin
# We have not specified how this particle interacts with other particles
# besides itself. Later on you must do this.

View File

@ -0,0 +1,87 @@
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between atoms in the
# protein and a chaperone provided in the supplemental materials section of:
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
# This is stored in a tabulated force field with a singularity at a distance R.
#
# To calculate the table for interaction between
# ...the chaperone and a hydrophobic bead (2004 PNAS paper), use this table:
# ./calc_chaperone_table.py 1.0 1.0 6.0 0.475 0.0 5.9 1181
# ...the chaperone and a hydrophilic bead (2004 PNAS paper), use this table:
# ./calc_chaperone_table.py 1.0 1.0 6.0 0.0 0.0 5.9 1181
# ...the chaperone and a hydrophobic bead (2006 JMB paper), use this table:
# ./calc_chaperone_table.py 1.0 1.0 3.0 0.60 3.1 8.0 981 True
# ...the chaperone and a hydrophilic bead (2006 JMB paper), use this table:
# ./calc_chaperone_table.py 1.0 1.0 3.0 0.0 3.1 8.0 981 True
from math import *
import sys
def U(r, eps, sigma, R, h):
#print('r='+str(r)+' eps='+str(eps)+' s='+str(sigma)+' R='+str(R)+' h='+str(h))
# Formula is undefined at r=0, but you can take the limit:
if r <= 0:
return 4.0*pi*R*R*4.0*eps*(pow((sigma/R), 12.0)
- h*pow((sigma/R), 6.0))
xp = sigma/(r+R)
xm = sigma/(r-R)
term10 = pow(xm, 10.0) - pow(xp, 10.0)
term4 = pow(xm, 4.0) - pow(xp, 4.0)
return 4.0*pi*eps*(R/r) * (0.2*term10 - 0.5*h*term4)
def F(r, eps, sigma, R, h):
# Formula is undefined at r=0, but you can take the limit:
if r <= 0:
return 0.0
product_term_a = U(r, eps, sigma, R, h) / r
ixp = (r+R)/sigma
ixm = (r-R)/sigma
dix_dr = 1.0/sigma
term10 = (10.0/sigma)*(pow(ixm, -11.0) - pow(ixp, -11.0))
term4 = (4.0/sigma)*(pow(ixm, -5.0) - pow(ixp, -5.0))
product_term_b = 4.0*eps*pi*(R/r) * (0.2*term10 - 0.5*h*term4)
return product_term_a + product_term_b
class InputError(Exception):
""" A generic exception object containing a string for error reporting.
"""
def __init__(self, err_msg):
self.err_msg = err_msg
def __str__(self):
return self.err_msg
def __repr__(self):
return str(self)
if len(sys.argv) < 8:
sys.stderr.write("Error: expected 7 arguments:\n"
"\n"
"Usage: "+sys.argv[0]+" epsilon sigma R h rmin rmax N\n\n")
sys.exit(-1)
epsilon = float(sys.argv[1])
sigma = float(sys.argv[2])
R = float(sys.argv[3])
h = float(sys.argv[4])
rmin = float(sys.argv[5])
rmax = float(sys.argv[6])
N = int(sys.argv[7])
subtract_Urcut = False
if len(sys.argv) == 9:
subtract_Urcut = True
rcut = rmax
for i in range(0,N):
r = rmin + i*(rmax-rmin)/(N-1)
U_r = U(r, epsilon, sigma, R, h)
F_r = F(r, epsilon, sigma, R, h)
if subtract_Urcut:
U_r -= U(rcut, epsilon, sigma, R, h)
if (r >= rcut) or (i==N-1):
U_r = 0.0
F_r = 0.0
print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))

View File

@ -0,0 +1,67 @@
#!/usr/bin/env python
# Calculate a table of dihedral angle interactions used in the alpha-helix
# and beta-sheet regions of the frustrated protein model described in
# provided in figure 8 of the supplemental materials section of:
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
# Note that the "A" and "B" parameters were incorrectly reported to be
# 5.4*epsilon and 6.0*epsilon. The values used were 5.6 and 6.0 epsilon.
# The phiA and phiB values were 57.29577951308232 degrees (1 rad)
# and 180 degrees, respectively. Both expA and expB were 6.0.
#
# To generate the table used for the alpha-helix (1 degree resolution) use this:
# ./calc_dihedral_table.py 6.0 57.29577951308232 6 5.6 180 6 0.0 359 360
# To generate the table used for the beta-sheets (1 degree resolution) use this:
# ./calc_dihedral_table.py 5.6 57.29577951308232 6 6.0 180 6 0.0 359 360
#
# (If you're curious as to why I set the location of the minima at phi_alpha
# to 1.0 radians (57.2957795 degrees), there was no particularly good reason.
# I think the correct value turns out to be something closer to 50 degrees.)
from math import *
import sys
# The previous version included the repulsive core term
def U(phi, A, phiA, expA, B, phiB, expB, use_radians=False):
conv_units = pi/180.0
if use_radians:
conv_units = 1.0
termA = pow(cos(0.5*(phi-phiA)*conv_units), expA)
termB = pow(cos(0.5*(phi-phiB)*conv_units), expB)
return -A*termA - B*termB
# The previous version included the repulsive core term
def F(phi, A, phiA, expA, B, phiB, expB, use_radians=False):
conv_units = pi/180.0
if use_radians:
conv_units = 1.0
termA = (0.5*sin(0.5*(phi-phiA)*conv_units) *
expA * pow(cos(0.5*(phi-phiA)*conv_units), expA-1.0))
termB = (0.5*sin(0.5*(phi-phiB)*conv_units) *
expB * pow(cos(0.5*(phi-phiB)*conv_units), expB-1.0))
return -conv_units*(A*termA + B*termB)
if len(sys.argv) != 10:
sys.stderr.write("Error: expected 9 arguments:\n"
"\n"
"Usage: "+sys.argv[0]+" A phiA expA B phiB expB phiMin phiMax N\n\n")
sys.exit(-1)
A = float(sys.argv[1])
phiA = float(sys.argv[2])
expA = float(sys.argv[3])
B = float(sys.argv[4])
phiB = float(sys.argv[5])
expB = float(sys.argv[6])
phi_min = float(sys.argv[7])
phi_max = float(sys.argv[8])
N = int(sys.argv[9])
for i in range(0,N):
phi = phi_min + i*(phi_max - phi_min)/(N-1)
U_phi = U(phi, A, phiA, expA, B, phiB, expB, use_radians=False)
F_phi = F(phi, A, phiA, expA, B, phiB, expB, use_radians=False)
print(str(i+1)+' '+str(phi)+' '+str(U_phi)+' '+str(F_phi))

View File

@ -0,0 +1,45 @@
write_once("Data Boundary") {
0.0 20.0 xlo xhi
0.0 20.0 ylo yhi
0.0 20.0 zlo zhi
}
import "1beadFrustrated_variants.lt"
import "chaperonin.lt"
protein = new 1beadMisfolded # (frustrated protein, misfolded conformation)
chaperinin = new Chaperonin # (hollow chaperonin cavity. usually immobile)
# ---- Now define interactions between the atoms in the protein ----
# ---- (named "B", "L", "N") and the atom which represents the ----
# ---- chaperone ("C"). These interactions are tabulated. ----
write_once("In Settings") {
pair_coeff @atom:Chaperonin/C @atom:1beadFrustrated/B table table_chaperonin_h=0.475.dat CH_H0.475
pair_coeff @atom:Chaperonin/C @atom:1beadFrustrated/L table table_chaperonin_h=0.dat CH_H0
pair_coeff @atom:Chaperonin/C @atom:1beadFrustrated/N table table_chaperonin_h=0.dat CH_H0
}
# Note: If you want to use a "hydrophilic" chaperone (with h=0, not h=0.475)
# then replace "table_chaperonin_h=0_475.dat CH_H0.475"
# with "table_chaperonin_h=0.dat CH_H0"
# LAMMPS has many available force field styles (and atom styles). Here we
# select the ones which work well for the full combine system. (This should
# override any settings made in "1beadFrustrated.lt" or "chaperonin.lt")
write_once("In Init") {
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 table spline 1181
pair_modify mix arithmetic
special_bonds lj 0.0 0.0 1.0 #(turn on "1-4" interactions)
}

View File

@ -0,0 +1,735 @@
# Table of the potential and its negative derivative for frustrated alpha helix
# (Note: Derivatives are in units of energy/radians, not energy/degrees.)
# ./calc_dihedral_table.py 6.0 57.29577951308232 6 5.6 180 6 0.0 359 360
FRUSTRATED_ALPHA
N 360 DEGREES
1 0.0 -2.74081145103 0.0783990792662
2 1.0 -2.81950869101 0.0789852583442
3 2.0 -2.89876136749 0.0795096391909
4 3.0 -2.97850675562 0.0799703813963
5 4.0 -3.05868032959 0.0803657243943
6 5.0 -3.13921584545 0.0806939935737
7 6.0 -3.22004543014 0.0809536062381
8 7.0 -3.30109967628 0.0811430773977
9 8.0 -3.38230774267 0.0812610253741
10 9.0 -3.46359746038 0.0813061772009
11 10.0 -3.54489544401 0.0812773738039
12 11.0 -3.62612720812 0.0811735749433
13 12.0 -3.70721728841 0.0809938639029
14 13.0 -3.78808936748 0.080737451911
15 14.0 -3.86866640485 0.0804036822781
16 15.0 -3.94887077101 0.0799920342374
17 16.0 -4.02862438516 0.0795021264757
18 17.0 -4.10784885622 0.0789337203415
19 18.0 -4.18646562704 0.0782867227197
20 19.0 -4.26439612115 0.0775611885609
21 20.0 -4.34156189202 0.0767573230567
22 21.0 -4.41788477419 0.0758754834523
23 22.0 -4.49328703609 0.0749161804868
24 23.0 -4.56769153408 0.0738800794563
25 24.0 -4.64102186743 0.0727680008923
26 25.0 -4.71320253365 0.0715809208518
27 26.0 -4.78415908407 0.0703199708131
28 27.0 -4.85381827903 0.0689864371778
29 28.0 -4.92210824234 0.067581760373
30 29.0 -4.98895861476 0.0661075335571
31 30.0 -5.05430070586 0.0645655009259
32 31.0 -5.11806764409 0.0629575556235
33 32.0 -5.18019452449 0.061285737258
34 33.0 -5.24061855376 0.0595522290273
35 34.0 -5.29927919225 0.0577593544584
36 35.0 -5.3561182925 0.0559095737673
37 36.0 -5.41108023395 0.0540054798439
38 37.0 -5.46411205346 0.0520497938726
39 38.0 -5.51516357127 0.0500453605949
40 39.0 -5.56418751203 0.0479951432253
41 40.0 -5.61113962059 0.0459022180302
42 41.0 -5.65597877221 0.0437697685824
43 42.0 -5.69866707689 0.0416010797029
44 43.0 -5.7391699774 0.0393995311046
45 44.0 -5.77745634094 0.0371685907508
46 45.0 -5.81349854393 0.034911807945
47 46.0 -5.84727254977 0.0326328061676
48 47.0 -5.87875797937 0.030335275675
49 48.0 -5.90793817411 0.0280229658805
50 49.0 -5.93480025113 0.0256996775336
51 50.0 -5.95933515063 0.0233692547166
52 51.0 -5.98153767519 0.0210355766777
53 52.0 -6.00140652074 0.0187025495211
54 53.0 -6.01894429926 0.016374097773
55 54.0 -6.03415755288 0.0140541558448
56 55.0 -6.04705675953 0.0117466594146
57 56.0 -6.05765632981 0.00945553674764
58 57.0 -6.06597459526 0.00718469997761
59 58.0 -6.07203378786 0.00493803637051
60 59.0 -6.07586001075 0.00271939959245
61 60.0 -6.07748320034 0.000532601003776
62 61.0 -6.07693707962 -0.00161859899905
63 62.0 -6.07425910291 -0.00373049957158
64 63.0 -6.06949039207 -0.00579946791801
65 64.0 -6.06267566421 -0.00782194767468
66 65.0 -6.05386315117 -0.00979446715893
67 66.0 -6.04310451074 -0.0117136474624
68 67.0 -6.03045472992 -0.0135762103679
69 68.0 -6.01597202036 -0.0153789860691
70 69.0 -5.99971770618 -0.0171189206741
71 70.0 -5.98175610439 -0.0187930834719
72 71.0 -5.9621543982 -0.0203986739443
73 72.0 -5.9409825034 -0.0219330285036
74 73.0 -5.91831292823 -0.0233936269399
75 74.0 -5.89422062685 -0.0247780985587
76 75.0 -5.86878284696 -0.0260842279959
77 76.0 -5.84207897162 -0.0273099606906
78 77.0 -5.81419035593 -0.0284534080045
79 78.0 -5.78520015867 -0.0295128519729
80 79.0 -5.7551931694 -0.0304867496727
81 80.0 -5.72425563141 -0.0313737371989
82 81.0 -5.6924750609 -0.0321726332348
83 82.0 -5.65994006273 -0.0328824422092
84 83.0 -5.62674014332 -0.0335023570292
85 84.0 -5.59296552097 -0.0340317613814
86 85.0 -5.55870693409 -0.0344702315961
87 86.0 -5.52405544786 -0.0348175380654
88 87.0 -5.48910225957 -0.0350736462148
89 88.0 -5.45393850338 -0.0352387170203
90 89.0 -5.41865505462 -0.0353131070729
91 90.0 -5.38334233438 -0.0352973681855
92 91.0 -5.34809011465 -0.0351922465446
93 92.0 -5.31298732458 -0.0349986814067
94 93.0 -5.27812185824 -0.034717803342
95 94.0 -5.24358038438 -0.0343509320285
96 95.0 -5.2094481586 -0.0338995736008
97 96.0 -5.17580883839 -0.0333654175598
98 97.0 -5.14274430152 -0.0327503332496
99 98.0 -5.11033446814 -0.0320563659092
100 99.0 -5.07865712698 -0.0312857323082
101 100.0 -5.04778776623 -0.0304408159764
102 101.0 -5.01779940929 -0.0295241620384
103 102.0 -4.98876245596 -0.0285384716647
104 103.0 -4.96074452928 -0.0274865961525
105 104.0 -4.93381032851 -0.0263715306507
106 105.0 -4.90802148862 -0.0251964075427
107 106.0 -4.88343644644 -0.0239644895038
108 107.0 -4.86011031397 -0.0226791622487
109 108.0 -4.83809475914 -0.0213439269874
110 109.0 -4.81743789414 -0.0199623926068
111 110.0 -4.79818417182 -0.0185382675969
112 111.0 -4.78037429015 -0.0170753517415
113 112.0 -4.76404510526 -0.0155775275918
114 113.0 -4.74922955293 -0.0140487517461
115 114.0 -4.73595657904 -0.0124930459538
116 115.0 -4.7242510789 -0.0109144880672
117 116.0 -4.71413384576 -0.00931720286182
118 117.0 -4.70562152846 -0.00770535274772
119 118.0 -4.69872659855 -0.00608312839491
120 119.0 -4.69345732669 -0.00445473929448
121 120.0 -4.6898177686 -0.00282440427898
122 121.0 -4.68780776044 -0.00119634202478
123 122.0 -4.68742292374 0.000425238440527
124 123.0 -4.68865467977 0.0020361472029
125 124.0 -4.69149027336 0.00363222287571
126 125.0 -4.69591280613 0.00520934194008
127 126.0 -4.70190127895 0.0067634279891
128 127.0 -4.70943064365 0.00829046085365
129 128.0 -4.71847186379 0.00978648558781
130 129.0 -4.72899198423 0.0112476212922
131 130.0 -4.74095420961 0.0126700697544
132 131.0 -4.7543179912 0.0140501238848
133 132.0 -4.76903912216 0.0153841759291
134 133.0 -4.78506984093 0.0166687254364
135 134.0 -4.80235894235 0.0179003869651
136 135.0 -4.82085189642 0.0190758975074
137 136.0 -4.84049097437 0.0201921236154
138 137.0 -4.86121538156 0.0212460682116
139 138.0 -4.88296139722 0.0222348770682
140 139.0 -4.90566252032 0.0231558449399
141 140.0 -4.9292496215 0.0240064213355
142 141.0 -4.95365110055 0.0247842159162
143 142.0 -4.97879304911 0.0254870035063
144 143.0 -5.00459941816 0.0261127287073
145 144.0 -5.03099218995 0.0266595101027
146 145.0 -5.05789155387 0.0271256440463
147 146.0 -5.08521608601 0.0275096080241
148 147.0 -5.11288293171 0.0278100635833
149 148.0 -5.14080799097 0.0280258588231
150 149.0 -5.16890610603 0.0281560304409
151 150.0 -5.19709125082 0.0281998053314
152 151.0 -5.22527672173 0.0281566017347
153 152.0 -5.25337532941 0.0280260299338
154 153.0 -5.28129959092 0.0278078924984
155 154.0 -5.30896192196 0.0275021840788
156 155.0 -5.33627482866 0.0271090907491
157 156.0 -5.36315109852 0.0266289889046
158 157.0 -5.38950398994 0.026062443717
159 158.0 -5.41524742011 0.0254102071518
160 159.0 -5.44029615055 0.0246732155563
161 160.0 -5.46456597019 0.0238525868232
162 161.0 -5.48797387528 0.0229496171403
163 162.0 -5.51043824587 0.0219657773349
164 163.0 -5.53187901853 0.0209027088232
165 164.0 -5.55221785468 0.0197622191769
166 165.0 -5.57137830441 0.0185462773191
167 166.0 -5.58928596528 0.0172570083629
168 167.0 -5.60586863576 0.0158966881068
169 168.0 -5.62105646307 0.0144677372016
170 169.0 -5.63478208493 0.0129727150063
171 170.0 -5.64698076513 0.0114143131467
172 171.0 -5.65759052241 0.00979534879707
173 172.0 -5.66655225257 0.00811875770075
174 173.0 -5.67380984344 0.00638758694863
175 174.0 -5.67931028251 0.00460498753534
176 175.0 -5.68300375706 0.00277420671195
177 176.0 -5.68484374646 0.000898580155594
178 177.0 -5.68478710669 -0.00101847602368
179 178.0 -5.68279414663 -0.00297347341791
180 179.0 -5.67882869631 -0.00496285957718
181 180.0 -5.67285816674 -0.00698302636509
182 181.0 -5.6648536014 -0.00903031839234
183 182.0 -5.65478971926 -0.0111010415069
184 183.0 -5.64264494925 -0.0131914713189
185 184.0 -5.62840145627 -0.0152978617389
186 185.0 -5.6120451586 -0.017416453508
187 186.0 -5.59356573683 -0.0195434826976
188 187.0 -5.57295663425 -0.0216751891584
189 188.0 -5.55021504898 -0.0238078248974
190 189.0 -5.52534191754 -0.0259376623617
191 190.0 -5.4983418904 -0.0280610026087
192 191.0 -5.46922329932 -0.0301741833429
193 192.0 -5.43799811672 -0.0322735868002
194 193.0 -5.40468190731 -0.0343556474589
195 194.0 -5.36929377207 -0.0364168595607
196 195.0 -5.33185628476 -0.0384537844225
197 196.0 -5.29239542138 -0.0404630575223
198 197.0 -5.25094048245 -0.0424413953416
199 198.0 -5.20752400881 -0.0443856019501
200 199.0 -5.16218169074 -0.0462925753151
201 200.0 -5.11495227114 -0.0481593133234
202 201.0 -5.06587744261 -0.0499829195012
203 202.0 -5.01500173918 -0.0517606084187
204 203.0 -4.96237242264 -0.0534897107689
205 204.0 -4.90803936404 -0.055167678109
206 205.0 -4.85205492059 -0.0567920872546
207 206.0 -4.79447380837 -0.0583606443179
208 207.0 -4.73535297113 -0.0598711883816
209 208.0 -4.6747514457 -0.0613216948024
210 209.0 -4.61273022413 -0.0627102781377
211 210.0 -4.54935211328 -0.0640351946902
212 211.0 -4.4846815919 -0.0652948446678
213 212.0 -4.41878466581 -0.0664877739558
214 213.0 -4.35172872155 -0.0676126754981
215 214.0 -4.28358237872 -0.0686683902899
216 215.0 -4.21441534165 -0.0696539079796
217 216.0 -4.14429825061 -0.070568367083
218 217.0 -4.07330253293 -0.0714110548116
219 218.0 -4.00150025463 -0.0721814065199
220 219.0 -3.92896397266 -0.072879004774
221 220.0 -3.85576658834 -0.0735035780505
222 221.0 -3.78198120223 -0.0740549990687
223 222.0 -3.70768097086 -0.0745332827669
224 223.0 -3.63293896573 -0.0749385839297
225 224.0 -3.5578280347 -0.0752711944755
226 225.0 -3.48242066643 -0.075531540416
227 226.0 -3.4067888579 -0.0757201784978
228 227.0 -3.33100398548 -0.0758377925383
229 228.0 -3.25513667985 -0.0758851894693
230 229.0 -3.17925670492 -0.0758632951011
231 230.0 -3.10343284123 -0.0757731496217
232 231.0 -3.02773277394 -0.0756159028468
233 232.0 -2.95222298559 -0.0753928092342
234 233.0 -2.87696865416 -0.0751052226812
235 234.0 -2.80203355622 -0.0747545911191
236 235.0 -2.72747997572 -0.0743424509249
237 236.0 -2.65336861841 -0.073870421164
238 237.0 -2.57975853208 -0.0733401976859
239 238.0 -2.50670703279 -0.0727535470871
240 239.0 -2.4342696372 -0.0721123005638
241 240.0 -2.36250000104 -0.0714183476691
242 241.0 -2.29144986396 -0.0706736299971
243 242.0 -2.22116900065 -0.0698801348102
244 243.0 -2.15170517837 -0.0690398886302
245 244.0 -2.0831041209 -0.0681549508121
246 245.0 -2.01540947892 -0.067227407119
247 246.0 -1.94866280684 -0.0662593633171
248 247.0 -1.88290354594 -0.0652529388105
249 248.0 -1.81816901389 -0.0642102603325
250 249.0 -1.7544944006 -0.0631334557138
251 250.0 -1.69191277013 -0.0620246477436
252 251.0 -1.6304550688 -0.0608859481423
253 252.0 -1.57015013921 -0.059719451663
254 253.0 -1.51102474011 -0.0585272303374
255 254.0 -1.45310357187 -0.0573113278834
256 255.0 -1.39640930762 -0.0560737542899
257 256.0 -1.34096262951 -0.054816480593
258 257.0 -1.28678227024 -0.0535414338587
259 258.0 -1.23388505944 -0.0522504923856
260 259.0 -1.18228597475 -0.0509454811405
261 260.0 -1.13199819729 -0.0496281674395
262 261.0 -1.08303317143 -0.0483002568854
263 262.0 -1.03540066834 -0.046963389572
264 263.0 -0.989108853377 -0.0456191365664
265 264.0 -0.944164356669 -0.0442689966762
266 265.0 -0.900572346917 -0.0429143935113
267 266.0 -0.858336607922 -0.0415566728462
268 267.0 -0.817459617608 -0.0401971002897
269 268.0 -0.777942629232 -0.0388368592669
270 269.0 -0.739785754436 -0.0374770493178
271 270.0 -0.702988047855 -0.0361186847156
272 271.0 -0.667547592939 -0.0347626934072
273 272.0 -0.633461588675 -0.0334099162773
274 273.0 -0.600726436882 -0.0320611067354
275 274.0 -0.569337829756 -0.0307169306269
276 275.0 -0.539290837348 -0.0293779664649
277 276.0 -0.510579994645 -0.0280447059807
278 277.0 -0.483199387947 -0.0267175549897
279 278.0 -0.457142740217 -0.0253968345674
280 279.0 -0.432403495111 -0.0240827825309
281 280.0 -0.408974899365 -0.0227755552188
282 281.0 -0.386850083265 -0.0214752295619
283 282.0 -0.366022138902 -0.020181805438
284 283.0 -0.346484195932 -0.0188952082997
285 284.0 -0.328229494574 -0.0176152920667
286 285.0 -0.311251455597 -0.0163418422722
287 286.0 -0.295543747024 -0.0150745794496
288 287.0 -0.28110034735 -0.0138131627512
289 288.0 -0.267915605017 -0.0125571937823
290 289.0 -0.255984293962 -0.011306220639
291 290.0 -0.245301665026 -0.0100597421363
292 291.0 -0.235863493049 -0.00881721220956
293 292.0 -0.22766611948 -0.00757804447631
294 293.0 -0.220706490355 -0.00634161694135
295 294.0 -0.214982189503 -0.00510727682957
296 295.0 -0.210491466861 -0.00387434552992
297 296.0 -0.207233261801 -0.00264212363344
298 297.0 -0.205207221373 -0.00140989604849
299 298.0 -0.204413713408 -0.00017693717569
300 299.0 -0.204853834414 0.0010574838751
301 300.0 -0.206529412255 0.00229409804323
302 301.0 -0.209443003569 0.00353363106913
303 302.0 -0.213597885954 0.00477679825726
304 303.0 -0.218998044922 0.00602429926791
305 304.0 -0.22564815567 0.00727681295572
306 305.0 -0.23355355972 0.00853499227222
307 306.0 -0.2427202365 0.00979945924997
308 307.0 -0.253154769958 0.0110708000854
309 308.0 -0.264864310313 0.0123495603372
310 309.0 -0.277856531075 0.0136362402565
311 310.0 -0.292139581459 0.0149312902659
312 311.0 -0.307722034364 0.0162351066015
313 312.0 -0.324612830087 0.0175480271349
314 313.0 -0.342821215943 0.0188703273888
315 314.0 -0.362356682012 0.0202022167596
316 315.0 -0.383228893218 0.0215438349636
317 316.0 -0.405447617967 0.0228952487148
318 317.0 -0.429022653586 0.0242564486517
319 318.0 -0.45396374882 0.0256273465206
320 319.0 -0.480280523637 0.0270077726275
321 320.0 -0.507982386639 0.0283974735696
322 321.0 -0.537078450328 0.029796110253
323 322.0 -0.567577444555 0.0312032562068
324 323.0 -0.59948762842 0.0326183962009
325 324.0 -0.632816700956 0.0340409251716
326 325.0 -0.667571710883 0.0354701474639
327 326.0 -0.703758965776 0.0369052763923
328 327.0 -0.741383940946 0.038345434125
329 328.0 -0.780451188376 0.0397896518935
330 329.0 -0.820964246018 0.0412368705304
331 330.0 -0.862925547807 0.042685941334
332 331.0 -0.906336334692 0.0441356272615
333 332.0 -0.951196567028 0.045584604448
334 333.0 -0.997504838648 0.0470314640498
335 334.0 -1.04525829294 0.048474714408
336 335.0 -1.09445254125 0.0499127835288
337 336.0 -1.1450815839 0.0513440218749
338 337.0 -1.1971377342 0.0527667054614
339 338.0 -1.25061154564 0.0541790392498
340 339.0 -1.30549174267 0.0555791608316
341 340.0 -1.36176515529 0.0569651443923
342 341.0 -1.41941665773 0.0583350049463
343 342.0 -1.47842911151 0.0596867028317
344 343.0 -1.53878331313 0.061018148454
345 344.0 -1.60045794659 0.0623272072653
346 345.0 -1.66342954101 0.0636117049668
347 346.0 -1.72767243359 0.0648694329207
348 347.0 -1.79315873807 0.0660981537565
349 348.0 -1.85985831882 0.0672956071568
350 349.0 -1.92773877092 0.0684595158069
351 350.0 -1.99676540616 0.0695875914917
352 351.0 -2.06690124527 0.0706775413231
353 352.0 -2.13810701636 0.0717270740805
354 353.0 -2.21034115987 0.0727339066469
355 354.0 -2.28355983986 0.0736957705223
356 355.0 -2.35771696194 0.0746104183955
357 356.0 -2.43276419776 0.0754756307561
358 357.0 -2.50865101613 0.0762892225281
359 358.0 -2.58532472075 0.0770490497051
360 359.0 -2.66273049463 0.0777530159679
# Table of the potential and its negative derivative for frustrated beta sheet
# (Note: Derivatives are in units of energy/radians, not energy/degrees.)
# ./calc_dihedral_table.py 5.6 57.29577951308232 6 6.0 180 6 0.0 359 360
FRUSTRATED_BETA
N 360 DEGREES
1 0.0 -2.55809068762 0.0731724739818
2 1.0 -2.63154144494 0.0737195744566
3 2.0 -2.70551060968 0.0742089966437
4 3.0 -2.77993963883 0.074639023134
5 4.0 -2.85476830901 0.0750080115297
6 5.0 -2.92993479441 0.0753144003899
7 6.0 -3.00537575069 0.0755567150326
8 7.0 -3.08102640456 0.0757335731758
9 8.0 -3.15682064892 0.0758436903983
10 9.0 -3.23269114341 0.075885885404
11 10.0 -3.30856942003 0.0758590850738
12 11.0 -3.38438599377 0.0757623292865
13 12.0 -3.46007047791 0.0755947754951
14 13.0 -3.53555170381 0.0753557030426
15 14.0 -3.61075784476 0.0750445172025
16 15.0 -3.68561654392 0.0746607529305
17 16.0 -3.76005504566 0.0742040783151
18 17.0 -3.83400033034 0.0736742977129
19 18.0 -3.907379252 0.0730713545594
20 19.0 -3.98011867868 0.0723953338429
21 20.0 -4.0521456351 0.0716464642332
22 21.0 -4.12338744726 0.0708251198546
23 22.0 -4.19377188857 0.0699318216967
24 23.0 -4.26322732737 0.0689672386556
25 24.0 -4.33168287509 0.0679321881993
26 25.0 -4.39906853508 0.0668276366524
27 26.0 -4.46531535141 0.0656546990963
28 27.0 -4.53035555742 0.0644146388823
29 28.0 -4.59412272358 0.0631088667546
30 29.0 -4.65655190431 0.061738939584
31 30.0 -4.71757978327 0.0603065587109
32 31.0 -4.77714481686 0.0588135679005
33 32.0 -4.83518737548 0.057261950911
34 33.0 -4.89164988211 0.0556538286799
35 34.0 -4.94647694795 0.0539914561312
36 35.0 -4.99961550465 0.0522772186102
37 36.0 -5.05101493277 0.0505136279528
38 37.0 -5.10062718621 0.048703318195
39 38.0 -5.14840691207 0.0468490409338
40 39.0 -5.19431156578 0.0449536603471
41 40.0 -5.23830152101 0.0430201478838
42 41.0 -5.28034017422 0.0410515766363
43 42.0 -5.3203940433 0.0390511154063
44 43.0 -5.35843286021 0.0370220224793
45 44.0 -5.39442965726 0.0349676391193
46 45.0 -5.4283608467 0.0328913828015
47 46.0 -5.46020629342 0.0307967401964
48 47.0 -5.48994938059 0.028687259923
49 48.0 -5.51757706789 0.0265665450883
50 49.0 -5.54307994213 0.0244382456298
51 50.0 -5.56645226024 0.0223060504811
52 51.0 -5.58769198425 0.0201736795783
53 52.0 -5.60680080825 0.0180448757265
54 53.0 -5.62378417713 0.0159233963481
55 54.0 -5.63865129702 0.0138130051308
56 55.0 -5.6514151374 0.0117174635982
57 56.0 -5.66209242462 0.00964052262251
58 57.0 -5.67070362704 0.00758591390103
59 58.0 -5.67727293157 0.00555734141841
60 59.0 -5.6818282117 0.00355847291538
61 60.0 -5.68440098698 0.00159293138608
62 61.0 -5.68502637408 -0.000335713374531
63 62.0 -5.68374302934 -0.00222395315148
64 63.0 -5.68059308309 -0.0040683495974
65 64.0 -5.67562206565 -0.00586554240548
66 65.0 -5.66887882528 -0.00761225734683
67 66.0 -5.66041543813 -0.00930531415106
68 67.0 -5.65028711044 -0.0109416342099
69 68.0 -5.63855207307 -0.0125182480831
70 69.0 -5.6252714687 -0.0140323027883
71 70.0 -5.61050923182 -0.0154810688529
72 71.0 -5.59433196178 -0.0168619471125
73 72.0 -5.57680878923 -0.0181724752358
74 73.0 -5.5580112361 -0.019410333958
75 74.0 -5.53801306959 -0.0205733530082
76 75.0 -5.51689015031 -0.0216595167121
77 76.0 -5.49472027505 -0.0226669692568
78 77.0 -5.47158301441 -0.0235940196022
79 78.0 -5.44755954575 -0.0244391460249
80 79.0 -5.42273248172 -0.0252010002837
81 80.0 -5.3971856949 -0.0258784113929
82 81.0 -5.37100413881 -0.0264703889936
83 82.0 -5.34427366574 -0.0269761263135
84 83.0 -5.31708084192 -0.0273950027051
85 84.0 -5.28951276022 -0.0277265857564
86 85.0 -5.26165685114 -0.0279706329651
87 86.0 -5.23360069216 -0.0281270929735
88 87.0 -5.20543181621 -0.0281961063563
89 88.0 -5.17723751951 -0.0281780059613
90 89.0 -5.14910466934 -0.0280733167983
91 90.0 -5.12111951208 -0.0278827554757
92 91.0 -5.09336748214 -0.0276072291861
93 92.0 -5.06593301201 -0.0272478342399
94 93.0 -5.0388993441 -0.026805854151
95 94.0 -5.01234834466 -0.0262827572773
96 95.0 -4.98636032033 -0.0256801940208
97 96.0 -4.96101383762 -0.0249999935924
98 97.0 -4.93638554598 -0.0242441603499
99 98.0 -4.91255000457 -0.0234148697145
100 99.0 -4.88957951348 -0.0225144636776
101 100.0 -4.86754394953 -0.0215454459053
102 101.0 -4.84651060724 -0.0205104764546
103 102.0 -4.8265440452 -0.01941236611
104 103.0 -4.80770593836 -0.0182540703564
105 104.0 -4.79005493648 -0.0170386830008
106 105.0 -4.77364652914 -0.0157694294583
107 106.0 -4.7585329176 -0.0144496597171
108 107.0 -4.74476289391 -0.0130828410011
109 108.0 -4.73238172744 -0.0116725501446
110 109.0 -4.72143105919 -0.0102224657007
111 110.0 -4.71194880414 -0.00873635979846
112 111.0 -4.70396906182 -0.0072180897712
113 112.0 -4.69752203541 -0.00567158957449
114 113.0 -4.69263395945 -0.00410086101469
115 114.0 -4.68932703648 -0.00250996480925
116 115.0 -4.68761938265 -0.000903011500147
117 116.0 -4.68752498248 0.00071584775762
118 117.0 -4.68905365291 0.00234243051027
119 118.0 -4.69221101668 0.00397253239976
120 119.0 -4.69699848518 0.00560193661579
121 120.0 -4.70341325069 0.00722642338265
122 121.0 -4.71144828821 0.00884177945771
123 122.0 -4.72109236669 0.0104438076188
124 123.0 -4.73233006984 0.0120283361174
125 124.0 -4.74514182625 0.0135912280748
126 125.0 -4.75950394898 0.0151283907985
127 126.0 -4.77538868431 0.0166357849963
128 127.0 -4.79276426974 0.0181094338658
129 128.0 -4.81159500092 0.0195454320375
130 129.0 -4.83184130754 0.0209399543498
131 130.0 -4.8534598378 0.0222892644342
132 131.0 -4.87640355143 0.0235897230915
133 132.0 -4.90062182095 0.0248377964369
134 133.0 -4.92606054096 0.0260300637961
135 134.0 -4.95266224518 0.0271632253326
136 135.0 -4.98036623096 0.028234109388
137 136.0 -5.00910869107 0.0292396795182
138 137.0 -5.03882285221 0.0301770412082
139 138.0 -5.06943912022 0.0310434482505
140 139.0 -5.10088523142 0.0318363087705
141 140.0 -5.13308640979 0.0325531908865
142 141.0 -5.16596552963 0.0331918279898
143 142.0 -5.19944328334 0.0337501236332
144 143.0 -5.23343835383 0.0342261560164
145 144.0 -5.26786759123 0.0346181820585
146 145.0 -5.30264619353 0.0349246410472
147 146.0 -5.33768789051 0.0351441578585
148 147.0 -5.37290513082 0.0352755457383
149 148.0 -5.40820927152 0.0353178086401
150 149.0 -5.4435107698 0.0352701431151
151 150.0 -5.4787193763 0.0351319397498
152 151.0 -5.51374432971 0.0349027841491
153 152.0 -5.54849455206 0.0345824574643
154 153.0 -5.58287884436 0.0341709364636
155 154.0 -5.61680608206 0.0336683931487
156 155.0 -5.65018540988 0.0330751939177
157 156.0 -5.68292643563 0.0323918982779
158 157.0 -5.71493942249 0.0316192571138
159 158.0 -5.74613547931 0.0307582105139
160 159.0 -5.77642674856 0.029809885165
161 160.0 -5.80572659147 0.0287755913197
162 161.0 -5.83394976986 0.0276568193473
163 162.0 -5.86101262442 0.0264552358763
164 163.0 -5.8868332488 0.025172679541
165 164.0 -5.91133165941 0.0238111563427
166 165.0 -5.93442996024 0.0223728346376
167 166.0 -5.95605250261 0.0208600397671
168 167.0 -5.97612603931 0.0192752483425
169 168.0 -5.99457987285 0.0176210822011
170 169.0 -6.01134599757 0.015900302049
171 170.0 -6.02635923519 0.014115800807
172 171.0 -6.03955736358 0.0122705966784
173 172.0 -6.05088123845 0.0103678259555
174 173.0 -6.0602749078 0.00841073558436
175 174.0 -6.06768571866 0.00640267550713
176 175.0 -6.0730644163 0.00434709080102
177 176.0 -6.07636523524 0.00224751363529
178 177.0 -6.07754598232 0.000107555066143
179 178.0 -6.07656811141 -0.00206910330914
180 179.0 -6.07339678973 -0.00427871781763
181 180.0 -6.06800095563 -0.00651749127408
182 181.0 -6.06035336781 -0.00878158162059
183 182.0 -6.05043064586 -0.0110671106207
184 183.0 -6.03821330204 -0.0133701725859
185 184.0 -6.02368576439 -0.0156868431131
186 185.0 -6.00683639108 -0.0180131878107
187 186.0 -5.98765747603 -0.0203452709919
188 187.0 -5.96614524589 -0.0226791643135
189 188.0 -5.94229984843 -0.025010955339
190 189.0 -5.91612533236 -0.0273367560054
191 190.0 -5.88762961878 -0.0296527109716
192 191.0 -5.85682446433 -0.0319550058299
193 192.0 -5.82372541626 -0.0342398751598
194 193.0 -5.78835175943 -0.0365036104045
195 194.0 -5.75072645562 -0.0387425675516
196 195.0 -5.71087607524 -0.0409531746008
197 196.0 -5.66883072166 -0.0431319387984
198 197.0 -5.62462394846 -0.0452754536249
199 198.0 -5.57829266983 -0.0473804055171
200 199.0 -5.5298770643 -0.0494435803104
201 200.0 -5.47942047235 -0.0514618693867
202 201.0 -5.42696928781 -0.0534322755136
203 202.0 -5.37257284377 -0.055351918363
204 203.0 -5.316283293 -0.0572180396955
205 204.0 -5.25815548345 -0.059028008202
206 205.0 -5.19824682901 -0.0607793239895
207 206.0 -5.13661717604 -0.0624696227052
208 207.0 -5.0733286659 -0.0640966792879
209 208.0 -5.00844559393 -0.0656584113417
210 209.0 -4.94203426529 -0.0671528821253
211 210.0 -4.87416284794 -0.0685783031513
212 211.0 -4.80490122327 -0.0699330363936
213 212.0 -4.7343208347 -0.0712155960973
214 213.0 -4.66249453466 -0.0724246501921
215 214.0 -4.58949643037 -0.0735590213066
216 215.0 -4.51540172879 -0.0746176873849
217 216.0 -4.44028658118 -0.0755997819067
218 217.0 -4.3642279276 -0.0765045937139
219 218.0 -4.28730334182 -0.0773315664459
220 219.0 -4.20959087694 -0.0780802975905
221 220.0 -4.13116891218 -0.0787505371538
222 221.0 -4.0521160012 -0.0793421859574
223 222.0 -3.97251072229 -0.0798552935693
224 223.0 -3.89243153076 -0.0802900558785
225 224.0 -3.81195661404 -0.0806468123209
226 225.0 -3.73116374964 -0.0809260427693
227 226.0 -3.65013016636 -0.0811283640964
228 227.0 -3.56893240921 -0.0812545264246
229 228.0 -3.48764620813 -0.0813054090744
230 229.0 -3.4063463509 -0.0812820162266
231 230.0 -3.32510656064 -0.0811854723104
232 231.0 -3.24399937793 -0.081017017134
233 232.0 -3.16309604794 -0.0807780007742
234 233.0 -3.08246641287 -0.0804698782381
235 234.0 -3.00217880976 -0.0800942039176
236 235.0 -2.92229997393 -0.079652625851
237 236.0 -2.84289494829 -0.0791468798106
238 237.0 -2.76402699866 -0.0785787832348
239 238.0 -2.68575753514 -0.0779502290223
240 239.0 -2.60814603984 -0.077263179207
241 240.0 -2.53125000097 -0.0765196585342
242 241.0 -2.4551248533 -0.0757217479546
243 242.0 -2.37982392531 -0.0748715780578
244 243.0 -2.30539839282 -0.073971322463
245 244.0 -2.23189723927 -0.0730231911866
246 245.0 -2.15936722267 -0.072029424007
247 246.0 -2.0878528491 -0.0709922838436
248 247.0 -2.01739635293 -0.0699140501714
249 248.0 -1.94803768347 -0.0687970124882
250 249.0 -1.87981449824 -0.0676434638537
251 250.0 -1.81276216256 -0.0664556945194
252 251.0 -1.74691375554 -0.0652359856651
253 252.0 -1.68230008218 -0.0639866032624
254 253.0 -1.61894969164 -0.0627097920793
255 254.0 -1.55688890134 -0.0614077698443
256 255.0 -1.49614182687 -0.0600827215855
257 256.0 -1.43673041741 -0.05873679416
258 257.0 -1.37867449659 -0.0573720909874
259 258.0 -1.32199180845 -0.0559906670036
260 259.0 -1.26669806833 -0.0545945238457
261 260.0 -1.21280701853 -0.0531856052829
262 261.0 -1.1603304883 -0.0517657929031
263 262.0 -1.1092784581 -0.0503369020679
264 263.0 -1.05965912771 -0.0489006781451
265 264.0 -1.01147898802 -0.0474587930279
266 265.0 -0.964742896092 -0.0460128419505
267 266.0 -0.919454153297 -0.0445643406057
268 267.0 -0.875614586172 -0.0431147225719
269 268.0 -0.833224629688 -0.0416653370554
270 269.0 -0.792283412613 -0.0402174469521
271 270.0 -0.752788844664 -0.038772227232
272 271.0 -0.714737705101 -0.0373307636499
273 272.0 -0.67812573245 -0.0358940517831
274 273.0 -0.642947715028 -0.0344629963972
275 274.0 -0.609197581934 -0.0330384111393
276 275.0 -0.576868494182 -0.0316210185584
277 276.0 -0.545952935658 -0.0302114504483
278 277.0 -0.51644280357 -0.0288102485125
279 278.0 -0.488329498068 -0.0274178653447
280 279.0 -0.461604010741 -0.0260346657211
281 280.0 -0.436257011655 -0.0246609281969
282 281.0 -0.412278934657 -0.023296847002
283 282.0 -0.389660060626 -0.0219425342253
284 283.0 -0.368390598407 -0.0205980222818
285 284.0 -0.348460763137 -0.01926326665
286 285.0 -0.329860851704 -0.0179381488715
287 286.0 -0.312581315078 -0.0166224797996
288 287.0 -0.296612827279 -0.015316003087
289 288.0 -0.281946350734 -0.0140183988977
290 289.0 -0.268573197826 -0.0127292878319
291 290.0 -0.256485088408 -0.0114482350481
292 291.0 -0.245674203109 -0.0101747545698
293 292.0 -0.236133232246 -0.00890831375923
294 293.0 -0.227855420178 -0.00764833794542
295 294.0 -0.220834604976 -0.00639421518813
296 295.0 -0.215065253253 -0.00514530116277
297 296.0 -0.210542490065 -0.00390092414876
298 297.0 -0.207262123775 -0.00266039010467
299 298.0 -0.205220665805 -0.00142298781263
300 299.0 -0.204415345223 -0.000187994074493
301 300.0 -0.204844118104 0.00104532105779
302 301.0 -0.206505671662 0.00227768903543
303 302.0 -0.209399423126 0.0035098375675
304 303.0 -0.213525513386 0.00474248539479
305 304.0 -0.218884795423 0.00597633710062
306 305.0 -0.225478817581 0.00721207797616
307 306.0 -0.233309801737 0.00845036895769
308 307.0 -0.242380616448 0.00969184165314
309 308.0 -0.252694745185 0.0109370934746
310 309.0 -0.264256249747 0.0121866828936
311 310.0 -0.277069729013 0.0134411248358
312 311.0 -0.291140273151 0.0147008862297
313 312.0 -0.306473413467 0.0159663817261
314 313.0 -0.323075068066 0.0172379696031
315 314.0 -0.340951483513 0.018515947869
316 315.0 -0.360109172702 0.0198005505798
317 316.0 -0.380554849155 0.0210919443819
318 317.0 -0.402295357987 0.0223902252933
319 318.0 -0.425337603767 0.0236954157356
320 319.0 -0.449688475549 0.0250074618263
321 320.0 -0.475354769327 0.0263262309427
322 321.0 -0.50234310819 0.0276515095659
323 322.0 -0.530659860472 0.0289830014145
324 323.0 -0.560311056174 0.0303203258736
325 324.0 -0.59130230198 0.0316630167284
326 325.0 -0.623638695141 0.0330105212056
327 326.0 -0.657324736579 0.0343621993296
328 327.0 -0.692364243488 0.0357173235955
329 328.0 -0.728760261774 0.0370750789637
330 329.0 -0.766514978659 0.0384345631765
331 330.0 -0.805629635748 0.0397947873984
332 331.0 -0.846104442913 0.04115467718
333 332.0 -0.887938493289 0.042513073745
334 333.0 -0.93112967973 0.0438687355968
335 334.0 -0.975674613021 0.0452203404434
336 335.0 -1.02156854218 0.0465664874361
337 336.0 -1.06880527714 0.0479056997168
338 337.0 -1.11737711415 0.0492364272675
339 338.0 -1.16727476416 0.0505570500574
340 339.0 -1.2184872845 0.051865881477
341 340.0 -1.27100201415 0.0531611720525
342 341.0 -1.32480451282 0.0544411134304
343 342.0 -1.37987850417 0.055703842622
344 343.0 -1.43620582346 0.0569474464963
345 344.0 -1.49376636966 0.0581699665097
346 345.0 -1.55253806258 0.05936940366
347 346.0 -1.61249680493 0.0605437236497
348 347.0 -1.67361644969 0.0616908622471
349 348.0 -1.73586877296 0.0628087308273
350 349.0 -1.79922345238 0.0638952220804
351 350.0 -1.86364805137 0.0649482158688
352 351.0 -1.92910800931 0.0659655852184
353 352.0 -1.9955666377 0.066945202426
354 353.0 -2.06298512258 0.0678849452658
355 354.0 -2.13132253309 0.0687827032771
356 355.0 -2.20053583647 0.0696363841147
357 356.0 -2.27057991931 0.0704439199439
358 357.0 -2.3414076153 0.0712032738621
359 358.0 -2.41296973939 0.0719124463259
360 359.0 -2.48521512832 0.072569481568

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