adding Pair::single() support to python pair style and examples

with the single function, python pair styles can be massively
sped up and made compatible to accelerators, as one can translate
the analytic force and energy functions through LAMMPS into suitable
tables and then simply use the on-the-fly tables for production runs
This commit is contained in:
Axel Kohlmeyer 2017-05-17 17:20:56 -04:00
parent 45becfb235
commit 43efe9e417
9 changed files with 9743 additions and 2 deletions

9029
examples/python/data.spce Normal file

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,45 @@
units real
atom_style full
read_data data.spce
pair_style hybrid/overlay python 12.0 coul/long 12.0
kspace_style pppm 1.0e-6
pair_coeff * * coul/long
pair_coeff * * python potentials.LJCutSPCE OW NULL
pair_modify table 0
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47
special_bonds lj/coul 0.0 0.0 1.0
neighbor 2.0 bin
fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300.0 300.0 100.0
# create combined lj/coul table for all atom types
# generate tabulated potential from python variant
pair_write 1 1 2000 rsq 0.1 12 spce.table OW-OW -0.8472 -0.8472
pair_write 1 2 2000 rsq 0.1 12 spce.table OW-HW -0.8472 0.4236
pair_write 2 2 2000 rsq 0.1 12 spce.table HW-HW 0.4236 0.4236
# switch to tabulated potential
pair_style table linear 2000 pppm
pair_coeff 1 1 spce.table OW-OW
pair_coeff 1 2 spce.table OW-HW
pair_coeff 2 2 spce.table HW-HW
thermo 10
run 100
shell rm spce.table

View File

@ -0,0 +1,32 @@
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass * 1.0
velocity all create 3.0 87287
pair_style python 2.5
pair_coeff * * potentials.LJCutMelt lj
# generate tabulated potential from python variant
pair_write 1 1 10000 rsq 0.01 2.5 lj_1_1.table LJ
pair_style table linear 10000
pair_coeff 1 1 lj_1_1.table LJ
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
thermo 50
run 250
shell rm lj_1_1.table

View File

@ -0,0 +1,178 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units real
atom_style full
read_data data.spce
orthogonal box = (0.02645 0.02645 0.02641) to (35.5328 35.5328 35.4736)
1 by 1 by 1 MPI processor grid
reading atoms ...
4500 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
3000 bonds
reading angles ...
1500 angles
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
pair_style hybrid/overlay python 12.0 coul/long 12.0
kspace_style pppm 1.0e-6
pair_coeff * * coul/long
pair_coeff * * python potentials.LJCutSPCE OW NULL
pair_modify table 0
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47
special_bonds lj/coul 0.0 0.0 1.0
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
2 = max # of special neighbors
neighbor 2.0 bin
fix 1 all shake 0.0001 20 0 b 1 a 1
0 = # of size 2 clusters
0 = # of size 3 clusters
0 = # of size 4 clusters
1500 = # of frozen angles
fix 2 all nvt temp 300.0 300.0 100.0
# create combined lj/coul table for all atom types
# generate tabulated potential from python variant
pair_write 1 1 2000 rsq 0.1 12 spce.table OW-OW -0.8472 -0.8472
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 103823 64000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair python, perpetual, skip from (2)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
pair_write 1 2 2000 rsq 0.1 12 spce.table OW-HW -0.8472 0.4236
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 103823 64000
pair_write 2 2 2000 rsq 0.1 12 spce.table HW-HW 0.4236 0.4236
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 103823 64000
# switch to tabulated potential
pair_style table linear 2000 pppm
pair_coeff 1 1 spce.table OW-OW
pair_coeff 1 2 spce.table OW-HW
pair_coeff 2 2 spce.table HW-HW
thermo 10
run 100
PPPM initialization ...
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394674
estimated relative force accuracy = 1.18855e-06
using double precision FFTs
3d grid and FFT values/proc = 103823 64000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 35.26 | 35.26 | 35.26 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -100272.97 0 -100272.97 -1282.0708
10 120.61568 -101350.63 0 -100272.39 -4077.5051
20 136.11379 -101465.43 0 -100248.65 -5136.5677
30 137.01602 -101455.3 0 -100230.46 -5347.8311
40 153.424 -101582.46 0 -100210.93 -5223.1676
50 167.73654 -101686.24 0 -100186.77 -4468.6687
60 163.11642 -101618.16 0 -100159.99 -3291.7815
70 169.64512 -101647.89 0 -100131.35 -2611.638
80 182.9979 -101737.01 0 -100101.11 -2390.6293
90 191.33873 -101778.71 0 -100068.24 -2239.386
100 194.7458 -101775.84 0 -100034.92 -1951.9128
Loop time of 7.60221 on 1 procs for 100 steps with 4500 atoms
Performance: 1.137 ns/day, 21.117 hours/ns, 13.154 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.7401 | 5.7401 | 5.7401 | 0.0 | 75.51
Bond | 0.00017881 | 0.00017881 | 0.00017881 | 0.0 | 0.00
Kspace | 1.5387 | 1.5387 | 1.5387 | 0.0 | 20.24
Neigh | 0.2299 | 0.2299 | 0.2299 | 0.0 | 3.02
Comm | 0.024311 | 0.024311 | 0.024311 | 0.0 | 0.32
Output | 0.00057936 | 0.00057936 | 0.00057936 | 0.0 | 0.01
Modify | 0.063158 | 0.063158 | 0.063158 | 0.0 | 0.83
Other | | 0.005243 | | | 0.07
Nlocal: 4500 ave 4500 max 4500 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 21216 ave 21216 max 21216 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 2.60177e+06 ave 2.60177e+06 max 2.60177e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 2601766
Ave neighs/atom = 578.17
Ave special neighs/atom = 2
Neighbor list builds = 3
Dangerous builds = 0
shell rm spce.table
Total wall time: 0:00:07

