jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #524 from martok/package-meamc 

Package USER-MEAMC
This commit is contained in:
sjplimp 2017-07-03 12:30:01 -06:00 committed by GitHub
parent 5c0c8bb4cd 22f3db4723
commit 6fc0a94e87
30 changed files with 4629 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Section_commands.txt
View File

@ -1039,6 +1039,7 @@ package"_Section_start.html#start_3.
"lj/sdk (gko)"_pair_sdk.html,
"lj/sdk/coul/long (go)"_pair_sdk.html,
"lj/sdk/coul/msm (o)"_pair_sdk.html,
"meam/c"_pair_meam.html,
"meam/spline (o)"_pair_meam_spline.html,
"meam/sw/spline"_pair_meam_sw_spline.html,
"mgpt"_pair_mgpt.html,

31
doc/src/Section_packages.txt
View File

@ -121,6 +121,7 @@ Package, Description, Doc page, Example, Library
"USER-INTEL"_#USER-INTEL, optimized Intel CPU and KNL styles,"Section 5.3.2"_accelerate_intel.html, WWW bench, -
"USER-LB"_#USER-LB, Lattice Boltzmann fluid,"fix lb/fluid"_fix_lb_fluid.html, USER/lb, -
"USER-MANIFOLD"_#USER-MANIFOLD, motion on 2d surfaces,"fix manifoldforce"_fix_manifoldforce.html, USER/manifold, -
"USER-MEAMC"_#USER-MEAMC, modified EAM potential (C++), "pair_style meam/c"_pair_meam.html, meam, -
"USER-MGPT"_#USER-MGPT, fast MGPT multi-ion potentials, "pair_style mgpt"_pair_mgpt.html, USER/mgpt, -
"USER-MISC"_#USER-MISC, single-file contributions, USER-MISC/README, USER/misc, -
"USER-MOLFILE"_#USER-MOLFILE, "VMD"_vmd_home molfile plug-ins,"dump molfile"_dump_molfile.html, -, ext
@ -2051,6 +2052,36 @@ http://lammps.sandia.gov/movies.html#manifold :ul
:line
USER-MEAMC package :link(USER-MEAMC),h4
[Contents:]
A pair style for the modified embedded atom (MEAM) potential
translated from the Fortran version in the "MEAM"_MEAM package
to plain C++. In contrast to the MEAM package, no library
needs to be compiled and the pair style can be instantiated
multiple times.
[Author:] Sebastian Huetter, (Otto-von-Guericke University Magdeburg)
based on the work of Greg Wagner (Northwestern U) while at Sandia.
[Install or un-install:]
make yes-user-meamc
make machine :pre
make no-user-meamc
make machine :pre
[Supporting info:]
src/USER-MEAMC: filenames -> commands
src/USER-MEAMC/README
"pair meam/c"_pair_meam.html
examples/meam :ul
:line
USER-MOLFILE package :link(USER-MOLFILE),h4
[Contents:]

18
doc/src/pair_meam.txt
View File

@ -7,10 +7,13 @@
:line
pair_style meam command :h3
pair_style meam/c command :h3
[Syntax:]
pair_style meam :pre
pair_style style :pre
style = {meam} or {meam/c}
[Examples:]
@ -30,7 +33,8 @@ using modified embedded-atom method (MEAM) potentials
"EAM potentials"_pair_eam.html which adds angular forces. It is
thus suitable for modeling metals and alloys with fcc, bcc, hcp and
diamond cubic structures, as well as covalently bonded materials like
silicon and carbon.
silicon and carbon. Style {meam/c} is a translation of the {meam} code
from (mostly) Fortran to C++. It is functionally equivalent to {meam}.
In the MEAM formulation, the total energy E of a system of atoms is
given by:
@ -331,10 +335,14 @@ This pair style can only be used via the {pair} keyword of the
[Restrictions:]
This style is part of the MEAM package. It is only enabled if LAMMPS
The {meam} style is part of the MEAM package. It is only enabled if LAMMPS
was built with that package, which also requires the MEAM library be
built and linked with LAMMPS. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
built and linked with LAMMPS.
The {meam/c} style is provided in the USER-MEAMC package. It is only enabled
if LAMMPS was built with that package. In contrast to the {meam} style,
{meam/c} does not require a separate library to be compiled and it can be
instantiated multiple times in a "hybrid"_pair_hybrid.html pair style.
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
[Related commands:]

30
examples/meam/in.meamc Normal file
View File

@ -0,0 +1,30 @@
# Test of MEAM potential for SiC system
units metal
boundary p p p
atom_style atomic
read_data data.meam
pair_style meam/c
pair_coeff * * library.meam Si C SiC.meam Si C
neighbor 0.3 bin
neigh_modify delay 10
fix 1 all nve
thermo 10
timestep 0.001
#dump 1 all atom 50 dump.meam
#dump 2 all image 10 image.*.jpg element element &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3 element Si C
#dump 3 all movie 10 movie.mpg element element &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3 element Si C
run 100

79
examples/meam/in.meamc.shear Normal file
View File

@ -0,0 +1,79 @@
# 3d metal shear simulation
units metal
boundary s s p
atom_style atomic
lattice fcc 3.52
region box block 0 16.0 0 10.0 0 2.828427
create_box 3 box
lattice fcc 3.52 orient x 1 0 0 orient y 0 1 1 orient z 0 -1 1 &
origin 0.5 0 0
create_atoms 1 box
pair_style meam/c
pair_coeff * * library.meam Ni4 Ni.meam Ni4 Ni4 Ni4
neighbor 0.3 bin
neigh_modify delay 5
region lower block INF INF INF 0.9 INF INF
region upper block INF INF 6.1 INF INF INF
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary
set group lower type 2
set group upper type 3
# void
#region void cylinder z 8 5 2.5 INF INF
#delete_atoms region void
# temp controllers
compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1
# equilibrate
velocity mobile create 300.0 5812775 temp new3d
fix 1 all nve
fix 2 boundary setforce 0.0 0.0 0.0
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new3d
thermo 25
thermo_modify temp new3d
timestep 0.001
run 100
# shear
velocity upper set 1.0 0 0
velocity mobile ramp vx 0.0 1.0 y 1.4 8.6 sum yes
unfix 3
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d
#dump 1 all atom 500 dump.meam.shear
#dump 2 all image 100 image.*.jpg type type &
# axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 2 pad 4
#dump 3 all movie 100 movie.mpg type type &
# axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 3 pad 4
thermo 100
thermo_modify temp new2d
reset_timestep 0
run 3000

37
examples/meam/log.5Oct16.meam.icc.1 → examples/meam/log.19May17.meam.g++.1
View File

@ -1,4 +1,5 @@
LAMMPS (5 Oct 2016)
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# Test of MEAM potential for SiC system
units metal
@ -34,13 +35,23 @@ timestep 0.001
run 100
Neighbor list info ...
2 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15 -> bins = 6 6 6
Memory usage per processor = 7.39054 Mbytes
binsize = 2.15, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.103 | 8.103 | 8.103 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -636.38121 0 -636.38121 -76571.819
10 1807.8862 -666.21959 0 -636.54126 -150571.49
@ -53,20 +64,20 @@ Step Temp E_pair E_mol TotEng Press
80 2126.0903 -671.43755 0 -636.53557 -97475.831
90 2065.755 -670.4349 0 -636.52338 -95858.837
100 2051.4553 -670.20799 0 -636.53122 -107068.9
Loop time of 0.094512 on 1 procs for 100 steps with 128 atoms
Loop time of 0.0864418 on 1 procs for 100 steps with 128 atoms
Performance: 91.417 ns/day, 0.263 hours/ns, 1058.067 timesteps/s
99.4% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 99.952 ns/day, 0.240 hours/ns, 1156.848 timesteps/s
99.6% 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.091268 | 0.091268 | 0.091268 | 0.0 | 96.57
Neigh | 0.0021861 | 0.0021861 | 0.0021861 | 0.0 | 2.31
Comm | 0.00059438 | 0.00059438 | 0.00059438 | 0.0 | 0.63
Output | 9.0837e-05 | 9.0837e-05 | 9.0837e-05 | 0.0 | 0.10
Modify | 0.00024438 | 0.00024438 | 0.00024438 | 0.0 | 0.26
Other | | 0.000128 | | | 0.14
Pair | 0.082592 | 0.082592 | 0.082592 | 0.0 | 95.55
Neigh | 0.0028124 | 0.0028124 | 0.0028124 | 0.0 | 3.25
Comm | 0.00049043 | 0.00049043 | 0.00049043 | 0.0 | 0.57
Output | 0.00014329 | 0.00014329 | 0.00014329 | 0.0 | 0.17
Modify | 0.00025415 | 0.00025415 | 0.00025415 | 0.0 | 0.29
Other | | 0.0001497 | | | 0.17
Nlocal: 128 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0

37
examples/meam/log.5Oct16.meam.icc.4 → examples/meam/log.19May17.meam.g++.4
View File

@ -1,4 +1,5 @@
LAMMPS (5 Oct 2016)
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# Test of MEAM potential for SiC system
units metal
@ -34,13 +35,23 @@ timestep 0.001
run 100
Neighbor list info ...
2 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15 -> bins = 6 6 6
Memory usage per processor = 7.319 Mbytes
binsize = 2.15, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.024 | 8.026 | 8.027 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -636.38121 0 -636.38121 -76571.819
10 1807.8862 -666.21959 0 -636.54126 -150571.49
@ -53,20 +64,20 @@ Step Temp E_pair E_mol TotEng Press
80 2126.0903 -671.43755 0 -636.53557 -97475.831
90 2065.755 -670.4349 0 -636.52338 -95858.837
100 2051.4553 -670.20799 0 -636.53122 -107068.9
Loop time of 0.0350628 on 4 procs for 100 steps with 128 atoms
Loop time of 0.0389078 on 4 procs for 100 steps with 128 atoms
Performance: 246.415 ns/day, 0.097 hours/ns, 2852.026 timesteps/s
98.4% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 222.063 ns/day, 0.108 hours/ns, 2570.177 timesteps/s
99.4% 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.030952 | 0.031776 | 0.032203 | 0.3 | 90.63
Neigh | 0.00058937 | 0.00061423 | 0.00063896 | 0.1 | 1.75
Comm | 0.0018125 | 0.0022421 | 0.0030777 | 1.1 | 6.39
Output | 0.00018525 | 0.00019765 | 0.00021911 | 0.1 | 0.56
Modify | 8.0585e-05 | 9.0539e-05 | 9.7752e-05 | 0.1 | 0.26
Other | | 0.0001422 | | | 0.41
Pair | 0.031958 | 0.033267 | 0.034691 | 0.6 | 85.50
Neigh | 0.00079155 | 0.00086409 | 0.00098801 | 0.0 | 2.22
Comm | 0.0025704 | 0.0041765 | 0.005573 | 1.9 | 10.73
Output | 0.0002749 | 0.00029588 | 0.00033569 | 0.0 | 0.76
Modify | 9.4891e-05 | 0.00010347 | 0.00011587 | 0.0 | 0.27
Other | | 0.000201 | | | 0.52
Nlocal: 32 ave 36 max 30 min
Histogram: 1 2 0 0 0 0 0 0 0 1

115
examples/meam/log.5Oct16.meam.shear.icc.1 → examples/meam/log.19May17.meam.shear.g++.1
View File

@ -1,4 +1,5 @@
LAMMPS (5 Oct 2016)
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d metal shear simulation
units metal
@ -62,38 +63,48 @@ fix_modify 3 temp new3d
thermo 25
thermo_modify temp new3d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:474)
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
timestep 0.001
run 100
Neighbor list info ...
2 neighbor list requests
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15 -> bins = 27 17 5
Memory usage per processor = 8.55725 Mbytes
binsize = 2.15, bins = 27 17 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.788 | 9.788 | 9.788 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -8232.7767 0 -8179.1466 1386.6643 19547.02
25 222.78953 -8188.1215 0 -8148.2941 9095.9008 19547.02
50 300 -8149.7654 0 -8096.1353 10633.141 19684.382
75 304.80657 -8163.4557 0 -8108.9665 7045.457 19759.745
100 300 -8173.6884 0 -8120.0584 5952.521 19886.589
Loop time of 1.72323 on 1 procs for 100 steps with 1912 atoms
Loop time of 1.58103 on 1 procs for 100 steps with 1912 atoms
Performance: 5.014 ns/day, 4.787 hours/ns, 58.031 timesteps/s
99.8% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 5.465 ns/day, 4.392 hours/ns, 63.250 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.7026 | 1.7026 | 1.7026 | 0.0 | 98.80
Neigh | 0.014496 | 0.014496 | 0.014496 | 0.0 | 0.84
Comm | 0.0015783 | 0.0015783 | 0.0015783 | 0.0 | 0.09
Output | 6.0081e-05 | 6.0081e-05 | 6.0081e-05 | 0.0 | 0.00
Modify | 0.0034628 | 0.0034628 | 0.0034628 | 0.0 | 0.20
Other | | 0.00101 | | | 0.06
Pair | 1.5561 | 1.5561 | 1.5561 | 0.0 | 98.42
Neigh | 0.018544 | 0.018544 | 0.018544 | 0.0 | 1.17
Comm | 0.0013864 | 0.0013864 | 0.0013864 | 0.0 | 0.09
Output | 0.00011396 | 0.00011396 | 0.00011396 | 0.0 | 0.01
Modify | 0.0038245 | 0.0038245 | 0.0038245 | 0.0 | 0.24
Other | | 0.001096 | | | 0.07
Nlocal: 1912 ave 1912 max 1912 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -128,11 +139,11 @@ fix_modify 3 temp new2d
thermo 100
thermo_modify temp new2d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:474)
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
reset_timestep 0
run 3000
Memory usage per processor = 8.73384 Mbytes
Per MPI rank memory allocation (min/avg/max) = 9.964 | 9.964 | 9.964 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300.39988 -8173.6884 0 -8137.8874 4992.9811 19894.297
100 292.06374 -8177.7096 0 -8142.9021 2568.3762 19871.53
@ -144,53 +155,53 @@ Step Temp E_pair E_mol TotEng Press Volume
700 300 -8148.1328 0 -8112.3794 3506.9234 20435.302
800 300 -8139.1821 0 -8103.4288 3628.3957 20509.519
900 305.03425 -8126.7734 0 -8090.4201 5316.2206 20638.992
1000 304.00321 -8112.1616 0 -8075.9311 7441.9639 20767.243
1100 304.14051 -8096.5041 0 -8060.2573 9646.698 20888.167
1200 302.78461 -8080.5931 0 -8044.5079 11516.21 20995.917
1300 308.67046 -8061.6724 0 -8024.8857 11496.487 21130.013
1400 309.83019 -8046.2701 0 -8009.3452 12926.847 21247.271
1500 300 -8035.0322 0 -7999.2789 15346.188 21370.637
1600 300 -8030.6678 0 -7994.9144 14802.342 21496.446
1700 300 -8024.5988 0 -7988.8454 13177.445 21611.262
1800 300 -8023.045 0 -7987.2916 10240.041 21740.735
1900 300 -8028.2797 0 -7992.5263 6912.1441 21866.544
2000 300 -8036.4487 0 -8000.6953 3561.8365 21977.695
2100 300 -8037.8249 0 -8002.0715 2879.2618 22109.611
2200 300 -8033.6682 0 -7997.9148 4936.3695 22224.427
2300 304.49349 -8033.4561 0 -7997.1673 5593.0915 22356.343
2400 300 -8033.2969 0 -7997.5436 7537.0891 22473.601
2500 300 -8033.1874 0 -7997.4341 11476.447 22598.189
2600 307.77395 -8026.9234 0 -7990.2436 15758.81 22720.333
2700 300 -8021.1736 0 -7985.4203 17948.896 22832.706
2800 300 -8017.0863 0 -7981.3329 17154.618 22957.293
2900 300 -8012.0514 0 -7976.298 13224.292 23089.209
3000 304.58031 -8008.1654 0 -7971.8661 8572.9227 23211.354
Loop time of 55.136 on 1 procs for 3000 steps with 1912 atoms
1000 304.00321 -8112.1616 0 -8075.9311 7441.9638 20767.243
1100 304.14047 -8096.5041 0 -8060.2573 9646.6976 20888.167
1200 302.78457 -8080.5931 0 -8044.5079 11516.209 20995.917
1300 308.67054 -8061.6724 0 -8024.8857 11496.479 21130.013
1400 309.8301 -8046.27 0 -8009.3452 12926.835 21247.271
1500 300 -8035.0321 0 -7999.2788 15346.455 21370.637
1600 300 -8030.6657 0 -7994.9123 14802.869 21496.446
1700 300 -8024.5251 0 -7988.7718 13176.923 21611.262
1800 300 -8022.9963 0 -7987.243 10260.665 21741.956
1900 300 -8028.0522 0 -7992.2988 6955.1367 21867.765
2000 300 -8037.2708 0 -8001.5175 3520.3749 21987.467
2100 300 -8035.9945 0 -8000.2412 3237.6121 22109.611
2200 300 -8036.1882 0 -8000.4348 4295.1386 22223.205
2300 300 -8032.7689 0 -7997.0155 6231.2346 22355.121
2400 300 -8034.6621 0 -7998.9088 8585.3799 22468.716
2500 300 -8031.7774 0 -7996.0241 11627.871 22587.196
2600 300 -8024.032 0 -7988.2786 16254.69 22716.669
2700 308.64215 -8017.407 0 -7980.6237 18440.226 22835.149
2800 300 -8015.7348 0 -7979.9814 17601.716 22957.293
2900 305.82558 -8013.7448 0 -7977.2972 14123.494 23089.209
3000 300 -8011.0289 0 -7975.2755 9957.2884 23205.247
Loop time of 51.9315 on 1 procs for 3000 steps with 1912 atoms
Performance: 4.701 ns/day, 5.105 hours/ns, 54.411 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 4.991 ns/day, 4.808 hours/ns, 57.768 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 | 54.317 | 54.317 | 54.317 | -nan | 98.51
Neigh | 0.63189 | 0.63189 | 0.63189 | 0.0 | 1.15
Comm | 0.051245 | 0.051245 | 0.051245 | 0.0 | 0.09
Output | 0.0005548 | 0.0005548 | 0.0005548 | 0.0 | 0.00
Modify | 0.10452 | 0.10452 | 0.10452 | 0.0 | 0.19
Other | | 0.03128 | | | 0.06
Pair | 50.918 | 50.918 | 50.918 | 0.0 | 98.05
Neigh | 0.81033 | 0.81033 | 0.81033 | 0.0 | 1.56
Comm | 0.047331 | 0.047331 | 0.047331 | 0.0 | 0.09
Output | 0.0011976 | 0.0011976 | 0.0011976 | 0.0 | 0.00
Modify | 0.11889 | 0.11889 | 0.11889 | 0.0 | 0.23
Other | | 0.03606 | | | 0.07
Nlocal: 1912 ave 1912 max 1912 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1667 ave 1667 max 1667 min
Nghost: 1672 ave 1672 max 1672 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 23365 ave 23365 max 23365 min
Neighs: 23557 ave 23557 max 23557 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 46730 ave 46730 max 46730 min
FullNghs: 47114 ave 47114 max 47114 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 46730
Ave neighs/atom = 24.4404
Total # of neighbors = 47114
Ave neighs/atom = 24.6412
Neighbor list builds = 221
Dangerous builds = 0
Total wall time: 0:00:56
Total wall time: 0:00:53

