Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85562598
generate_simple.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Sep 29, 23:20
Size
9 KB
Mime Type
text/x-python
Expires
Tue, Oct 1, 23:20 (2 d)
Engine
blob
Format
Raw Data
Handle
21199734
Attached To
rLAMMPS lammps
generate_simple.py
View Options
# Setup tool for oxDNA input in LAMMPS format.
import
math
,
numpy
as
np
,
sys
,
os
# system size
lxmin
=
-
115.0
lxmax
=
+
115.0
lymin
=
-
115.0
lymax
=
+
115.0
lzmin
=
-
115.0
lzmax
=
+
115.0
# rise in z-direction
r0
=
0.7
# definition of single untwisted strand
def
single
():
strand
=
inp
[
1
]
.
split
(
':'
)
com_start
=
strand
[
0
]
.
split
(
','
)
posx
=
float
(
com_start
[
0
])
posy
=
float
(
com_start
[
1
])
posz
=
float
(
com_start
[
2
])
risex
=
0
risey
=
0
risez
=
r0
strandstart
=
len
(
nucleotide
)
+
1
for
letter
in
strand
[
1
]:
temp
=
[]
temp
.
append
(
nt2num
[
letter
])
temp
.
append
([
posx
,
posy
,
posz
])
vel
=
[
0
,
0
,
0
,
0
,
0
,
0
]
temp
.
append
(
vel
)
temp
.
append
(
shape
)
quat
=
[
1
,
0
,
0
,
0
]
temp
.
append
(
quat
)
posx
=
posx
+
risex
posy
=
posy
+
risey
posz
=
posz
+
risez
if
(
len
(
nucleotide
)
+
1
>
strandstart
):
topology
.
append
([
1
,
len
(
nucleotide
),
len
(
nucleotide
)
+
1
])
nucleotide
.
append
(
temp
)
return
# definition of single twisted strand
def
single_helix
():
strand
=
inp
[
1
]
.
split
(
':'
)
com_start
=
strand
[
0
]
.
split
(
','
)
twist
=
0.6
posx
=
float
(
com_start
[
0
])
posy
=
float
(
com_start
[
1
])
posz
=
float
(
com_start
[
2
])
risex
=
0
risey
=
0
risez
=
math
.
sqrt
(
r0
**
2
-
4.0
*
math
.
sin
(
0.5
*
twist
)
**
2
)
dcomh
=
0.76
axisx
=
dcomh
+
posx
axisy
=
posy
strandstart
=
len
(
nucleotide
)
+
1
quat
=
[
1
,
0
,
0
,
0
]
qrot0
=
math
.
cos
(
0.5
*
twist
)
qrot1
=
0
qrot2
=
0
qrot3
=
math
.
sin
(
0.5
*
twist
)
for
letter
in
strand
[
1
]:
temp
=
[]
temp
.
append
(
nt2num
[
letter
])
temp
.
append
([
posx
,
posy
,
posz
])
vel
=
[
0
,
0
,
0
,
0
,
0
,
0
]
temp
.
append
(
vel
)
temp
.
append
(
shape
)
temp
.
append
(
quat
)
quat0
=
quat
[
0
]
*
qrot0
-
quat
[
1
]
*
qrot1
-
quat
[
2
]
*
qrot2
-
quat
[
3
]
*
qrot3
quat1
=
quat
[
0
]
*
qrot1
+
quat
[
1
]
*
qrot0
+
quat
[
2
]
*
qrot3
-
quat
[
3
]
*
qrot2
quat2
=
quat
[
0
]
*
qrot2
+
quat
[
2
]
*
qrot0
+
quat
[
3
]
*
qrot1
-
quat
[
1
]
*
qrot3
quat3
=
quat
[
0
]
*
qrot3
+
quat
[
3
]
*
qrot0
+
quat
[
1
]
*
qrot2
+
quat
[
2
]
*
qrot1
quat
=
[
quat0
,
quat1
,
quat2
,
quat3
]
posx
=
axisx
-
dcomh
*
(
quat
[
0
]
**
2
+
quat
[
1
]
**
2
-
quat
[
2
]
**
2
-
quat
[
3
]
**
2
)
posy
=
axisy
-
dcomh
*
(
2
*
(
quat
[
1
]
*
quat
[
2
]
+
quat
[
0
]
*
quat
[
3
]))
posz
=
posz
+
risez
if
(
len
(
nucleotide
)
+
1
>
strandstart
):
topology
.
append
([
1
,
len
(
nucleotide
),
len
(
nucleotide
)
+
1
])
nucleotide
.
append
(
temp
)
return
# definition of twisted duplex
def
duplex
():
strand
=
inp
[
1
]
.
split
(
':'
)
com_start
=
strand
[
0
]
.
split
(
','
)
twist
=
0.6
compstrand
=
[]
comptopo
=
[]
posx1
=
float
(
com_start
[
0
])
posy1
=
float
(
com_start
[
1
])
posz1
=
float
(
com_start
[
2
])
risex
=
0
risey
=
0
risez
=
math
.
sqrt
(
r0
**
2
-
4.0
*
math
.
sin
(
0.5
*
twist
)
**
2
)
dcomh
=
0.76
axisx
=
dcomh
+
posx1
axisy
=
posy1
posx2
=
axisx
+
dcomh
posy2
=
posy1
posz2
=
posz1
strandstart
=
len
(
nucleotide
)
+
1
quat1
=
[
1
,
0
,
0
,
0
]
quat2
=
[
0
,
0
,
-
1
,
0
]
qrot0
=
math
.
cos
(
0.5
*
twist
)
qrot1
=
0
qrot2
=
0
qrot3
=
math
.
sin
(
0.5
*
twist
)
for
letter
in
strand
[
1
]:
temp1
=
[]
temp2
=
[]
temp1
.
append
(
nt2num
[
letter
])
temp2
.
append
(
compnt2num
[
letter
])
temp1
.
append
([
posx1
,
posy1
,
posz1
])
temp2
.
append
([
posx2
,
posy2
,
posz2
])
vel
=
[
0
,
0
,
0
,
0
,
0
,
0
]
temp1
.
append
(
vel
)
temp2
.
append
(
vel
)
temp1
.
append
(
shape
)
temp2
.
append
(
shape
)
temp1
.
append
(
quat1
)
temp2
.
append
(
quat2
)
quat1_0
=
quat1
[
0
]
*
qrot0
-
quat1
[
1
]
*
qrot1
-
quat1
[
2
]
*
qrot2
-
quat1
[
3
]
*
qrot3
quat1_1
=
quat1
[
0
]
*
qrot1
+
quat1
[
1
]
*
qrot0
+
quat1
[
2
]
*
qrot3
-
quat1
[
3
]
*
qrot2
quat1_2
=
quat1
[
0
]
*
qrot2
+
quat1
[
2
]
*
qrot0
+
quat1
[
3
]
*
qrot1
-
quat1
[
1
]
*
qrot3
quat1_3
=
quat1
[
0
]
*
qrot3
+
quat1
[
3
]
*
qrot0
+
quat1
[
1
]
*
qrot2
+
quat1
[
2
]
*
qrot1
quat1
=
[
quat1_0
,
quat1_1
,
quat1_2
,
quat1_3
]
posx1
=
axisx
-
dcomh
*
(
quat1
[
0
]
**
2
+
quat1
[
1
]
**
2
-
quat1
[
2
]
**
2
-
quat1
[
3
]
**
2
)
posy1
=
axisy
-
dcomh
*
(
2
*
(
quat1
[
1
]
*
quat1
[
2
]
+
quat1
[
0
]
*
quat1
[
3
]))
posz1
=
posz1
+
risez
quat2_0
=
quat2
[
0
]
*
qrot0
-
quat2
[
1
]
*
qrot1
-
quat2
[
2
]
*
qrot2
+
quat2
[
3
]
*
qrot3
quat2_1
=
quat2
[
0
]
*
qrot1
+
quat2
[
1
]
*
qrot0
-
quat2
[
2
]
*
qrot3
-
quat2
[
3
]
*
qrot2
quat2_2
=
quat2
[
0
]
*
qrot2
+
quat2
[
2
]
*
qrot0
+
quat2
[
3
]
*
qrot1
+
quat2
[
1
]
*
qrot3
quat2_3
=-
quat2
[
0
]
*
qrot3
+
quat2
[
3
]
*
qrot0
+
quat2
[
1
]
*
qrot2
+
quat2
[
2
]
*
qrot1
quat2
=
[
quat2_0
,
quat2_1
,
quat2_2
,
quat2_3
]
posx2
=
axisx
+
dcomh
*
(
quat1
[
0
]
**
2
+
quat1
[
1
]
**
2
-
quat1
[
2
]
**
2
-
quat1
[
3
]
**
2
)
posy2
=
axisy
+
dcomh
*
(
2
*
(
quat1
[
1
]
*
quat1
[
2
]
+
quat1
[
0
]
*
quat1
[
3
]))
posz2
=
posz1
if
(
len
(
nucleotide
)
+
1
>
strandstart
):
topology
.
append
([
1
,
len
(
nucleotide
),
len
(
nucleotide
)
+
1
])
comptopo
.
append
([
1
,
len
(
nucleotide
)
+
len
(
strand
[
1
]),
len
(
nucleotide
)
+
len
(
strand
[
1
])
+
1
])
nucleotide
.
append
(
temp1
)
compstrand
.
append
(
temp2
)
for
ib
in
range
(
len
(
compstrand
)):
nucleotide
.
append
(
compstrand
[
len
(
compstrand
)
-
1
-
ib
])
for
ib
in
range
(
len
(
comptopo
)):
topology
.
append
(
comptopo
[
ib
])
return
# definition of array of duplexes
def
duplex_array
():
strand
=
inp
[
1
]
.
split
(
':'
)
number
=
strand
[
0
]
.
split
(
','
)
posz1_0
=
float
(
strand
[
1
])
twist
=
0.6
nx
=
int
(
number
[
0
])
ny
=
int
(
number
[
1
])
dx
=
(
lxmax
-
lxmin
)
/
nx
dy
=
(
lymax
-
lymin
)
/
ny
risex
=
0
risey
=
0
risez
=
math
.
sqrt
(
r0
**
2
-
4.0
*
math
.
sin
(
0.5
*
twist
)
**
2
)
dcomh
=
0.76
for
ix
in
range
(
nx
):
axisx
=
lxmin
+
dx
/
2
+
ix
*
dx
for
iy
in
range
(
ny
):
axisy
=
lymin
+
dy
/
2
+
iy
*
dy
compstrand
=
[]
comptopo
=
[]
posx1
=
axisx
-
dcomh
posy1
=
axisy
posz1
=
posz1_0
posx2
=
axisx
+
dcomh
posy2
=
posy1
posz2
=
posz1
strandstart
=
len
(
nucleotide
)
+
1
quat1
=
[
1
,
0
,
0
,
0
]
quat2
=
[
0
,
0
,
-
1
,
0
]
qrot0
=
math
.
cos
(
0.5
*
twist
)
qrot1
=
0
qrot2
=
0
qrot3
=
math
.
sin
(
0.5
*
twist
)
for
letter
in
strand
[
2
]:
temp1
=
[]
temp2
=
[]
temp1
.
append
(
nt2num
[
letter
])
temp2
.
append
(
compnt2num
[
letter
])
temp1
.
append
([
posx1
,
posy1
,
posz1
])
temp2
.
append
([
posx2
,
posy2
,
posz2
])
vel
=
[
0
,
0
,
0
,
0
,
0
,
0
]
temp1
.
append
(
vel
)
temp2
.
append
(
vel
)
temp1
.
append
(
shape
)
temp2
.
append
(
shape
)
temp1
.
append
(
quat1
)
temp2
.
append
(
quat2
)
quat1_0
=
quat1
[
0
]
*
qrot0
-
quat1
[
1
]
*
qrot1
-
quat1
[
2
]
*
qrot2
-
quat1
[
3
]
*
qrot3
quat1_1
=
quat1
[
0
]
*
qrot1
+
quat1
[
1
]
*
qrot0
+
quat1
[
2
]
*
qrot3
-
quat1
[
3
]
*
qrot2
quat1_2
=
quat1
[
0
]
*
qrot2
+
quat1
[
2
]
*
qrot0
+
quat1
[
3
]
*
qrot1
-
quat1
[
1
]
*
qrot3
quat1_3
=
quat1
[
0
]
*
qrot3
+
quat1
[
3
]
*
qrot0
+
quat1
[
1
]
*
qrot2
+
quat1
[
2
]
*
qrot1
quat1
=
[
quat1_0
,
quat1_1
,
quat1_2
,
quat1_3
]
posx1
=
axisx
-
dcomh
*
(
quat1
[
0
]
**
2
+
quat1
[
1
]
**
2
-
quat1
[
2
]
**
2
-
quat1
[
3
]
**
2
)
posy1
=
axisy
-
dcomh
*
(
2
*
(
quat1
[
1
]
*
quat1
[
2
]
+
quat1
[
0
]
*
quat1
[
3
]))
posz1
=
posz1
+
risez
quat2_0
=
quat2
[
0
]
*
qrot0
-
quat2
[
1
]
*
qrot1
-
quat2
[
2
]
*
qrot2
+
quat2
[
3
]
*
qrot3
quat2_1
=
quat2
[
0
]
*
qrot1
+
quat2
[
1
]
*
qrot0
-
quat2
[
2
]
*
qrot3
-
quat2
[
3
]
*
qrot2
quat2_2
=
quat2
[
0
]
*
qrot2
+
quat2
[
2
]
*
qrot0
+
quat2
[
3
]
*
qrot1
+
quat2
[
1
]
*
qrot3
quat2_3
=-
quat2
[
0
]
*
qrot3
+
quat2
[
3
]
*
qrot0
+
quat2
[
1
]
*
qrot2
+
quat2
[
2
]
*
qrot1
quat2
=
[
quat2_0
,
quat2_1
,
quat2_2
,
quat2_3
]
posx2
=
axisx
+
dcomh
*
(
quat1
[
0
]
**
2
+
quat1
[
1
]
**
2
-
quat1
[
2
]
**
2
-
quat1
[
3
]
**
2
)
posy2
=
axisy
+
dcomh
*
(
2
*
(
quat1
[
1
]
*
quat1
[
2
]
+
quat1
[
0
]
*
quat1
[
3
]))
posz2
=
posz1
if
(
len
(
nucleotide
)
+
1
>
strandstart
):
topology
.
append
([
1
,
len
(
nucleotide
),
len
(
nucleotide
)
+
1
])
comptopo
.
append
([
1
,
len
(
nucleotide
)
+
len
(
strand
[
2
]),
len
(
nucleotide
)
+
len
(
strand
[
2
])
+
1
])
nucleotide
.
append
(
temp1
)
compstrand
.
append
(
temp2
)
for
ib
in
range
(
len
(
compstrand
)):
nucleotide
.
append
(
compstrand
[
len
(
compstrand
)
-
1
-
ib
])
for
ib
in
range
(
len
(
comptopo
)):
topology
.
append
(
comptopo
[
ib
])
return
# main part
nt2num
=
{
'A'
:
1
,
'C'
:
2
,
'G'
:
3
,
'T'
:
4
}
compnt2num
=
{
'T'
:
1
,
'G'
:
2
,
'C'
:
3
,
'A'
:
4
}
shape
=
[
1.1739845031423408
,
1.1739845031423408
,
1.1739845031423408
]
nucleotide
=
[]
topology
=
[]
seqfile
=
open
(
sys
.
argv
[
1
],
'r'
)
# process sequence file line by line
for
line
in
seqfile
:
inp
=
line
.
split
()
if
inp
[
0
]
==
'single'
:
single
()
if
inp
[
0
]
==
'single_helix'
:
single_helix
()
if
inp
[
0
]
==
'duplex'
:
duplex
()
if
inp
[
0
]
==
'duplex_array'
:
duplex_array
()
# output atom data in LAMMPS format
out
=
open
(
sys
.
argv
[
2
],
'w'
)
out
.
write
(
'# LAMMPS data file
\n
'
)
out
.
write
(
'
%d
atoms
\n
'
%
len
(
nucleotide
))
out
.
write
(
'
%d
ellipsoids
\n
'
%
len
(
nucleotide
))
out
.
write
(
'
%d
bonds
\n
'
%
len
(
topology
))
out
.
write
(
'
\n
'
)
out
.
write
(
'4 atom types
\n
'
)
out
.
write
(
'1 bond types
\n
'
)
out
.
write
(
'
\n
'
)
out
.
write
(
'# System size
\n
'
)
out
.
write
(
'
%f
%f
xlo xhi
\n
'
%
(
lxmin
,
lxmax
))
out
.
write
(
'
%f
%f
ylo yhi
\n
'
%
(
lymin
,
lymax
))
out
.
write
(
'
%f
%f
zlo zhi
\n
'
%
(
lzmin
,
lzmax
))
out
.
write
(
'
\n
'
)
out
.
write
(
'Masses
\n
'
)
out
.
write
(
'
\n
'
)
out
.
write
(
'1 3.1575
\n
'
)
out
.
write
(
'2 3.1575
\n
'
)
out
.
write
(
'3 3.1575
\n
'
)
out
.
write
(
'4 3.1575
\n
'
)
out
.
write
(
'
\n
'
)
out
.
write
(
'# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
\n
'
)
out
.
write
(
'Atoms
\n
'
)
out
.
write
(
'
\n
'
)
for
ib
in
range
(
len
(
nucleotide
)):
out
.
write
(
"
%d
%d
%22.16le
%22.16le
%22.16le
1 1 1
\n
"
%
(
ib
+
1
,
nucleotide
[
ib
][
0
],
nucleotide
[
ib
][
1
][
0
],
nucleotide
[
ib
][
1
][
1
],
nucleotide
[
ib
][
1
][
2
]))
out
.
write
(
'
\n
'
)
out
.
write
(
'# Atom-ID, translational, rotational velocity
\n
'
)
out
.
write
(
'Velocities
\n
'
)
out
.
write
(
'
\n
'
)
for
ib
in
range
(
len
(
nucleotide
)):
out
.
write
(
"
%d
%22.16le
%22.16le
%22.16le
%22.16le
%22.16le
%22.16le
\n
"
%
(
ib
+
1
,
nucleotide
[
ib
][
2
][
0
],
nucleotide
[
ib
][
2
][
1
],
nucleotide
[
ib
][
2
][
2
],
nucleotide
[
ib
][
2
][
3
],
nucleotide
[
ib
][
2
][
4
],
nucleotide
[
ib
][
2
][
5
]))
out
.
write
(
'
\n
'
)
out
.
write
(
'# Atom-ID, shape, quaternion
\n
'
)
out
.
write
(
'Ellipsoids
\n
'
)
out
.
write
(
'
\n
'
)
for
ib
in
range
(
len
(
nucleotide
)):
out
.
write
(
"
%d
%22.16le
%22.16le
%22.16le
%22.16le
%22.16le
%22.16le
%22.16le
\n
"
%
(
ib
+
1
,
nucleotide
[
ib
][
3
][
0
],
nucleotide
[
ib
][
3
][
1
],
nucleotide
[
ib
][
3
][
2
],
nucleotide
[
ib
][
4
][
0
],
nucleotide
[
ib
][
4
][
1
],
nucleotide
[
ib
][
4
][
2
],
nucleotide
[
ib
][
4
][
3
]))
out
.
write
(
'
\n
'
)
out
.
write
(
'# Bond topology
\n
'
)
out
.
write
(
'Bonds
\n
'
)
out
.
write
(
'
\n
'
)
for
ib
in
range
(
len
(
topology
)):
out
.
write
(
"
%d
%d
%d
%d
\n
"
%
(
ib
+
1
,
topology
[
ib
][
0
],
topology
[
ib
][
1
],
topology
[
ib
][
2
]))
out
.
close
()
seqfile
.
close
()
sys
.
exit
(
0
)
Event Timeline
Log In to Comment