View File

@ -0,0 +1,178 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units real
atom_style full
read_data data.spce
orthogonal box = (0.02645 0.02645 0.02641) to (35.5328 35.5328 35.4736)
2 by 2 by 1 MPI processor grid
reading atoms ...
4500 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
3000 bonds
reading angles ...
1500 angles
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
pair_style hybrid/overlay python 12.0 coul/long 12.0
kspace_style pppm 1.0e-6
pair_coeff * * coul/long
pair_coeff * * python potentials.LJCutSPCE OW NULL
pair_modify table 0
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47
special_bonds lj/coul 0.0 0.0 1.0
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
2 = max # of special neighbors
neighbor 2.0 bin
fix 1 all shake 0.0001 20 0 b 1 a 1
0 = # of size 2 clusters
0 = # of size 3 clusters
0 = # of size 4 clusters
1500 = # of frozen angles
fix 2 all nvt temp 300.0 300.0 100.0
# create combined lj/coul table for all atom types
# generate tabulated potential from python variant
pair_write 1 1 2000 rsq 0.1 12 spce.table OW-OW -0.8472 -0.8472
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 34263 16000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair python, perpetual, skip from (2)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
pair_write 1 2 2000 rsq 0.1 12 spce.table OW-HW -0.8472 0.4236
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 34263 16000
pair_write 2 2 2000 rsq 0.1 12 spce.table HW-HW 0.4236 0.4236
PPPM initialization ...
WARNING: Using polynomial approximation for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394206
estimated relative force accuracy = 1.18714e-06
using double precision FFTs
3d grid and FFT values/proc = 34263 16000
# switch to tabulated potential
pair_style table linear 2000 pppm
pair_coeff 1 1 spce.table OW-OW
pair_coeff 1 2 spce.table OW-HW
pair_coeff 2 2 spce.table HW-HW
thermo 10
run 100
PPPM initialization ...
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.279652
grid = 40 40 40
stencil order = 5
estimated absolute RMS force accuracy = 0.000394674
estimated relative force accuracy = 1.18855e-06
using double precision FFTs
3d grid and FFT values/proc = 34263 16000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 13.05 | 13.05 | 13.05 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -100272.97 0 -100272.97 -1282.0708
10 120.61568 -101350.63 0 -100272.39 -4077.5051
20 136.11379 -101465.43 0 -100248.65 -5136.5677
30 137.01602 -101455.3 0 -100230.46 -5347.8311
40 153.424 -101582.46 0 -100210.93 -5223.1676
50 167.73654 -101686.24 0 -100186.77 -4468.6687
60 163.11642 -101618.16 0 -100159.99 -3291.7815
70 169.64512 -101647.89 0 -100131.35 -2611.638
80 182.9979 -101737.01 0 -100101.11 -2390.6293
90 191.33873 -101778.71 0 -100068.24 -2239.386
100 194.7458 -101775.84 0 -100034.92 -1951.9128
Loop time of 2.38392 on 4 procs for 100 steps with 4500 atoms
Performance: 3.624 ns/day, 6.622 hours/ns, 41.948 timesteps/s
99.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.4377 | 1.5465 | 1.6848 | 7.3 | 64.87
Bond | 0.00010276 | 0.00012648 | 0.0001452 | 0.0 | 0.01
Kspace | 0.53311 | 0.66842 | 0.77484 | 10.9 | 28.04
Neigh | 0.066 | 0.066074 | 0.066101 | 0.0 | 2.77
Comm | 0.045355 | 0.048344 | 0.050747 | 1.0 | 2.03
Output | 0.00042391 | 0.00044996 | 0.00052667 | 0.0 | 0.02
Modify | 0.049891 | 0.050191 | 0.050336 | 0.1 | 2.11
Other | | 0.003771 | | | 0.16
Nlocal: 1125 ave 1154 max 1092 min
Histogram: 1 0 0 0 1 0 0 1 0 1
Nghost: 12256.2 ave 12296 max 12213 min
Histogram: 1 0 1 0 0 0 0 0 1 1
Neighs: 650442 ave 678824 max 626375 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Total # of neighbors = 2601766
Ave neighs/atom = 578.17
Ave special neighs/atom = 2
Neighbor list builds = 3
Dangerous builds = 0
shell rm spce.table
Total wall time: 0:00:02