123
examples/meam/log.5Oct16.meam.shear.icc.4 → examples/meam/log.19May17.meam.shear.g++.4
View File

@ -1,4 +1,5 @@
LAMMPS (5 Oct 2016)
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d metal shear simulation
units metal
@ -62,38 +63,48 @@ fix_modify 3 temp new3d
thermo 25
thermo_modify temp new3d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:474)
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
timestep 0.001
run 100
Neighbor list info ...
2 neighbor list requests
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15 -> bins = 27 17 5
Memory usage per processor = 7.74146 Mbytes
binsize = 2.15, bins = 27 17 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.954 | 8.957 | 8.959 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -8232.7767 0 -8179.1466 1386.6643 19547.02
25 221.59546 -8187.6813 0 -8148.0673 9100.4509 19547.02
50 300 -8150.0685 0 -8096.4384 10317.407 19685.743
75 307.76021 -8164.6669 0 -8109.6496 6289.7138 19757.814
100 300 -8176.5141 0 -8122.884 4162.2559 19873.327
Loop time of 0.469502 on 4 procs for 100 steps with 1912 atoms
Loop time of 0.482293 on 4 procs for 100 steps with 1912 atoms
Performance: 18.402 ns/day, 1.304 hours/ns, 212.992 timesteps/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 17.914 ns/day, 1.340 hours/ns, 207.343 timesteps/s
98.7% 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.44052 | 0.45213 | 0.45813 | 1.0 | 96.30
Neigh | 0.0036478 | 0.0037832 | 0.003854 | 0.1 | 0.81
Comm | 0.0055377 | 0.011533 | 0.02316 | 6.5 | 2.46
Output | 9.0837e-05 | 9.8228e-05 | 0.00011325 | 0.1 | 0.02
Modify | 0.00098062 | 0.0010158 | 0.0010564 | 0.1 | 0.22
Other | | 0.0009408 | | | 0.20
Pair | 0.44374 | 0.45604 | 0.46922 | 1.4 | 94.56
Neigh | 0.0047338 | 0.0049097 | 0.0051899 | 0.2 | 1.02
Comm | 0.0054841 | 0.019044 | 0.031472 | 6.9 | 3.95
Output | 0.00012755 | 0.00013644 | 0.00015831 | 0.0 | 0.03
Modify | 0.0011139 | 0.0011852 | 0.0012643 | 0.2 | 0.25
Other | | 0.0009753 | | | 0.20
Nlocal: 478 ave 492 max 465 min
Histogram: 2 0 0 0 0 0 0 0 1 1
@ -128,11 +139,11 @@ fix_modify 3 temp new2d
thermo 100
thermo_modify temp new2d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:474)
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
reset_timestep 0
run 3000
Memory usage per processor = 7.78572 Mbytes
Per MPI rank memory allocation (min/avg/max) = 8.999 | 9.002 | 9.005 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 295.32113 -8176.5141 0 -8141.3183 3169.3113 19886.93
100 292.00251 -8176.5358 0 -8141.7356 -825.04802 19918.765
@ -144,53 +155,53 @@ Step Temp E_pair E_mol TotEng Press Volume
700 296.44479 -8149.7914 0 -8114.4618 1985.4155 20421.046
800 307.75738 -8139.1649 0 -8102.487 4319.078 20513.183
900 304.61422 -8126.9246 0 -8090.6214 6654.0963 20640.213
1000 300 -8113.8464 0 -8078.0931 7760.1242 20768.465
1100 300.17874 -8097.7469 0 -8061.9722 8438.1263 20874.731
1200 306.01444 -8083.3367 0 -8046.8665 10835.585 20994.432
1000 300 -8113.8464 0 -8078.0931 7760.1239 20768.465
1100 300.17873 -8097.7469 0 -8061.9722 8438.126 20874.731
1200 306.01441 -8083.3367 0 -8046.8665 10835.586 20994.432
1300 300 -8067.022 0 -8031.2686 11216.067 21126.348
1400 300 -8053.223 0 -8017.4697 10570.21 21253.378
1500 300 -8043.4848 0 -8007.7314 11360.829 21375.523
1600 300 -8034.6216 0 -7998.8683 11371.642 21498.889
1700 300 -8028.6774 0 -7992.924 9595.8772 21613.705
1800 300 -8033.0808 0 -7997.3274 8767.6261 21743.178
1900 303.30302 -8035.1958 0 -7999.0488 8059.5152 21859.215
2000 300 -8025.0857 0 -7989.3323 9308.9938 21980.138
2100 300 -8041.5796 0 -8005.8263 6656.0066 22108.39
2200 300 -8039.6315 0 -8003.8781 7532.9687 22226.87
2300 300 -8053.203 0 -8017.4497 8466.9094 22356.343
2400 300 -8050.9154 0 -8015.162 11784.136 22467.494
2500 300 -8037.6394 0 -8001.886 16464.786 22588.417
2600 300 -8030.9221 0 -7995.1688 16807.326 22719.112
2700 300 -8025.2102 0 -7989.4569 13729.61 22837.592
2800 300 -8014.5312 0 -7978.7779 6784.6283 22953.629
2900 300 -8007.4768 0 -7971.7234 1362.3131 23084.324
3000 300 -7994.5614 0 -7958.808 -1726.5273 23194.254
Loop time of 14.8108 on 4 procs for 3000 steps with 1912 atoms
1400 300 -8053.223 0 -8017.4697 10570.206 21253.378
1500 300 -8043.4849 0 -8007.7315 11360.766 21375.523
1600 300 -8034.621 0 -7998.8676 11371.584 21498.889
1700 300 -8028.6783 0 -7992.925 9596.524 21613.705
1800 300 -8033.0818 0 -7997.3285 8767.2651 21743.178
1900 303.18912 -8035.194 0 -7999.0606 8059.9558 21859.215
2000 300 -8025.0327 0 -7989.2794 9305.7521 21980.138
2100 300 -8041.4626 0 -8005.7092 6623.8789 22108.39
2200 300 -8040.3133 0 -8004.5599 7512.9368 22225.648
2300 300 -8055.6567 0 -8019.9033 8281.354 22344.128
2400 304.05922 -8050.289 0 -8014.0518 11964.826 22476.044
2500 305.75646 -8037.0481 0 -8000.6087 16594.032 22595.746
2600 307.71105 -8031.2253 0 -7994.5529 18381.745 22708.119
2700 307.397 -8026.5338 0 -7989.8988 13944.653 22829.042
2800 309.3455 -8020.2305 0 -7983.3634 7037.4046 22954.851
2900 301.2859 -8010.4731 0 -7974.5665 3843.8972 23072.109
3000 303.29908 -8000.0395 0 -7963.8929 364.90172 23207.69
Loop time of 14.5278 on 4 procs for 3000 steps with 1912 atoms
Performance: 17.501 ns/day, 1.371 hours/ns, 202.555 timesteps/s
99.8% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 17.842 ns/day, 1.345 hours/ns, 206.500 timesteps/s
99.4% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 14.05 | 14.237 | 14.332 | 2.9 | 96.12
Neigh | 0.1592 | 0.16414 | 0.1671 | 0.8 | 1.11
Comm | 0.26002 | 0.35589 | 0.54696 | 18.8 | 2.40
Output | 0.00061679 | 0.00065172 | 0.0007441 | 0.2 | 0.00
Modify | 0.02895 | 0.030174 | 0.03104 | 0.5 | 0.20
Other | | 0.02338 | | | 0.16
Pair | 13.872 | 13.929 | 13.998 | 1.4 | 95.88
Neigh | 0.20891 | 0.21114 | 0.21272 | 0.3 | 1.45
Comm | 0.25364 | 0.32377 | 0.37706 | 8.9 | 2.23
Output | 0.0011427 | 0.0012097 | 0.0013931 | 0.3 | 0.01
Modify | 0.033687 | 0.033991 | 0.034694 | 0.2 | 0.23
Other | | 0.02871 | | | 0.20
Nlocal: 478 ave 509 max 448 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 799.25 ave 844 max 756 min
Histogram: 1 1 0 0 0 0 0 1 0 1
Neighs: 5813.25 ave 6081 max 5602 min
Histogram: 2 0 0 0 0 0 1 0 0 1
FullNghs: 11626.5 ave 12151 max 11205 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Nlocal: 478 ave 509 max 445 min
Histogram: 1 1 0 0 0 0 0 0 1 1
Nghost: 804 ave 845 max 759 min
Histogram: 1 1 0 0 0 0 0 0 1 1
Neighs: 5827 ave 6177 max 5496 min
Histogram: 1 0 0 1 0 1 0 0 0 1
FullNghs: 11654 ave 12330 max 11039 min
Histogram: 1 0 0 1 0 1 0 0 0 1
Total # of neighbors = 46506
Ave neighs/atom = 24.3232
Neighbor list builds = 225
Total # of neighbors = 46616
Ave neighs/atom = 24.3808
Neighbor list builds = 223
Dangerous builds = 0
Total wall time: 0:00:15

95
examples/meam/log.19May17.meamc.g++.1 Normal file
View File

@ -0,0 +1,95 @@
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# Test of MEAM potential for SiC system
units metal
boundary p p p
atom_style atomic
read_data data.meam
orthogonal box = (-6 -6 -6) to (5.97232 5.97232 5.97232)
1 by 1 by 1 MPI processor grid
reading atoms ...
128 atoms
pair_style meam/c
pair_coeff * * library.meam Si C SiC.meam Si C
Reading potential file library.meam with DATE: 2012-06-29
Reading potential file SiC.meam with DATE: 2007-06-11
neighbor 0.3 bin
neigh_modify delay 10
fix 1 all nve
thermo 10
timestep 0.001
#dump 1 all atom 50 dump.meam
#dump 2 all image 10 image.*.jpg element element # axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3 element Si C
#dump 3 all movie 10 movie.mpg element element # axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3 element Si C
run 100
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/c, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/c, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.103 | 8.103 | 8.103 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -636.38121 0 -636.38121 -76571.819
10 1807.8862 -666.21959 0 -636.54126 -150571.49
20 1932.4467 -668.2581 0 -636.53498 -120223.52
30 1951.3652 -668.58139 0 -636.54771 -100508.4
40 2172.5974 -672.22715 0 -636.5617 -110753.34
50 2056.9149 -670.33108 0 -636.56468 -105418.07
60 1947.9564 -668.52788 0 -636.55015 -111413.04
70 1994.7712 -669.28849 0 -636.54225 -109645.76
80 2126.0903 -671.43755 0 -636.53557 -97475.832
90 2065.7549 -670.4349 0 -636.52338 -95858.836
100 2051.4553 -670.20799 0 -636.53122 -107068.89
Loop time of 0.0778153 on 1 procs for 100 steps with 128 atoms
Performance: 111.032 ns/day, 0.216 hours/ns, 1285.094 timesteps/s
99.5% 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.073801 | 0.073801 | 0.073801 | 0.0 | 94.84
Neigh | 0.0029731 | 0.0029731 | 0.0029731 | 0.0 | 3.82
Comm | 0.00047708 | 0.00047708 | 0.00047708 | 0.0 | 0.61
Output | 0.00015664 | 0.00015664 | 0.00015664 | 0.0 | 0.20
Modify | 0.00025702 | 0.00025702 | 0.00025702 | 0.0 | 0.33
Other | | 0.0001504 | | | 0.19
Nlocal: 128 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 543 ave 543 max 543 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1526 ave 1526 max 1526 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 3052 ave 3052 max 3052 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 3052
Ave neighs/atom = 23.8438
Neighbor list builds = 10
Dangerous builds = 10
Total wall time: 0:00:00

95
examples/meam/log.19May17.meamc.g++.4 Normal file
View File

@ -0,0 +1,95 @@
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# Test of MEAM potential for SiC system
units metal
boundary p p p
atom_style atomic
read_data data.meam
orthogonal box = (-6 -6 -6) to (5.97232 5.97232 5.97232)
1 by 2 by 2 MPI processor grid
reading atoms ...
128 atoms
pair_style meam/c
pair_coeff * * library.meam Si C SiC.meam Si C
Reading potential file library.meam with DATE: 2012-06-29
Reading potential file SiC.meam with DATE: 2007-06-11
neighbor 0.3 bin
neigh_modify delay 10
fix 1 all nve
thermo 10
timestep 0.001
#dump 1 all atom 50 dump.meam
#dump 2 all image 10 image.*.jpg element element # axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3 element Si C
#dump 3 all movie 10 movie.mpg element element # axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3 element Si C
run 100
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/c, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/c, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.024 | 8.026 | 8.027 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -636.38121 0 -636.38121 -76571.819
10 1807.8862 -666.21959 0 -636.54126 -150571.49
20 1932.4467 -668.2581 0 -636.53498 -120223.52
30 1951.3652 -668.58139 0 -636.54771 -100508.4
40 2172.5974 -672.22715 0 -636.5617 -110753.34
50 2056.9149 -670.33108 0 -636.56468 -105418.07
60 1947.9564 -668.52788 0 -636.55015 -111413.04
70 1994.7712 -669.28849 0 -636.54225 -109645.76
80 2126.0903 -671.43755 0 -636.53557 -97475.832
90 2065.7549 -670.4349 0 -636.52338 -95858.836
100 2051.4553 -670.20799 0 -636.53122 -107068.89
Loop time of 0.037066 on 4 procs for 100 steps with 128 atoms
Performance: 233.097 ns/day, 0.103 hours/ns, 2697.887 timesteps/s
97.4% 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.029985 | 0.031596 | 0.033021 | 0.8 | 85.24
Neigh | 0.0007906 | 0.00085384 | 0.00088596 | 0.0 | 2.30
Comm | 0.0025654 | 0.0040313 | 0.0057514 | 2.2 | 10.88
Output | 0.00027013 | 0.00029153 | 0.00033426 | 0.0 | 0.79
Modify | 9.5367e-05 | 0.00010639 | 0.00012016 | 0.0 | 0.29
Other | | 0.0001866 | | | 0.50
Nlocal: 32 ave 36 max 30 min
Histogram: 1 2 0 0 0 0 0 0 0 1
Nghost: 293.75 ave 305 max 285 min
Histogram: 2 0 0 0 0 0 0 1 0 1
Neighs: 381.5 ave 413 max 334 min
Histogram: 1 0 0 0 1 0 0 0 0 2
FullNghs: 763 ave 866 max 678 min
Histogram: 1 0 1 0 0 1 0 0 0 1
Total # of neighbors = 3052
Ave neighs/atom = 23.8438
Neighbor list builds = 10
Dangerous builds = 10
Total wall time: 0:00:00

207
examples/meam/log.19May17.meamc.shear.g++.1 Normal file
View File

@ -0,0 +1,207 @@
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d metal shear simulation
units metal
boundary s s p
atom_style atomic
lattice fcc 3.52
Lattice spacing in x,y,z = 3.52 3.52 3.52
region box block 0 16.0 0 10.0 0 2.828427
create_box 3 box
Created orthogonal box = (0 0 0) to (56.32 35.2 9.95606)
1 by 1 by 1 MPI processor grid
lattice fcc 3.52 orient x 1 0 0 orient y 0 1 1 orient z 0 -1 1 origin 0.5 0 0
Lattice spacing in x,y,z = 3.52 4.97803 4.97803
create_atoms 1 box
Created 1912 atoms
pair_style meam/c
pair_coeff * * library.meam Ni4 Ni.meam Ni4 Ni4 Ni4
Reading potential file library.meam with DATE: 2012-06-29
Reading potential file Ni.meam with DATE: 2007-06-11
neighbor 0.3 bin
neigh_modify delay 5
region lower block INF INF INF 0.9 INF INF
region upper block INF INF 6.1 INF INF INF
group lower region lower
264 atoms in group lower
group upper region upper
264 atoms in group upper
group boundary union lower upper
528 atoms in group boundary
group mobile subtract all boundary
1384 atoms in group mobile
set group lower type 2
264 settings made for type
set group upper type 3
264 settings made for type
# void
#region void cylinder z 8 5 2.5 INF INF
#delete_atoms region void
# temp controllers
compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1
# equilibrate
velocity mobile create 300.0 5812775 temp new3d
fix 1 all nve
fix 2 boundary setforce 0.0 0.0 0.0
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new3d
thermo 25
thermo_modify temp new3d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
timestep 0.001
run 100
Neighbor list info ...
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15, bins = 27 17 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/c, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/c, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.788 | 9.788 | 9.788 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -8232.7767 0 -8179.1466 1386.6643 19547.02
25 222.78953 -8188.1215 0 -8148.2941 9095.9003 19547.02
50 300 -8149.7654 0 -8096.1353 10633.139 19684.382
75 304.80657 -8163.4557 0 -8108.9665 7045.4555 19759.745
100 300 -8173.6884 0 -8120.0584 5952.5197 19886.589
Loop time of 1.46986 on 1 procs for 100 steps with 1912 atoms
Performance: 5.878 ns/day, 4.083 hours/ns, 68.034 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.445 | 1.445 | 1.445 | 0.0 | 98.31
Neigh | 0.018564 | 0.018564 | 0.018564 | 0.0 | 1.26
Comm | 0.0012956 | 0.0012956 | 0.0012956 | 0.0 | 0.09
Output | 0.00012612 | 0.00012612 | 0.00012612 | 0.0 | 0.01
Modify | 0.0038197 | 0.0038197 | 0.0038197 | 0.0 | 0.26
Other | | 0.001095 | | | 0.07
Nlocal: 1912 ave 1912 max 1912 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1672 ave 1672 max 1672 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 23806 ave 23806 max 23806 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 47612 ave 47612 max 47612 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 47612
Ave neighs/atom = 24.9017
Neighbor list builds = 5
Dangerous builds = 0
# shear
velocity upper set 1.0 0 0
velocity mobile ramp vx 0.0 1.0 y 1.4 8.6 sum yes
unfix 3
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d
#dump 1 all atom 500 dump.meam.shear
#dump 2 all image 100 image.*.jpg type type # axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 2 pad 4
#dump 3 all movie 100 movie.mpg type type # axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 3 pad 4
thermo 100
thermo_modify temp new2d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
reset_timestep 0
run 3000
Per MPI rank memory allocation (min/avg/max) = 9.964 | 9.964 | 9.964 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300.39988 -8173.6884 0 -8137.8874 4992.9799 19894.297
100 292.06374 -8177.7096 0 -8142.9021 2568.3756 19871.53
200 306.69894 -8177.1357 0 -8140.584 874.24118 20047.24
300 295.68216 -8172.9213 0 -8137.6825 -1049.0799 20091.759
400 308.99955 -8169.6355 0 -8132.8096 -1785.9554 20121.698
500 303.85688 -8163.9842 0 -8127.7712 -150.60892 20183.813
600 300 -8157.7627 0 -8122.0093 1492.8514 20279.887
700 300 -8148.1314 0 -8112.3781 3507.1949 20435.297
800 300 -8139.1805 0 -8103.4272 3628.5908 20509.519
900 305.03217 -8126.7741 0 -8090.421 5313.7881 20638.992
1000 303.7648 -8112.1574 0 -8075.9554 7433.3181 20767.243
1100 302.39719 -8096.1399 0 -8060.1009 9681.7685 20888.167
1200 304.04919 -8080.7022 0 -8044.4663 11621.974 21011.532
1300 303.56395 -8062.0984 0 -8025.9203 11410.793 21125.127
1400 309.92338 -8046.0008 0 -8009.0648 12408.158 21246.05
1500 300 -8034.7094 0 -7998.956 14845.312 21363.308
1600 300 -8028.4585 0 -7992.7051 15120.908 21489.117
1700 308.23904 -8015.9618 0 -7979.2265 14714.73 21612.483
1800 300 -8013.5458 0 -7977.7924 11955.065 21737.07
1900 300 -8012.2984 0 -7976.545 6667.1353 21854.329
2000 300 -8025.6019 0 -7989.8485 2006.6545 21981.359
2100 300 -8027.6823 0 -7991.9289 16.47633 22109.611
2200 300 -8029.6905 0 -7993.9372 -603.98293 22224.427
2300 300 -8033.2663 0 -7997.513 -464.68645 22351.457
2400 300 -8040.6863 0 -8004.9329 -640.54641 22467.494
2500 300 -8037.0332 0 -8001.2799 1504.2444 22587.196
2600 300 -8036.0909 0 -8000.3375 4190.2052 22708.119
2700 308.97892 -8028.5269 0 -7991.7035 5755.7418 22832.706
2800 305.27189 -8023.8286 0 -7987.4469 2611.7551 22962.179
2900 301.94251 -8017.4523 0 -7981.4675 358.11928 23078.216
3000 305.67682 -8009.853 0 -7973.4231 -2345.487 23197.918
Loop time of 48.351 on 1 procs for 3000 steps with 1912 atoms
Performance: 5.361 ns/day, 4.477 hours/ns, 62.046 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 | 47.356 | 47.356 | 47.356 | 0.0 | 97.94
Neigh | 0.79977 | 0.79977 | 0.79977 | 0.0 | 1.65
Comm | 0.043133 | 0.043133 | 0.043133 | 0.0 | 0.09
Output | 0.0011899 | 0.0011899 | 0.0011899 | 0.0 | 0.00
Modify | 0.11648 | 0.11648 | 0.11648 | 0.0 | 0.24
Other | | 0.03404 | | | 0.07
Nlocal: 1912 ave 1912 max 1912 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1662 ave 1662 max 1662 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 23143 ave 23143 max 23143 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 46286 ave 46286 max 46286 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 46286
Ave neighs/atom = 24.2082
Neighbor list builds = 220
Dangerous builds = 0
Total wall time: 0:00:49

