Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86582727
lmpsdata.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
Mon, Oct 7, 08:47
Size
81 KB
Mime Type
text/x-python
Expires
Wed, Oct 9, 08:47 (2 d)
Engine
blob
Format
Raw Data
Handle
21444206
Attached To
rLAMMPS lammps
lmpsdata.py
View Options
#
#
# lmpsdata.py
#
# For reading, writing and manipulating lammps data files
# For calculation of certain properities using lammps data files
# For creating VMD input text files using lammps data files
# All x,y,z calculations assume the information includes image flags
class
Lmpsdata
:
def
__init__
(
self
,
file
,
atomtype
):
"""initiates lammps data structures"""
self
.
atomtype
=
atomtype
self
.
keywords
=
[]
self
.
atoms
=
[]
self
.
angles
=
[]
self
.
bonds
=
[]
self
.
dihedrals
=
[]
self
.
dipoles
=
[]
self
.
impropers
=
[]
self
.
masses
=
[]
self
.
shapes
=
[]
self
.
velocities
=
[]
self
.
anglecoef
=
[]
self
.
bondcoef
=
[]
self
.
dihedralcoef
=
[]
self
.
impropercoef
=
[]
self
.
paircoef
=
[]
self
.
read
(
file
)
def
read
(
self
,
file
):
"""Reads in lammps data file
Skips the first line of the file (Comment line)
First reads the header portion of the file (sectflag=1)
blank lines are skipped
header keywords delineate an assignment
if no header keyword is found, body portion begins
Second reads the body portion of the file (sectflag=2)
first line of a section has only a keyword
next line is skipped
remaining lines contain values
a blank line signifies the end of that section
if no value is listed on a line than a keyword must be used
File is read until the end
The keywords read in are stored in self.keywords"""
if
file
==
''
:
print
'no file is given. Will have to build keywords and class structures manually'
return
sectflag
=
1
f
=
open
(
file
,
'r'
)
f
.
readline
()
for
line
in
f
:
row
=
line
.
split
()
if
sectflag
==
1
:
if
len
(
row
)
==
0
:
#skip line the line is blank
pass
elif
len
(
row
)
==
1
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
len
(
row
)
==
2
:
if
row
[
1
]
==
'atoms'
:
self
.
atomnum
=
row
[
0
]
self
.
keywords
.
append
(
'atoms'
)
elif
row
[
1
]
==
'bonds'
:
self
.
bondnum
=
row
[
0
]
self
.
keywords
.
append
(
'bonds'
)
elif
row
[
1
]
==
'angles'
:
self
.
anglenum
=
row
[
0
]
self
.
keywords
.
append
(
'angles'
)
elif
row
[
1
]
==
'dihedrals'
:
self
.
dihedralnum
=
row
[
0
]
self
.
keywords
.
append
(
'dihedrals'
)
elif
row
[
1
]
==
'impropers'
:
self
.
impropernum
=
row
[
0
]
self
.
keywords
.
append
(
'impropers'
)
else
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
len
(
row
)
==
3
:
if
row
[
1
]
==
'atom'
and
row
[
2
]
==
'types'
:
self
.
atomtypenum
=
row
[
0
]
self
.
keywords
.
append
(
'atom types'
)
elif
row
[
1
]
==
'bond'
and
row
[
2
]
==
'types'
:
self
.
bondtypenum
=
row
[
0
]
self
.
keywords
.
append
(
'bond types'
)
elif
row
[
1
]
==
'angle'
and
row
[
2
]
==
'types'
:
self
.
angletypenum
=
row
[
0
]
self
.
keywords
.
append
(
'angle types'
)
elif
row
[
1
]
==
'dihedral'
and
row
[
2
]
==
'types'
:
self
.
dihedraltypenum
=
row
[
0
]
self
.
keywords
.
append
(
'dihedral types'
)
elif
row
[
1
]
==
'improper'
and
row
[
2
]
==
'types'
:
self
.
impropertypenum
=
row
[
0
]
self
.
keywords
.
append
(
'improper types'
)
else
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
len
(
row
)
==
4
:
if
row
[
2
]
==
'xlo'
and
row
[
3
]
==
'xhi'
:
self
.
xdimension
=
[
row
[
0
],
row
[
1
]]
self
.
keywords
.
append
(
'xlo xhi'
)
elif
row
[
2
]
==
'ylo'
and
row
[
3
]
==
'yhi'
:
self
.
ydimension
=
[
row
[
0
],
row
[
1
]]
self
.
keywords
.
append
(
'ylo yhi'
)
elif
row
[
2
]
==
'zlo'
and
row
[
3
]
==
'zhi'
:
self
.
zdimension
=
[
row
[
0
],
row
[
1
]]
self
.
keywords
.
append
(
'zlo zhi'
)
else
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
len
(
row
)
==
5
:
if
row
[
1
]
==
'extra'
and
row
[
2
]
==
'bond'
and
row
[
3
]
==
'per'
and
row
[
4
]
==
'atom'
:
self
.
extrabonds
=
row
[
0
]
self
.
keywords
.
append
(
'extra bond per atom'
)
else
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
len
(
row
)
==
6
:
if
row
[
3
]
==
'xy'
and
row
[
4
]
==
'xz'
and
row
[
5
]
==
'yz'
:
self
.
tilt
=
[
row
[
0
],
row
[
1
],
row
[
2
]]
self
.
keywords
.
append
(
'xy xz yz'
)
else
:
# Set Sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
else
:
# set sectflag to 2, assume line is a keyword
sectflag
=
2
checkkey
=
1
#ensures keyword will be checked in body portion
keyword
=
row
elif
sectflag
==
2
:
if
checkkey
==
1
:
if
len
(
keyword
)
==
1
:
if
keyword
[
0
]
==
'Atoms'
or
keyword
[
0
]
==
'Velocities'
or
keyword
[
0
]
==
'Masses'
or
\
keyword
[
0
]
==
'Shapes'
or
keyword
[
0
]
==
'Dipoles'
or
keyword
[
0
]
==
'Bonds'
or
\
keyword
[
0
]
==
'Angles'
or
keyword
[
0
]
==
'Dihedrals'
or
keyword
[
0
]
==
'Impropers'
:
bodyflag
=
1
blanknum
=
0
self
.
keywords
.
append
(
keyword
[
0
])
checkkey
=
0
else
:
bodyflag
=
0
checkkey
=
0
elif
len
(
keyword
)
==
2
:
if
row
[
1
]
==
'Coeffs'
and
(
row
[
0
]
==
'Pair'
or
row
[
0
]
==
'Bond'
or
row
[
0
]
==
'Angle'
or
\
row
[
0
]
==
'Dihedral'
or
row
[
0
]
==
'Improper'
):
#class 2 force field keywords not included
bodyflag
=
1
blanknum
=
0
self
.
keywords
.
append
(
'{0} {1}'
.
format
(
keyword
[
0
],
keyword
[
1
]))
checkkey
=
0
else
:
bodyflag
=
0
checkkey
=
0
else
:
#egnore line and set bodyflag=0
bodyflag
=
0
checkkey
=
0
if
bodyflag
==
0
:
#bodyflag 0 means no body keyword has been found
if
len
(
row
)
==
1
:
if
row
[
0
]
==
'Atoms'
or
row
[
0
]
==
'Velocities'
or
row
[
0
]
==
'Masses'
or
\
row
[
0
]
==
'Shapes'
or
row
[
0
]
==
'Dipoles'
or
row
[
0
]
==
'Bonds'
or
\
row
[
0
]
==
'Angles'
or
row
[
0
]
==
'Dihedrals'
or
row
[
0
]
==
'Impropers'
:
bodyflag
=
1
blanknum
=
0
keyword
=
row
self
.
keywords
.
append
(
keyword
[
0
])
else
:
#egnore line
pass
elif
len
(
row
)
==
2
:
if
row
[
1
]
==
'Coeffs'
and
(
row
[
0
]
==
'Pair'
or
row
[
0
]
==
'Bond'
or
row
[
0
]
==
'Angle'
or
\
row
[
0
]
==
'Dihedral'
or
row
[
0
]
==
'Improper'
):
#class 2 force field keywords not included
bodyflag
=
1
blanknum
=
0
keyword
=
row
self
.
keywords
.
append
(
'{0} {1}'
.
format
(
keyword
[
0
],
keyword
[
1
]))
else
:
#egnore line
pass
else
:
#egnore line
pass
elif
bodyflag
==
1
:
#currently assumes 1 or more blank lines are between body data keywords
if
len
(
row
)
==
0
:
blanknum
+=
1
if
blanknum
>
1
:
bodyflag
=
0
elif
len
(
keyword
)
==
1
:
if
keyword
[
0
]
==
'Atoms'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
atoms
.
append
(
row
)
elif
keyword
[
0
]
==
'Velocities'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
velocities
.
append
(
row
)
elif
keyword
[
0
]
==
'Masses'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
masses
.
append
(
row
)
elif
keyword
[
0
]
==
'Shapes'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
shapes
.
append
(
row
)
elif
keyword
[
0
]
==
'Dipoles'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
dipoles
.
append
(
row
)
elif
keyword
[
0
]
==
'Bonds'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
bonds
.
append
(
row
)
elif
keyword
[
0
]
==
'Angles'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
angles
.
append
(
row
)
elif
keyword
[
0
]
==
'Dihedrals'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
dihedrals
.
append
(
row
)
elif
keyword
[
0
]
==
'Impropers'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
impropers
.
append
(
row
)
else
:
#egnore line and change bodyflag to 0
bodyflag
=
0
elif
len
(
keyword
)
==
2
:
if
keyword
[
0
]
==
'Pair'
and
keyword
[
1
]
==
'Coeffs'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
paircoef
.
append
(
row
)
elif
keyword
[
0
]
==
'Bond'
and
keyword
[
1
]
==
'Coeffs'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
bondcoef
.
append
(
row
)
elif
keyword
[
0
]
==
'Angle'
and
keyword
[
1
]
==
'Coeffs'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
anglecoef
.
append
(
row
)
elif
keyword
[
0
]
==
'Dihedral'
and
keyword
[
1
]
==
'Coeffs'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
dihedralcoef
.
append
(
row
)
elif
keyword
[
0
]
==
'Improper'
and
keyword
[
1
]
==
'Coeffs'
:
try
:
int
(
row
[
0
])
except
ValueError
:
keyword
=
row
checkkey
=
1
else
:
self
.
impropercoef
.
append
(
row
)
else
:
#egnore line and change bodyflag to 0
bodyflag
=
0
else
:
#egnore line and change bodyflag to 0
bodyflag
=
0
f
.
close
()
def
write
(
self
,
file
,
modflag
):
"""Write lammps data files using the lammps keywords and lammpsdata structures
writes first line of the file as a Comment line
if no modifications to any of the lammpsdata structures (modflag=0)
Use Keywords to write lammpsdata structures directly
if modifications to any of the lammpsdata structures (modflag=1)
Key Lammpsdata structures like atom numbers, coeficient numbers
need to be modified to match the other modified lammpsdata structures
This section will still use the keywords to write lammpsdata structures.
For all modflags, the code will write data to the file until all of the
keyword's data structures have been finished writing.
The keywords are stored in self.keywords"""
f
=
open
(
file
,
'w'
)
f
.
write
(
'polymer data file
\n
'
)
for
row
in
self
.
keywords
:
if
row
==
'atoms'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
atomnum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
atoms
),
row
))
elif
row
==
'bonds'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
bondnum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
bonds
),
row
))
elif
row
==
'angles'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
anglenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
angles
),
row
))
elif
row
==
'dihedrals'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
dihedralnum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
dihedrals
),
row
))
elif
row
==
'impropers'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
impropernum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
impropers
),
row
))
elif
row
==
'atom types'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
atomtypenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
masses
),
row
))
elif
row
==
'bond types'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
bondtypenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
bondcoef
),
row
))
elif
row
==
'angle types'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
angletypenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
anglecoef
),
row
))
elif
row
==
'dihedral types'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
dihedraltypenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
dihedralcoef
),
row
))
elif
row
==
'improper types'
:
if
modflag
==
0
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
impropertypenum
,
row
))
elif
modflag
==
1
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
len
(
self
.
impropercoef
),
row
))
elif
row
==
'xlo xhi'
:
f
.
write
(
'{0} {1} {2}
\n
'
.
format
(
self
.
xdimension
[
0
],
self
.
xdimension
[
1
],
row
))
elif
row
==
'ylo yhi'
:
f
.
write
(
'{0} {1} {2}
\n
'
.
format
(
self
.
ydimension
[
0
],
self
.
ydimension
[
1
],
row
))
elif
row
==
'zlo zhi'
:
f
.
write
(
'{0} {1} {2}
\n
'
.
format
(
self
.
zdimension
[
0
],
self
.
zdimension
[
1
],
row
))
elif
row
==
'extra bond per atom'
:
f
.
write
(
'{0} {1}
\n
'
.
format
(
self
.
extrabonds
,
row
))
elif
row
==
'xy xz yz'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
self
.
tilt
[
0
],
self
.
tilt
[
1
],
self
.
tilt
[
2
],
row
))
elif
row
==
'Atoms'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
atoms
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Velocities'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
velocities
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Masses'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
masses
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Shapes'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
shapes
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Dipoles'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
dipoles
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Bonds'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
bonds
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Angles'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
angles
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Dihedrals'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
dihedrals
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
elif
row
==
'Impropers'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
impropers
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Pair Coeffs'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
paircoef
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Bond Coeffs'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
bondcoef
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Angle Coeffs'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
anglecoef
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Dihedral Coeffs'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
dihedralcoef
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
elif
row
==
'Improper Coeffs'
:
f
.
write
(
'
\n
{0}'
.
format
(
row
))
#new line between header and body portion or two body keywords
f
.
write
(
'
\n
'
)
#new line between body keyword and body data
for
line
in
self
.
impropercoef
:
f
.
write
(
'
\n
'
)
#creates a new line for adding body data
for
item
in
line
:
f
.
write
(
' {0}'
.
format
(
item
))
#adds in each peice of body data with space imbetween
f
.
write
(
'
\n
'
)
#allows space to be added between the end of body data and a new keyword
else
:
pass
f
.
close
()
def
atomorder
(
self
):
"""Takes self.atoms and organizes the atom id from least to greatest.
If the atom ids are allready ordered this algorithm will do nothing."""
current
=
range
(
len
(
self
.
atoms
[
0
]))
# initialize current [assumes self.atoms coloumn #'s does not change]
for
i
in
range
(
1
,
len
(
self
.
atoms
)):
#when i=0, self.atoms will not change; therefore its skipped
for
k
in
range
(
len
(
self
.
atoms
[
i
])):
current
[
k
]
=
self
.
atoms
[
i
][
k
]
for
j
in
range
(
i
-
1
,
-
1
,
-
1
):
for
k
in
range
(
len
(
self
.
atoms
[
j
])):
self
.
atoms
[
j
+
1
][
k
]
=
self
.
atoms
[
j
][
k
]
if
int
(
current
[
0
])
>
int
(
self
.
atoms
[
j
][
0
]):
for
k
in
range
(
len
(
current
)):
self
.
atoms
[
j
+
1
][
k
]
=
current
[
k
]
break
elif
j
==
0
:
for
k
in
range
(
len
(
current
)):
self
.
atoms
[
j
][
k
]
=
current
[
k
]
#dont need a break here because this is the last j value in the for loop
def
addatoms
(
self
,
atoms
,
retlist
=
False
):
"""Appends atoms to self.atoms. Assumes atoms are written in the correct atomtype format.
Change added atom numbers in self.atoms so self.atoms will be in increasing sequential order.
If retlist is True return a list of the modified atom numbers otherwise exit."""
initpos
=
len
(
self
.
atoms
)
#Store index of first appended atoms
for
item
in
atoms
:
self
.
atoms
.
append
(
range
(
len
(
item
)))
#initializing spots for the new atoms to go
numberchange
=
[]
#initiate number change where the changed atom numbers will be stored
count
=
0
for
i
in
range
(
initpos
,
len
(
self
.
atoms
)):
for
j
in
range
(
len
(
self
.
atoms
[
i
])):
if
j
==
0
:
self
.
atoms
[
i
][
0
]
=
str
(
i
+
1
)
numberchange
.
append
(
str
(
i
+
1
))
else
:
self
.
atoms
[
i
][
j
]
=
atoms
[
count
][
j
]
count
+=
1
if
retlist
:
return
numberchange
else
:
return
#need to redo this algorithm so their is no referencing going on.
def
adddata
(
self
,
data
,
keyword
):
"""Adds data to a keyword's existing data structure
All body keywords are viable except Atoms which has its own method
All header keywords will do nothing in this algrorithm because they have their own algorithm
changes the body id number so id numbers will be in sequential order"""
if
keyword
==
'Velocities'
:
pos
=
len
(
self
.
velocities
)
for
item
in
data
:
self
.
velocities
.
append
(
item
)
#adds data to existing data structure
self
.
velocities
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Masses'
:
pos
=
len
(
self
.
masses
)
for
item
in
data
:
self
.
masses
.
append
(
item
)
#adds data to existing data structure
self
.
masses
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Shapes'
:
pos
=
len
(
self
.
shapes
)
for
item
in
data
:
self
.
shapes
.
append
(
item
)
#adds data to existing data structure
self
.
shapes
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Dipoles'
:
pos
=
len
(
self
.
dipoles
)
for
item
in
data
:
self
.
dipoles
.
append
(
item
)
#adds data to existing data structure
self
.
dipoles
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Bonds'
:
pos
=
len
(
self
.
bonds
)
for
item
in
data
:
self
.
bonds
.
append
(
item
)
#adds data to existing data structure
self
.
bonds
[
pos
][
0
]
=
str
(
pos
+
1
)
pos
+=
1
elif
keyword
==
'Angles'
:
pos
=
len
(
self
.
angles
)
for
item
in
data
:
self
.
angles
.
append
(
item
)
#adds data to existing data structure
self
.
angles
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Dihedrals'
:
pos
=
len
(
self
.
dihedrals
)
for
item
in
data
:
self
.
dihedrals
.
append
(
item
)
#adds data to existing data structure
self
.
dihedrals
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Impropers'
:
pos
=
len
(
self
.
impropers
)
for
item
in
data
:
self
.
impropers
.
append
(
item
)
#adds data to existing data structure
self
.
impropers
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Pair Coeffs'
:
pos
=
len
(
self
.
paircoef
)
for
item
in
data
:
self
.
paircoef
.
append
(
item
)
#adds data to existing data structure
self
.
paircoef
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Bond Coeffs'
:
pos
=
len
(
self
.
bondcoef
)
for
item
in
data
:
self
.
bondcoef
.
append
(
item
)
#adds data to existing data structure
self
.
bondcoef
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Angle Coeffs'
:
pos
=
len
(
self
.
anglecoef
)
for
item
in
data
:
self
.
anglecoef
.
append
(
item
)
#adds data to existing data structure
self
.
anglecoef
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Dihedral Coeffs'
:
pos
=
len
(
self
.
dihedralcoef
)
for
item
in
data
:
self
.
dihedralcoef
.
append
(
item
)
#adds data to existing data structure
self
.
dihedralcoef
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
elif
keyword
==
'Improper Coeffs'
:
pos
=
len
(
self
.
impropercoef
)
for
item
in
data
:
self
.
impropercoef
.
append
(
item
)
#adds data to existing data structure
self
.
impropercoef
[
pos
][
0
]
=
str
(
pos
+
1
)
#changing body id number
pos
+=
1
# def modifiydata(data,keyword): Will not implement
# """for modifying header data"""
def
changeatomnum
(
self
,
data
,
originalatoms
,
atomchanges
,
vflag
=
False
):
"""Takes data which contains atom ids
and alters those ids from the original atom id system in originalatoms
to the new atom id system in atomchanges."""
#set up boolean array to match the size of data and with all True values
# print atomchanges
array
=
booleanarray
(
len
(
data
),
len
(
data
[
0
]),
True
)
# print originalatoms
if
vflag
:
# for changing atomnumbers for velocities
for
i
in
range
(
len
(
originalatoms
)):
#len of originalatoms should match len of atomchanges
if
originalatoms
[
i
][
0
]
==
atomchanges
[
i
]:
continue
for
j
in
range
(
len
(
data
)):
for
k
in
range
(
1
):
if
data
[
j
][
k
]
==
originalatoms
[
i
][
0
]:
#checks if the atom id in data
#matches the atom id in origianalatoms
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
#set the boolean array to False and
#change the atom id in data to atomchanges
array
.
setelement
(
j
,
k
,
False
)
data
[
j
][
k
]
=
atomchanges
[
i
]
break
#the change of boolean array to False insures
#the atom id in data will only be changed once.
else
:
# for changing atom numbers for everything else
for
i
in
range
(
len
(
originalatoms
)):
#len of originalatoms should match len of atomchanges
if
originalatoms
[
i
][
0
]
==
atomchanges
[
i
]:
continue
for
j
in
range
(
len
(
data
)):
for
k
in
range
(
2
,
len
(
data
[
j
])):
if
data
[
j
][
k
]
==
originalatoms
[
i
][
0
]:
#checks if the atom id in data
#matches the atom id in origianalatoms
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
#set the boolean array to False and
#change the atom id in data to atomchanges
array
.
setelement
(
j
,
k
,
False
)
data
[
j
][
k
]
=
atomchanges
[
i
]
break
#the change of boolean array to False insures
#the atom id in data will only be changed once.
# print data
return
data
def
deletebodydata
(
self
,
keyword
):
"""Changes a keyword's class data structure to an empty structure"""
if
keyword
==
'Velocities'
:
self
.
velocities
=
[]
elif
keyword
==
'Masses'
:
self
.
masses
=
[]
elif
keyword
==
'Shapes'
:
self
.
shapes
=
[]
elif
keyword
==
'Dipoles'
:
self
.
dipoles
=
[]
elif
keyword
==
'Bonds'
:
self
.
bonds
=
[]
elif
keyword
==
'Angles'
:
self
.
angles
=
[]
elif
keyword
==
'Dihedrals'
:
self
.
dihedrals
=
[]
elif
keyword
==
'Impropers'
:
self
.
impropers
=
[]
elif
keyword
==
'Pair Coeffs'
:
self
.
paircoef
=
[]
elif
keyword
==
'Bond Coeffs'
:
self
.
bondcoef
=
[]
elif
keyword
==
'Angle Coeffs'
:
self
.
anglecoef
=
[]
elif
keyword
==
'Dihedral Coeffs'
:
self
.
dihedralcoef
=
[]
elif
keyword
==
'Improper Coeffs'
:
self
.
impropercoef
=
[]
elif
keyword
==
'Atoms'
:
self
.
atoms
=
[]
def
extractmolecules
(
self
,
molecule
):
"""Takes the variable molecule and
extracts the individual molecules' data back into lmpsdata.
This extraction takes place through a 4 step process.
Step 1: Use a molecule's keywords to alter the lmpsdata data structures to empty.
To acomplish this procedure use the lmpsdata method deletebodydata.
Step 2: Add the molecules' atoms to lmpsdata's atoms using the method addatoms.
Return a list of atom id changes for each molecule.
Step 3: Utilize each molecules list of atom id changes to change their data's atom id numbers.
Uses changeatomnum and returns the altered molecule's data.
Step 4: Add the altered molecules' data to lmpsdata's data using the method adddata"""
#Use molecule index 0's keyword to change the equivalent lmpsdata structures to empty
print
'extracting the molecules back to data'
for
keyword
in
molecule
[
0
]
.
keywords
:
# step 1
self
.
deletebodydata
(
keyword
)
atomchanges
=
[]
#initializing step 2
for
i
in
range
(
len
(
molecule
)):
# step 2
atomchanges
.
append
(
self
.
addatoms
(
molecule
[
i
]
.
atoms
,
True
))
for
i
in
range
(
len
(
molecule
)):
# step 3
for
keyword
in
molecule
[
i
]
.
keywords
:
if
keyword
==
'Angles'
:
molecule
[
i
]
.
angles
=
self
.
changeatomnum
(
molecule
[
i
]
.
angles
,
molecule
[
i
]
.
atoms
,
atomchanges
[
i
])
if
keyword
==
'Bonds'
:
molecule
[
i
]
.
bonds
=
self
.
changeatomnum
(
molecule
[
i
]
.
bonds
,
molecule
[
i
]
.
atoms
,
atomchanges
[
i
])
elif
keyword
==
'Dihedrals'
:
molecule
[
i
]
.
dihedrals
=
self
.
changeatomnum
(
molecule
[
i
]
.
dihedrals
,
molecule
[
i
]
.
atoms
,
atomchanges
[
i
])
elif
keyword
==
'Velocities'
:
molecule
[
i
]
.
velocities
=
self
.
changeatomnum
(
molecule
[
i
]
.
velocities
,
molecule
[
i
]
.
atoms
,
atomchanges
[
i
],
True
)
elif
keyword
==
'Impropers'
:
molecule
[
i
]
.
impropers
=
self
.
changeatomnum
(
molecule
[
i
]
.
impropers
,
molecule
[
i
]
.
atoms
,
atomchanges
[
i
])
for
i
in
range
(
len
(
molecule
)):
# step 4
for
keyword
in
molecule
[
i
]
.
keywords
:
if
keyword
==
'Angles'
:
self
.
adddata
(
molecule
[
i
]
.
angles
,
keyword
)
elif
keyword
==
'Bonds'
:
self
.
adddata
(
molecule
[
i
]
.
bonds
,
keyword
)
elif
keyword
==
'Dihedrals'
:
self
.
adddata
(
molecule
[
i
]
.
dihedrals
,
keyword
)
elif
keyword
==
'Velocities'
:
self
.
adddata
(
molecule
[
i
]
.
velocities
,
keyword
)
elif
keyword
==
'Impropers'
:
self
.
adddata
(
molecule
[
i
]
.
impropers
,
keyword
)
def
density
(
self
,
ringsize
,
init
,
final
,
file
=
' '
):
"""Create spherical shells from the initial radius to the final radius
The spherical shell's thickness is defined by ringsize
Calculate the mass of particles within each spherical shell
Calculate the volume of each spherical shell
Then calculate the density in each spherical shell
Current code is written with the assumption the spheres are centered around the origin
If a file is listed place the results in the variable file otherwise return the results
The distance calculations are written assuming image flags are in the data file"""
rho
=
[]
r
=
[
init
]
#initial radius to examine
radius
=
init
+
ringsize
while
radius
<
final
:
#finding the rest of the radii to examine
r
.
append
(
radius
)
radius
+=
ringsize
for
radius
in
r
:
volume
=
4.0
/
3.0
*
3.1416
*
((
radius
+
ringsize
)
**
3
-
radius
**
3
)
mass
=
0.0
for
row
in
self
.
atoms
:
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'bond'
or
\
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'full'
or
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'molecular'
or
\
self
.
atomtype
==
'peri'
:
l
=
len
(
row
)
-
1
dist
=
float
(
row
[
l
-
3
])
**
2
+
float
(
row
[
l
-
4
])
**
2
+
float
(
row
[
l
-
5
])
**
2
if
dist
<
(
radius
+
ringsize
)
**
2
and
dist
>=
radius
**
2
:
if
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'charge'
or
\
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'peri'
:
type
=
int
(
row
[
1
])
-
1
else
:
type
=
int
(
row
[
2
])
-
1
mass
+=
float
(
self
.
masses
[
type
][
1
])
#assumes self.masses is written in atomtype order
elif
self
.
atomtype
==
'dipole'
:
l
=
len
(
row
)
-
1
dist
=
float
(
row
[
l
-
6
])
**
2
+
float
(
row
[
l
-
7
])
**
2
+
float
(
row
[
l
-
8
])
**
2
if
dist
<
(
radius
+
ringsize
)
**
2
and
dist
>=
radius
**
2
:
type
=
int
(
row
[
1
])
-
1
mass
+=
float
(
self
.
masses
[
type
][
1
])
#assumes self.masses is written in atomtype order
elif
self
.
atomtype
==
'ellipsoid'
:
l
=
len
(
row
)
-
1
dist
=
float
(
row
[
l
-
7
])
**
2
+
float
(
row
[
l
-
8
])
**
2
+
float
(
row
[
l
-
9
])
**
2
if
dist
<
(
radius
+
ringsize
)
**
2
and
dist
>=
radius
**
2
:
type
=
int
(
row
[
1
])
mass
+=
float
(
self
.
masses
[
type
][
1
])
#assumes self.masses is written in atomtype order
elif
self
.
atomtype
==
'hybrid'
:
dist
=
float
(
row
[
2
])
**
2
+
float
(
row
[
3
])
**
2
+
float
(
row
[
4
])
**
2
if
dist
<
(
radius
+
ringsize
)
**
2
and
dist
>=
radius
**
2
:
type
=
int
(
row
[
1
])
mass
+=
float
(
self
.
masses
[
type
][
1
])
#assumes self.masses is written in atomtype order
rho
.
append
(
mass
/
volume
)
if
file
==
' '
:
return
r
,
rho
else
:
f
=
open
(
file
,
'w'
)
for
i
in
range
(
len
(
r
)):
f
.
write
(
'{0} {1}
\n
'
.
format
(
r
[
i
],
rho
[
i
]))
f
.
close
()
return
def
createxyz
(
self
,
file
,
routine
=
'mass'
,
values
=
None
):
"""Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f
=
open
(
file
,
'w'
)
f
.
write
(
'{0}
\n
'
.
format
(
len
(
self
.
atoms
)))
f
.
write
(
'atoms
\n
'
)
if
routine
==
'mass'
:
for
line
in
self
.
atoms
:
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
type
=
int
(
line
[
2
])
-
1
#3rd position
else
:
type
=
int
(
line
[
1
])
-
1
#2nd position
mass
=
self
.
masses
[
type
][
1
]
if
mass
==
'12.0107'
:
elementnum
=
6
elif
mass
==
'15.9994'
:
elementnum
=
8
elif
mass
==
'26.9815'
:
elementnum
=
13
else
:
print
'no matching mass value. The method will exit'
return
l
=
len
(
line
)
-
1
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'bond'
or
\
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'full'
or
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'molecular'
or
\
self
.
atomtype
==
'peri'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
elif
self
.
atomtype
==
'dipole'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
8
],
line
[
l
-
7
],
line
[
l
-
6
]))
elif
self
.
atomtype
==
'ellipsoid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
9
],
line
[
l
-
8
],
line
[
l
-
7
]))
elif
self
.
atomtype
==
'hybrid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
2
],
line
[
3
],
line
[
4
]))
f
.
close
()
elif
routine
==
'atomtype'
:
for
line
in
self
.
atoms
:
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
type
=
int
(
line
[
2
])
-
1
#3rd position
else
:
type
=
int
(
line
[
1
])
-
1
#2nd position
elementnum
=
values
[
type
]
l
=
len
(
line
)
-
1
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'bond'
or
\
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'full'
or
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'molecular'
or
\
self
.
atomtype
==
'peri'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
elif
self
.
atomtype
==
'dipole'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
8
],
line
[
l
-
7
],
line
[
l
-
6
]))
elif
self
.
atomtype
==
'ellipsoid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
9
],
line
[
l
-
8
],
line
[
l
-
7
]))
elif
self
.
atomtype
==
'hybrid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
2
],
line
[
3
],
line
[
4
]))
f
.
close
()
class
booleanarray
:
"""A class that stores boolean values in a list of lists."""
def
__init__
(
self
,
rownum
,
colnum
,
initval
):
""" initialise a list of lists (array) with
rownum correspondinig to the number of lists in the list and
colnum corresponding to the number of elements in the list's list.
initval is the value the list of lists will be initialized with.
initval should be a boolean value."""
# initializing self.array
self
.
array
=
[]
for
i
in
range
(
rownum
):
self
.
array
.
append
(
range
(
colnum
))
# setting elements of self.array to initval
for
i
in
range
(
rownum
):
for
j
in
range
(
colnum
):
self
.
setelement
(
i
,
j
,
initval
)
def
setelement
(
self
,
rownum
,
colnum
,
value
):
"""Assigns value to the list of lists (array) element at rownum and colnum."""
self
.
array
[
rownum
][
colnum
]
=
value
def
getelement
(
self
,
rownum
,
colnum
):
"""Returns element in the list of lists (array) at rownum and colnum."""
return
self
.
array
[
rownum
][
colnum
]
def
atomdistance
(
a
,
b
,
atomtype
):
"""Returns the distance between atom a and b.
Atomtype is the style the atom data is written in.
All atom data is assumed to have image flags in the data.
Atom a and atom b are assumed to be the same atomtype."""
from
math
import
sqrt
if
atomtype
==
'angle'
or
atomtype
==
'atomic'
or
atomtype
==
'bond'
or
\
atomtype
==
'charge'
or
atomtype
==
'colloid'
or
atomtype
==
'electron'
or
\
atomtype
==
'full'
or
atomtype
==
'granular'
or
atomtype
==
'molecular'
or
\
atomtype
==
'peri'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
3
])
-
float
(
b
[
l
-
3
]))
**
2
\
+
(
float
(
a
[
l
-
4
])
-
float
(
b
[
l
-
4
]))
**
2
\
+
(
float
(
a
[
l
-
5
])
-
float
(
b
[
l
-
5
]))
**
2
)
elif
atomtype
==
'dipole'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
6
])
-
float
(
b
[
l
-
6
]))
**
2
\
+
(
float
(
a
[
l
-
7
])
-
float
(
b
[
l
-
7
]))
**
2
\
+
(
float
(
a
[
l
-
8
])
-
float
(
b
[
l
-
8
]))
**
2
)
elif
atomtype
==
'ellipsoid'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
7
])
-
float
(
b
[
l
-
7
]))
**
2
\
+
(
float
(
a
[
l
-
8
])
-
float
(
b
[
l
-
8
]))
**
2
\
+
(
float
(
a
[
l
-
9
])
-
float
(
b
[
l
-
9
]))
**
2
)
elif
atomtype
==
'hybrid'
:
dist
=
sqrt
((
float
(
a
[
2
])
-
float
(
b
[
2
]))
**
2
\
+
(
float
(
a
[
3
])
-
float
(
b
[
3
]))
**
2
\
+
(
float
(
a
[
4
])
-
float
(
b
[
4
]))
**
2
)
return
dist
def
distance
(
a
,
coord
,
atomtype
):
"""Returns the distance between atom a and coord.
Atomtype is the style the atom data is written in.
All atom data is assumed to have image flags in the data."""
from
math
import
sqrt
#need to alter this slightly
if
atomtype
==
'angle'
or
atomtype
==
'atomic'
or
atomtype
==
'bond'
or
\
atomtype
==
'charge'
or
atomtype
==
'colloid'
or
atomtype
==
'electron'
or
\
atomtype
==
'full'
or
atomtype
==
'granular'
or
atomtype
==
'molecular'
or
\
atomtype
==
'peri'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
3
])
-
coord
[
2
])
**
2
\
+
(
float
(
a
[
l
-
4
])
-
coord
[
1
])
**
2
\
+
(
float
(
a
[
l
-
5
])
-
coord
[
0
])
**
2
)
elif
atomtype
==
'dipole'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
6
])
-
coord
[
2
])
**
2
\
+
(
float
(
a
[
l
-
7
])
-
coord
[
1
])
**
2
\
+
(
float
(
a
[
l
-
8
])
-
coord
[
0
])
**
2
)
elif
atomtype
==
'ellipsoid'
:
l
=
len
(
a
)
-
1
dist
=
sqrt
((
float
(
a
[
l
-
7
])
-
coord
[
2
])
**
2
\
+
(
float
(
a
[
l
-
8
])
-
coord
[
1
])
**
2
\
+
(
float
(
a
[
l
-
9
])
-
coord
[
0
])
**
2
)
elif
atomtype
==
'hybrid'
:
dist
=
sqrt
((
float
(
a
[
2
])
-
coord
[
0
])
**
2
\
+
(
float
(
a
[
3
])
-
coord
[
1
])
**
2
\
+
(
float
(
a
[
4
])
-
coord
[
0
])
**
2
)
return
dist
class
particlesurface
:
def
__init__
(
self
,
particle
,
cutoff
,
atomid
,
atomtype
,
shape
=
'sphere'
):
"""Builds a particle surface with a specific shape from a particle
The atoms choosen from the surface will have the specific atomid
atomid will be given in terms of an integer and not a string."""
self
.
particle
=
particle
.
atoms
self
.
atomtype
=
atomtype
self
.
cutoff
=
cutoff
self
.
surface
=
[]
self
.
particlelocator
=
[]
#stores surface particle locations with respect to the particle
self
.
deletedatoms
=
[]
if
shape
==
'sphere'
:
self
.
createspheresurf
(
atomid
)
def
createspheresurf
(
self
,
atomid
):
"""Assumes sphere is centered around 0,0,0.
Finds maximum distance particle with atomid.
Than all atoms within cutoff distance of max distance
will be included into self.surface
Also will build a list of where the surface particles are
in relationship to the particle.
Will create a booleanarray to keep track of any changes made to the surface."""
particledist
=
[]
for
row
in
self
.
particle
:
#Stores all of the distances of the particle in particledist
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
type
=
int
(
row
[
2
])
#3rd position
else
:
type
=
int
(
row
[
1
])
#2nd position
if
type
==
atomid
:
#only calculates the distance if the atom has the correct atomid number
particledist
.
append
(
distance
(
row
,[
0
,
0
,
0
],
self
.
atomtype
))
else
:
particledist
.
append
(
0.0
)
self
.
maxdist
=
max
(
particledist
)
# finds the max dist.
for
i
in
range
(
len
(
particledist
)):
if
particledist
[
i
]
>=
(
self
.
maxdist
-
self
.
cutoff
):
self
.
particlelocator
.
append
(
i
)
self
.
surface
.
append
(
self
.
particle
[
i
])
self
.
bool
=
booleanarray
(
len
(
self
.
surface
),
1
,
True
)
def
getsurfatom
(
self
,
rownum
):
return
self
.
surface
[
rownum
]
def
getbool
(
self
,
rownum
):
return
self
.
bool
.
getelement
(
rownum
,
0
)
def
setbool
(
self
,
rownum
,
value
):
self
.
bool
.
setelement
(
rownum
,
0
,
value
)
def
removesurfatom
(
self
,
rownum
):
"""note this does not actually remove the surface atom.
Instead this adds a value to the list deletedatoms.
This list is later used in extractparticle to actually delete those atoms"""
self
.
deletedatoms
.
append
(
self
.
particlelocator
[
rownum
])
def
extractparticle
(
self
):
""""deletesatoms from particle from the list of deletedatoms
and than returns the particle."""
self
.
particle
=
self
.
deleteatoms
(
self
.
particle
,
self
.
deletedatoms
)
return
self
.
particle
def
deleteatoms
(
self
,
structure
,
rows
):
"""delete atoms from particle and shifts the structure down"""
new
=
[]
#mulitple copying of b to the rows being replaced.
if
rows
==
[]:
for
line
in
structure
:
new
.
append
(
line
)
#if no rows need replacing copy structure
return
new
for
i
in
range
(
rows
[
0
]):
new
.
append
(
structure
[
i
])
#copy structure to new until the first replaced row
count
=
0
for
i
in
range
(
rows
[
0
],
len
(
structure
)
-
len
(
rows
)):
# to replace rows and shift undeleted rows over
for
j
in
range
(
i
+
1
+
count
,
len
(
structure
)):
for
val
in
rows
:
if
val
==
j
:
count
+=
1
break
if
val
==
j
:
continue
else
:
new
.
append
(
structure
[
j
])
break
return
new
def
addatom
(
self
,
atomtype
,
charge
=
None
,
moleculenum
=
None
):
"""Adds an atom to the particle surface between maxdist and the cutoff distance.
This atom is stored in self.particle and self.surface.
In the current coding the particle is centered at (0,0,0).
This method has a required input of atomtype and optional inputs of charge and moleculenum.
This method currently does not support correctly self.atomtype of dipole, electron, ellipsoid, granular, peri or hybrid
Will use the image flags 0,0,0.
Note: The atomtype here is different from the atomtype attribute as part of this class
Note: If the added atom will be required for bonding later than you will need to extract the particle data
and rebuild the surface, because this method doesn't include updates to the required variables due to some coding issues"""
import
math
,
random
# Checks to make sure self.atomtype is not set to an unsupported value
if
self
.
atomtype
==
'dipole'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'ellipsoid'
or
self
.
atomtype
==
'peri'
or
self
.
atomtype
==
'hybrid'
:
print
self
.
atomtype
,
'is not currently supported. But, you can add support by adding the needed functionality and inputs'
return
print
'adding an atom to the surface'
# Randomly assigns the position of a new atom in a spherical shell between self.cutoff and self.maxdist
foundpoint
=
False
while
(
foundpoint
==
False
):
x
=
random
.
uniform
(
-
self
.
maxdist
,
self
.
maxdist
)
y
=
random
.
uniform
(
-
self
.
maxdist
,
self
.
maxdist
)
try
:
z
=
random
.
uniform
(
math
.
sqrt
(
self
.
cutoff
**
2
-
x
**
2
-
y
**
2
),
math
.
sqrt
(
self
.
maxdist
**
2
-
x
**
2
-
y
**
2
))
except
ValueError
:
continue
distance
=
math
.
sqrt
(
x
**
2
+
y
**
2
+
z
**
2
)
if
distance
<
self
.
maxdist
and
distance
>=
self
.
maxdist
-
self
.
cutoff
:
foundpoint
=
True
# initialize new row in self.particle and self.surface
self
.
particle
.
append
([])
self
.
surface
.
append
([])
pposition
=
len
(
self
.
particle
)
-
1
#particle postion
sposition
=
len
(
self
.
surface
)
-
1
#surface position
# Adds the atom id number to the new row
self
.
particle
[
pposition
]
.
append
(
str
(
pposition
+
1
))
self
.
surface
[
sposition
]
.
append
(
str
(
pposition
+
1
))
# Adds the atomtype and the moleculenum if needed to the new row
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
self
.
particle
[
pposition
]
.
append
(
str
(
moleculenum
))
self
.
surface
[
sposition
]
.
append
(
str
(
moleculenum
))
self
.
particle
[
pposition
]
.
append
(
str
(
atomtype
))
self
.
surface
[
sposition
]
.
append
(
str
(
atomtype
))
else
:
self
.
particle
[
pposition
]
.
append
(
str
(
atomtype
))
self
.
surface
[
sposition
]
.
append
(
str
(
atomtype
))
# Adds the charge if needed to the new row
if
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'full'
:
self
.
particle
[
pposition
]
.
append
(
str
(
charge
))
self
.
surface
[
sposition
]
.
append
(
str
(
charge
))
# Adds the atom's position to the new row
self
.
particle
[
pposition
]
.
append
(
str
(
x
))
self
.
surface
[
sposition
]
.
append
(
str
(
x
))
self
.
particle
[
pposition
]
.
append
(
str
(
y
))
self
.
surface
[
sposition
]
.
append
(
str
(
y
))
self
.
particle
[
pposition
]
.
append
(
str
(
z
))
self
.
surface
[
sposition
]
.
append
(
str
(
z
))
# Adds the atom's image flags to the new row
self
.
particle
[
pposition
]
.
append
(
'0'
)
self
.
surface
[
sposition
]
.
append
(
'0'
)
self
.
particle
[
pposition
]
.
append
(
'0'
)
self
.
surface
[
sposition
]
.
append
(
'0'
)
self
.
particle
[
pposition
]
.
append
(
'0'
)
self
.
surface
[
sposition
]
.
append
(
'0'
)
def
createxyz
(
self
,
file
,
data
,
routine
=
'mass'
,
values
=
None
):
"""This shows the particle surface. To show the particle surface after bonding has occured,
you will need to extract the particle than reinsert the particle into the class and use createxyz.
Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f
=
open
(
file
,
'w'
)
f
.
write
(
'{0}
\n
'
.
format
(
len
(
self
.
surface
)))
f
.
write
(
'atoms
\n
'
)
if
routine
==
'mass'
:
for
line
in
self
.
surface
:
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
type
=
int
(
line
[
2
])
-
1
#3rd position
else
:
type
=
int
(
line
[
1
])
-
1
#2nd position
mass
=
data
.
masses
[
type
][
1
]
if
mass
==
'12.0107'
:
elementnum
=
6
elif
mass
==
'15.9994'
:
elementnum
=
8
elif
mass
==
'26.981539'
:
elementnum
=
13
else
:
print
'no matching mass value. The method will exit'
return
l
=
len
(
line
)
-
1
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'bond'
or
\
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'full'
or
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'molecular'
or
\
self
.
atomtype
==
'peri'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
elif
self
.
atomtype
==
'dipole'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
8
],
line
[
l
-
7
],
line
[
l
-
6
]))
elif
self
.
atomtype
==
'ellipsoid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
9
],
line
[
l
-
8
],
line
[
l
-
7
]))
elif
self
.
atomtype
==
'hybrid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
2
],
line
[
3
],
line
[
4
]))
f
.
close
()
elif
routine
==
'atomtype'
:
for
line
in
self
.
surface
:
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'bond'
or
self
.
atomtype
==
'full'
or
\
self
.
atomtype
==
'molecular'
:
type
=
int
(
line
[
2
])
-
1
#3rd position
else
:
type
=
int
(
line
[
1
])
-
1
#2nd position
elementnum
=
values
[
type
]
l
=
len
(
line
)
-
1
if
self
.
atomtype
==
'angle'
or
self
.
atomtype
==
'atomic'
or
self
.
atomtype
==
'bond'
or
\
self
.
atomtype
==
'charge'
or
self
.
atomtype
==
'colloid'
or
self
.
atomtype
==
'electron'
or
\
self
.
atomtype
==
'full'
or
self
.
atomtype
==
'granular'
or
self
.
atomtype
==
'molecular'
or
\
self
.
atomtype
==
'peri'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
elif
self
.
atomtype
==
'dipole'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
8
],
line
[
l
-
7
],
line
[
l
-
6
]))
elif
self
.
atomtype
==
'ellipsoid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
9
],
line
[
l
-
8
],
line
[
l
-
7
]))
elif
self
.
atomtype
==
'hybrid'
:
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
2
],
line
[
3
],
line
[
4
]))
f
.
close
()
def
molecules
(
data
,
init
,
final
,
processors
,
method
=
'all'
):
""" There are two ways to initialize the class Lmpsmolecule.
Method controls which way is chosen.
Acceptable values for Method are 'all', 'atom'.
The default method is 'all'"""
molecule
=
[]
from
multiprocessing
import
Pool
p
=
Pool
(
processors
)
for
i
in
range
(
init
,
final
+
1
):
molecule
.
append
(
p
.
apply_async
(
Lmpsmolecule
,(
i
,
data
,
method
,)))
for
i
in
range
(
len
(
molecule
)):
molecule
[
i
]
=
molecule
[
i
]
.
get
()
p
.
close
()
p
.
join
()
return
molecule
class
Lmpsmolecule
:
#Technically should be a meta class but written as a seperate class for easier coding.
def
__init__
(
self
,
moleculenum
,
data
,
method
):
"""initiates lammps molecule structures
and than extract the appropriate molecular structures from the base class data"""
# ***** Lammps Molecule Structures
self
.
keywords
=
[
'Atoms'
]
#include atoms here since this will be in every instance of this class
self
.
atoms
=
[]
#extract
self
.
angles
=
[]
#extract if keyword is there
self
.
bonds
=
[]
#extract if keyword is there
self
.
dihedrals
=
[]
#extract if keyword is there
self
.
impropers
=
[]
#extract if keyword is there
self
.
velocities
=
[]
#extract if keyword is there
# ***** Molecule number
self
.
moleculenum
=
moleculenum
# ***** Extracting the moleculue's atom information from data
self
.
extract
(
data
,
method
)
# ***** If methods is 'all' then self.extract will extract all of the molecule's information from data
def
extract
(
self
,
data
,
method
):
"""uses base class data to extract the molecules atoms
and any other molecule data from molecule moleculenum. This extraction is done in a two step process.
Step 1: extract data.atoms with moleculenum and renumber self.atoms beginning from 1
store extracted atom numbers from data.atoms in a list called temp
Step 2: pass temp and data to method changeatomnum
which extracts the rest of the Lammps Molecule structures from data
and changes the atomnumbers in those structures to match the ones already in self.atoms"""
#checking to make sure atomtype is a valid molecule
#if not a valid molecule print error and exit method
if
data
.
atomtype
==
'full'
or
data
.
atomtype
==
'molecular'
:
# ***** Sets the molecule's atomtype
self
.
atomtype
=
data
.
atomtype
else
:
print
"not a valid molecular structure"
return
#extract the molecule self.moleculenum from data.atom to self.atom
atomnum
=
1
temp
=
[]
for
i
in
range
(
len
(
data
.
atoms
)):
if
int
(
data
.
atoms
[
i
][
1
])
==
self
.
moleculenum
:
temp
.
append
(
data
.
atoms
[
i
][
0
])
#store extracted atomnumbers into temp
self
.
atoms
.
append
(
data
.
atoms
[
i
])
#store extracted atom into self.atom
self
.
atoms
[
atomnum
-
1
][
0
]
=
str
(
atomnum
)
#change extracted atom to the correct atomnumber
atomnum
+=
1
#update atomnumber
#extract the rest of the Lammps Molecule structures from data using changeatomnum
if
method
==
'atom'
:
return
else
:
self
.
changeatomnum
(
temp
,
data
)
def
changeatomnum
(
self
,
temp
,
data
=
' '
):
"""changes the atomnumbers in Lmpsmolecule structures angles, bonds, dihedrals, impropers,
and velocities in a three step process.
If data is defined extract the rest of the keywords used for Lmpsmolecule.
Step 1A:If data is defined copy from data one of the above structures.
If data is not defined skip this step.
Step 1B:If data is defined remove the rows of copy which do not contain any values from temp
Step 2: If data is defined extract from copy to an above Lmpsmolecule structure and remove the extracted rows
If data is not defined skip this step.
Step 3: Use temp and self.atoms to convert the atomnumbers in Lmpsmolecule structures
to the correct values.
temp and self.atoms line up, so temp[i] and self.atoms[i][0]
correspond to the current Lmpsmolecule structures atom numbers
and correct values
will use booleanarray class to ensure that lmpsmolecule's structure values are altered only once."""
#Step 1 and Step 2
if
data
!=
' '
:
#extract molecular keywords from data.keywords
for
keywords
in
data
.
keywords
:
if
keywords
==
'Angles'
or
keywords
==
'Bonds'
or
keywords
==
'Dihedrals'
or
\
keywords
==
'Impropers'
or
keywords
==
'Velocities'
:
self
.
keywords
.
append
(
keywords
)
for
keywords
in
self
.
keywords
:
if
keywords
!=
'Atoms'
:
print
'extracting the data from'
,
keywords
if
keywords
==
'Angles'
:
copy
=
self
.
copy
(
data
.
angles
)
#Step 1A
dnr
=
[]
#dnr means do not remove a 0 means remove and a 1 means keep
for
j
in
range
(
len
(
copy
)):
#Step 1B
dnr
.
append
(
0
)
#adds 0 to all list elements of dnr
for
item
in
temp
:
#finds copied structure that has temp values
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
dnr
[
j
]
=
1
#changes jth list element of dnr to 1
break
remove
=
[]
for
j
in
range
(
len
(
dnr
)):
# finds dnr values that are still 0
if
dnr
[
j
]
==
0
:
remove
.
append
(
j
)
#and appends their index value to the list remove
copy
=
self
.
deleterows
(
copy
,
remove
)
#removes all unneeded rows from copy
structnum
=
1
print
'the length of data is'
,
len
(
copy
)
for
item
in
temp
:
#Step 2
found
=
[]
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
found
.
append
(
j
)
self
.
angles
.
append
(
copy
[
j
])
l
=
len
(
self
.
angles
)
-
1
self
.
angles
[
l
][
0
]
=
str
(
structnum
)
structnum
+=
1
break
# print 'the item is', item
# print 'found is', found
copy
=
self
.
deleterows
(
copy
,
found
)
elif
keywords
==
'Bonds'
:
copy
=
self
.
copy
(
data
.
bonds
)
#Step 1A
dnr
=
[]
#dnr means do not remove a 0 means remove and a 1 means keep
for
j
in
range
(
len
(
copy
)):
#Step 1B
dnr
.
append
(
0
)
#adds 0 to all list elements of dnr
for
item
in
temp
:
#finds copied structure that has temp values
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
dnr
[
j
]
=
1
#changes jth list element of dnr to 1
break
remove
=
[]
for
j
in
range
(
len
(
dnr
)):
# finds dnr values that are still 0
if
dnr
[
j
]
==
0
:
remove
.
append
(
j
)
#and appends their index value to the list remove
copy
=
self
.
deleterows
(
copy
,
remove
)
#removes all unneeded rows from copy
structnum
=
1
print
'the length of data is'
,
len
(
copy
)
for
item
in
temp
:
#Step 2
found
=
[]
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
found
.
append
(
j
)
self
.
bonds
.
append
(
copy
[
j
])
l
=
len
(
self
.
bonds
)
-
1
self
.
bonds
[
l
][
0
]
=
str
(
structnum
)
structnum
+=
1
break
# print 'the item is', item
# print 'found is', found
copy
=
self
.
deleterows
(
copy
,
found
)
elif
keywords
==
'Dihedrals'
:
copy
=
self
.
copy
(
data
.
dihedrals
)
#Step 1A
dnr
=
[]
#dnr means do not remove a 0 means remove and a 1 means keep
for
j
in
range
(
len
(
copy
)):
#Step 1B
dnr
.
append
(
0
)
#adds 0 to all list elements of dnr
for
item
in
temp
:
#finds copied structure that has temp values
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
dnr
[
j
]
=
1
#changes jth list element of dnr to 1
break
remove
=
[]
for
j
in
range
(
len
(
dnr
)):
# finds dnr values that are still 0
if
dnr
[
j
]
==
0
:
remove
.
append
(
j
)
#and appends their index value to the list remove
copy
=
self
.
deleterows
(
copy
,
remove
)
#removes all unneeded rows from copy
structnum
=
1
print
'the length of data is'
,
len
(
copy
)
structnum
=
1
for
item
in
temp
:
#Step 2
found
=
[]
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
found
.
append
(
j
)
self
.
dihedrals
.
append
(
copy
[
j
])
l
=
len
(
self
.
dihedrals
)
-
1
self
.
dihedrals
[
l
][
0
]
=
str
(
structnum
)
structnum
+=
1
break
# print 'the item is', item
# print 'found is', found
copy
=
self
.
deleterows
(
copy
,
found
)
elif
keywords
==
'Impropers'
:
copy
=
self
.
copy
(
data
.
impropers
)
#Step 1B
dnr
=
[]
#dnr means do not remove a 0 means remove and a 1 means keep
for
j
in
range
(
len
(
copy
)):
#Step 1B
dnr
.
append
(
0
)
#adds 0 to all list elements of dnr
for
item
in
temp
:
#finds copied structure that has temp values
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
dnr
[
j
]
=
1
#changes jth list element of dnr to 1
break
remove
=
[]
for
j
in
range
(
len
(
dnr
)):
# finds dnr values that are still 0
if
dnr
[
j
]
==
0
:
remove
.
append
(
j
)
#and appends their index value to the list remove
copy
=
self
.
deleterows
(
copy
,
remove
)
#removes all unneeded rows from copy
structnum
=
1
print
'the length of data is'
,
len
(
copy
)
for
item
in
temp
:
#Step 2
found
=
[]
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
i
]
==
item
:
found
.
append
(
j
)
self
.
impropers
.
append
(
copy
[
j
])
l
=
len
(
self
.
impropers
)
-
1
self
.
impropers
[
l
][
0
]
=
str
(
structnum
)
structnum
+=
1
break
# print 'the item is', item
# print 'found is', found
copy
=
self
.
deleterows
(
copy
,
found
)
elif
keywords
==
'Velocities'
:
copy
=
self
.
copy
(
data
.
velocities
)
#Step 1
dnr
=
[]
#dnr means do not remove a 0 means remove and a 1 means keep
for
j
in
range
(
len
(
copy
)):
#Step 1B
dnr
.
append
(
0
)
#adds 0 to all list elements of dnr
for
item
in
temp
:
#finds copied structure that has temp values
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
1
):
if
copy
[
j
][
i
]
==
item
:
dnr
[
j
]
=
1
#changes jth list element of dnr to 1
break
remove
=
[]
for
j
in
range
(
len
(
dnr
)):
# finds dnr values that are still 0
if
dnr
[
j
]
==
0
:
remove
.
append
(
j
)
#and appends their index value to the list remove
copy
=
self
.
deleterows
(
copy
,
remove
)
#removes all unneeded rows from copy
structnum
=
1
print
'the length of data is'
,
len
(
copy
)
for
item
in
temp
:
#Step 2
found
=
[]
for
j
in
range
(
len
(
copy
)):
for
i
in
range
(
1
):
if
copy
[
j
][
i
]
==
item
:
found
.
append
(
j
)
self
.
velocities
.
append
(
copy
[
j
])
l
=
len
(
self
.
velocities
)
-
1
self
.
velocities
[
l
][
0
]
=
str
(
structnum
)
structnum
+=
1
break
# print 'the item is', item
# print 'found is', found
copy
=
self
.
deleterows
(
copy
,
found
)
#Step 3
for
keywords
in
self
.
keywords
:
if
keywords
!=
'Atoms'
:
print
'altering data structure values for'
,
keywords
if
keywords
==
'Angles'
:
copy
=
self
.
copy
(
self
.
angles
)
array
=
booleanarray
(
len
(
copy
),
len
(
copy
[
0
]),
True
)
#creating a booleanarray with true values
for
i
in
range
(
len
(
temp
)):
if
temp
[
i
]
==
self
.
atoms
[
i
][
0
]:
continue
for
j
in
range
(
len
(
copy
)):
for
k
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
k
]
==
temp
[
i
]:
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
self
.
angles
[
j
][
k
]
=
self
.
atoms
[
i
][
0
]
array
.
setelement
(
j
,
k
,
False
)
break
elif
keywords
==
'Bonds'
:
copy
=
self
.
copy
(
self
.
bonds
)
array
=
booleanarray
(
len
(
copy
),
len
(
copy
[
0
]),
True
)
#creating a booleanarray with true values
for
i
in
range
(
len
(
temp
)):
if
temp
[
i
]
==
self
.
atoms
[
i
][
0
]:
continue
for
j
in
range
(
len
(
copy
)):
for
k
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
k
]
==
temp
[
i
]:
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
self
.
bonds
[
j
][
k
]
=
self
.
atoms
[
i
][
0
]
array
.
setelement
(
j
,
k
,
False
)
break
elif
keywords
==
'Dihedrals'
:
copy
=
self
.
copy
(
self
.
dihedrals
)
array
=
booleanarray
(
len
(
copy
),
len
(
copy
[
0
]),
True
)
#creating a booleanarray with true values
for
i
in
range
(
len
(
temp
)):
if
temp
[
i
]
==
self
.
atoms
[
i
][
0
]:
continue
for
j
in
range
(
len
(
copy
)):
for
k
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
k
]
==
temp
[
i
]:
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
self
.
dihedrals
[
j
][
k
]
=
self
.
atoms
[
i
][
0
]
array
.
setelement
(
j
,
k
,
False
)
break
elif
keywords
==
'Impropers'
:
copy
=
self
.
copy
(
self
.
impropers
)
array
=
booleanarray
(
len
(
copy
),
len
(
copy
[
0
]),
True
)
#creating a booleanarray with true values
for
i
in
range
(
len
(
temp
)):
if
temp
[
i
]
==
self
.
atoms
[
i
][
0
]:
continue
for
j
in
range
(
len
(
copy
)):
for
k
in
range
(
2
,
len
(
copy
[
j
])):
if
copy
[
j
][
k
]
==
temp
[
i
]:
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
self
.
impropers
[
j
][
k
]
=
self
.
atoms
[
i
][
0
]
array
.
setelement
(
j
,
k
,
False
)
break
elif
keywords
==
'Velocities'
:
copy
=
self
.
copy
(
self
.
velocities
)
array
=
booleanarray
(
len
(
copy
),
len
(
copy
[
0
]),
True
)
#creating a booleanarray with true values
for
i
in
range
(
len
(
temp
)):
if
temp
[
i
]
==
self
.
atoms
[
i
][
0
]:
continue
for
j
in
range
(
len
(
copy
)):
for
k
in
range
(
1
):
if
copy
[
j
][
k
]
==
temp
[
i
]:
if
array
.
getelement
(
j
,
k
):
#if the boolean array is true
self
.
velocities
[
j
][
k
]
=
self
.
atoms
[
i
][
0
]
array
.
setelement
(
j
,
k
,
False
)
break
def
copy
(
self
,
structure
):
"""copies structure to same and returns same."""
same
=
[]
for
row
in
structure
:
same
.
append
(
row
)
return
same
def
deleterows
(
self
,
structure
,
rows
):
# run through this {structure is list of strings and rows is list
#of numbers}
"""delete rows in a structure and shifts the structure up
rows must be in increasing order for this algorithm to work correctly"""
new
=
[]
#mulitple copying of b to the rows being replaced.
if
rows
==
[]:
for
line
in
structure
:
new
.
append
(
line
)
#if no rows need replacing copy structure
return
new
for
i
in
range
(
rows
[
0
]):
new
.
append
(
structure
[
i
])
#copy structure to new until the first replaced row
count
=
0
for
i
in
range
(
rows
[
0
],
len
(
structure
)
-
len
(
rows
)):
# to replace rows and shift undeleted rows over
for
j
in
range
(
i
+
1
+
count
,
len
(
structure
)):
for
val
in
rows
:
if
val
==
j
:
count
+=
1
break
if
val
==
j
:
continue
else
:
new
.
append
(
structure
[
j
])
break
return
new
def
modifyatom
(
self
,
atomnumber
,
column
,
newvalue
):
"""modifies self.atom[atomnumber-1][column] to newvalue.
*note: newvalue needs to be a string
if column is 0 than changeatomnum method might need runnning for all atoms
that have had their atomnumber(column=0) modified.
changeatomnum method is not ran in this method when column=0"""
self
.
atoms
[
atomnumber
-
1
][
column
]
=
newvalue
def
deleteatoms
(
self
,
atomnumbers
,
atomid
):
"""Algorithm to find all atoms bonded to atomnumbers in the direction of atomid.
Atomid can be a list or a single integer. The single integer corresponds to atom's atomtype value.
The list corresponds to the atom's atomid values. The only difference between both cases is in the top portion of code.
When all bonded atoms have been found; ie: (the modified return values from findbonds yields an empty list);
delete those bonded atoms and all molecule structures which contain those atoms.
Convert the list of bonded atoms into rows and delete those rows from molecule.atoms"""
print
'finding atoms to delete'
bondedatoms
=
[]
# 1st iteration
nextatoms
=
self
.
findbonds
(
atomnumbers
)
try
:
# tests whether atomid is a list or a single integer
atomid
[
0
]
except
TypeError
:
#need to remove atoms from nextatoms which dont have the proper atom id
testflag
=
True
i
=
0
while
testflag
:
if
i
>
len
(
nextatoms
)
-
1
:
break
#For the case where i becomes greater than the len of nextatoms break loop
if
int
(
self
.
atoms
[
nextatoms
[
i
]
-
1
][
2
])
!=
atomid
:
#uses the atom id for the row in atoms and than checks atom id
del
nextatoms
[
i
]
#delets the atom at i
i
-=
1
#and than decreases i by 1 so next atom in the list will line up when i is increased
i
+=
1
#increase i by 1
else
:
#need to remove atoms from nextatoms which dont have the proper atom id
testflag
=
True
i
=
0
while
testflag
:
if
i
>
len
(
nextatoms
)
-
1
:
break
#For the case where i becomes greater than the len of nextatoms break loop
keep
=
False
for
id
in
atomid
:
if
int
(
self
.
atoms
[
nextatoms
[
i
]
-
1
][
0
])
==
id
:
#checking if atomid is in next atom
keep
=
True
#keep this atom
break
if
not
keep
:
del
nextatoms
[
i
]
#delets the atom at i
i
-=
1
#and than decreases i by 1 so next atom in the list will line up when i is increased
i
+=
1
#increase i by 1
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms
=
[]
for
atom
in
nextatoms
:
bondedatoms
.
append
(
atom
)
prevatoms
.
append
(
atom
)
#2nd iteration
if
prevatoms
==
[]:
print
'no bonds were found in first iteration that had atomid criteria'
return
nextatoms
=
self
.
findbonds
(
prevatoms
)
#need to remove atoms from nextatoms which are in atomnumbers
for
atom
in
atomnumbers
:
for
i
in
range
(
len
(
nextatoms
)):
if
nextatoms
[
i
]
==
atom
:
#checking if atom is in next atom
del
nextatoms
[
i
]
#delete the atom at i
break
if
nextatoms
==
[]:
break
#all bonds from find bonds have allready been added to bondedatoms
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms
=
[]
for
atom
in
nextatoms
:
bondedatoms
.
append
(
atom
)
prevatoms
.
append
(
atom
)
#iterative proccess for finding the rest of the atoms bonded to atomnumbers in the direction of atomid.
while
prevatoms
!=
[]:
nextatoms
=
self
.
findbonds
(
prevatoms
)
#need to remove atoms from nextatoms which are in the prevatoms
for
atom
in
bondedatoms
:
for
i
in
range
(
len
(
nextatoms
)):
if
nextatoms
[
i
]
==
atom
:
#checking if atom is in next atom
del
nextatoms
[
i
]
#delete the atom at i
break
if
nextatoms
==
[]:
break
#all bonds from find bonds have allready been added to bondedatoms
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms
=
[]
for
atom
in
nextatoms
:
bondedatoms
.
append
(
atom
)
prevatoms
.
append
(
atom
)
print
'the atoms to delete are'
,
bondedatoms
print
'deleting atoms from structures'
#delete bonded atoms from the molecule's structures except the atom structure
for
i
in
range
(
1
,
len
(
self
.
keywords
)):
#goes through all keywords except the atom keyword
print
self
.
keywords
[
i
]
if
self
.
keywords
[
i
]
==
'Angles'
:
rows
=
self
.
findatomnumbers
(
bondedatoms
,
self
.
angles
,
False
)
rows
=
self
.
listorder
(
rows
)
#to order the rows in increasing order
self
.
angles
=
self
.
deleterows
(
self
.
angles
,
rows
)
#requires that rows be in increasing order
elif
self
.
keywords
[
i
]
==
'Bonds'
:
rows
=
self
.
findatomnumbers
(
bondedatoms
,
self
.
bonds
,
False
)
rows
=
self
.
listorder
(
rows
)
#to order the rows in increasing order
self
.
bonds
=
self
.
deleterows
(
self
.
bonds
,
rows
)
#requires that rows be in increasing order
elif
self
.
keywords
[
i
]
==
'Dihedrals'
:
rows
=
self
.
findatomnumbers
(
bondedatoms
,
self
.
dihedrals
,
False
)
rows
=
self
.
listorder
(
rows
)
#to order the rows in increasing order
self
.
dihedrals
=
self
.
deleterows
(
self
.
dihedrals
,
rows
)
#requires that rows be in increasing order
elif
self
.
keywords
[
i
]
==
'Impropers'
:
rows
=
self
.
findatomnumbers
(
bondedatoms
,
self
.
impropers
,
False
)
rows
=
self
.
listorder
(
rows
)
#to order the rows in increasing order
self
.
impropers
=
self
.
deleterows
(
self
.
impropers
,
rows
)
#requires that rows be in increasing order
elif
self
.
keywords
[
i
]
==
'Velocities'
:
rows
=
self
.
findatomnumbers
(
bondedatoms
,
self
.
velocities
,
True
)
rows
=
self
.
listorder
(
rows
)
#to order the rows in increasing order
self
.
velocities
=
self
.
deleterows
(
self
.
velocities
,
rows
)
#requires that rows be in increasing order
print
'Atoms'
#convert bondedatoms from atom numbers to row numbers
for
i
in
range
(
len
(
bondedatoms
)):
bondedatoms
[
i
]
-=
1
bondedatoms
=
self
.
listorder
(
bondedatoms
)
#to order the row numbers in increasing order
#delete bonded atoms (row numbers) from the atom structure
self
.
atoms
=
self
.
deleterows
(
self
.
atoms
,
bondedatoms
)
#requires that row numbers be in increasing order
def
findatomnumbers
(
self
,
atomnumbers
,
structure
,
vflag
):
#need to read through this algorithm
"""Algorithm to find atomnumbers in a molecule structure except.
Atoms structure is not handled in here.
Returns a list of the rows in which the atomnumbers are contained in the molecule structure"""
rows
=
[]
if
vflag
:
#for handling velocity structure
for
atom
in
atomnumbers
:
for
i
in
range
(
len
(
structure
)):
if
int
(
structure
[
i
][
0
])
==
atom
:
#duplicate rows for vflag=True are not possible.
rows
.
append
(
i
)
break
else
:
#for handling all other structures except atoms
for
atom
in
atomnumbers
:
for
i
in
range
(
len
(
structure
)):
for
j
in
range
(
2
,
len
(
structure
[
i
])):
if
int
(
structure
[
i
][
j
])
==
atom
:
#need to make sure duplicate value of rows are not being added
if
rows
==
[]:
rows
.
append
(
i
)
#appends the row number
else
:
#checking for duplicate values of rows
duplicateflag
=
0
for
k
in
range
(
len
(
rows
)):
if
rows
[
k
]
==
i
:
duplicateflag
=
1
break
if
duplicateflag
==
0
:
#if no duplicates adds bond number
rows
.
append
(
i
)
#appends the row number
break
print
'the finished row is'
,
rows
return
rows
def
listorder
(
self
,
struct
):
"""Takes struct and organizes the list from least to greatest.
If the list is allready ordered this algorithm will do nothing."""
if
len
(
struct
)
==
1
:
return
struct
#with the length at 1; thier is only one element and theirfore nothing to order
for
i
in
range
(
1
,
len
(
struct
)):
#when i=0, struct will not change; therefore its skipped
copy
=
struct
[
i
]
for
j
in
range
(
i
-
1
,
-
1
,
-
1
):
struct
[
j
+
1
]
=
struct
[
j
]
if
copy
>
struct
[
j
]:
struct
[
j
+
1
]
=
copy
break
elif
j
==
0
:
struct
[
j
]
=
copy
#dont need a break here because this is the last j value in the for loop
print
'the organized row is'
,
struct
return
struct
def
findbonds
(
self
,
atomnumbers
):
"""Algorithm to find all atoms bonded to atomnumbers.
Returns a list of atomnumbers"""
#finds the bonds in which the atomnumbers are located
bondids
=
[]
for
i
in
range
(
len
(
atomnumbers
)):
for
j
in
range
(
len
(
self
.
bonds
)):
for
k
in
range
(
2
,
len
(
self
.
bonds
[
j
])):
if
int
(
self
.
bonds
[
j
][
k
])
==
atomnumbers
[
i
]:
if
bondids
==
[]:
bondids
.
append
(
int
(
self
.
bonds
[
j
][
0
]))
#appends the bond number
else
:
#checking for duplicates of bondids
duplicateflag
=
0
for
l
in
range
(
len
(
bondids
)):
if
bondids
[
l
]
==
int
(
self
.
bonds
[
j
][
0
]):
duplicateflag
=
1
break
if
duplicateflag
==
0
:
#if no duplicates adds bond number
bondids
.
append
(
int
(
self
.
bonds
[
j
][
0
]))
break
#Using bondids find the atoms bonded to atomnumbers
bondedatoms
=
[]
from
math
import
fabs
for
id
in
bondids
:
for
atoms
in
atomnumbers
:
found
=
False
for
i
in
range
(
2
,
len
(
self
.
bonds
[
id
-
1
])):
if
int
(
self
.
bonds
[
id
-
1
][
i
])
==
atoms
:
j
=
int
(
fabs
(
i
-
5
))
#switches the index from the atomnumber location to the other location
bondedatoms
.
append
(
int
(
self
.
bonds
[
id
-
1
][
j
]))
#appends the atomnuber at the other location
found
=
True
break
if
found
==
True
:
break
return
bondedatoms
def
findparticlebondingpoints
(
self
,
particle
,
atomid
,
cutoffdistance
,
bondnumber
):
"""Particle is a particlesurfaceobject, atomid is an int.
Finds atoms with atomid in the molecule which are less than the cutoffdistance from atoms in the particle.
The found atoms and particles locations are stored in a 3 dimensional list called possiblebonds
The first index corresponds with the molecule's atom
The second index correspons with the molecule/particle combination
The third index corresponds with whether the value is the molecule or the particle
Always bonds the two ends of the molecule that meet cutoff requirement.
All other possible bonds are randomly choosen until the required bondnumbers are met.
After every bond is choosen, the particle object's boolean list is updated, and possiblebonds is updated.
The update to possiblebonds involves removing the row from which the bonded molecule's atom is located
and also removing the particle atom and it's corresponding bonded atom from other rows of possiblebonds.
The final bonds are all stored in self.bonding as a 2 dimensional list"""
possiblebonds
=
[]
for
i
in
range
(
len
(
self
.
atoms
)):
if
int
(
self
.
atoms
[
i
][
2
])
==
atomid
:
row
=
[]
for
j
in
range
(
len
(
particle
.
surface
)):
#assumes molecule and particle have same atomtype if not true than atomdistance will give bad value
if
atomdistance
(
self
.
atoms
[
i
],
particle
.
surface
[
j
],
self
.
atomtype
)
<=
cutoffdistance
:
if
particle
.
getbool
(
j
):
row
.
append
([
i
,
j
])
#if not bonded than can form a possible bond
if
row
!=
[]:
possiblebonds
.
append
(
row
)
#need to correct this....
#initiate section which assigns bonds into bonding information
bonds
=
0
self
.
bondinginformation
=
[]
#initiate new class member
#Checks to see if no bonds are possible [2 possible cases]
if
possiblebonds
==
[]:
print
'no possible bonds can be formed'
return
if
bondnumber
==
0
:
print
'bondnumber is 0; so, no possible bonds can form'
return
#section which assigns a bond to the first molecule atom which can be bonded.
self
.
bondinginformation
.
append
(
self
.
particlebondinginfo
(
particle
,
possiblebonds
,
0
))
bonds
+=
1
if
bonds
==
bondnumber
:
return
del
possiblebonds
[
0
]
#deletes possible bonds to the first molecule which can be bonded
l
=
len
(
self
.
bondinginformation
)
-
1
possiblebonds
=
self
.
updatepossiblebonds
(
possiblebonds
,
self
.
bondinginformation
[
l
][
1
])
#updates possiblebonds
#by removing any bonds which contain the newly bonded particle.
if
possiblebonds
==
[]:
return
#section which finds the last molecule atom which can be bonded
l
=
len
(
possiblebonds
)
-
1
while
possiblebonds
[
l
]
==
[]:
# to find the last molecule atom which can be bonded
if
possiblebonds
==
[]:
return
del
possiblebonds
[
l
]
#since there are no more possible bonds in this location delete
l
-=
1
#go to next possible spot where the last molecule atom could be bonded
#section which assigns a bond to the last molecule atom which can be bonded.
self
.
bondinginformation
.
append
(
self
.
particlebondinginfo
(
particle
,
possiblebonds
,
l
))
bonds
+=
1
if
bonds
==
bondnumber
:
return
del
possiblebonds
[
l
]
#deletes possible bonds to the last molecule which can be bonded
l
=
len
(
self
.
bondinginformation
)
-
1
possiblebonds
=
self
.
updatepossiblebonds
(
possiblebonds
,
self
.
bondinginformation
[
l
][
1
])
if
bondnumber
-
bonds
>=
len
(
possiblebonds
):
#the rest of the bonds are assigned in order
while
possiblebonds
!=
[]:
if
possiblebonds
[
0
]
==
[]:
#if row of possible bonds is empty then delete the row
del
possiblebonds
[
0
]
else
:
#else find a bond in the row of possible bonds
self
.
bondinginformation
.
append
(
self
.
particlebondinginfo
(
particle
,
possiblebonds
,
0
))
#dont need to update bonds since possiblebonds will become empty at the same time or before bonds
#is equal to bondnumber
del
possiblebonds
[
0
]
l
=
len
(
self
.
bondinginformation
)
-
1
possiblebonds
=
self
.
updatepossiblebonds
(
possiblebonds
,
self
.
bondinginformation
[
l
][
1
])
return
else
:
#the rest of the bonds are assigned randomly
from
random
import
randint
#use to randomly choose an index to bond.
while
bonds
<
bondnumber
:
if
possiblebonds
==
[]:
break
#exits while loop when possiblebonds has become an empty set.
#this ensures that in the case there are not egnough viable bonds from possiblebonds to
#reach the bondnumber. Than, the while loop will not become infinite.
l
=
len
(
possiblebonds
)
-
1
i
=
randint
(
0
,
l
)
if
possiblebonds
[
i
]
==
[]:
#if row of possible bonds is empty then delete the row
del
possiblebonds
[
i
]
else
:
#else find a bond in the row of possible bonds
self
.
bondinginformation
.
append
(
self
.
particlebondinginfo
(
particle
,
possiblebonds
,
i
))
bonds
+=
1
del
possiblebonds
[
i
]
l
=
len
(
self
.
bondinginformation
)
-
1
possiblebonds
=
self
.
updatepossiblebonds
(
possiblebonds
,
self
.
bondinginformation
[
l
][
1
])
return
def
particlebondinginfo
(
self
,
particle
,
possiblebonds
,
index1
):
"""Takes list of possiblebonds and assigns a bond from possiblebonds[index1].
Then assigns false to particle.setbool at the bonding location
to ensure no more bonds can form with that particle.
Returns the assigned bond."""
from
random
import
randint
#use to randomly choose 0 and 1 where 1 is keep.
for
i
in
range
(
len
(
possiblebonds
[
index1
])):
if
i
==
len
(
possiblebonds
[
index1
])
-
1
:
#automattically bonds under this condition
break
else
:
#bond has random chance of forming
if
randint
(
0
,
1
)
==
1
:
#bond forms
break
else
:
#no bond forms
continue
particle
.
setbool
(
possiblebonds
[
index1
][
i
][
1
],
False
)
#makes sure no more bonds can form
return
possiblebonds
[
index1
][
i
]
def
updatepossiblebonds
(
self
,
possiblebonds
,
particlenum
):
"""finds the particle number in the remaining possible bonds than deletes that bonding information.
This algorithm assumes that particlenum can exist only once in each row of possiblebonds."""
for
i
in
range
(
len
(
possiblebonds
)):
for
j
in
range
(
len
(
possiblebonds
[
i
])):
if
possiblebonds
[
i
][
j
][
1
]
==
particlenum
:
del
possiblebonds
[
i
][
j
]
break
#go to next row of possiblebonds
return
possiblebonds
def
bondtoparticle
(
self
,
particle
,
atomid
,
newid
,
newcharge
):
"""Bonds molecule to particle. Particle is a particlesurface object.
Moves the molecule's atoms bonded to the particle to the particle's position.
Alters the atomtype of the molecule's atoms bonded to the particle to newid.
Alters the charge of the molecule's atoms bonded to the particle to newcharge
Because bonds have formed between the molecule's atoms and the particle's atoms,
atoms with atomid on the molecule need to be removed otherwise the molecule's atoms will have to many bonds.
self.deleteatoms will take care of deleting the atoms atomid and
atoms down the polymer chain in the direction away from the molecule/particle bond.
Now remove the bonded particle atoms from the particle surface becaused the bonded molecule atoms
have replaced those particle atoms and their bonds with the rest of the particle.
Newid must be a string. Atomid can now be a list or an integer.
The list is a list of atom's id values and the single integer is atom's atomtype"""
#moves the molecule's atoms bonded to the particle to the particle's position
#need to do this step first before the atom id's and row indexes get out of sync
#which will occur after atoms get deleted.
print
'modifying the bonded molecules atom information'
for
i
in
range
(
len
(
self
.
bondinginformation
)):
#patom=particle.getsurfatom(self.bondinginformation[i][1])
#if self.atomtype=='molecular':#3-5->x,y,z[molecular]
# for j in range(3,6):
# self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here
#elif self.atomtype=='full':#4-6->x,y,z[full]
# for j in range(4,7):
# self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here
#alters the atomtype to newid of the molecule's atoms bonded to the particle.
self
.
modifyatom
(
self
.
bondinginformation
[
i
][
0
]
+
1
,
2
,
newid
)
#uses the atomid here
if
self
.
atomtype
==
'full'
:
# Alters the charge of the molecule's atoms bonded to the particle to newcharge
self
.
modifyatom
(
self
.
bondinginformation
[
i
][
0
]
+
1
,
3
,
newcharge
)
#This is a seperate loop so atom id's and row indexes for previous steps wont get out of sync
#create atomnumbers to begin deletion process.
atomnumbers
=
[]
for
i
in
range
(
len
(
self
.
bondinginformation
)):
atomnumbers
.
append
(
self
.
bondinginformation
[
i
][
0
]
+
1
)
#using actual atomnumbers rather than the row index
#Call algorithm to find all atoms bonded to atomnumbers in the direction of atomid.
#Than delete those atoms and all molecule structures which contain those atoms.
self
.
deleteatoms
(
atomnumbers
,
atomid
)
print
'begining deletion process of the surface atoms for which the molecule atoms have replaced'
#Goes through the bondinginformation and superficially removes the surfaceatom
for
i
in
range
(
len
(
self
.
bondinginformation
)):
particle
.
removesurfatom
(
self
.
bondinginformation
[
i
][
1
])
#uses the row number
#Allows the particle extract method used outside this class to remove this atom.
def
createxyz
(
self
,
file
,
data
,
routine
=
'mass'
,
values
=
None
):
"""Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f
=
open
(
file
,
'w'
)
f
.
write
(
'{0}
\n
'
.
format
(
len
(
self
.
atoms
)))
f
.
write
(
'atoms
\n
'
)
if
routine
==
'mass'
:
for
line
in
self
.
atoms
:
type
=
int
(
line
[
2
])
-
1
mass
=
data
.
masses
[
type
][
1
]
if
mass
==
'12.0107'
:
elementnum
=
6
elif
mass
==
'15.9994'
:
elementnum
=
8
elif
mass
==
'26.9815'
:
elementnum
=
13
else
:
print
'no matching mass value. The method will exit'
return
l
=
len
(
line
)
-
1
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
f
.
close
()
elif
routine
==
'atomtype'
:
for
line
in
self
.
atoms
:
type
=
int
(
line
[
2
])
-
1
elementnum
=
values
[
type
]
l
=
len
(
line
)
-
1
f
.
write
(
'{0} {1} {2} {3}
\n
'
.
format
(
elementnum
,
line
[
l
-
5
],
line
[
l
-
4
],
line
[
l
-
3
]))
f
.
close
()
Event Timeline
Log In to Comment