diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index 0fbab732c8..e17645e6d0 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -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, diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt index 24506379c3..18030162cf 100644 --- a/doc/src/Section_packages.txt +++ b/doc/src/Section_packages.txt @@ -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:] diff --git a/doc/src/pair_meam.txt b/doc/src/pair_meam.txt index 4fcb7a2e6c..62fa59f406 100644 --- a/doc/src/pair_meam.txt +++ b/doc/src/pair_meam.txt @@ -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:] diff --git a/examples/meam/in.meamc b/examples/meam/in.meamc new file mode 100644 index 0000000000..f6815cd7d4 --- /dev/null +++ b/examples/meam/in.meamc @@ -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 diff --git a/examples/meam/in.meamc.shear b/examples/meam/in.meamc.shear new file mode 100644 index 0000000000..e4584d9744 --- /dev/null +++ b/examples/meam/in.meamc.shear @@ -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 diff --git a/examples/meam/log.5Oct16.meam.icc.1 b/examples/meam/log.19May17.meam.g++.1 similarity index 68% rename from examples/meam/log.5Oct16.meam.icc.1 rename to examples/meam/log.19May17.meam.g++.1 index 200da68e06..a98de97f8e 100644 --- a/examples/meam/log.5Oct16.meam.icc.1 +++ b/examples/meam/log.19May17.meam.g++.1 @@ -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 diff --git a/examples/meam/log.5Oct16.meam.icc.4 b/examples/meam/log.19May17.meam.g++.4 similarity index 68% rename from examples/meam/log.5Oct16.meam.icc.4 rename to examples/meam/log.19May17.meam.g++.4 index 51a6619e3b..adc34d08aa 100644 --- a/examples/meam/log.5Oct16.meam.icc.4 +++ b/examples/meam/log.19May17.meam.g++.4 @@ -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 diff --git a/examples/meam/log.5Oct16.meam.shear.icc.1 b/examples/meam/log.19May17.meam.shear.g++.1 similarity index 55% rename from examples/meam/log.5Oct16.meam.shear.icc.1 rename to examples/meam/log.19May17.meam.shear.g++.1 index 57f48d5ee2..77b9688455 100644 --- a/examples/meam/log.5Oct16.meam.shear.icc.1 +++ b/examples/meam/log.19May17.meam.shear.g++.1 @@ -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 diff --git a/examples/meam/log.5Oct16.meam.shear.icc.4 b/examples/meam/log.19May17.meam.shear.g++.4 similarity index 55% rename from examples/meam/log.5Oct16.meam.shear.icc.4 rename to examples/meam/log.19May17.meam.shear.g++.4 index 2f197de920..84cb94f4b9 100644 --- a/examples/meam/log.5Oct16.meam.shear.icc.4 +++ b/examples/meam/log.19May17.meam.shear.g++.4 @@ -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 diff --git a/examples/meam/log.19May17.meamc.g++.1 b/examples/meam/log.19May17.meamc.g++.1 new file mode 100644 index 0000000000..4aa1f7bcc3 --- /dev/null +++ b/examples/meam/log.19May17.meamc.g++.1 @@ -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 diff --git a/examples/meam/log.19May17.meamc.g++.4 b/examples/meam/log.19May17.meamc.g++.4 new file mode 100644 index 0000000000..3701fb80d8 --- /dev/null +++ b/examples/meam/log.19May17.meamc.g++.4 @@ -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 diff --git a/examples/meam/log.19May17.meamc.shear.g++.1 b/examples/meam/log.19May17.meamc.shear.g++.1 new file mode 100644 index 0000000000..5eb06592d8 --- /dev/null +++ b/examples/meam/log.19May17.meamc.shear.g++.1 @@ -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 diff --git a/examples/meam/log.19May17.meamc.shear.g++.4 b/examples/meam/log.19May17.meamc.shear.g++.4 new file mode 100644 index 0000000000..136459a44a --- /dev/null +++ b/examples/meam/log.19May17.meamc.shear.g++.4 @@ -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 diff --git a/src/.gitignore b/src/.gitignore index 3bdf64a4cb..ae5b3d957c 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/Makefile b/src/Makefile index 92a430a747..c7b20dcb13 100644 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/USER-MEAMC/README b/src/USER-MEAMC/README new file mode 100644 index 0000000000..c1faf7c0c4 --- /dev/null +++ b/src/USER-MEAMC/README @@ -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. + diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h new file mode 100644 index 0000000000..944b3fdf28 --- /dev/null +++ b/src/USER-MEAMC/meam.h @@ -0,0 +1,252 @@ +#ifndef LMP_MEAM_H +#define LMP_MEAM_H + +#include "memory.h" +#include +#include + +#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 +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 +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 diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp new file mode 100644 index 0000000000..2adf4000f1 --- /dev/null +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -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; + } + } + } + } +} + diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp new file mode 100644 index 0000000000..e1a7509ab3 --- /dev/null +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -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 0cutforcesq) { + 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; + } + } + } + } + } + } +} + diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp new file mode 100644 index 0000000000..cbed31fd5b --- /dev/null +++ b/src/USER-MEAMC/meam_force.cpp @@ -0,0 +1,531 @@ +#include "meam.h" +#include "math_special.h" +#include + +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 + } +} diff --git a/src/USER-MEAMC/meam_funcs.cpp b/src/USER-MEAMC/meam_funcs.cpp new file mode 100644 index 0000000000..3c0dcb9d0b --- /dev/null +++ b/src/USER-MEAMC/meam_funcs.cpp @@ -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 + +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; +} diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp new file mode 100644 index 0000000000..25c53785e4 --- /dev/null +++ b/src/USER-MEAMC/meam_impl.cpp @@ -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); +} diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp new file mode 100644 index 0000000000..5ef6595253 --- /dev/null +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -0,0 +1,792 @@ +#include "meam.h" +#include "math_special.h" +#include +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) { + 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 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 ineltypes; 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; +} diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp new file mode 100644 index 0000000000..7062d5ebca --- /dev/null +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -0,0 +1,70 @@ +#include "meam.h" +#include +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; +} diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp new file mode 100644 index 0000000000..585b5b5712 --- /dev/null +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -0,0 +1,209 @@ +#include "meam.h" +#include +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; + } +} diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp new file mode 100644 index 0000000000..e35c54352e --- /dev/null +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -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 +#include +#include +#include +#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; + } +} diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h new file mode 100644 index 0000000000..476a70dd04 --- /dev/null +++ b/src/USER-MEAMC/pair_meamc.h @@ -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. + +*/ diff --git a/src/USER-OMP/pair_agni_omp.cpp b/src/USER-OMP/pair_agni_omp.cpp index 20bb3987bd..72e8331718 100644 --- a/src/USER-OMP/pair_agni_omp.cpp +++ b/src/USER-OMP/pair_agni_omp.cpp @@ -29,127 +29,6 @@ using namespace LAMMPS_NS; using namespace MathSpecial; -/* - Copyright (c) 2012,2013 Axel Kohlmeyer - 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 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 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) : diff --git a/src/math_special.cpp b/src/math_special.cpp index 39487dd386..6777123ac8 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -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); +} + diff --git a/src/math_special.h b/src/math_special.h index e4b4998d54..8cd328f5fc 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -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)