207
examples/meam/log.19May17.meamc.shear.g++.4 Normal file
View File

@ -0,0 +1,207 @@
LAMMPS (19 May 2017)
using 1 OpenMP thread(s) per MPI task
# 3d metal shear simulation
units metal
boundary s s p
atom_style atomic
lattice fcc 3.52
Lattice spacing in x,y,z = 3.52 3.52 3.52
region box block 0 16.0 0 10.0 0 2.828427
create_box 3 box
Created orthogonal box = (0 0 0) to (56.32 35.2 9.95606)
2 by 2 by 1 MPI processor grid
lattice fcc 3.52 orient x 1 0 0 orient y 0 1 1 orient z 0 -1 1 origin 0.5 0 0
Lattice spacing in x,y,z = 3.52 4.97803 4.97803
create_atoms 1 box
Created 1912 atoms
pair_style meam/c
pair_coeff * * library.meam Ni4 Ni.meam Ni4 Ni4 Ni4
Reading potential file library.meam with DATE: 2012-06-29
Reading potential file Ni.meam with DATE: 2007-06-11
neighbor 0.3 bin
neigh_modify delay 5
region lower block INF INF INF 0.9 INF INF
region upper block INF INF 6.1 INF INF INF
group lower region lower
264 atoms in group lower
group upper region upper
264 atoms in group upper
group boundary union lower upper
528 atoms in group boundary
group mobile subtract all boundary
1384 atoms in group mobile
set group lower type 2
264 settings made for type
set group upper type 3
264 settings made for type
# void
#region void cylinder z 8 5 2.5 INF INF
#delete_atoms region void
# temp controllers
compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1
# equilibrate
velocity mobile create 300.0 5812775 temp new3d
fix 1 all nve
fix 2 boundary setforce 0.0 0.0 0.0
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new3d
thermo 25
thermo_modify temp new3d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
timestep 0.001
run 100
Neighbor list info ...
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15, bins = 27 17 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/c, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/c, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.954 | 8.957 | 8.959 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -8232.7767 0 -8179.1466 1386.6643 19547.02
25 221.59546 -8187.6813 0 -8148.0673 9100.4505 19547.02
50 300 -8150.0685 0 -8096.4384 10317.406 19685.743
75 307.76021 -8164.6669 0 -8109.6496 6289.7123 19757.814
100 300 -8176.5141 0 -8122.884 4162.2548 19873.327
Loop time of 0.435874 on 4 procs for 100 steps with 1912 atoms
Performance: 19.822 ns/day, 1.211 hours/ns, 229.424 timesteps/s
98.8% 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.4092 | 0.41563 | 0.42184 | 0.7 | 95.35
Neigh | 0.0048575 | 0.004932 | 0.0049984 | 0.1 | 1.13
Comm | 0.0069079 | 0.013151 | 0.019513 | 4.2 | 3.02
Output | 0.00012398 | 0.00013405 | 0.00015688 | 0.0 | 0.03
Modify | 0.0011275 | 0.0011509 | 0.0011735 | 0.1 | 0.26
Other | | 0.0008795 | | | 0.20
Nlocal: 478 ave 492 max 465 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Nghost: 809 ave 822 max 795 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Neighs: 5916 ave 6133 max 5658 min
Histogram: 1 0 0 1 0 0 0 0 1 1
FullNghs: 11832 ave 12277 max 11299 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Total # of neighbors = 47328
Ave neighs/atom = 24.7531
Neighbor list builds = 5
Dangerous builds = 0
# shear
velocity upper set 1.0 0 0
velocity mobile ramp vx 0.0 1.0 y 1.4 8.6 sum yes
unfix 3
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d
#dump 1 all atom 500 dump.meam.shear
#dump 2 all image 100 image.*.jpg type type # axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 2 pad 4
#dump 3 all movie 100 movie.mpg type type # axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
#dump_modify 3 pad 4
thermo 100
thermo_modify temp new2d
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:489)
reset_timestep 0
run 3000
Per MPI rank memory allocation (min/avg/max) = 8.999 | 9.002 | 9.005 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 295.32113 -8176.5141 0 -8141.3183 3169.3102 19886.93
100 292.0025 -8176.5358 0 -8141.7356 -825.04852 19918.765
200 306.11682 -8176.7718 0 -8140.2895 -1370.6881 19948.878
300 300 -8172.6262 0 -8136.8729 -1735.9794 20085.715
400 306.88504 -8168.4351 0 -8131.8612 -933.05532 20117.008
500 308.99111 -8166.2909 0 -8129.466 -1049.3442 20198.267
600 304.22555 -8158.0946 0 -8121.8377 583.69142 20328.753
700 296.41606 -8149.7777 0 -8114.4515 1986.7982 20421.032
800 307.88628 -8139.1709 0 -8102.4776 4311.4142 20513.183
900 303.37209 -8126.9382 0 -8090.7829 6712.7316 20640.213
1000 300 -8113.7973 0 -8078.044 7630.2594 20750.143
1100 300.07815 -8098.1383 0 -8062.3756 8423.7063 20879.616
1200 300 -8083.3163 0 -8047.563 10772.917 21000.539
1300 300 -8066.6741 0 -8030.9208 10834.336 21128.791
1400 300 -8050.8799 0 -8015.1265 10842.382 21257.043
1500 300 -8040.3206 0 -8004.5673 11852.589 21362.087
1600 300 -8033.2471 0 -7997.4937 11298.745 21492.782
1700 300 -8030.6375 0 -7994.8842 10850.43 21610.04
1800 300 -8035.1631 0 -7999.4097 9985.6107 21734.628
1900 308.56562 -8040.1954 0 -8003.4213 9865.7107 21859.215
2000 300 -8030.1651 0 -7994.4117 11817.502 21980.138
2100 300 -8031.6147 0 -7995.8613 12791.357 22101.061
2200 300 -8033.7793 0 -7998.0259 13823.043 22234.198
2300 300 -8040.5964 0 -8004.8431 16204.549 22350.236
2400 309.42867 -8045.9414 0 -8009.0644 18506.386 22465.051
2500 300 -8054.2629 0 -8018.5095 20099.612 22593.303
2600 300 -8054.9562 0 -8019.2028 20036.318 22721.555
2700 300 -8051.4788 0 -8015.7254 16993.437 22844.921
2800 300 -8050.6908 0 -8014.9374 9048.5896 22964.622
2900 309.90783 -8044.3096 0 -8007.3754 5017.0198 23080.659
3000 300 -8035.8165 0 -8000.0632 2084.8999 23207.69
Loop time of 13.4901 on 4 procs for 3000 steps with 1912 atoms
Performance: 19.214 ns/day, 1.249 hours/ns, 222.386 timesteps/s
99.4% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 12.851 | 12.919 | 12.945 | 1.1 | 95.76
Neigh | 0.20449 | 0.20777 | 0.21169 | 0.6 | 1.54
Comm | 0.27479 | 0.30264 | 0.36667 | 6.8 | 2.24
Output | 0.0010171 | 0.0010971 | 0.0012388 | 0.3 | 0.01
Modify | 0.033248 | 0.033878 | 0.034567 | 0.3 | 0.25
Other | | 0.02594 | | | 0.19
Nlocal: 478 ave 506 max 450 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 790.5 ave 822 max 751 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Neighs: 5829.5 ave 6014 max 5672 min
Histogram: 2 0 0 0 0 0 0 0 1 1
FullNghs: 11659 ave 12027 max 11320 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Total # of neighbors = 46636
Ave neighs/atom = 24.3912
Neighbor list builds = 220
Dangerous builds = 0
Total wall time: 0:00:13

5
src/.gitignore vendored
View File

@ -31,6 +31,11 @@
/fix_*manifold*.cpp
/fix_*manifold*.h
/meam*.h
/meam*.cpp
/pair_meamc.cpp
/pair_meamc.h
/fix_qeq*.cpp
/fix_qeq*.h

2
src/Makefile
View File

@ -59,7 +59,7 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
PACKUSER = user-atc user-awpmd user-cgdna user-cgsdk user-colvars \
user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
user-intel user-lb user-manifold user-mgpt user-misc user-molfile \
user-intel user-lb user-manifold user-meamc user-mgpt user-misc user-molfile \
user-netcdf user-omp user-phonon user-qmmm user-qtb \
user-quip user-reaxc user-smd user-smtbq user-sph user-tally \
user-vtk

26
src/USER-MEAMC/README Normal file
View File

@ -0,0 +1,26 @@
This package implements the MEAM/C potential as a LAMMPS pair style.
==============================================================================
This package is a translation of the MEAM package to native C++. See
that package as well as the Fortran code distributed in lib/meam for
the original sources and information.
Translation by
Sebastian Hütter, sebastian.huetter@ovgu.de
Institute of Materials and Joining Technology
Otto-von-Guericke University Magdeburg, Germany
The original Fortran implementation was created by
Greg Wagner (while at Sandia, now at Northwestern U).
==============================================================================
Use "make yes-user-meamc" to enable this package when building LAMMPS.
In your LAMMPS input script, specifiy
pair_style meam/c
to enable the use of this implementation. All parameters, input files and
outputs are exactly identical to these used with pair_style meam.

252
src/USER-MEAMC/meam.h Normal file
View File

@ -0,0 +1,252 @@
#ifndef LMP_MEAM_H
#define LMP_MEAM_H
#include "memory.h"
#include <math.h>
#include <stdlib.h>
#define maxelt 5
namespace LAMMPS_NS {
typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;
class MEAM
{
public:
MEAM(Memory* mem);
~MEAM();
private:
Memory* memory;
// cutforce = force cutoff
// cutforcesq = force cutoff squared
double cutforce, cutforcesq;
// Ec_meam = cohesive energy
// re_meam = nearest-neighbor distance
// Omega_meam = atomic volume
// B_meam = bulk modulus
// Z_meam = number of first neighbors for reference structure
// ielt_meam = atomic number of element
// A_meam = adjustable parameter
// alpha_meam = sqrt(9*Omega*B/Ec)
// rho0_meam = density scaling parameter
// delta_meam = heat of formation for alloys
// beta[0-3]_meam = electron density constants
// t[0-3]_meam = coefficients on densities in Gamma computation
// rho_ref_meam = background density for reference structure
// ibar_meam(i) = selection parameter for Gamma function for elt i,
// lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
// neltypes = maximum number of element type defined
// eltind = index number of pair (similar to Voigt notation; ij = ji)
// phir = pair potential function array
// phirar[1-6] = spline coeffs
// attrac_meam = attraction parameter in Rose energy
// repuls_meam = repulsion parameter in Rose energy
// nn2_meam = 1 if second nearest neighbors are to be computed, else 0
// zbl_meam = 1 if zbl potential for small r to be use, else 0
// emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0
// bkgd_dyn = 1 if reference densities follows Dynamo, else 0
// Cmin_meam, Cmax_meam = min and max values in screening cutoff
// rc_meam = cutoff distance for meam
// delr_meam = cutoff region for meam
// ebound_meam = factor giving maximum boundary of sceen fcn ellipse
// augt1 = flag for whether t1 coefficient should be augmented
// ialloy = flag for newer alloy formulation (as in dynamo code)
// mix_ref_t = flag to recover "old" way of computing t in reference config
// erose_form = selection parameter for form of E_rose function
// gsmooth_factor = factor determining length of G smoothing region
// vind[23]D = Voight notation index maps for 2 and 3D
// v2D,v3D = array of factors to apply for Voight notation
// nr,dr = pair function discretization parameters
// nrar,rdrar = spline coeff array parameters
double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
double Omega_meam[maxelt], Z_meam[maxelt];
double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];
double delta_meam[maxelt][maxelt];
double beta0_meam[maxelt], beta1_meam[maxelt];
double beta2_meam[maxelt], beta3_meam[maxelt];
double t0_meam[maxelt], t1_meam[maxelt];
double t2_meam[maxelt], t3_meam[maxelt];
double rho_ref_meam[maxelt];
int ibar_meam[maxelt], ielt_meam[maxelt];
lattice_t lattce_meam[maxelt][maxelt];
int nn2_meam[maxelt][maxelt];
int zbl_meam[maxelt][maxelt];
int eltind[maxelt][maxelt];
int neltypes;
double** phir;
double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, **phirar6;
double attrac_meam[maxelt][maxelt], repuls_meam[maxelt][maxelt];
double Cmin_meam[maxelt][maxelt][maxelt];
double Cmax_meam[maxelt][maxelt][maxelt];
double rc_meam, delr_meam, ebound_meam[maxelt][maxelt];
int augt1, ialloy, mix_ref_t, erose_form;
int emb_lin_neg, bkgd_dyn;
double gsmooth_factor;
int vind2D[3][3], vind3D[3][3][3];
int v2D[6], v3D[10];
int nr, nrar;
double dr, rdrar;
public:
int nmax;
double *rho, *rho0, *rho1, *rho2, *rho3, *frhop;
double *gamma, *dgamma1, *dgamma2, *dgamma3, *arho2b;
double **arho1, **arho2, **arho3, **arho3b, **t_ave, **tsq_ave;
int maxneigh;
double *scrfcn, *dscrfcn, *fcpair;
protected:
// meam_funcs.cpp
//-----------------------------------------------------------------------------
// Cutoff function
//
static double fcut(const double xi) {
double a;
if (xi >= 1.0)
return 1.0;
else if (xi <= 0.0)
return 0.0;
else {
a = 1.0 - xi;
a *= a; a *= a;
a = 1.0 - a;
return a * a;
}
}
//-----------------------------------------------------------------------------
// Cutoff function and derivative
//
static double dfcut(const double xi, double& dfc) {
double a, a3, a4, a1m4;
if (xi >= 1.0) {
dfc = 0.0;
return 1.0;
} else if (xi <= 0.0) {
dfc = 0.0;
return 0.0;
} else {
a = 1.0 - xi;
a3 = a * a * a;
a4 = a * a3;
a1m4 = 1.0-a4;
dfc = 8 * a1m4 * a3;
return a1m4*a1m4;
}
}
//-----------------------------------------------------------------------------
// Derivative of Cikj w.r.t. rij
// Inputs: rij,rij2,rik2,rjk2
//
static double dCfunc(const double rij2, const double rik2, const double rjk2) {
double rij4, a, asq, b,denom;
rij4 = rij2 * rij2;
a = rik2 - rjk2;
b = rik2 + rjk2;
asq = a*a;
denom = rij4 - asq;
denom = denom * denom;
return -4 * (-2 * rij2 * asq + rij4 * b + asq * b) / denom;
}
//-----------------------------------------------------------------------------
// Derivative of Cikj w.r.t. rik and rjk
// Inputs: rij,rij2,rik2,rjk2
//
static void dCfunc2(const double rij2, const double rik2, const double rjk2,
double& dCikj1, double& dCikj2) {
double rij4, rik4, rjk4, a, denom;
rij4 = rij2 * rij2;
rik4 = rik2 * rik2;
rjk4 = rjk2 * rjk2;
a = rik2 - rjk2;
denom = rij4 - a * a;
denom = denom * denom;
dCikj1 = 4 * rij2 * (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;
dCikj2 = 4 * rij2 * (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;
}
double G_gam(const double gamma, const int ibar, int &errorflag) const;
double dG_gam(const double gamma, const int ibar, double &dG) const;
static double zbl(const double r, const int z1, const int z2);
static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);
static void get_shpfcn(const lattice_t latt, double (&s)[3]);
static int get_Zij(const lattice_t latt);
static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);
protected:
void meam_checkindex(int, int, int, int*, int*);
void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
int* firstneigh, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap);
void calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* scrfcn, double* fcpair);
void alloyparams();
void compute_pair_meam();
double phi_meam(double, int, int);
void compute_reference_density();
void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double,
double, double, double, int, int, lattice_t);
void get_sijk(double, int, int, int, double*);
void get_densref(double, int, int, double*, double*, double*, double*, double*, double*, double*, double*);
void interpolate_meam(int);
double compute_phi(double, int, int);
public:
void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar);
void meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag);
void meam_setup_done(double* cutmax);
void meam_dens_setup(int atom_nmax, int nall, int n_neigh);
void meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset);
void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, int& errorflag);
void meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom);
};
// Functions we need for compat
static inline bool iszero(const double f) {
return fabs(f) < 1e-20;
}
template <typename TYPE, size_t maxi, size_t maxj>
static inline void setall2d(TYPE (&arr)[maxi][maxj], const TYPE v) {
for (size_t i = 0; i < maxi; i++)
for (size_t j = 0; j < maxj; j++)
arr[i][j] = v;
}
template <typename TYPE, size_t maxi, size_t maxj, size_t maxk>
static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
for (size_t i = 0; i < maxi; i++)
for (size_t j = 0; j < maxj; j++)
for (size_t k = 0; k < maxk; k++)
arr[i][j][k] = v;
}
};
#endif