View File

@ -0,0 +1,99 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6796 1.6796 1.6796
region box block 0 10 0 10 0 10
create_box 1 box
Created orthogonal box = (0 0 0) to (16.796 16.796 16.796)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 4000 atoms
mass * 1.0
velocity all create 3.0 87287
pair_style python 2.5
pair_coeff * * potentials.LJCutMelt lj
# generate tabulated potential from python variant
pair_write 1 1 10000 rsq 0.01 2.5 lj_1_1.table LJ
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair python, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
pair_style table linear 10000
pair_coeff 1 1 lj_1_1.table LJ
WARNING: 1 of 1000 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (../pair_table.cpp:476)
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
thermo 50
run 250
Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.184 | 3.184 | 3.184 Mbytes
Step Temp E_pair E_mol TotEng Press
0 3 -6.7733675 0 -2.2744925 -3.7033435
50 1.6758875 -4.7951764 0 -2.2819736 5.6705794
100 1.6458266 -4.7488945 0 -2.2807717 5.8696895
150 1.6324439 -4.7283321 0 -2.2802784 5.9594952
200 1.6630547 -4.7746809 0 -2.2807225 5.7372657
250 1.6278968 -4.7226363 0 -2.2814016 5.9559236
Loop time of 1.0498 on 1 procs for 250 steps with 4000 atoms
Performance: 102877.190 tau/day, 238.142 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.93242 | 0.93242 | 0.93242 | 0.0 | 88.82
Neigh | 0.088495 | 0.088495 | 0.088495 | 0.0 | 8.43
Comm | 0.012153 | 0.012153 | 0.012153 | 0.0 | 1.16
Output | 0.00013924 | 0.00013924 | 0.00013924 | 0.0 | 0.01
Modify | 0.013729 | 0.013729 | 0.013729 | 0.0 | 1.31
Other | | 0.002855 | | | 0.27
Nlocal: 4000 ave 4000 max 4000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 5504 ave 5504 max 5504 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 151497 ave 151497 max 151497 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 151497
Ave neighs/atom = 37.8743
Neighbor list builds = 12
Dangerous builds not checked
shell rm lj_1_1.table
Total wall time: 0:00:01

View File

