Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87867954
analyse.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
Tue, Oct 15, 11:33
Size
6 KB
Mime Type
text/x-python
Expires
Thu, Oct 17, 11:33 (2 d)
Engine
blob
Format
Raw Data
Handle
21670985
Attached To
rLAMMPS lammps
analyse.py
View Options
#!/usr/bin/env python3
import
numpy
as
np
import
subprocess
as
sub
import
shlex
import
matplotlib.pyplot
as
plt
sceleton
=
"""
# 3d Lennard-Jones melt
units metal
atom_style atomic
atom_modify map array
timestep 0
region box block 0 12 0 12 0 12
boundary p p p
create_box 3 box
create_atoms 1 single {pos1} units box
create_atoms 2 single {pos2} units box
create_atoms 3 single {pos3} units box
create_atoms 3 single {pos4} units box
mass 1 1.0
mass 2 1.0
mass 3 1.0
pair_style {pot}neuronet
pair_coeff * * in.const.NN_short in.params.NN_short
neighbor 0.3 bin
compute 1 all pe
variable pppp equal c_1
variable f1x equal "fx[1]"
variable f1y equal "fy[1]"
variable f1z equal "fz[1]"
variable f2x equal "fx[2]"
variable f2y equal "fy[2]"
variable f2z equal "fz[2]"
variable r1x equal "x[1]"
variable r1y equal "y[1]"
variable r1z equal "z[1]"
variable r2x equal "x[2]"
variable r2y equal "y[2]"
variable r2z equal "z[2]"
fix 1 all nve
fix 2 all print 1 "POTENTIAL ENERGY DIRECT ${pppp}"
fix 3 all print 1 "FORCE VECTOR DIRECT ${{f1x}} ${{f1y}} ${{f1z}} ${{f2x}} ${{f2y}} ${{f2z}} "
fix 4 all print 1 "POSITION VECTOR DIRECT ${{r1x}} ${{r1y}} ${{r1z}} ${{r2x}} ${{r2y}} ${{r2z}} "
dump 1 all custom 1 dump.out id xs ys zs fx fy fz
run 2
"""
def
get_sceleton
(
del_x
=
0.
,
pot
=
''
,
shift_x
=
0
,
shift_y
=
0
,
shift_z
=
0
):
box_siz
=
12.
pos1
=
"{} {} {}"
.
format
((
0
+
shift_x
)
%
box_siz
,
(
0
+
shift_y
)
%
box_siz
,
(
del_x
+
shift_z
)
%
box_siz
)
pos2
=
"{} {} {}"
.
format
((
3
+
shift_x
)
%
box_siz
,
(
3
+
del_x
+
shift_y
)
%
box_siz
,
(
0
+
shift_z
)
%
box_siz
)
pos3
=
"{} {} {}"
.
format
((
3
+
shift_x
)
%
box_siz
,
(
0
+
shift_y
)
%
box_siz
,
(
3
+
shift_z
)
%
box_siz
)
pos4
=
"{} {} {}"
.
format
((
3
+
shift_x
)
%
box_siz
,
(
3
+
shift_y
)
%
box_siz
,
(
3
+
shift_z
)
%
box_siz
)
return
sceleton
.
format
(
pos1
=
pos1
,
pos2
=
pos2
,
pos3
=
pos3
,
pos4
=
pos4
,
pot
=
pot
,
pppp
=
'{pppp}'
)
def
nextline
(
fh
,
do_print
=
False
):
line
=
fh
.
__next__
()
if
do_print
:
print
(
line
)
return
line
def
get_force_del_pos
(
fname
=
"dump.out"
):
do_print
=
False
with
open
(
fname
,
'r'
)
as
fh
:
line
=
''
while
not
line
.
startswith
(
"ITEM: NUMBER OF ATOMS"
):
line
=
nextline
(
fh
,
do_print
)
nb_atoms
=
int
(
nextline
(
fh
,
do_print
))
while
not
line
.
startswith
(
"ITEM: ATOMS id xs ys zs fx fy fz"
):
line
=
nextline
(
fh
,
do_print
)
force
=
[]
pos_init
=
[]
for
i
in
range
(
nb_atoms
):
vals
=
[
float
(
val
)
for
val
in
nextline
(
fh
,
do_print
)
.
split
()[
-
6
:]]
pos_init
+=
vals
[:
3
]
force
+=
vals
[
-
3
:]
pos_fin
=
[]
line
=
''
while
not
line
.
startswith
(
"ITEM: ATOMS id xs ys zs fx fy fz"
):
line
=
nextline
(
fh
,
do_print
)
for
i
in
range
(
nb_atoms
):
vals
=
[
float
(
val
)
for
val
in
nextline
(
fh
,
do_print
)
.
split
()[
-
6
:]]
pos_fin
+=
vals
[:
3
]
force
=
np
.
array
(
force
)
pos_init
=
12
*
np
.
array
(
pos_init
)
pos_fin
=
12
*
np
.
array
(
pos_fin
)
del_pos
=
pos_fin
-
pos_init
print
(
"force = {}"
.
format
(
force
))
print
(
"pos = {}"
.
format
(
pos_init
))
return
force
,
pos_init
def
get_energy
(
fname
=
"log.lammps"
):
do_print
=
False
with
open
(
fname
,
'r'
)
as
fh
:
line
=
""
while
not
line
.
startswith
(
"potential energy"
):
line
=
nextline
(
fh
,
do_print
)
epot_init
=
float
(
line
.
split
()[
-
1
])
print
(
"E_pot = {}"
.
format
(
epot_init
))
return
epot_init
def
get_everything
(
output
):
do_print
=
False
def
nextliner
(
do_print
):
for
line
in
output
.
decode
()
.
split
(
'
\n
'
):
if
do_print
:
print
(
line
)
yield
line
gen
=
nextliner
(
do_print
)
def
nextline
(
do_print
):
return
gen
.
__next__
()
line
=
""
while
not
line
.
startswith
(
"POTENTIAL ENERGY DIRECT"
):
line
=
nextline
(
do_print
)
e_pot
=
float
(
line
.
split
()[
-
1
])
force
=
np
.
array
([
float
(
item
)
for
item
in
nextline
(
do_print
)
.
split
()[
3
:]])
pos
=
np
.
array
([
float
(
item
)
for
item
in
nextline
(
do_print
)
.
split
()[
3
:]])
return
e_pot
,
force
,
pos
def
eval_derivative
(
del_pos
,
force
,
del_epot
):
estim
=
np
.
dot
(
-
force
,
del_pos
)
print
(
"Δx = {}"
.
format
(
del_pos
))
print
(
"ΔE = {}, Δx.-f = {}, err = {}"
.
format
(
del_epot
,
estim
,
1
-
estim
/
del_epot
))
def
run
(
del_x
,
pot
,
shift_x
=
0
,
shift_y
=
0
,
shift_z
=
0
):
command
=
"../../src/lmp_serial"
with
sub
.
Popen
(
shlex
.
split
(
command
),
stdin
=
sub
.
PIPE
,
stdout
=
sub
.
PIPE
)
as
proc
:
input
=
get_sceleton
(
del_x
,
pot
,
shift_x
,
shift_y
,
shift_z
)
#print(input)
out
,
err
=
proc
.
communicate
(
input
=
input
.
encode
())
proc
.
wait
()
#print(out.decode())
return
out
def
analyse_case
(
name
,
del_x
,
pot
):
print
(
"For {}:"
.
format
(
name
))
output
=
run
(
0.
,
pot
)
e_i
,
f
,
x_i
=
get_everything
(
output
)
print
(
"x_i = {}"
.
format
(
x_i
))
output
=
run
(
del_x
,
pot
)
e_f
,
dummy
,
x_f
=
get_everything
(
output
)
print
(
"x_f = {}"
.
format
(
x_f
))
print
(
"E_i = {}, E_f = {}"
.
format
(
e_i
,
e_f
))
eval_derivative
(
x_f
-
x_i
,
f
,
e_f
-
e_i
)
def
shift_case
(
name
,
del_x
,
pot
,
dir
,
nb_step
):
print
(
"For {}:"
.
format
(
name
))
e_pots
=
np
.
zeros
(
nb_step
,
dtype
=
float
)
errs
=
np
.
zeros_like
(
e_pots
)
estims
=
np
.
zeros_like
(
e_pots
)
for
i
in
range
(
nb_step
):
shifts
=
tuple
((
d
*
i
for
d
in
dir
))
output
=
run
(
0
,
pot
,
*
shifts
)
e_pots
[
i
],
f
,
x_i
=
get_everything
(
output
)
output
=
run
(
del_x
,
pot
,
*
shifts
)
e_f
,
dummy
,
x_f
=
get_everything
(
output
)
estims
[
i
]
=
np
.
dot
(
-
f
,
x_f
-
x_i
)
errs
[
i
]
=
1
-
estims
[
i
]
/
(
e_f
-
e_pots
[
i
])
del_epots
=
e_pots
-
e_pots
[
0
]
print
(
"E_pot0 = {}"
.
format
(
e_pots
[
0
]))
print
(
"ΔE_pot = {}"
.
format
(
del_epots
))
print
(
"estimated = {}"
.
format
(
estims
))
print
(
"force_err = {}"
.
format
(
errs
))
return
del_epots
def
main
():
del_x
=
1.e-8
cases
=
[(
"Ryo"
,
del_x
,
''
),
(
"Me"
,
del_x
,
'c'
)]
for
case
in
cases
:
analyse_case
(
*
case
)
print
(
"
\n\n
Shifting solid motion
\n
"
)
dir
=
[
3
,
4
,
6.48
]
for
case
in
cases
:
shift_case
(
*
case
,
dir
,
4
)
if
__name__
==
"__main__"
:
main
()
Event Timeline
Log In to Comment