148
src/USER-MEAMC/meam_dens_final.cpp Normal file
View File

@ -0,0 +1,148 @@
#include "meam.h"
#include "math_special.h"
using namespace LAMMPS_NS;
void
MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, int& errorflag)
{
int i, elti;
int m;
double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;
double B, denom, rho_bkgd;
// Complete the calculation of density
for (i = 0; i < nlocal; i++) {
elti = fmap[type[i]];
if (elti >= 0) {
rho1[i] = 0.0;
rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];
rho3[i] = 0.0;
for (m = 0; m < 3; m++) {
rho1[i] = rho1[i] + arho1[i][m] * arho1[i][m];
rho3[i] = rho3[i] - 3.0 / 5.0 * arho3b[i][m] * arho3b[i][m];
}
for (m = 0; m < 6; m++) {
rho2[i] = rho2[i] + this->v2D[m] * arho2[i][m] * arho2[i][m];
}
for (m = 0; m < 10; m++) {
rho3[i] = rho3[i] + this->v3D[m] * arho3[i][m] * arho3[i][m];
}
if (rho0[i] > 0.0) {
if (this->ialloy == 1) {
t_ave[i][0] = t_ave[i][0] / tsq_ave[i][0];
t_ave[i][1] = t_ave[i][1] / tsq_ave[i][1];
t_ave[i][2] = t_ave[i][2] / tsq_ave[i][2];
} else if (this->ialloy == 2) {
t_ave[i][0] = this->t1_meam[elti];
t_ave[i][1] = this->t2_meam[elti];
t_ave[i][2] = this->t3_meam[elti];
} else {
t_ave[i][0] = t_ave[i][0] / rho0[i];
t_ave[i][1] = t_ave[i][1] / rho0[i];
t_ave[i][2] = t_ave[i][2] / rho0[i];
}
}
gamma[i] = t_ave[i][0] * rho1[i] + t_ave[i][1] * rho2[i] + t_ave[i][2] * rho3[i];
if (rho0[i] > 0.0) {
gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
}
Z = this->Z_meam[elti];
G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
if (errorflag != 0)
return;
get_shpfcn(this->lattce_meam[elti][elti], shp);
if (this->ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
if (this->mix_ref_t == 1) {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
} else {
gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /
(Z * Z);
}
Gbar = G_gam(gam, this->ibar_meam[elti], errorflag);
}
rho[i] = rho0[i] * G;
if (this->mix_ref_t == 1) {
if (this->ibar_meam[elti] <= 0) {
Gbar = 1.0;
dGbar = 0.0;
} else {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
Gbar = dG_gam(gam, this->ibar_meam[elti], dGbar);
}
rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd = this->rho0_meam[elti] * Z;
} else {
rho_bkgd = this->rho_ref_meam[elti];
}
}
rhob = rho[i] / rho_bkgd;
denom = 1.0 / rho_bkgd;
G = dG_gam(gamma[i], this->ibar_meam[elti], dG);
dgamma1[i] = (G - 2 * dG * gamma[i]) * denom;
if (!iszero(rho0[i])) {
dgamma2[i] = (dG / rho0[i]) * denom;
} else {
dgamma2[i] = 0.0;
}
// dgamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (this->mix_ref_t == 1) {
dgamma3[i] = rho0[i] * G * dGbar / (Gbar * Z * Z) * denom;
} else {
dgamma3[i] = 0.0;
}
B = this->A_meam[elti] * this->Ec_meam[elti][elti];
if (!iszero(rhob)) {
if (this->emb_lin_neg == 1 && rhob <= 0) {
frhop[i] = -B;
} else {
frhop[i] = B * (log(rhob) + 1.0);
}
if (eflag_either != 0) {
if (eflag_global != 0) {
if (this->emb_lin_neg == 1 && rhob <= 0) {
*eng_vdwl = *eng_vdwl - B * rhob;
} else {
*eng_vdwl = *eng_vdwl + B * rhob * log(rhob);
}
}
if (eflag_atom != 0) {
if (this->emb_lin_neg == 1 && rhob <= 0) {
eatom[i] = eatom[i] - B * rhob;
} else {
eatom[i] = eatom[i] + B * rhob * log(rhob);
}
}
}
} else {
if (this->emb_lin_neg == 1) {
frhop[i] = -B;
} else {
frhop[i] = B;
}
}
}
}
}

355
src/USER-MEAMC/meam_dens_init.cpp Normal file
View File

@ -0,0 +1,355 @@
#include "meam.h"
#include "math_special.h"
using namespace LAMMPS_NS;
void
MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
{
int i, j;
// grow local arrays if necessary
if (atom_nmax > nmax) {
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
nmax = atom_nmax;
memory->create(rho, nmax, "pair:rho");
memory->create(rho0, nmax, "pair:rho0");
memory->create(rho1, nmax, "pair:rho1");
memory->create(rho2, nmax, "pair:rho2");
memory->create(rho3, nmax, "pair:rho3");
memory->create(frhop, nmax, "pair:frhop");
memory->create(gamma, nmax, "pair:gamma");
memory->create(dgamma1, nmax, "pair:dgamma1");
memory->create(dgamma2, nmax, "pair:dgamma2");
memory->create(dgamma3, nmax, "pair:dgamma3");
memory->create(arho2b, nmax, "pair:arho2b");
memory->create(arho1, nmax, 3, "pair:arho1");
memory->create(arho2, nmax, 6, "pair:arho2");
memory->create(arho3, nmax, 10, "pair:arho3");
memory->create(arho3b, nmax, 3, "pair:arho3b");
memory->create(t_ave, nmax, 3, "pair:t_ave");
memory->create(tsq_ave, nmax, 3, "pair:tsq_ave");
}
if (n_neigh > maxneigh) {
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
maxneigh = n_neigh;
memory->create(scrfcn, maxneigh, "pair:scrfcn");
memory->create(dscrfcn, maxneigh, "pair:dscrfcn");
memory->create(fcpair, maxneigh, "pair:fcpair");
}
// zero out local arrays
for (i = 0; i < nall; i++) {
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
for (j = 0; j < 6; j++)
arho2[i][j] = 0.0;
for (j = 0; j < 10; j++)
arho3[i][j] = 0.0;
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
}
void
MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x,
int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset)
{
// Compute screening function and derivatives
getscreen(i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, numneigh, firstneigh,
numneigh_full, firstneigh_full, ntype, type, fmap);
// Calculate intermediate density terms to be communicated
calc_rho1(i, ntype, type, fmap, x, numneigh, firstneigh, &scrfcn[fnoffset], &fcpair[fnoffset]);
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
int* firstneigh, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap)
{
int jn, j, kn, k;
int elti, eltj, eltk;
double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij;
double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/;
double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/;
double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj;
double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2;
double dCikj;
double rnorm, fc, dfc, drinv;
drinv = 1.0 / this->delr_meam;
elti = fmap[type[i]];
if (elti < 0) return;
xitmp = x[i][0];
yitmp = x[i][1];
zitmp = x[i][2];
for (jn = 0; jn < numneigh; jn++) {
j = firstneigh[jn];
eltj = fmap[type[j]];
if (eltj < 0) continue;
// First compute screening function itself, sij
xjtmp = x[j][0];
yjtmp = x[j][1];
zjtmp = x[j][2];
delxij = xjtmp - xitmp;
delyij = yjtmp - yitmp;
delzij = zjtmp - zitmp;
rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
rij = sqrt(rij2);
if (rij > this->rc_meam) {
fcij = 0.0;
dfcij = 0.0;
sij = 0.0;
} else {
rnorm = (this->rc_meam - rij) * drinv;
sij = 1.0;
// if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
const double rbound = this->ebound_meam[elti][eltj] * rij2;
for (kn = 0; kn < numneigh_full; kn++) {
k = firstneigh_full[kn];
eltk = fmap[type[k]];
if (eltk < 0) continue;
if (k == j) continue;
delxjk = x[k][0] - xjtmp;
delyjk = x[k][1] - yjtmp;
delzjk = x[k][2] - zjtmp;
rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
if (rjk2 > rbound) continue;
delxik = x[k][0] - xitmp;
delyik = x[k][1] - yitmp;
delzik = x[k][2] - zitmp;
rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
if (rik2 > rbound) continue;
xik = rik2 / rij2;
xjk = rjk2 / rij2;
a = 1 - (xik - xjk) * (xik - xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a <= 0.0) continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = this->Cmax_meam[elti][eltj][eltk];
Cmin = this->Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax) continue;
// note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case here
// (other negative cikj cases were handled by the test on "a" above)
else if (cikj <= Cmin) {
sij = 0.0;
break;
} else {
delc = Cmax - Cmin;
cikj = (cikj - Cmin) / delc;
sikj = fcut(cikj);
}
sij *= sikj;
}
fc = dfcut(rnorm, dfc);
fcij = fc;
dfcij = dfc * drinv;
}
// Now compute derivatives
dscrfcn[jn] = 0.0;
sfcij = sij * fcij;
if (iszero(sfcij) || iszero(sfcij - 1.0))
goto LABEL_100;
rbound = this->ebound_meam[elti][eltj] * rij2;
for (kn = 0; kn < numneigh_full; kn++) {
k = firstneigh_full[kn];
if (k == j) continue;
eltk = fmap[type[k]];
if (eltk < 0) continue;
xktmp = x[k][0];
yktmp = x[k][1];
zktmp = x[k][2];
delxjk = xktmp - xjtmp;
delyjk = yktmp - yjtmp;
delzjk = zktmp - zjtmp;
rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
if (rjk2 > rbound) continue;
delxik = xktmp - xitmp;
delyik = yktmp - yitmp;
delzik = zktmp - zitmp;
rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
if (rik2 > rbound) continue;
xik = rik2 / rij2;
xjk = rjk2 / rij2;
a = 1 - (xik - xjk) * (xik - xjk);
// if a < 0, then ellipse equation doesn't describe this case and
// atom k can't possibly screen i-j
if (a <= 0.0) continue;
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
Cmax = this->Cmax_meam[elti][eltj][eltk];
Cmin = this->Cmin_meam[elti][eltj][eltk];
if (cikj >= Cmax) {
continue;
// Note that cikj may be slightly negative (within numerical
// tolerance) if atoms are colinear, so don't reject that case
// here
// (other negative cikj cases were handled by the test on "a"
// above)
// Note that we never have 0<cikj<Cmin here, else sij=0
// (rejected above)
} else {
delc = Cmax - Cmin;
cikj = (cikj - Cmin) / delc;
sikj = dfcut(cikj, dfikj);
coef1 = dfikj / (delc * sikj);
dCikj = dCfunc(rij2, rik2, rjk2);
dscrfcn[jn] = dscrfcn[jn] + coef1 * dCikj;
}
}
coef1 = sfcij;
coef2 = sij * dfcij / rij;
dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2;
LABEL_100:
scrfcn[jn] = sij;
fcpair[jn] = fcij;
}
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* scrfcn, double* fcpair)
{
int jn, j, m, n, p, elti, eltj;
int nv2, nv3;
double xtmp, ytmp, ztmp, delij[3], rij2, rij, sij;
double ai, aj, rhoa0j, rhoa1j, rhoa2j, rhoa3j, A1j, A2j, A3j;
// double G,Gbar,gam,shp[3+1];
double ro0i, ro0j;
double rhoa0i, rhoa1i, rhoa2i, rhoa3i, A1i, A2i, A3i;
elti = fmap[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
for (jn = 0; jn < numneigh; jn++) {
if (!iszero(scrfcn[jn])) {
j = firstneigh[jn];
sij = scrfcn[jn] * fcpair[jn];
delij[0] = x[j][0] - xtmp;
delij[1] = x[j][1] - ytmp;
delij[2] = x[j][2] - ztmp;
rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
if (rij2 < this->cutforcesq) {
eltj = fmap[type[j]];
rij = sqrt(rij2);
ai = rij / this->re_meam[elti][elti] - 1.0;
aj = rij / this->re_meam[eltj][eltj] - 1.0;
ro0i = this->rho0_meam[elti];
ro0j = this->rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj) * sij;
rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj) * sij;
rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj) * sij;
rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj) * sij;
rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai) * sij;
rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai) * sij;
rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai) * sij;
rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai) * sij;
if (this->ialloy == 1) {
rhoa1j = rhoa1j * this->t1_meam[eltj];
rhoa2j = rhoa2j * this->t2_meam[eltj];
rhoa3j = rhoa3j * this->t3_meam[eltj];
rhoa1i = rhoa1i * this->t1_meam[elti];
rhoa2i = rhoa2i * this->t2_meam[elti];
rhoa3i = rhoa3i * this->t3_meam[elti];
}
rho0[i] = rho0[i] + rhoa0j;
rho0[j] = rho0[j] + rhoa0i;
// For ialloy = 2, use single-element value (not average)
if (this->ialloy != 2) {
t_ave[i][0] = t_ave[i][0] + this->t1_meam[eltj] * rhoa0j;
t_ave[i][1] = t_ave[i][1] + this->t2_meam[eltj] * rhoa0j;
t_ave[i][2] = t_ave[i][2] + this->t3_meam[eltj] * rhoa0j;
t_ave[j][0] = t_ave[j][0] + this->t1_meam[elti] * rhoa0i;
t_ave[j][1] = t_ave[j][1] + this->t2_meam[elti] * rhoa0i;
t_ave[j][2] = t_ave[j][2] + this->t3_meam[elti] * rhoa0i;
}
if (this->ialloy == 1) {
tsq_ave[i][0] = tsq_ave[i][0] + this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j;
tsq_ave[i][1] = tsq_ave[i][1] + this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j;
tsq_ave[i][2] = tsq_ave[i][2] + this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j;
tsq_ave[j][0] = tsq_ave[j][0] + this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i;
tsq_ave[j][1] = tsq_ave[j][1] + this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i;
tsq_ave[j][2] = tsq_ave[j][2] + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i;
}
arho2b[i] = arho2b[i] + rhoa2j;
arho2b[j] = arho2b[j] + rhoa2i;
A1j = rhoa1j / rij;
A2j = rhoa2j / rij2;
A3j = rhoa3j / (rij2 * rij);
A1i = rhoa1i / rij;
A2i = rhoa2i / rij2;
A3i = rhoa3i / (rij2 * rij);
nv2 = 0;
nv3 = 0;
for (m = 0; m < 3; m++) {
arho1[i][m] = arho1[i][m] + A1j * delij[m];
arho1[j][m] = arho1[j][m] - A1i * delij[m];
arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij;
arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij;
for (n = m; n < 3; n++) {
arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n];
arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n];
nv2 = nv2 + 1;
for (p = n; p < 3; p++) {
arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p];
arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p];
nv3 = nv3 + 1;
}
}
}
}
}
}
}

531
src/USER-MEAMC/meam_force.cpp Normal file
View File