@ -0,0 +1,99 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6796 1.6796 1.6796
region box block 0 10 0 10 0 10
create_box 1 box
Created orthogonal box = (0 0 0) to (16.796 16.796 16.796)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 4000 atoms
mass * 1.0
velocity all create 3.0 87287
pair_style python 2.5
pair_coeff * * potentials.LJCutMelt lj
# generate tabulated potential from python variant
pair_write 1 1 10000 rsq 0.01 2.5 lj_1_1.table LJ
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair python, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
pair_style table linear 10000
pair_coeff 1 1 lj_1_1.table LJ
WARNING: 1 of 10000 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (../pair_table.cpp:476)
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
thermo 50
run 250
Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.69 | 2.69 | 2.69 Mbytes
Step Temp E_pair E_mol TotEng Press
0 3 -6.7733675 0 -2.2744925 -3.7033435
50 1.6754092 -4.794723 0 -2.2822376 5.6616601
100 1.6503295 -4.7559815 0 -2.2811061 5.8051261
150 1.6596603 -4.7699379 0 -2.2810699 5.7830168
200 1.6371948 -4.7365549 0 -2.2813766 5.9245585
250 1.6321199 -4.7288017 0 -2.2812339 5.9776124
Loop time of 0.313548 on 4 procs for 250 steps with 4000 atoms
Performance: 344444.576 tau/day, 797.325 timesteps/s
99.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.24963 | 0.25443 | 0.2632 | 1.1 | 81.15
Neigh | 0.023249 | 0.023735 | 0.024497 | 0.3 | 7.57
Comm | 0.020689 | 0.030402 | 0.035249 | 3.4 | 9.70
Output | 0.00020766 | 0.00021476 | 0.00023031 | 0.0 | 0.07
Modify | 0.0034959 | 0.0035564 | 0.0036762 | 0.1 | 1.13
Other | | 0.001206 | | | 0.38
Nlocal: 1000 ave 1010 max 982 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Nghost: 2703.75 ave 2713 max 2689 min
Histogram: 1 0 0 0 0 0 0 2 0 1
Neighs: 37915.5 ave 39231 max 36202 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Total # of neighbors = 151662
Ave neighs/atom = 37.9155
Neighbor list builds = 12
Dangerous builds not checked
shell rm lj_1_1.table
Total wall time: 0:00:00

View File

@ -35,7 +35,7 @@ using namespace LAMMPS_NS;
PairPython::PairPython(LAMMPS *lmp) : Pair(lmp) {
respa_enable = 0;
single_enable = 0;
single_enable = 1;
writedata = 0;
restartinfo = 0;
one_coeff = 1;
@ -46,7 +46,7 @@ PairPython::PairPython(LAMMPS *lmp) : Pair(lmp) {
py_potential = NULL;
// add current directory to PYTHONPATH
PyObject * py_path = PySys_GetObject("path");
PyObject * py_path = PySys_GetObject((char *)"path");
PyList_Append(py_path, PY_STRING_FROM_STRING("."));
// if LAMMPS_POTENTIALS environment variable is set, add it to PYTHONPATH as well
@ -353,3 +353,83 @@ double PairPython::init_one(int, int)
return cut_global;
}
/* ---------------------------------------------------------------------- */
double PairPython::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
// prepare access to compute_force and compute_energy functions
PyGILState_STATE gstate = PyGILState_Ensure();
PyObject *py_pair_instance = (PyObject *) py_potential;
PyObject *py_compute_force
= PyObject_GetAttrString(py_pair_instance,"compute_force");
if (!py_compute_force) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Could not find 'compute_force' method'");
}
if (!PyCallable_Check(py_compute_force)) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Python 'compute_force' is not callable");
}
PyObject *py_compute_energy
= PyObject_GetAttrString(py_pair_instance,"compute_energy");
if (!py_compute_energy) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Could not find 'compute_energy' method'");
}
if (!PyCallable_Check(py_compute_energy)) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Python 'compute_energy' is not callable");
}
PyObject *py_rsq, *py_itype, *py_jtype, *py_value;
PyObject *py_compute_args = PyTuple_New(3);
if (!py_compute_args) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Could not create tuple for 'compute' function arguments");
}
py_itype = PY_INT_FROM_LONG(itype);
PyTuple_SetItem(py_compute_args,1,py_itype);
py_jtype = PY_INT_FROM_LONG(jtype);
PyTuple_SetItem(py_compute_args,2,py_jtype);
py_rsq = PyFloat_FromDouble(rsq);
PyTuple_SetItem(py_compute_args,0,py_rsq);
py_value = PyObject_CallObject(py_compute_force,py_compute_args);
if (!py_value) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Calling 'compute_force' function failed");
}
fforce = factor_lj*PyFloat_AsDouble(py_value)/rsq;
py_value = PyObject_CallObject(py_compute_energy,py_compute_args);
if (!py_value) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Calling 'compute_energy' function failed");
}
double evdwl = factor_lj*PyFloat_AsDouble(py_value);
Py_DECREF(py_compute_args);
PyGILState_Release(gstate);
return evdwl;
}

View File

@ -42,6 +42,7 @@ class PairPython : public Pair {
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
double single(int, int, int, int, double, double, double, double &);
protected:
double cut_global;