@ -0,0 +1,531 @@
#include "meam.h"
#include "math_special.h"
#include <algorithm>
using namespace LAMMPS_NS;
void
MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,
double* eatom, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom)
{
int j, jn, k, kn, kk, m, n, p, q;
int nv2, nv3, elti, eltj, eltk, ind;
double xitmp, yitmp, zitmp, delij[3], rij2, rij, rij3;
double v[6], fi[3], fj[3];
double third, sixth;
double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem;
double r, recip, phi, phip;
double sij;
double a1, a1i, a1j, a2, a2i, a2j;
double a3i, a3j;
double shpi[3], shpj[3];
double ai, aj, ro0i, ro0j, invrei, invrej;
double rhoa0j, drhoa0j, rhoa0i, drhoa0i;
double rhoa1j, drhoa1j, rhoa1i, drhoa1i;
double rhoa2j, drhoa2j, rhoa2i, drhoa2i;
double a3, a3a, rhoa3j, drhoa3j, rhoa3i, drhoa3i;
double drho0dr1, drho0dr2, drho0ds1, drho0ds2;
double drho1dr1, drho1dr2, drho1ds1, drho1ds2;
double drho1drm1[3], drho1drm2[3];
double drho2dr1, drho2dr2, drho2ds1, drho2ds2;
double drho2drm1[3], drho2drm2[3];
double drho3dr1, drho3dr2, drho3ds1, drho3ds2;
double drho3drm1[3], drho3drm2[3];
double dt1dr1, dt1dr2, dt1ds1, dt1ds2;
double dt2dr1, dt2dr2, dt2ds1, dt2ds2;
double dt3dr1, dt3dr2, dt3ds1, dt3ds2;
double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3], drhodrm2[3];
double arg;
double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;
double dsij1, dsij2, force1, force2;
double t1i, t2i, t3i, t1j, t2j, t3j;
third = 1.0 / 3.0;
sixth = 1.0 / 6.0;
// Compute forces atom i
elti = fmap[type[i]];
if (elti < 0) return;
xitmp = x[i][0];
yitmp = x[i][1];
zitmp = x[i][2];
// Treat each pair
for (jn = 0; jn < numneigh; jn++) {
j = firstneigh[jn];
eltj = fmap[type[j]];
if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {
sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn];
delij[0] = x[j][0] - xitmp;
delij[1] = x[j][1] - yitmp;
delij[2] = x[j][2] - zitmp;
rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
if (rij2 < this->cutforcesq) {
rij = sqrt(rij2);
r = rij;
// Compute phi and phip
ind = this->eltind[elti][eltj];
pp = rij * this->rdrar;
kk = (int)pp;
kk = std::min(kk, this->nrar - 2);
pp = pp - kk;
pp = std::min(pp, 1.0);
phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
this->phirar[ind][kk];
phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
recip = 1.0 / r;
if (eflag_either != 0) {
if (eflag_global != 0)
*eng_vdwl = *eng_vdwl + phi * sij;
if (eflag_atom != 0) {
eatom[i] = eatom[i] + 0.5 * phi * sij;
eatom[j] = eatom[j] + 0.5 * phi * sij;
}
}
// write(1,*) "force_meamf: phi: ",phi
// write(1,*) "force_meamf: phip: ",phip
// Compute pair densities and derivatives
invrei = 1.0 / this->re_meam[elti][elti];
ai = rij * invrei - 1.0;
ro0i = this->rho0_meam[elti];
rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai);
drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i;
rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai);
drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i;
rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai);
drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i;
rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai);
drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i;
if (elti != eltj) {
invrej = 1.0 / this->re_meam[eltj][eltj];
aj = rij * invrej - 1.0;
ro0j = this->rho0_meam[eltj];
rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj);
drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j;
rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj);
drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j;
rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj);
drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j;
rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj);
drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j;
} else {
rhoa0j = rhoa0i;
drhoa0j = drhoa0i;
rhoa1j = rhoa1i;
drhoa1j = drhoa1i;
rhoa2j = rhoa2i;
drhoa2j = drhoa2i;
rhoa3j = rhoa3i;
drhoa3j = drhoa3i;
}
const double t1mi = this->t1_meam[elti];
const double t2mi = this->t2_meam[elti];
const double t3mi = this->t3_meam[elti];
const double t1mj = this->t1_meam[eltj];
const double t2mj = this->t2_meam[eltj];
const double t3mj = this->t3_meam[eltj];
if (this->ialloy == 1) {
rhoa1j *= t1mj;
rhoa2j *= t2mj;
rhoa3j *= t3mj;
rhoa1i *= t1mi;
rhoa2i *= t2mi;
rhoa3i *= t3mi;
drhoa1j *= t1mj;
drhoa2j *= t2mj;
drhoa3j *= t3mj;
drhoa1i *= t1mi;
drhoa2i *= t2mi;
drhoa3i *= t3mi;
}
nv2 = 0;
nv3 = 0;
arg1i1 = 0.0;
arg1j1 = 0.0;
arg1i2 = 0.0;
arg1j2 = 0.0;
arg1i3 = 0.0;
arg1j3 = 0.0;
arg3i3 = 0.0;
arg3j3 = 0.0;
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
for (q = p; q < 3; q++) {
arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
arg1i3 = arg1i3 + arho3[i][nv3] * arg;
arg1j3 = arg1j3 - arho3[j][nv3] * arg;
nv3 = nv3 + 1;
}
arg = delij[n] * delij[p] * this->v2D[nv2];
arg1i2 = arg1i2 + arho2[i][nv2] * arg;
arg1j2 = arg1j2 + arho2[j][nv2] * arg;
nv2 = nv2 + 1;
}
arg1i1 = arg1i1 + arho1[i][n] * delij[n];
arg1j1 = arg1j1 - arho1[j][n] * delij[n];
arg3i3 = arg3i3 + arho3b[i][n] * delij[n];
arg3j3 = arg3j3 - arho3b[j][n] * delij[n];
}
// rho0 terms
drho0dr1 = drhoa0j * sij;
drho0dr2 = drhoa0i * sij;
// rho1 terms
a1 = 2 * sij / rij;
drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1;
drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1;
a1 = 2.0 * sij / rij;
for (m = 0; m < 3; m++) {
drho1drm1[m] = a1 * rhoa1j * arho1[i][m];
drho1drm2[m] = -a1 * rhoa1i * arho1[j][m];
}
// rho2 terms
a2 = 2 * sij / rij2;
drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arho2b[i] * drhoa2j * sij;
drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arho2b[j] * drhoa2i * sij;
a2 = 4 * sij / rij2;
for (m = 0; m < 3; m++) {
drho2drm1[m] = 0.0;
drho2drm2[m] = 0.0;
for (n = 0; n < 3; n++) {
drho2drm1[m] = drho2drm1[m] + arho2[i][this->vind2D[m][n]] * delij[n];
drho2drm2[m] = drho2drm2[m] - arho2[j][this->vind2D[m][n]] * delij[n];
}
drho2drm1[m] = a2 * rhoa2j * drho2drm1[m];
drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m];
}
// rho3 terms
rij3 = rij * rij2;
a3 = 2 * sij / rij3;
a3a = 6.0 / 5.0 * sij / rij;
drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3;
drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3;
a3 = 6 * sij / rij3;
a3a = 6 * sij / (5 * rij);
for (m = 0; m < 3; m++) {
drho3drm1[m] = 0.0;
drho3drm2[m] = 0.0;
nv2 = 0;
for (n = 0; n < 3; n++) {
for (p = n; p < 3; p++) {
arg = delij[n] * delij[p] * this->v2D[nv2];
drho3drm1[m] = drho3drm1[m] + arho3[i][this->vind3D[m][n][p]] * arg;
drho3drm2[m] = drho3drm2[m] + arho3[j][this->vind3D[m][n][p]] * arg;
nv2 = nv2 + 1;
}
}
drho3drm1[m] = (a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j;
drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
}
// Compute derivatives of weighting functions t wrt rij
t1i = t_ave[i][0];
t2i = t_ave[i][1];
t3i = t_ave[i][2];
t1j = t_ave[j][0];
t2j = t_ave[j][1];
t3j = t_ave[j][2];
if (this->ialloy == 1) {
a1i = 0.0;
a1j = 0.0;
a2i = 0.0;
a2j = 0.0;
a3i = 0.0;
a3j = 0.0;
if (!iszero(tsq_ave[i][0]))
a1i = drhoa0j * sij / tsq_ave[i][0];
if (!iszero(tsq_ave[j][0]))
a1j = drhoa0i * sij / tsq_ave[j][0];
if (!iszero(tsq_ave[i][1]))
a2i = drhoa0j * sij / tsq_ave[i][1];
if (!iszero(tsq_ave[j][1]))
a2j = drhoa0i * sij / tsq_ave[j][1];
if (!iszero(tsq_ave[i][2]))
a3i = drhoa0j * sij / tsq_ave[i][2];
if (!iszero(tsq_ave[j][2]))
a3j = drhoa0i * sij / tsq_ave[j][2];
dt1dr1 = a1i * (t1mj - t1i * t1mj*t1mj);
dt1dr2 = a1j * (t1mi - t1j * t1mi*t1mi);
dt2dr1 = a2i * (t2mj - t2i * t2mj*t2mj);
dt2dr2 = a2j * (t2mi - t2j * t2mi*t2mi);
dt3dr1 = a3i * (t3mj - t3i * t3mj*t3mj);
dt3dr2 = a3j * (t3mi - t3j * t3mi*t3mi);
} else if (this->ialloy == 2) {
dt1dr1 = 0.0;
dt1dr2 = 0.0;
dt2dr1 = 0.0;
dt2dr2 = 0.0;
dt3dr1 = 0.0;
dt3dr2 = 0.0;
} else {
ai = 0.0;
if (!iszero(rho0[i]))
ai = drhoa0j * sij / rho0[i];
aj = 0.0;
if (!iszero(rho0[j]))
aj = drhoa0i * sij / rho0[j];
dt1dr1 = ai * (t1mj - t1i);
dt1dr2 = aj * (t1mi - t1j);
dt2dr1 = ai * (t2mj - t2i);
dt2dr2 = aj * (t2mi - t2j);
dt3dr1 = ai * (t3mj - t3i);
dt3dr2 = aj * (t3mi - t3j);
}
// Compute derivatives of total density wrt rij, sij and rij(3)
get_shpfcn(this->lattce_meam[elti][elti], shpi);
get_shpfcn(this->lattce_meam[eltj][eltj], shpj);
drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
dt3dr1 * rho3[i] + t3i * drho3dr1) -
dgamma3[i] * (shpi[0] * dt1dr1 + shpi[1] * dt2dr1 + shpi[2] * dt3dr1);
drhodr2 = dgamma1[j] * drho0dr2 +
dgamma2[j] * (dt1dr2 * rho1[j] + t1j * drho1dr2 + dt2dr2 * rho2[j] + t2j * drho2dr2 +
dt3dr2 * rho3[j] + t3j * drho3dr2) -
dgamma3[j] * (shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2);
for (m = 0; m < 3; m++) {
drhodrm1[m] = 0.0;
drhodrm2[m] = 0.0;
drhodrm1[m] = dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]);
drhodrm2[m] = dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]);
}
// Compute derivatives wrt sij, but only if necessary
if (!iszero(dscrfcn[fnoffset + jn])) {
drho0ds1 = rhoa0j;
drho0ds2 = rhoa0i;
a1 = 2.0 / rij;
drho1ds1 = a1 * rhoa1j * arg1i1;
drho1ds2 = a1 * rhoa1i * arg1j1;
a2 = 2.0 / rij2;
drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j;
drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i;
a3 = 2.0 / rij3;
a3a = 6.0 / (5.0 * rij);
drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
if (this->ialloy == 1) {
a1i = 0.0;
a1j = 0.0;
a2i = 0.0;
a2j = 0.0;
a3i = 0.0;
a3j = 0.0;
if (!iszero(tsq_ave[i][0]))
a1i = rhoa0j / tsq_ave[i][0];
if (!iszero(tsq_ave[j][0]))
a1j = rhoa0i / tsq_ave[j][0];
if (!iszero(tsq_ave[i][1]))
a2i = rhoa0j / tsq_ave[i][1];
if (!iszero(tsq_ave[j][1]))
a2j = rhoa0i / tsq_ave[j][1];
if (!iszero(tsq_ave[i][2]))
a3i = rhoa0j / tsq_ave[i][2];
if (!iszero(tsq_ave[j][2]))
a3j = rhoa0i / tsq_ave[j][2];
dt1ds1 = a1i * (t1mj - t1i * pow(t1mj, 2));
dt1ds2 = a1j * (t1mi - t1j * pow(t1mi, 2));
dt2ds1 = a2i * (t2mj - t2i * pow(t2mj, 2));
dt2ds2 = a2j * (t2mi - t2j * pow(t2mi, 2));
dt3ds1 = a3i * (t3mj - t3i * pow(t3mj, 2));
dt3ds2 = a3j * (t3mi - t3j * pow(t3mi, 2));
} else if (this->ialloy == 2) {
dt1ds1 = 0.0;
dt1ds2 = 0.0;
dt2ds1 = 0.0;
dt2ds2 = 0.0;
dt3ds1 = 0.0;
dt3ds2 = 0.0;
} else {
ai = 0.0;
if (!iszero(rho0[i]))
ai = rhoa0j / rho0[i];
aj = 0.0;
if (!iszero(rho0[j]))
aj = rhoa0i / rho0[j];
dt1ds1 = ai * (t1mj - t1i);
dt1ds2 = aj * (t1mi - t1j);
dt2ds1 = ai * (t2mj - t2i);
dt2ds2 = aj * (t2mi - t2j);
dt3ds1 = ai * (t3mj - t3i);
dt3ds2 = aj * (t3mi - t3j);
}
drhods1 = dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 +
dt3ds1 * rho3[i] + t3i * drho3ds1) -
dgamma3[i] * (shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
drhods2 = dgamma1[j] * drho0ds2 +
dgamma2[j] * (dt1ds2 * rho1[j] + t1j * drho1ds2 + dt2ds2 * rho2[j] + t2j * drho2ds2 +
dt3ds2 * rho3[j] + t3j * drho3ds2) -
dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
}
// Compute derivatives of energy wrt rij, sij and rij[3]
dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2;
dUdsij = 0.0;
if (!iszero(dscrfcn[fnoffset + jn])) {
dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2;
}
for (m = 0; m < 3; m++) {
dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
}
// Add the part of the force due to dUdrij and dUdsij
force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn];
for (m = 0; m < 3; m++) {
forcem = delij[m] * force + dUdrijm[m];
f[i][m] = f[i][m] + forcem;
f[j][m] = f[j][m] - forcem;
}
// Tabulate per-atom virial as symmetrized stress tensor
if (vflag_atom != 0) {
fi[0] = delij[0] * force + dUdrijm[0];
fi[1] = delij[1] * force + dUdrijm[1];
fi[2] = delij[2] * force + dUdrijm[2];
v[0] = -0.5 * (delij[0] * fi[0]);
v[1] = -0.5 * (delij[1] * fi[1]);
v[2] = -0.5 * (delij[2] * fi[2]);
v[3] = -0.25 * (delij[0] * fi[1] + delij[1] * fi[0]);
v[4] = -0.25 * (delij[0] * fi[2] + delij[2] * fi[0]);
v[5] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]);
for (m = 0; m < 6; m++) {
vatom[i][m] = vatom[i][m] + v[m];
vatom[j][m] = vatom[j][m] + v[m];
}
}
// Now compute forces on other atoms k due to change in sij
if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop
double dxik(0), dyik(0), dzik(0);
double dxjk(0), dyjk(0), dzjk(0);
for (kn = 0; kn < numneigh_full; kn++) {
k = firstneigh_full[kn];
eltk = fmap[type[k]];
if (k != j && eltk >= 0) {
double xik, xjk, cikj, sikj, dfc, a;
double dCikj1, dCikj2;
double delc, rik2, rjk2;
sij = scrfcn[jn+fnoffset] * fcpair[jn+fnoffset];
const double Cmax = this->Cmax_meam[elti][eltj][eltk];
const double Cmin = this->Cmin_meam[elti][eltj][eltk];
dsij1 = 0.0;
dsij2 = 0.0;
if (!iszero(sij) && !iszero(sij - 1.0)) {
const double rbound = rij2 * this->ebound_meam[elti][eltj];
delc = Cmax - Cmin;
dxjk = x[k][0] - x[j][0];
dyjk = x[k][1] - x[j][1];
dzjk = x[k][2] - x[j][2];
rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
if (rjk2 <= rbound) {
dxik = x[k][0] - x[i][0];
dyik = x[k][1] - x[i][1];
dzik = x[k][2] - x[i][2];
rik2 = dxik * dxik + dyik * dyik + dzik * dzik;
if (rik2 <= rbound) {
xik = rik2 / rij2;
xjk = rjk2 / rij2;
a = 1 - (xik - xjk) * (xik - xjk);
if (!iszero(a)) {
cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
if (cikj >= Cmin && cikj <= Cmax) {
cikj = (cikj - Cmin) / delc;
sikj = dfcut(cikj, dfc);
dCfunc2(rij2, rik2, rjk2, dCikj1, dCikj2);
a = sij / delc * dfc / sikj;
dsij1 = a * dCikj1;
dsij2 = a * dCikj2;
}
}
}
}
}
if (!iszero(dsij1) || !iszero(dsij2)) {
force1 = dUdsij * dsij1;
force2 = dUdsij * dsij2;
f[i][0] += force1 * dxik;
f[i][1] += force1 * dyik;
f[i][2] += force1 * dzik;
f[j][0] += force2 * dxjk;
f[j][1] += force2 * dyjk;
f[j][2] += force2 * dzjk;
f[k][0] -= force1 * dxik + force2 * dxjk;
f[k][1] -= force1 * dyik + force2 * dyjk;
f[k][2] -= force1 * dzik + force2 * dzjk;
// Tabulate per-atom virial as symmetrized stress tensor
if (vflag_atom != 0) {
fi[0] = force1 * dxik;
fi[1] = force1 * dyik;
fi[2] = force1 * dzik;
fj[0] = force2 * dxjk;
fj[1] = force2 * dyjk;
fj[2] = force2 * dzjk;
v[0] = -third * (dxik * fi[0] + dxjk * fj[0]);
v[1] = -third * (dyik * fi[1] + dyjk * fj[1]);
v[2] = -third * (dzik * fi[2] + dzjk * fj[2]);
v[3] = -sixth * (dxik * fi[1] + dxjk * fj[1] + dyik * fi[0] + dyjk * fj[0]);
v[4] = -sixth * (dxik * fi[2] + dxjk * fj[2] + dzik * fi[0] + dzjk * fj[0]);
v[5] = -sixth * (dyik * fi[2] + dyjk * fj[2] + dzik * fi[1] + dzjk * fj[1]);
for (m = 0; m < 6; m++) {
vatom[i][m] = vatom[i][m] + v[m];
vatom[j][m] = vatom[j][m] + v[m];
vatom[k][m] = vatom[k][m] + v[m];
}
}
}
}
// end of k loop
}
}
}
// end of j loop
}
}

321
src/USER-MEAMC/meam_funcs.cpp Normal file
View File

@ -0,0 +1,321 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Sebastian Hütter (OvGU)
------------------------------------------------------------------------- */
#include "meam.h"
#include "math_special.h"
#include <math.h>
using namespace LAMMPS_NS;
//-----------------------------------------------------------------------------
// Compute G(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+gamma)
// 1 => G = exp(gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+exp(-gamma))
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
//
double
MEAM::G_gam(const double gamma, const int ibar, int& errorflag) const
{
double gsmooth_switchpoint;
switch (ibar) {
case 0:
case 4:
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/gamma)**99
double G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
return sqrt(G);
} else {
return sqrt(1.0 + gamma);
}
case 1:
return MathSpecial::fm_exp(gamma / 2.0);
case 3:
return 2.0 / (1.0 + exp(-gamma));
case -5:
if ((1.0 + gamma) >= 0) {
return sqrt(1.0 + gamma);
} else {
return -sqrt(-1.0 - gamma);
}
}
errorflag = 1;
return 0.0;
}
//-----------------------------------------------------------------------------
// Compute G(gamma and dG(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+gamma)
// 1 => G = exp(gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+exp(-gamma))
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
//
double
MEAM::dG_gam(const double gamma, const int ibar, double& dG) const
{
double gsmooth_switchpoint;
double G;
switch (ibar) {
case 0:
case 4:
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/gamma)**99
double G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
G = sqrt(G);
dG = -gsmooth_factor * G / (2.0 * gamma);
return G;
} else {
G = sqrt(1.0 + gamma);
dG = 1.0 / (2.0 * G);
return G;
}
case 1:
G = MathSpecial::fm_exp(gamma / 2.0);
dG = G / 2.0;
return G;
case 3:
G = 2.0 / (1.0 + MathSpecial::fm_exp(-gamma));
dG = G * (2.0 - G) / 2;
return G;
case -5:
if ((1.0 + gamma) >= 0) {
G = sqrt(1.0 + gamma);
dG = 1.0 / (2.0 * G);
return G;
} else {
G = -sqrt(-1.0 - gamma);
dG = -1.0 / (2.0 * G);
return G;
}
}
dG = 1.0;
return 0.0;
}
//-----------------------------------------------------------------------------
// Compute ZBL potential
//
double
MEAM::zbl(const double r, const int z1, const int z2)
{
int i;
const double c[] = { 0.028171, 0.28022, 0.50986, 0.18175 };
const double d[] = { 0.20162, 0.40290, 0.94229, 3.1998 };
const double azero = 0.4685;
const double cc = 14.3997;
double a, x;
// azero = (9pi^2/128)^1/3 (0.529) Angstroms
a = azero / (pow(z1, 0.23) + pow(z2, 0.23));
double result = 0.0;
x = r / a;
for (i = 0; i <= 3; i++) {
result = result + c[i] * MathSpecial::fm_exp(-d[i] * x);
}
if (r > 0.0)
result = result * z1 * z2 / r * cc;
return result;
}
//-----------------------------------------------------------------------------
// Compute Rose energy function, I.16
//
double
MEAM::erose(const double r, const double re, const double alpha, const double Ec, const double repuls,
const double attrac, const int form)
{
double astar, a3;
double result = 0.0;
if (r > 0.0) {
astar = alpha * (r / re - 1.0);
a3 = 0.0;
if (astar >= 0)
a3 = attrac;
else if (astar < 0)
a3 = repuls;
if (form == 1)
result = -Ec * (1 + astar + (-attrac + repuls / r) * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
else if (form == 2)
result = -Ec * (1 + astar + a3 * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
else
result = -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * MathSpecial::fm_exp(-astar);
}
return result;
}
//-----------------------------------------------------------------------------
// Shape factors for various configurations
//
void
MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
{
switch (latt) {
case FCC:
case BCC:
case B1:
case B2:
s[0] = 0.0;
s[1] = 0.0;
s[2] = 0.0;
break;
case HCP:
s[0] = 0.0;
s[1] = 0.0;
s[2] = 1.0 / 3.0;
break;
case DIA:
s[0] = 0.0;
s[1] = 0.0;
s[2] = 32.0 / 9.0;
break;
case DIM:
s[0] = 1.0;
s[1] = 2.0 / 3.0;
// s(3) = 1.d0
s[2] = 0.40;
break;
default:
s[0] = 0.0;
// call error('Lattice not defined in get_shpfcn.')
}
}
//-----------------------------------------------------------------------------
// Number of neighbors for the reference structure
//
int
MEAM::get_Zij(const lattice_t latt)
{
switch (latt) {
case FCC:
return 12;
case BCC:
return 8;
case HCP:
return 12;
case B1:
return 6;
case DIA:
return 4;
case DIM:
return 1;
case C11:
return 10;
case L12:
return 12;
case B2:
return 8;
// call error('Lattice not defined in get_Zij.')
}
return 0;
}
//-----------------------------------------------------------------------------
// Number of second neighbors for the reference structure
// a = distance ratio R1/R2
// S = second neighbor screening function
//
int
MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, double& a, double& S)
{
double C, x, sijk;
int Zij2 = 0, numscr = 0;
switch (latt) {
case FCC:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case BCC:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case HCP:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case B1:
Zij2 = 12;
a = sqrt(2.0);
numscr = 2;
break;
case DIA:
Zij2 = 0;
a = sqrt(8.0 / 3.0);
numscr = 4;
if (cmin < 0.500001) {
// call error('can not do 2NN MEAM for dia')
}
break;
case DIM:
// this really shouldn't be allowed; make sure screening is zero
a = 1.0;
S = 0.0;
return 0;
case L12:
Zij2 = 6;
a = sqrt(2.0);
numscr = 4;
break;
case B2:
Zij2 = 6;
a = 2.0 / sqrt(3.0);
numscr = 4;
break;
case C11:
// unsupported lattice flag C11 in get_Zij
break;
default:
// unknown lattic flag in get Zij
// call error('Lattice not defined in get_Zij.')
break;
}
// Compute screening for each first neighbor
C = 4.0 / (a * a) - 1.0;
x = (C - cmin) / (cmax - cmin);
sijk = fcut(x);
// There are numscr first neighbors screening the second neighbors
S = MathSpecial::powint(sijk, numscr);
return Zij2;
}

72
src/USER-MEAMC/meam_impl.cpp Normal file
View File

@ -0,0 +1,72 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Sebastian Hütter (OvGU)
------------------------------------------------------------------------- */
#include "meam.h"
#include "memory.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MEAM::MEAM(Memory* mem)
: memory(mem)
{
phir = phirar = phirar1 = phirar2 = phirar3 = phirar4 = phirar5 = phirar6 = NULL;
nmax = 0;
rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
maxneigh = 0;
scrfcn = dscrfcn = fcpair = NULL;
}
MEAM::~MEAM()
{
memory->destroy(this->phirar6);
memory->destroy(this->phirar5);
memory->destroy(this->phirar4);
memory->destroy(this->phirar3);
memory->destroy(this->phirar2);
memory->destroy(this->phirar1);
memory->destroy(this->phirar);
memory->destroy(this->phir);
memory->destroy(this->rho);
memory->destroy(this->rho0);
memory->destroy(this->rho1);
memory->destroy(this->rho2);
memory->destroy(this->rho3);
memory->destroy(this->frhop);
memory->destroy(this->gamma);
memory->destroy(this->dgamma1);
memory->destroy(this->dgamma2);
memory->destroy(this->dgamma3);
memory->destroy(this->arho2b);
memory->destroy(this->arho1);
memory->destroy(this->arho2);
memory->destroy(this->arho3);
memory->destroy(this->arho3b);
memory->destroy(this->t_ave);
memory->destroy(this->tsq_ave);
memory->destroy(this->scrfcn);
memory->destroy(this->dscrfcn);
memory->destroy(this->fcpair);
}

792
src/USER-MEAMC/meam_setup_done.cpp Normal file
View File

@ -0,0 +1,792 @@
#include "meam.h"
#include "math_special.h"
#include <algorithm>
using namespace LAMMPS_NS;
void
MEAM::meam_setup_done(double* cutmax)
{
int nv2, nv3, m, n, p;
// Force cutoff
this->cutforce = this->rc_meam;
this->cutforcesq = this->cutforce * this->cutforce;
// Pass cutoff back to calling program
*cutmax = this->cutforce;
// Augment t1 term
for (int i = 0; i < maxelt; i++)
this->t1_meam[i] = this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i];
// Compute off-diagonal alloy parameters
alloyparams();
// indices and factors for Voight notation
nv2 = 0;
nv3 = 0;
for (m = 0; m < 3; m++) {
for (n = m; n < 3; n++) {
this->vind2D[m][n] = nv2;
this->vind2D[n][m] = nv2;
nv2 = nv2 + 1;
for (p = n; p < 3; p++) {
this->vind3D[m][n][p] = nv3;
this->vind3D[m][p][n] = nv3;
this->vind3D[n][m][p] = nv3;
this->vind3D[n][p][m] = nv3;
this->vind3D[p][m][n] = nv3;
this->vind3D[p][n][m] = nv3;
nv3 = nv3 + 1;
}
}
}
this->v2D[0] = 1;
this->v2D[1] = 2;
this->v2D[2] = 2;
this->v2D[3] = 1;
this->v2D[4] = 2;
this->v2D[5] = 1;
this->v3D[0] = 1;
this->v3D[1] = 3;
this->v3D[2] = 3;
this->v3D[3] = 3;
this->v3D[4] = 6;
this->v3D[5] = 3;
this->v3D[6] = 1;
this->v3D[7] = 3;
this->v3D[8] = 3;
this->v3D[9] = 1;
nv2 = 0;
for (m = 0; m < this->neltypes; m++) {
for (n = m; n < this->neltypes; n++) {
this->eltind[m][n] = nv2;
this->eltind[n][m] = nv2;
nv2 = nv2 + 1;
}
}
// Compute background densities for reference structure
compute_reference_density();
// Compute pair potentials and setup arrays for interpolation
this->nr = 1000;
this->dr = 1.1 * this->rc_meam / this->nr;
compute_pair_meam();
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// Fill off-diagonal alloy parameters
void
MEAM::alloyparams(void)
{
int i, j, k;
double eb;
// Loop over pairs
for (i = 0; i < this->neltypes; i++) {
for (j = 0; j < this->neltypes; j++) {
// Treat off-diagonal pairs
// If i>j, set all equal to i<j case (which has aready been set,
// here or in the input file)
if (i > j) {
this->re_meam[i][j] = this->re_meam[j][i];
this->Ec_meam[i][j] = this->Ec_meam[j][i];
this->alpha_meam[i][j] = this->alpha_meam[j][i];
this->lattce_meam[i][j] = this->lattce_meam[j][i];
this->nn2_meam[i][j] = this->nn2_meam[j][i];
// If i<j and term is unset, use default values (e.g. mean of i-i and
// j-j)
} else if (j > i) {
if (iszero(this->Ec_meam[i][j])) {
if (this->lattce_meam[i][j] == L12)
this->Ec_meam[i][j] =
(3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 - this->delta_meam[i][j];
else if (this->lattce_meam[i][j] == C11) {
if (this->lattce_meam[i][i] == DIA)
this->Ec_meam[i][j] =
(2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
else
this->Ec_meam[i][j] =
(this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
} else
this->Ec_meam[i][j] = (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - this->delta_meam[i][j];
}
if (iszero(this->alpha_meam[i][j]))
this->alpha_meam[i][j] = (this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0;
if (iszero(this->re_meam[i][j]))
this->re_meam[i][j] = (this->re_meam[i][i] + this->re_meam[j][j]) / 2.0;
}
}
}
// Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets
// where i>j, set equal to the i<j element. Likewise for Cmax.
for (i = 1; i < this->neltypes; i++) {
for (j = 0; j < i; j++) {
for (k = 0; k < this->neltypes; k++) {
this->Cmin_meam[i][j][k] = this->Cmin_meam[j][i][k];
this->Cmax_meam[i][j][k] = this->Cmax_meam[j][i][k];
}
}
}
// ebound gives the squared distance such that, for rik2 or rjk2>ebound,
// atom k definitely lies outside the screening function ellipse (so
// there is no need to calculate its effects). Here, compute it for all
// triplets [i][j][k] so that ebound[i][j] is the maximized over k
for (i = 0; i < this->neltypes; i++) {
for (j = 0; j < this->neltypes; j++) {
for (k = 0; k < this->neltypes; k++) {
eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0));
this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb);
}
}
}
}
//-----------------------------------------------------------------------
// compute MEAM pair potential for each pair of element types
//
void
MEAM::compute_pair_meam(void)
{
double r /*ununsed:, temp*/;
int j, a, b, nv2;
double astar, frac, phizbl;
int n, nmax, Z1, Z2;
double arat, rarat, scrn, scrn2;
double phiaa, phibb /*unused:,phitmp*/;
double C, s111, s112, s221, S11, S22;
// check for previously allocated arrays and free them
if (this->phir != NULL)
memory->destroy(this->phir);
if (this->phirar != NULL)
memory->destroy(this->phirar);
if (this->phirar1 != NULL)
memory->destroy(this->phirar1);
if (this->phirar2 != NULL)
memory->destroy(this->phirar2);
if (this->phirar3 != NULL)
memory->destroy(this->phirar3);
if (this->phirar4 != NULL)
memory->destroy(this->phirar4);
if (this->phirar5 != NULL)
memory->destroy(this->phirar5);
if (this->phirar6 != NULL)
memory->destroy(this->phirar6);
// allocate memory for array that defines the potential
memory->create(this->phir, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phir");
// allocate coeff memory
memory->create(this->phirar, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar");
memory->create(this->phirar1, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar1");
memory->create(this->phirar2, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar2");
memory->create(this->phirar3, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar3");
memory->create(this->phirar4, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar4");
memory->create(this->phirar5, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar5");
memory->create(this->phirar6, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar6");
// loop over pairs of element types
nv2 = 0;
for (a = 0; a < this->neltypes; a++) {
for (b = a; b < this->neltypes; b++) {
// loop over r values and compute
for (j = 0; j < this->nr; j++) {
r = j * this->dr;
this->phir[nv2][j] = phi_meam(r, a, b);
// if using second-nearest neighbor, solve recursive problem
// (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
if (this->nn2_meam[a][b] == 1) {
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b], arat, scrn);
// The B1, B2, and L12 cases with NN2 have a trick to them; we
// need to
// compute the contributions from second nearest neighbors, like
// a-a
// pairs, but need to include NN2 contributions to those pairs as
// well.
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
this->lattce_meam[a][b] == L12) {
rarat = r * arat;
// phi_aa
phiaa = phi_meam(rarat, a, a);
Z1 = get_Zij(this->lattce_meam[a][a]);
Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);
}
}
// phi_bb
phibb = phi_meam(rarat, b, b);
Z1 = get_Zij(this->lattce_meam[b][b]);
Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
this->Cmax_meam[b][b][b], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);
}
}
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2) {
// Add contributions to the B1 or B2 potential
Z1 = get_Zij(this->lattce_meam[a][b]);
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b], arat, scrn);
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a],
this->Cmax_meam[b][b][a], arat, scrn2);
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
} else if (this->lattce_meam[a][b] == L12) {
// The L12 case has one last trick; we have to be careful to
// compute
// the correct screening between 2nd-neighbor pairs. 1-1
// second-neighbor pairs are screened by 2 type 1 atoms and
// two type
// 2 atoms. 2-2 second-neighbor pairs are screened by 4 type
// 1
// atoms.
C = 1.0;
get_sijk(C, a, a, a, &s111);
get_sijk(C, a, a, b, &s112);
get_sijk(C, b, b, a, &s221);
S11 = s111 * s111 * s112 * s112;
S22 = pow(s221, 4);
this->phir[nv2][j] = this->phir[nv2][j] - 0.75 * S11 * phiaa - 0.25 * S22 * phibb;
}
} else {
nmax = 10;
for (n = 1; n <= nmax; n++) {
this->phir[nv2][j] =
this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);
}
}
}
// For Zbl potential:
// if astar <= -3
// potential is zbl potential
// else if -3 < astar < -1
// potential is linear combination with zbl potential
// endif
if (this->zbl_meam[a][b] == 1) {
astar = this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0);
if (astar <= -3.0)
this->phir[nv2][j] = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
else if (astar > -3.0 && astar < -1.0) {
frac = fcut(1 - (astar + 1.0) / (-3.0 + 1.0));
phizbl = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
this->phir[nv2][j] = frac * this->phir[nv2][j] + (1 - frac) * phizbl;
}
}
}
// call interpolation
interpolate_meam(nv2);
nv2 = nv2 + 1;
}
}
}
//----------------------------------------------------------------------c
// Compute MEAM pair potential for distance r, element types a and b
//
double
MEAM::phi_meam(double r, int a, int b)
{
/*unused:double a1,a2,a12;*/
double t11av, t21av, t31av, t12av, t22av, t32av;
double G1, G2, s1[3], s2[3], rho0_1, rho0_2;
double Gam1, Gam2, Z1, Z2;
double rhobar1, rhobar2, F1, F2;
double rho01, rho11, rho21, rho31;
double rho02, rho12, rho22, rho32;
double scalfac, phiaa, phibb;
double Eu;
double arat, scrn /*unused:,scrn2*/;
int Z12, errorflag;
int n, nmax, Z1nn, Z2nn;
lattice_t latta /*unused:,lattb*/;
double rho_bkgd1, rho_bkgd2;
double phi_m = 0.0;
// Equation numbers below refer to:
// I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
// get number of neighbors in the reference structure
// Nref[i][j] = # of i's neighbors of type j
Z12 = get_Zij(this->lattce_meam[a][b]);
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32);
// if densities are too small, numerical problems may result; just return zero
if (rho01 <= 1e-14 && rho02 <= 1e-14)
return 0.0;
// calculate average weighting factors for the reference structure
if (this->lattce_meam[a][b] == C11) {
if (this->ialloy == 2) {
t11av = this->t1_meam[a];
t12av = this->t1_meam[b];
t21av = this->t2_meam[a];
t22av = this->t2_meam[b];
t31av = this->t3_meam[a];
t32av = this->t3_meam[b];
} else {
scalfac = 1.0 / (rho01 + rho02);
t11av = scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02);
t12av = t11av;
t21av = scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02);
t22av = t21av;
t31av = scalfac * (this->t3_meam[a] * rho01 + this->t3_meam[b] * rho02);
t32av = t31av;
}
} else {
// average weighting factors for the reference structure, eqn. I.8
get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, this->t1_meam[a], this->t2_meam[a],
this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b,
this->lattce_meam[a][b]);
}
// for c11b structure, calculate background electron densities
if (this->lattce_meam[a][b] == C11) {
latta = this->lattce_meam[a][a];
if (latta == DIA) {
rhobar1 = pow(((Z12 / 2) * (rho02 + rho01)), 2) + t11av * pow((rho12 - rho11), 2) +
t21av / 6.0 * pow(rho22 + rho21, 2) + 121.0 / 40.0 * t31av * pow((rho32 - rho31), 2);
rhobar1 = sqrt(rhobar1);
rhobar2 = pow(Z12 * rho01, 2) + 2.0 / 3.0 * t21av * pow(rho21, 2);
rhobar2 = sqrt(rhobar2);
} else {
rhobar2 = pow(((Z12 / 2) * (rho01 + rho02)), 2) + t12av * pow((rho11 - rho12), 2) +
t22av / 6.0 * pow(rho21 + rho22, 2) + 121.0 / 40.0 * t32av * pow((rho31 - rho32), 2);
rhobar2 = sqrt(rhobar2);
rhobar1 = pow(Z12 * rho02, 2) + 2.0 / 3.0 * t22av * pow(rho22, 2);
rhobar1 = sqrt(rhobar1);
}
} else {
// for other structures, use formalism developed in Huang's paper
//
// composition-dependent scaling, equation I.7
// If using mixing rule for t, apply to reference structure; else
// use precomputed values
if (this->mix_ref_t == 1) {
Z1 = this->Z_meam[a];
Z2 = this->Z_meam[b];
if (this->ibar_meam[a] <= 0)
G1 = 1.0;
else {
get_shpfcn(this->lattce_meam[a][a], s1);
Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);
G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
}
if (this->ibar_meam[b] <= 0)
G2 = 1.0;
else {
get_shpfcn(this->lattce_meam[b][b], s2);
Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);
G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
}
rho0_1 = this->rho0_meam[a] * Z1 * G1;
rho0_2 = this->rho0_meam[b] * Z2 * G2;
}
Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31);
if (rho01 < 1.0e-14)
Gam1 = 0.0;
else
Gam1 = Gam1 / (rho01 * rho01);
Gam2 = (t12av * rho12 + t22av * rho22 + t32av * rho32);
if (rho02 < 1.0e-14)
Gam2 = 0.0;
else
Gam2 = Gam2 / (rho02 * rho02);
G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
if (this->mix_ref_t == 1) {
rho_bkgd1 = rho0_1;
rho_bkgd2 = rho0_2;
} else {
if (this->bkgd_dyn == 1) {
rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a];
rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b];
} else {
rho_bkgd1 = this->rho_ref_meam[a];
rho_bkgd2 = this->rho_ref_meam[b];
}
}
rhobar1 = rho01 / rho_bkgd1 * G1;
rhobar2 = rho02 / rho_bkgd2 * G2;
}
// compute embedding functions, eqn I.5
if (iszero(rhobar1))
F1 = 0.0;
else {
if (this->emb_lin_neg == 1 && rhobar1 <= 0)
F1 = -this->A_meam[a] * this->Ec_meam[a][a] * rhobar1;
else
F1 = this->A_meam[a] * this->Ec_meam[a][a] * rhobar1 * log(rhobar1);
}
if (iszero(rhobar2))
F2 = 0.0;
else {
if (this->emb_lin_neg == 1 && rhobar2 <= 0)
F2 = -this->A_meam[b] * this->Ec_meam[b][b] * rhobar2;
else
F2 = this->A_meam[b] * this->Ec_meam[b][b] * rhobar2 * log(rhobar2);
}
// compute Rose function, I.16
Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],
this->attrac_meam[a][b], this->erose_form);
// calculate the pair energy
if (this->lattce_meam[a][b] == C11) {
latta = this->lattce_meam[a][a];
if (latta == DIA) {
phiaa = phi_meam(r, a, a);
phi_m = (3 * Eu - F2 - 2 * F1 - 5 * phiaa) / Z12;
} else {
phibb = phi_meam(r, b, b);
phi_m = (3 * Eu - F1 - 2 * F2 - 5 * phibb) / Z12;
}
} else if (this->lattce_meam[a][b] == L12) {
phiaa = phi_meam(r, a, a);
// account for second neighbor a-a potential here...
Z1nn = get_Zij(this->lattce_meam[a][a]);
Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
nmax = 10;
if (scrn > 0.0) {
for (n = 1; n <= nmax; n++) {
phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a);
}
}
phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
} else {
//
// potential is computed from Rose function and embedding energy
phi_m = (2 * Eu - F1 - F2) / Z12;
//
}
// if r = 0, just return 0
if (iszero(r)) {
phi_m = 0.0;
}
return phi_m;
}
//----------------------------------------------------------------------c
// Compute background density for reference structure of each element
void
MEAM::compute_reference_density(void)
{
int a, Z, Z2, errorflag;
double gam, Gbar, shp[3];
double rho0, rho0_2nn, arat, scrn;
// loop over element types
for (a = 0; a < this->neltypes; a++) {
Z = (int)this->Z_meam[a];
if (this->ibar_meam[a] <= 0)
Gbar = 1.0;
else {
get_shpfcn(this->lattce_meam[a][a], shp);
gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);
Gbar = G_gam(gam, this->ibar_meam[a], errorflag);
}
// The zeroth order density in the reference structure, with
// equilibrium spacing, is just the number of first neighbors times
// the rho0_meam coefficient...
rho0 = this->rho0_meam[a] * Z;
// ...unless we have unscreened second neighbors, in which case we
// add on the contribution from those (accounting for partial
// screening)
if (this->nn2_meam[a][a] == 1) {
Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a], arat, scrn);
rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
rho0 = rho0 + Z2 * rho0_2nn * scrn;
}
this->rho_ref_meam[a] = rho0 * Gbar;
}
}
//------------------------------------------------------------------------------c
// Average weighting factors for the reference structure
void
MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, double* t22av, double* t32av,
double t11, double t21, double t31, double t12, double t22, double t32, double r, int a,
int b, lattice_t latt)
{
double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/;
// For ialloy = 2, no averaging is done
if (this->ialloy == 2) {
*t11av = t11;
*t21av = t21;
*t31av = t31;
*t12av = t12;
*t22av = t22;
*t32av = t32;
} else {
if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || latt == B1 || latt == DIM || latt == B2) {
// all neighbors are of the opposite type
*t11av = t12;
*t21av = t22;
*t31av = t32;
*t12av = t11;
*t22av = t21;
*t32av = t31;
} else {
a1 = r / this->re_meam[a][a] - 1.0;
a2 = r / this->re_meam[b][b] - 1.0;
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
if (latt == L12) {
rho01 = 8 * rhoa01 + 4 * rhoa02;
*t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
*t12av = t11;
*t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
*t22av = t21;
*t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
*t32av = t31;
} else {
// call error('Lattice not defined in get_tavref.')
}
}
}
}
//------------------------------------------------------------------------------c
void
MEAM::get_sijk(double C, int i, int j, int k, double* sijk)
{
double x;
x = (C - this->Cmin_meam[i][j][k]) / (this->Cmax_meam[i][j][k] - this->Cmin_meam[i][j][k]);
*sijk = fcut(x);
}
//------------------------------------------------------------------------------c
// Calculate density functions, assuming reference configuration
void
MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, double* rho31,
double* rho02, double* rho12, double* rho22, double* rho32)
{
double a1, a2;
double s[3];
lattice_t lat;
int Zij2nn;
double rhoa01nn, rhoa02nn;
double rhoa01, rhoa11, rhoa21, rhoa31;
double rhoa02, rhoa12, rhoa22, rhoa32;
double arat, scrn, denom;
double C, s111, s112, s221, S11, S22;
a1 = r / this->re_meam[a][a] - 1.0;
a2 = r / this->re_meam[b][b] - 1.0;
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1);
rhoa21 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1);
rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1);
rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
rhoa12 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2);
rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2);
rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2);
lat = this->lattce_meam[a][b];
*rho11 = 0.0;
*rho21 = 0.0;
*rho31 = 0.0;
*rho12 = 0.0;
*rho22 = 0.0;
*rho32 = 0.0;
if (lat == FCC) {
*rho01 = 12.0 * rhoa02;
*rho02 = 12.0 * rhoa01;
} else if (lat == BCC) {
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
} else if (lat == B1) {
*rho01 = 6.0 * rhoa02;
*rho02 = 6.0 * rhoa01;
} else if (lat == DIA) {
*rho01 = 4.0 * rhoa02;
*rho02 = 4.0 * rhoa01;
*rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
*rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
} else if (lat == HCP) {
*rho01 = 12 * rhoa02;
*rho02 = 12 * rhoa01;
*rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
*rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
} else if (lat == DIM) {
get_shpfcn(DIM, s);
*rho01 = rhoa02;
*rho02 = rhoa01;
*rho11 = s[0] * rhoa12 * rhoa12;
*rho12 = s[0] * rhoa11 * rhoa11;
*rho21 = s[1] * rhoa22 * rhoa22;
*rho22 = s[1] * rhoa21 * rhoa21;
*rho31 = s[2] * rhoa32 * rhoa32;
*rho32 = s[2] * rhoa31 * rhoa31;
} else if (lat == C11) {
*rho01 = rhoa01;
*rho02 = rhoa02;
*rho11 = rhoa11;
*rho12 = rhoa12;
*rho21 = rhoa21;
*rho22 = rhoa22;
*rho31 = rhoa31;
*rho32 = rhoa32;
} else if (lat == L12) {
*rho01 = 8 * rhoa01 + 4 * rhoa02;
*rho02 = 12 * rhoa01;
if (this->ialloy == 1) {
*rho21 = 8. / 3. * pow(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b], 2);
denom = 8 * rhoa01 * pow(this->t2_meam[a], 2) + 4 * rhoa02 * pow(this->t2_meam[b], 2);
if (denom > 0.)
*rho21 = *rho21 / denom * *rho01;
} else
*rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
} else if (lat == B2) {
*rho01 = 8.0 * rhoa02;
*rho02 = 8.0 * rhoa01;
} else {
// call error('Lattice not defined in get_densref.')
}
if (this->nn2_meam[a][b] == 1) {
Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], arat, scrn);
a1 = arat * r / this->re_meam[a][a] - 1.0;
a2 = arat * r / this->re_meam[b][b] - 1.0;
rhoa01nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
rhoa02nn = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
if (lat == L12) {
// As usual, L12 thinks it's special; we need to be careful computing
// the screening functions
C = 1.0;
get_sijk(C, a, a, a, &s111);
get_sijk(C, a, a, b, &s112);
get_sijk(C, b, b, a, &s221);
S11 = s111 * s111 * s112 * s112;
S22 = pow(s221, 4);
*rho01 = *rho01 + 6 * S11 * rhoa01nn;
*rho02 = *rho02 + 6 * S22 * rhoa02nn;
} else {
// For other cases, assume that second neighbor is of same type,
// first neighbor may be of different type
*rho01 = *rho01 + Zij2nn * scrn * rhoa01nn;
// Assume Zij2nn and arat don't depend on order, but scrn might
Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], arat, scrn);
*rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
}
}
}
void
MEAM::interpolate_meam(int ind)
{
int j;
double drar;
// map to coefficient space
this->nrar = this->nr;
drar = this->dr;
this->rdrar = 1.0 / drar;
// phir interp
for (j = 0; j < this->nrar; j++) {
this->phirar[ind][j] = this->phir[ind][j];
}
this->phirar1[ind][0] = this->phirar[ind][1] - this->phirar[ind][0];
this->phirar1[ind][1] = 0.5 * (this->phirar[ind][2] - this->phirar[ind][0]);
this->phirar1[ind][this->nrar - 2] =
0.5 * (this->phirar[ind][this->nrar - 1] - this->phirar[ind][this->nrar - 3]);
this->phirar1[ind][this->nrar - 1] = 0.0;
for (j = 2; j < this->nrar - 2; j++) {
this->phirar1[ind][j] = ((this->phirar[ind][j - 2] - this->phirar[ind][j + 2]) +
8.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j - 1])) /
12.;
}
for (j = 0; j < this->nrar - 1; j++) {
this->phirar2[ind][j] = 3.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]) -
2.0 * this->phirar1[ind][j] - this->phirar1[ind][j + 1];
this->phirar3[ind][j] = this->phirar1[ind][j] + this->phirar1[ind][j + 1] -
2.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]);
}
this->phirar2[ind][this->nrar - 1] = 0.0;
this->phirar3[ind][this->nrar - 1] = 0.0;
for (j = 0; j < this->nrar; j++) {
this->phirar4[ind][j] = this->phirar1[ind][j] / drar;
this->phirar5[ind][j] = 2.0 * this->phirar2[ind][j] / drar;
this->phirar6[ind][j] = 3.0 * this->phirar3[ind][j] / drar;
}
}
//---------------------------------------------------------------------
// Compute Rose energy function, I.16
//
double
MEAM::compute_phi(double rij, int elti, int eltj)
{
double pp;
int ind, kk;
ind = this->eltind[elti][eltj];
pp = rij * this->rdrar;
kk = (int)pp;
kk = std::min(kk, this->nrar - 2);
pp = pp - kk;
pp = std::min(pp, 1.0);
double result =
((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
this->phirar[ind][kk];
return result;
}

70
src/USER-MEAMC/meam_setup_global.cpp Normal file
View File

@ -0,0 +1,70 @@
#include "meam.h"
#include <math.h>
using namespace LAMMPS_NS;
void
MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar)
{
int i;
double tmplat[maxelt];
this->neltypes = nelt;
for (i = 0; i < nelt; i++) {
this->lattce_meam[i][i] = lat[i];
this->Z_meam[i] = z[i];
this->ielt_meam[i] = ielement[i];
this->alpha_meam[i][i] = alpha[i];
this->beta0_meam[i] = b0[i];
this->beta1_meam[i] = b1[i];
this->beta2_meam[i] = b2[i];
this->beta3_meam[i] = b3[i];
tmplat[i] = alat[i];
this->Ec_meam[i][i] = esub[i];
this->A_meam[i] = asub[i];
this->t0_meam[i] = t0[i];
this->t1_meam[i] = t1[i];
this->t2_meam[i] = t2[i];
this->t3_meam[i] = t3[i];
this->rho0_meam[i] = rozero[i];
this->ibar_meam[i] = ibar[i];
if (this->lattce_meam[i][i] == FCC)
this->re_meam[i][i] = tmplat[i] / sqrt(2.0);
else if (this->lattce_meam[i][i] == BCC)
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
else if (this->lattce_meam[i][i] == HCP)
this->re_meam[i][i] = tmplat[i];
else if (this->lattce_meam[i][i] == DIM)
this->re_meam[i][i] = tmplat[i];
else if (this->lattce_meam[i][i] == DIA)
this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
else {
// error
}
}
// Set some defaults
this->rc_meam = 4.0;
this->delr_meam = 0.1;
setall2d(this->attrac_meam, 0.0);
setall2d(this->repuls_meam, 0.0);
setall3d(this->Cmax_meam, 2.8);
setall3d(this->Cmin_meam, 2.0);
setall2d(this->ebound_meam, pow(2.8, 2) / (4.0 * (2.8 - 1.0)));
setall2d(this->delta_meam, 0.0);
setall2d(this->nn2_meam, 0);
setall2d(this->zbl_meam, 1);
this->gsmooth_factor = 99.0;
this->augt1 = 1;
this->ialloy = 0;
this->mix_ref_t = 0;
this->emb_lin_neg = 0;
this->bkgd_dyn = 0;
this->erose_form = 0;
}

209
src/USER-MEAMC/meam_setup_param.cpp Normal file
View File

@ -0,0 +1,209 @@
#include "meam.h"
#include <algorithm>
using namespace LAMMPS_NS;
//
// do a sanity check on index parameters
void
MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr)
{
//: idx[0..2]
*ierr = 0;
if (nidx < num) {
*ierr = 2;
return;
}
for (int i = 0; i < num; i++) {
if ((idx[i] < 0) || (idx[i] >= lim)) {
*ierr = 3;
return;
}
}
}
// The "which" argument corresponds to the index of the "keyword" array
// in pair_meam.cpp:
//
// 0 = Ec_meam
// 1 = alpha_meam
// 2 = rho0_meam
// 3 = delta_meam
// 4 = lattce_meam
// 5 = attrac_meam
// 6 = repuls_meam
// 7 = nn2_meam
// 8 = Cmin_meam
// 9 = Cmax_meam
// 10 = rc_meam
// 11 = delr_meam
// 12 = augt1
// 13 = gsmooth_factor
// 14 = re_meam
// 15 = ialloy
// 16 = mixture_ref_t
// 17 = erose_form
// 18 = zbl_meam
// 19 = emb_lin_neg
// 20 = bkgd_dyn
void
MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag)
{
//: index[0..2]
int i1, i2;
lattice_t vlat;
*errorflag = 0;
switch (which) {
// 0 = Ec_meam
case 0:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Ec_meam[index[0]][index[1]] = value;
break;
// 1 = alpha_meam
case 1:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->alpha_meam[index[0]][index[1]] = value;
break;
// 2 = rho0_meam
case 2:
meam_checkindex(1, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->rho0_meam[index[0]] = value;
break;
// 3 = delta_meam
case 3:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->delta_meam[index[0]][index[1]] = value;
break;
// 4 = lattce_meam
case 4:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
vlat = (lattice_t)value;
this->lattce_meam[index[0]][index[1]] = vlat;
break;
// 5 = attrac_meam
case 5:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->attrac_meam[index[0]][index[1]] = value;
break;
// 6 = repuls_meam
case 6:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->repuls_meam[index[0]][index[1]] = value;
break;
// 7 = nn2_meam
case 7:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
i1 = std::min(index[0], index[1]);
i2 = std::max(index[0], index[1]);
this->nn2_meam[i1][i2] = (int)value;
break;
// 8 = Cmin_meam
case 8:
meam_checkindex(3, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Cmin_meam[index[0]][index[1]][index[2]] = value;
break;
// 9 = Cmax_meam
case 9:
meam_checkindex(3, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->Cmax_meam[index[0]][index[1]][index[2]] = value;
break;
// 10 = rc_meam
case 10:
this->rc_meam = value;
break;
// 11 = delr_meam
case 11:
this->delr_meam = value;
break;
// 12 = augt1
case 12:
this->augt1 = (int)value;
break;
// 13 = gsmooth
case 13:
this->gsmooth_factor = value;
break;
// 14 = re_meam
case 14:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
this->re_meam[index[0]][index[1]] = value;
break;
// 15 = ialloy
case 15:
this->ialloy = (int)value;
break;
// 16 = mixture_ref_t
case 16:
this->mix_ref_t = (int)value;
break;
// 17 = erose_form
case 17:
this->erose_form = (int)value;
break;
// 18 = zbl_meam
case 18:
meam_checkindex(2, neltypes, nindex, index, errorflag);
if (*errorflag != 0)
return;
i1 = std::min(index[0], index[1]);
i2 = std::max(index[0], index[1]);
this->zbl_meam[i1][i2] = (int)value;
break;
// 19 = emb_lin_neg
case 19:
this->emb_lin_neg = (int)value;
break;
// 20 = bkgd_dyn
case 20:
this->bkgd_dyn = (int)value;
break;
default:
*errorflag = 1;
}
}

782
src/USER-MEAMC/pair_meamc.cpp Normal file
View File

@ -0,0 +1,782 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Greg Wagner (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "meam.h"
#include "pair_meamc.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
static const int nkeywords = 21;
static const char *keywords[] = {
"Ec","alpha","rho0","delta","lattce",
"attrac","repuls","nn2","Cmin","Cmax","rc","delr",
"augt1","gsmooth_factor","re","ialloy",
"mixture_ref_t","erose_form","zbl",
"emb_lin_neg","bkgd_dyn"};
/* ---------------------------------------------------------------------- */
PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
allocated = 0;
nelements = 0;
elements = NULL;
mass = NULL;
meam_inst = new MEAM(memory);
// set comm size needed by this Pair
comm_forward = 38;
comm_reverse = 30;
}
/* ----------------------------------------------------------------------
free all arrays
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairMEAMC::~PairMEAMC()
{
delete meam_inst;
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
}
}
/* ---------------------------------------------------------------------- */
void PairMEAMC::compute(int eflag, int vflag)
{
int i,ii,n,inum_half,errorflag;
int *ilist_half,*numneigh_half,**firstneigh_half;
int *numneigh_full,**firstneigh_full;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
// neighbor list info
inum_half = listhalf->inum;
ilist_half = listhalf->ilist;
numneigh_half = listhalf->numneigh;
firstneigh_half = listhalf->firstneigh;
numneigh_full = listfull->numneigh;
firstneigh_full = listfull->firstneigh;
// strip neighbor lists of any special bond flags before using with MEAM
// necessary before doing neigh_f2c and neigh_c2f conversions each step
if (neighbor->ago == 0) {
neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
}
// check size of scrfcn based on half neighbor list
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
n = 0;
for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
meam_inst->meam_dens_setup(atom->nmax, nall, n);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int ntype = atom->ntypes;
// 3 stages of MEAM calculation
// loop over my atoms followed by communication
int offset = 0;
errorflag = 0;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
meam_inst->meam_dens_init(i,ntype,type,map,x,
numneigh_half[i],firstneigh_half[i],
numneigh_full[i],firstneigh_full[i],
offset);
offset += numneigh_half[i];
}
comm->reverse_comm_pair(this);
meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom,
&eng_vdwl,eatom,ntype,type,map,errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
comm->forward_comm_pair(this);
offset = 0;
// vptr is first value in vatom if it will be used by meam_force()
// else vatom may not exist, so pass dummy ptr
double **vptr;
if (vflag_atom) vptr = vatom;
else vptr = NULL;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom,
vflag_atom,&eng_vdwl,eatom,ntype,type,map,x,
numneigh_half[i],firstneigh_half[i],
numneigh_full[i],firstneigh_full[i],
offset,f,vptr);
offset += numneigh_half[i];
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairMEAMC::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairMEAMC::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairMEAMC::coeff(int narg, char **arg)
{
int i,j,m,n;
if (!allocated) allocate();
if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read MEAM element names between 2 filenames
// nelements = # of MEAM elements
// elements = list of unique element names
if (nelements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
}
nelements = narg - 4 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
mass = new double[nelements];
for (i = 0; i < nelements; i++) {
n = strlen(arg[i+3]) + 1;
elements[i] = new char[n];
strcpy(elements[i],arg[i+3]);
}
// read MEAM library and parameter files
// pass all parameters to MEAM package
// tell MEAM package that setup is done
read_files(arg[2],arg[2+nelements+1]);
meam_inst->meam_setup_done(&cutmax);
// read args that map atom types to MEAM elements
// map[i] = which element the Ith atom type is, -1 if not mapped
for (i = 4 + nelements; i < narg; i++) {
m = i - (4+nelements) + 1;
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
if (j < nelements) map[m] = j;
else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
else error->all(FLERR,"Incorrect args for pair coefficients");
}
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass for i,i in atom class
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairMEAMC::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style MEAM requires newton pair on");
// need full and half neighbor list
int irequest_full = neighbor->request(this,instance_me);
neighbor->requests[irequest_full]->id = 1;
neighbor->requests[irequest_full]->half = 0;
neighbor->requests[irequest_full]->full = 1;
int irequest_half = neighbor->request(this,instance_me);
neighbor->requests[irequest_half]->id = 2;
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
half or full
------------------------------------------------------------------------- */
void PairMEAMC::init_list(int id, NeighList *ptr)
{
if (id == 1) listfull = ptr;
else if (id == 2) listhalf = ptr;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMEAMC::init_one(int i, int j)
{
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairMEAMC::read_files(char *globalfile, char *userfile)
{
// open global meamf file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(globalfile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open MEAM potential file %s",globalfile);
error->one(FLERR,str);
}
}
// allocate parameter arrays
int params_per_line = 19;
lattice_t *lat = new lattice_t[nelements];
int *ielement = new int[nelements];
int *ibar = new int[nelements];
double *z = new double[nelements];
double *atwt = new double[nelements];
double *alpha = new double[nelements];
double *b0 = new double[nelements];
double *b1 = new double[nelements];
double *b2 = new double[nelements];
double *b3 = new double[nelements];
double *alat = new double[nelements];
double *esub = new double[nelements];
double *asub = new double[nelements];
double *t0 = new double[nelements];
double *t1 = new double[nelements];
double *t2 = new double[nelements];
double *t3 = new double[nelements];
double *rozero = new double[nelements];
bool *found = new bool[nelements];
for (int i = 0; i < nelements; i++) found[i] = false;
// read each set of params from global MEAM file
// one set of params can span multiple lines
// store params if element name is in element list
// if element name appears multiple times, only store 1st entry
int i,n,nwords;
char **words = new char*[params_per_line+1];
char line[MAXLINE],*ptr;
int eof = 0;
int nset = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in MEAM potential file");
// words = ptrs to all words in line
// strip single and double quotes from words
nwords = 0;
words[nwords++] = strtok(line,"' \t\n\r\f");
while ((words[nwords++] = strtok(NULL,"' \t\n\r\f"))) continue;
// skip if element name isn't in element list
for (i = 0; i < nelements; i++)
if (strcmp(words[0],elements[i]) == 0) break;
if (i >= nelements) continue;
// skip if element already appeared
if (found[i] == true) continue;
found[i] = true;
// map lat string to an integer
if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
else if (strcmp(words[1],"dia") == 0) lat[i] = DIA;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
// store parameters
z[i] = atof(words[2]);
ielement[i] = atoi(words[3]);
atwt[i] = atof(words[4]);
alpha[i] = atof(words[5]);
b0[i] = atof(words[6]);
b1[i] = atof(words[7]);
b2[i] = atof(words[8]);
b3[i] = atof(words[9]);
alat[i] = atof(words[10]);
esub[i] = atof(words[11]);
asub[i] = atof(words[12]);
t0[i] = atof(words[13]);
t1[i] = atof(words[14]);
t2[i] = atof(words[15]);
t3[i] = atof(words[16]);
rozero[i] = atof(words[17]);
ibar[i] = atoi(words[18]);
nset++;
}
// error if didn't find all elements in file
if (nset != nelements)
error->all(FLERR,"Did not find all elements in MEAM library file");
// pass element parameters to MEAM package
meam_inst->meam_setup_global(nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
// set element masses
for (i = 0; i < nelements; i++) mass[i] = atwt[i];
// clean-up memory
delete [] words;
delete [] lat;
delete [] ielement;
delete [] ibar;
delete [] z;
delete [] atwt;
delete [] alpha;
delete [] b0;
delete [] b1;
delete [] b2;
delete [] b3;
delete [] alat;
delete [] esub;
delete [] asub;
delete [] t0;
delete [] t1;
delete [] t2;
delete [] t3;
delete [] rozero;
delete [] found;
// done if user param file is NULL
if (strcmp(userfile,"NULL") == 0) return;
// open user param file on proc 0
if (comm->me == 0) {
fp = force->open_potential(userfile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open MEAM potential file %s",userfile);
error->one(FLERR,str);
}
}
// read settings
// pass them one at a time to MEAM package
// match strings to list of corresponding ints
int which;
double value;
int nindex,index[3];
int maxparams = 6;
char **params = new char*[maxparams];
int nparams;
eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nparams = atom->count_words(line);
if (nparams == 0) continue;
// words = ptrs to all words in line
nparams = 0;
params[nparams++] = strtok(line,"=(), '\t\n\r\f");
while (nparams < maxparams &&
(params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
continue;
nparams--;
for (which = 0; which < nkeywords; which++)
if (strcmp(params[0],keywords[which]) == 0) break;
if (which == nkeywords) {
char str[128];
sprintf(str,"Keyword %s in MEAM parameter file not recognized",
params[0]);
error->all(FLERR,str);
}
nindex = nparams - 2;
for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]) - 1;
// map lattce_meam value to an integer
if (which == 4) {
if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
else if (strcmp(params[nparams-1],"dia") == 0) value = DIA;
else if (strcmp(params[nparams-1],"b1") == 0) value = B1;
else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
else if (strcmp(params[nparams-1],"b2") == 0) value = B2;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
}
else value = atof(params[nparams-1]);
// pass single setting to MEAM package
int errorflag = 0;
meam_inst->meam_setup_param(which,value,nindex,index,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->all(FLERR,str);
}
}
delete [] params;
}
/* ---------------------------------------------------------------------- */
int PairMEAMC::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,k,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = meam_inst->rho0[j];
buf[m++] = meam_inst->rho1[j];
buf[m++] = meam_inst->rho2[j];
buf[m++] = meam_inst->rho3[j];
buf[m++] = meam_inst->frhop[j];
buf[m++] = meam_inst->gamma[j];
buf[m++] = meam_inst->dgamma1[j];
buf[m++] = meam_inst->dgamma2[j];
buf[m++] = meam_inst->dgamma3[j];
buf[m++] = meam_inst->arho2b[j];
buf[m++] = meam_inst->arho1[j][0];
buf[m++] = meam_inst->arho1[j][1];
buf[m++] = meam_inst->arho1[j][2];
buf[m++] = meam_inst->arho2[j][0];
buf[m++] = meam_inst->arho2[j][1];
buf[m++] = meam_inst->arho2[j][2];
buf[m++] = meam_inst->arho2[j][3];
buf[m++] = meam_inst->arho2[j][4];
buf[m++] = meam_inst->arho2[j][5];
for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[j][k];
buf[m++] = meam_inst->arho3b[j][0];
buf[m++] = meam_inst->arho3b[j][1];
buf[m++] = meam_inst->arho3b[j][2];
buf[m++] = meam_inst->t_ave[j][0];
buf[m++] = meam_inst->t_ave[j][1];
buf[m++] = meam_inst->t_ave[j][2];
buf[m++] = meam_inst->tsq_ave[j][0];
buf[m++] = meam_inst->tsq_ave[j][1];
buf[m++] = meam_inst->tsq_ave[j][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairMEAMC::unpack_forward_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
meam_inst->rho0[i] = buf[m++];
meam_inst->rho1[i] = buf[m++];
meam_inst->rho2[i] = buf[m++];
meam_inst->rho3[i] = buf[m++];
meam_inst->frhop[i] = buf[m++];
meam_inst->gamma[i] = buf[m++];
meam_inst->dgamma1[i] = buf[m++];
meam_inst->dgamma2[i] = buf[m++];
meam_inst->dgamma3[i] = buf[m++];
meam_inst->arho2b[i] = buf[m++];
meam_inst->arho1[i][0] = buf[m++];
meam_inst->arho1[i][1] = buf[m++];
meam_inst->arho1[i][2] = buf[m++];
meam_inst->arho2[i][0] = buf[m++];
meam_inst->arho2[i][1] = buf[m++];
meam_inst->arho2[i][2] = buf[m++];
meam_inst->arho2[i][3] = buf[m++];
meam_inst->arho2[i][4] = buf[m++];
meam_inst->arho2[i][5] = buf[m++];
for (k = 0; k < 10; k++) meam_inst->arho3[i][k] = buf[m++];
meam_inst->arho3b[i][0] = buf[m++];
meam_inst->arho3b[i][1] = buf[m++];
meam_inst->arho3b[i][2] = buf[m++];
meam_inst->t_ave[i][0] = buf[m++];
meam_inst->t_ave[i][1] = buf[m++];
meam_inst->t_ave[i][2] = buf[m++];
meam_inst->tsq_ave[i][0] = buf[m++];
meam_inst->tsq_ave[i][1] = buf[m++];
meam_inst->tsq_ave[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairMEAMC::pack_reverse_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = meam_inst->rho0[i];
buf[m++] = meam_inst->arho2b[i];
buf[m++] = meam_inst->arho1[i][0];
buf[m++] = meam_inst->arho1[i][1];
buf[m++] = meam_inst->arho1[i][2];
buf[m++] = meam_inst->arho2[i][0];
buf[m++] = meam_inst->arho2[i][1];
buf[m++] = meam_inst->arho2[i][2];
buf[m++] = meam_inst->arho2[i][3];
buf[m++] = meam_inst->arho2[i][4];
buf[m++] = meam_inst->arho2[i][5];
for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[i][k];
buf[m++] = meam_inst->arho3b[i][0];
buf[m++] = meam_inst->arho3b[i][1];
buf[m++] = meam_inst->arho3b[i][2];
buf[m++] = meam_inst->t_ave[i][0];
buf[m++] = meam_inst->t_ave[i][1];
buf[m++] = meam_inst->t_ave[i][2];
buf[m++] = meam_inst->tsq_ave[i][0];
buf[m++] = meam_inst->tsq_ave[i][1];
buf[m++] = meam_inst->tsq_ave[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,k,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
meam_inst->rho0[j] += buf[m++];
meam_inst->arho2b[j] += buf[m++];
meam_inst->arho1[j][0] += buf[m++];
meam_inst->arho1[j][1] += buf[m++];
meam_inst->arho1[j][2] += buf[m++];
meam_inst->arho2[j][0] += buf[m++];
meam_inst->arho2[j][1] += buf[m++];
meam_inst->arho2[j][2] += buf[m++];
meam_inst->arho2[j][3] += buf[m++];
meam_inst->arho2[j][4] += buf[m++];
meam_inst->arho2[j][5] += buf[m++];
for (k = 0; k < 10; k++) meam_inst->arho3[j][k] += buf[m++];
meam_inst->arho3b[j][0] += buf[m++];
meam_inst->arho3b[j][1] += buf[m++];
meam_inst->arho3b[j][2] += buf[m++];
meam_inst->t_ave[j][0] += buf[m++];
meam_inst->t_ave[j][1] += buf[m++];
meam_inst->t_ave[j][2] += buf[m++];
meam_inst->tsq_ave[j][0] += buf[m++];
meam_inst->tsq_ave[j][1] += buf[m++];
meam_inst->tsq_ave[j][2] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairMEAMC::memory_usage()
{
double bytes = 11 * meam_inst->nmax * sizeof(double);
bytes += (3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double);
bytes += 3 * meam_inst->maxneigh * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
strip special bond flags from neighbor list entries
are not used with MEAM
need to do here so Fortran lib doesn't see them
done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
------------------------------------------------------------------------- */
void PairMEAMC::neigh_strip(int inum, int *ilist,
int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
}
}

112
src/USER-MEAMC/pair_meamc.h Normal file
View File

@ -0,0 +1,112 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(meam/c,PairMEAMC)
#else
#ifndef LMP_PAIR_MEAMC_H
#define LMP_PAIR_MEAMC_H
#include "pair.h"
namespace LAMMPS_NS {
class MEAM;
class PairMEAMC : public Pair {
public:
PairMEAMC(class LAMMPS *);
~PairMEAMC();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
private:
class MEAM *meam_inst;
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
double *mass; // mass of each element
int *map; // mapping from atom types (1-indexed) to elements (1-indexed)
void allocate();
void read_files(char *, char *);
void neigh_strip(int, int *, int *, int **);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: MEAM library error %d
A call to the MEAM Fortran library returned an error.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style MEAM requires newton pair on
See the newton command. This is a restriction to use the MEAM
potential.
E: Cannot open MEAM potential file %s
The specified MEAM potential file cannot be opened. Check that the
path and name are correct.
E: Incorrect format in MEAM potential file
Incorrect number of words per line in the potential file.
E: Unrecognized lattice type in MEAM file 1
The lattice type in an entry of the MEAM library file is not
valid.
E: Did not find all elements in MEAM library file
The requested elements were not found in the MEAM file.
E: Keyword %s in MEAM parameter file not recognized
Self-explanatory.
E: Unrecognized lattice type in MEAM file 2
The lattice type in an entry of the MEAM parameter file is not
valid.
*/

121
src/USER-OMP/pair_agni_omp.cpp
View File

@ -29,127 +29,6 @@
using namespace LAMMPS_NS;
using namespace MathSpecial;
/*
Copyright (c) 2012,2013 Axel Kohlmeyer <akohlmey@gmail.com>
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the <organization> nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/* faster versions of 2**x, e**x, and 10**x in single and double precision.
*
* Based on the Cephes math library 2.8
*/
/* internal definitions for the fastermath library */
/* IEEE 754 double precision floating point data manipulation */
typedef union
{
double f;
uint64_t u;
struct {int32_t i0,i1;};
} udi_t;
#define FM_DOUBLE_BIAS 1023
#define FM_DOUBLE_EMASK 2146435072
#define FM_DOUBLE_MBITS 20
#define FM_DOUBLE_MMASK 1048575
#define FM_DOUBLE_EZERO 1072693248
/* generate 2**num in floating point by bitshifting */
#define FM_DOUBLE_INIT_EXP(var,num) \
var.i0 = 0; \
var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20
/* double precision constants */
#define FM_DOUBLE_LOG2OFE 1.4426950408889634074
#define FM_DOUBLE_LOGEOF2 6.9314718055994530942e-1
#define FM_DOUBLE_LOG2OF10 3.32192809488736234789
#define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1
#define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1
#define FM_DOUBLE_SQRT2 1.41421356237309504880
#define FM_DOUBLE_SQRTH 0.70710678118654752440
/* optimizer friendly implementation of exp2(x).
*
* strategy:
*
* split argument into an integer part and a fraction:
* ipart = floor(x+0.5);
* fpart = x - ipart;
*
* compute exp2(ipart) from setting the ieee754 exponent
* compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[
*
* the result becomes: exp2(x) = exp2(ipart) * exp2(fpart)
*/
static const double fm_exp2_q[] = {
/* 1.00000000000000000000e0, */
2.33184211722314911771e2,
4.36821166879210612817e3
};
static const double fm_exp2_p[] = {
2.30933477057345225087e-2,
2.02020656693165307700e1,
1.51390680115615096133e3
};
static double fm_exp2(double x)
{
double ipart, fpart, px, qx;
udi_t epart;
ipart = floor(x+0.5);
fpart = x - ipart;
FM_DOUBLE_INIT_EXP(epart,ipart);
x = fpart*fpart;
px = fm_exp2_p[0];
px = px*x + fm_exp2_p[1];
qx = x + fm_exp2_q[0];
px = px*x + fm_exp2_p[2];
qx = qx*x + fm_exp2_q[1];
px = px * fpart;
x = 1.0 + 2.0*(px/(qx-px));
return epart.f*x;
}
static double fm_exp(double x)
{
#if defined(__BYTE_ORDER__)
#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
return fm_exp2(FM_DOUBLE_LOG2OFE * (x));
#endif
#endif
return exp(x);
}
/* ---------------------------------------------------------------------- */
PairAGNIOMP::PairAGNIOMP(LAMMPS *lmp) :

14
src/math_special.cpp
View File

@ -508,6 +508,9 @@ static const double fm_exp2_p[] = {
1.51390680115615096133e3
};
/* double precision constants */
#define FM_DOUBLE_LOG2OFE 1.4426950408889634074
double MathSpecial::exp2_x86(double x)
{
double ipart, fpart, px, qx;
@ -531,3 +534,14 @@ double MathSpecial::exp2_x86(double x)
x = 1.0 + 2.0*(px/(qx-px));
return epart.f*x;
}
double MathSpecial::fm_exp(double x)
{
#if defined(__BYTE_ORDER__)
#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
return exp2_x86(FM_DOUBLE_LOG2OFE * x);
#endif
#endif
return ::exp(x);
}

3
src/math_special.h
View File

@ -27,6 +27,9 @@ namespace MathSpecial {
// fast 2**x function without argument checks for little endian CPUs
extern double exp2_x86(double x);
// fast e**x function for little endian CPUs, falls back to libc on other platforms
extern double fm_exp(double x);
// scaled error function complement exp(x*x)*erfc(x) for coul/long styles
static inline double my_erfcx(const double x)