lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #466 from DallasTrinkle/meam-spline-multicomponent 

Meam spline multicomponent
This commit is contained in:
sjplimp 2017-05-11 09:22:25 -06:00 committed by GitHub
parent a6e2d5b5f7 d0da0639f0
commit deff6ffaac
17 changed files with 1725 additions and 458 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/Eqs/pair_meam_spline.jpg
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 10 KiB

After

Width:  |  Height:  |  Size: 21 KiB

5
doc/src/Eqs/pair_meam_spline.tex
View File

@ -1,13 +1,14 @@
\documentclass[12pt]{article}
\usepackage{amsmath}
\begin{document}
$$
E=\sum_{ij}\phi(r_{ij})+\sum_{i}U(\rho_{i}),
E=\sum_{i<j}\phi(r_{ij})+\sum_{i}U(n_{i}),
$$
$$
\rho_{i}=\sum_{j}\rho(r_{ij})+\sum_{jk}f(r_{ij})f(r_{ik})g[\cos(\theta_{jik})]
n_{i}=\sum_{j}\rho(r_{ij})+\sum_{\substack{j<k,\\j,k\neq i}}f(r_{ij})f(r_{ik})g[\cos(\theta_{jik})]
$$
\end{document}

BIN
doc/src/Eqs/pair_meam_spline_multicomponent.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 22 KiB

14
doc/src/Eqs/pair_meam_spline_multicomponent.tex Normal file
View File

@ -0,0 +1,14 @@
\documentclass[12pt]{article}
\usepackage{amsmath}
\begin{document}
$$
E=\sum_{i<j}\phi_{ij}(r_{ij})+\sum_{i}U_i(n_{i}),
$$
$$
n_{i}=\sum_{j\ne i}\rho_j(r_{ij})+\sum_{\substack{j<k,\\j,k\neq i}}f_{j}(r_{ij})f_{k}(r_{ik})g_{jk}[\cos(\theta_{jik})]
$$
\end{document}

40
doc/src/pair_meam_spline.txt
View File

@ -23,7 +23,8 @@ pair_coeff * * Ti.meam.spline Ti Ti Ti :pre
The {meam/spline} style computes pairwise interactions for metals
using a variant of modified embedded-atom method (MEAM) potentials
"(Lenosky)"_#Lenosky1. The total energy E is given by
"(Lenosky)"_#Lenosky1. For a single species ("old-style") MEAM,
the total energy E is given by
:c,image(Eqs/pair_meam_spline.jpg)
@ -31,6 +32,20 @@ where rho_i is the density at atom I, theta_jik is the angle between
atoms J, I, and K centered on atom I. The five functions Phi, U, rho,
f, and g are represented by cubic splines.
The {meam/spline} style also supports a new style multicomponent
modified embedded-atom method (MEAM) potential "(Zhang)"_#Zhang1, where
the total energy E is given by
:c,image(Eqs/pair_meam_spline_multicomponent.jpg)
where the five functions Phi, U, rho, f, and g depend on the chemistry
of the atoms in the interaction. In particular, if there are N different
chemistries, there are N different U, rho, and f functions, while there
are N(N+1)/2 different Phi and g functions. The new style multicomponent
MEAM potential files are indicated by the second line in the file starts
with "meam/spline" followed by the number of elements and the name of each
element.
The cutoffs and the coefficients for these spline functions are listed
in a parameter file which is specified by the
"pair_coeff"_pair_coeff.html command. Parameter files for different
@ -59,7 +74,7 @@ N element names = mapping of spline-based MEAM elements to atom types :ul
See the "pair_coeff"_pair_coeff.html doc page for alternate ways
to specify the path for the potential file.
As an example, imagine the Ti.meam.spline file has values for Ti. If
As an example, imagine the Ti.meam.spline file has values for Ti (old style). If
your LAMMPS simulation has 3 atoms types and they are all to be
treated with this potentials, you would use the following pair_coeff
command:
@ -72,10 +87,19 @@ in the potential file. If a mapping value is specified as NULL, the
mapping is not performed. This can be used when a {meam/spline}
potential is used as part of the {hybrid} pair style. The NULL values
are placeholders for atom types that will be used with other
potentials.
potentials. The old-style potential maps any non-NULL species named
on the command line to that single type.
NOTE: The {meam/spline} style currently supports only single-element
MEAM potentials. It may be extended for alloy systems in the future.
An example with a two component spline (new style) is TiO.meam.spline, where
the command
pair_coeff * * TiO.meam.spline Ti O :pre
will map the 1st atom type to Ti and the second atom type to O. Note
in this case that the species names need to match exactly with the
names of the elements in the TiO.meam.spline file; otherwise an
error will be raised. This behavior is different than the old style
MEAM files.
:line
@ -104,9 +128,6 @@ more instructions on how to use the accelerated styles effectively.
[Mixing, shift, table, tail correction, restart, rRESPA info]:
The current version of this pair style does not support multiple
element types or mixing. It has been designed for pure elements only.
This pair style does not support the "pair_modify"_pair_modify.html
shift, table, and tail options.
@ -142,3 +163,6 @@ for more info.
[(Lenosky)] Lenosky, Sadigh, Alonso, Bulatov, de la Rubia, Kim, Voter,
Kress, Modelling Simulation Materials Science Engineering, 8, 825
(2000).
:link(Zhang1)
[(Zhang)] Zhang and Trinkle, Computational Materials Science, 124, 204-210 (2016).

63
examples/USER/misc/meam_spline/Si_1.meam.spline Normal file
View File

@ -0,0 +1,63 @@
DATE: 2012-02-01 CONTRIBUTOR: Alexander Stukowski, stukowski@mm.tu-darmstadt.de CITATION: Lenosky, Sadigh, Alonso, Bulatov, de la Rubia, Kim, Voter and Kress, Modell Simul Mater Sci Eng, 8, 825 (2000) COMMENT: Spline-based MEAM potential for Si. Reference: T. J. Lenosky, B. Sadigh, E. Alonso, V. V. Bulatov, T. D. de la Rubia, J. Kim, A. F. Voter, and J. D. Kress, Modell. Simul. Mater. Sci. Eng. 8, 825 (2000)
10
-4.266966781858503300e+01 0.000000000000000000e+00
1 0 1 0
1.500000000000000000e+00 6.929943430771341000e+00 1.653321602557917600e+02
1.833333333333333300e+00 -4.399503747408950400e-01 3.941543472528634600e+01
2.166666666666666500e+00 -1.701233725061446700e+00 6.871065423413908100e+00
2.500000000000000000e+00 -1.624732919215791800e+00 5.340648014033163800e+00
2.833333333333333000e+00 -9.969641728342462100e-01 1.534811309391571000e+00
3.166666666666667000e+00 -2.739141845072665100e-01 -6.334706186546093900e+00
3.500000000000000000e+00 -2.499156963774082700e-02 -1.798864729909626500e+00
3.833333333333333500e+00 -1.784331481529976400e-02 4.743496636420091500e-01
4.166666666666666100e+00 -9.612303290166881000e-03 -4.006506271304824400e-02
4.500000000000000000e+00 0.000000000000000000e+00 -2.394996574779807200e-01
11
-1.000000000000000000e+00 0.000000000000000000e+00
1 0 0 0
1.500000000000000000e+00 1.374674212682983900e-01 -3.227795813279568500e+00
1.700000000000000000e+00 -1.483141815327918000e-01 -6.411648793604404900e+00
1.899999999999999900e+00 -5.597204896096039700e-01 1.003068519633888300e+01
2.100000000000000100e+00 -7.310964379372824100e-01 2.293461970618954700e+00
2.299999999999999800e+00 -7.628287071954063000e-01 1.742018781618444500e+00
2.500000000000000000e+00 -7.291769685066557000e-01 5.460640949384478700e-01
2.700000000000000200e+00 -6.662022220044453400e-01 4.721760106467195500e-01
2.899999999999999900e+00 -5.732830582550895200e-01 2.056894449546524200e+00
3.100000000000000100e+00 -4.069014309729406300e-01 2.319615721086100800e+00
3.299999999999999800e+00 -1.666155295956388300e-01 -2.497162196179187900e-01
3.500000000000000000e+00 0.000000000000000000e+00 -1.237130660986393100e+01
8
7.351364478015182100e-01 6.165217237728655200e-01
1 1 1 1
-1.770934559908718700e+00 -1.074925682941420000e+00 -1.482768170233858500e-01
-3.881557649503457600e-01 -2.004503493658201000e-01 -1.492100354067345500e-01
9.946230300080272100e-01 4.142241371345077300e-01 -7.012475119623896900e-02
2.377401824966400000e+00 8.793892953828742500e-01 -3.944355024164965900e-02
3.760180619924772900e+00 1.266888024536562100e+00 -1.581431192239436000e-02
5.142959414883146800e+00 1.629979548834614900e+00 2.611224310900800400e-02
6.525738209841518900e+00 1.977379549636293600e+00 -1.378738550324104500e-01
7.908517004799891800e+00 2.396177220616657200e+00 7.494253977092666400e-01
10
-3.618936018538757300e+00 0.000000000000000000e+00
1 0 1 0
1.500000000000000000e+00 1.250311510312851300e+00 2.790400588857243500e+01
1.722222222222222300e+00 8.682060369372680600e-01 -4.522554291731776900e+00
1.944444444444444400e+00 6.084604017544847900e-01 5.052931618779816800e+00
2.166666666666666500e+00 4.875624808097850400e-01 1.180825096539679600e+00
2.388888888888888800e+00 4.416345603457190700e-01 -6.673769465415171400e-01
2.611111111111111200e+00 3.760976313325982700e-01 -8.938118490837722000e-01
2.833333333333333000e+00 2.714524157414608400e-01 -5.090324763524399800e-01
3.055555555555555400e+00 1.481440300150710900e-01 6.623665830603995300e-01
3.277777777777777700e+00 4.854596610856590900e-02 7.403702452268122700e-01
3.500000000000000000e+00 0.000000000000000000e+00 2.578982318481970500e+00
8
-1.395041572145673000e+01 1.134616739799360700e+00
1 1 1 1
-1.000000000000000900e+00 5.254163992149617700e+00 1.582685381253900500e+01
-7.428367052748285900e-01 2.359149452448745100e+00 3.117611233789983400e+01
-4.856734105496561800e-01 1.195946960915646100e+00 1.658962813584905800e+01
-2.285101158244838800e-01 1.229952028074150000e+00 1.108360928564026400e+01
2.865317890068852500e-02 2.035650777568434500e+00 9.088861456447702400e+00
2.858164736258610400e-01 3.424741418405580000e+00 5.489943377538379500e+00
5.429797683510331200e-01 4.948585892304984100e+00 -1.882291580187675700e+01
8.001430630762056400e-01 5.617988713941801200e+00 -7.718625571850646200e+00

130
examples/USER/misc/meam_spline/TiO.meam.spline Normal file
View File

@ -0,0 +1,130 @@
# Ti-O cubic spline potential where O is in the dilute limit. DATE: 2016-06-05 CONTRIBUTOR: Pinchao Zhang, Dallas R. Trinkle
meam/spline 2 Ti O
spline3eq
13
-20 0
1.742692837 3.744277175966 99.4865081627958
2.05580176725 0.910839730906 10.8702523265355
2.3689106975 0.388045896634 -1.55322418749562
2.68201962775 -0.018840906533 2.43630041329215
2.995128558 -0.248098929639 2.67912713976835
3.30823748825 -0.264489550297 -0.125056384603077
3.6213464185 -0.227196189283 1.10662555360438
3.93445534875 -0.129293090176 -0.592053676745914
4.247564279 -0.059685366933 -0.470123414607672
4.56067320925 -0.031100025561 -0.0380739973059663
4.8737821395 -0.013847363202 -0.0711547960695406
5.18689106975 -0.003203412728 -0.081768292420175
5.5 0 -0.0571422964883619
spline3eq
5
0.155001355787331 0
1.9 0.533321679606674 0
2.8 0.456402081843862 -1.60311717015859
3.7 -0.324281383502201 1.19940299483249
4.6 -0.474029826906675 1.47909794595154
5.5 0 -2.49521499855605
spline3eq
13
0 0
1.742692837 0 0
2.05580176725 0 0
2.3689106975 0 0
2.68201962775 0 0
2.995128558 0 0
3.30823748825 0 0
3.6213464185 0 0
3.93445534875 0 0
4.247564279 0 0
4.56067320925 0 0
4.8737821395 0 0
5.18689106975 0 0
5.5 0 0
spline3eq
11
-1 0
2.055801767 1.7475279661 -525.869786904802
2.2912215903 -5.8677963945 252.796316927755
2.5266414136 -8.3376288737 71.7318388721015
2.7620612369 -5.8398712842 -1.93587742753693
2.9974810602 -3.1140648231 -39.2999192667503
3.2329008835 -1.7257245065 14.3424136002004
3.4683207068 -0.4428977017 -29.4925534559498
3.7037405301 -0.1466643003 -3.18010534572236
3.9391603534 -0.2095507945 3.33490838803603
4.1745801767 -0.1442384563 3.71918691359508
4.41 0 -9.66717019857564
spline3eq
5
-61.9827585211652 0
1.9 11.2293641315584 0
2.8 -27.9976343076148 122.648031332411
3.7 -8.32979773113248 -54.3340881766381
4.6 -1.00863195297399 3.23150064581724
5.5 0 -5.3514242228123
spline3eq
4
0.00776934946045395 0.105197706160344
-55.14233165 -0.29745568008 0.00152870603877451
-44.7409899033333 -0.15449458722 0.00038933722543571
-34.3396481566667 0.05098657168 0.00038124926922248
-23.93830641 0.57342694704 0.0156639264890892
spline3eq
5
-0.00676745157022662 -0.0159520381982146
-23.9928 0.297607384684645 0
-15.9241175 0.216691597077105 -0.0024248755353942
-7.855435 0.0637598673719069 0.00306245895013358
0.213247499999998 -0.00183450621970427 -0.00177588407633909
8.28193 -0.111277018874367 0
spline3eq
10
2.77327511656661 0
2.055801767 -0.1485215264 72.2010867146919
2.31737934844444 1.6845304918 -47.2744689053404
2.57895692988889 2.0113365977 -15.1859578405326
2.84053451133333 1.1444092747 3.33978204841873
3.10211209277778 0.2861606803 2.587867603808
3.36368967422222 -0.3459281126 6.14070694084556
3.62526725566667 -0.6257480601 3.7397696717154
3.88684483711111 -0.6119510826 4.64749084871402
4.14842241855556 -0.3112059651 2.83275746415936
4.41 0 -15.0612086827734
spline3eq
5
12.3315547862781 0
1.9 2.62105440156724 0
2.8 10.2850803058354 -25.439802988016
3.7 3.23933763743897 -7.20203673434025
4.6 -5.79049355858613 39.5509978688682
5.5 0 -41.221771373642
spline3eq
8
8.33642274810572 -60.4024574736564
-1 0.07651409193 -110.652321293778
-0.724509054371429 0.14155824541 44.8853405500508
-0.449018108742857 0.75788697341 -25.3065115342002
-0.173527163114286 0.63011570378 -2.48510144915082
0.101963782514286 0.09049597305 2.68769386908235
0.377454728142857 -0.35741586657 -1.01558570129633
0.652945673771428 -0.65293217647 13.4224786001212
0.9284366194 -6.00912190653 -452.752542694929
spline3eq
5
0.137191606537625 -1.55094230968985
-1 0.0513843442016519 0
-0.5 0.0179024412245673 -2.44986494990154
0 -0.260650876879273 3.91774583656401
0.5 -0.190163791764901 -4.84414871911743
1 -0.763795416646599 0
spline3eq
8
0 0
-1 0 0
-0.724509054371429 0 0
-0.449018108742857 0 0
-0.173527163114286 0 0
0.101963782514286 0 0
0.377454728142857 0 0
0.652945673771428 0 0
0.9284366194 0 0

22
examples/USER/misc/meam_spline/in.meam-spline.Si Normal file
View File

@ -0,0 +1,22 @@
# Si fcc phase
units metal
boundary p p p
atom_style atomic
lattice fcc 3.98
region box block 0 5 0 5 0 5
create_box 1 box
create_atoms 1 box
pair_style meam/spline
pair_coeff * * Si_1.meam.spline Si
mass * 28.085
velocity all create 500.0 44226611
fix 1 all nvt temp 500.0 500.0 1.0
thermo 50
run 500

92
examples/USER/misc/meam_spline/in.meam-spline.TiO2 Normal file
View File

@ -0,0 +1,92 @@
#
variable T_depart equal 300
variable dt equal 0.0002
variable a equal 4.5937
variable c equal 2.9587
variable ca equal ${c}/${a}
variable nx equal 6
variable ny equal 6
variable nz equal 11
variable bx equal ${a}*${nx}
variable by equal ${a}*${ny}
variable bz equal ${c}*${nz}
# =======================================================================
units metal
atom_style atomic
dimension 3
boundary p p p
lattice sc 1.0
region box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
create_box 2 box_vide
#lattice sc 1.0
#region box_TiO2 block 0 ${bx} 0 ${by} 0 ${bz}
# titanium atoms
lattice custom ${a} origin 0.0 0.0 0.0 &
orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 ${ca} &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.5
create_atoms 2 region box_vide
# Oxygen atoms
lattice custom ${a} origin 0.0 0.0 0.0 &
orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 ${ca} &
basis 0.30478 0.30478 0.0 &
basis 0.69522 0.69522 0.0 &
basis 0.19522 0.80478 0.5 &
basis 0.80478 0.19522 0.5
create_atoms 1 region box_vide
mass 1 16.00
group Oxy type 1
mass 2 47.867
group Ti type 2
velocity all create ${T_depart} 277387
pair_style meam/spline
pair_coeff * * TiO.meam.spline O Ti
neighbor 0.5 bin
neigh_modify every 2 delay 0 check yes
timestep ${dt}
thermo_style custom step temp press pe ke etotal lx ly lz vol
thermo 10
#dump 5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z
fix 3 all nve
run 100
unfix 3
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-3 1.0e-5 1000 10000
unfix 1
reset_timestep 0
thermo 50
fix 3 all npt temp 300 300 0.1 aniso 1.0 1.0 1.0
run 500

88
examples/USER/misc/meam_spline/log.4May2017.meam-spline.Si.g++.1 Normal file
View File

@ -0,0 +1,88 @@
LAMMPS (13 Apr 2017)
using 1 OpenMP thread(s) per MPI task
# Si fcc phase
units metal
boundary p p p
atom_style atomic
lattice fcc 3.98
Lattice spacing in x,y,z = 3.98 3.98 3.98
region box block 0 5 0 5 0 5
create_box 1 box
Created orthogonal box = (0 0 0) to (19.9 19.9 19.9)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 500 atoms
pair_style meam/spline
pair_coeff * * Si_1.meam.spline Si
Reading potential file Si_1.meam.spline with DATE: 2012-02-01
mass * 28.085
velocity all create 500.0 44226611
fix 1 all nvt temp 500.0 500.0 1.0
thermo 50
run 500
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.5
ghost atom cutoff = 6.5
binsize = 3.25, bins = 7 7 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/spline, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/spline, 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) = 3.892 | 3.892 | 3.892 Mbytes
Step Temp E_pair E_mol TotEng Press
0 500 -1847.729 0 -1815.4786 1813162.7
50 1934.0932 -1940.8016 0 -1816.051 -48657.676
100 2570.1286 -1984.8725 0 -1819.0971 8002.4248
150 2566.7917 -1990.2724 0 -1824.7123 16819.447
200 2555.1319 -1995.2233 0 -1830.4152 5891.5313
250 2487.2881 -1995.8302 0 -1835.3981 -4339.7172
300 2381.4836 -1994.2492 0 -1840.6415 16508.04
350 2330.8663 -1996.6588 0 -1846.3161 24194.447
400 2212.6035 -1994.9278 0 -1852.2131 -9856.3709
450 2257.7531 -2003.8187 0 -1858.1918 -8029.6019
500 2211.4385 -2006.9846 0 -1864.345 4152.4867
Loop time of 5.13837 on 1 procs for 500 steps with 500 atoms
Performance: 8.407 ns/day, 2.855 hours/ns, 97.307 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.0952 | 5.0952 | 5.0952 | 0.0 | 99.16
Neigh | 0.026447 | 0.026447 | 0.026447 | 0.0 | 0.51
Comm | 0.0063307 | 0.0063307 | 0.0063307 | 0.0 | 0.12
Output | 0.0001905 | 0.0001905 | 0.0001905 | 0.0 | 0.00
Modify | 0.0082877 | 0.0082877 | 0.0082877 | 0.0 | 0.16
Other | | 0.00187 | | | 0.04
Nlocal: 500 ave 500 max 500 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1767 ave 1767 max 1767 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18059 ave 18059 max 18059 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 36118 ave 36118 max 36118 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 36118
Ave neighs/atom = 72.236
Neighbor list builds = 14
Dangerous builds = 0
Total wall time: 0:00:05

88
examples/USER/misc/meam_spline/log.4May2017.meam-spline.Si.g++.4 Normal file
View File

@ -0,0 +1,88 @@
LAMMPS (13 Apr 2017)
using 1 OpenMP thread(s) per MPI task
# Si fcc phase
units metal
boundary p p p
atom_style atomic
lattice fcc 3.98
Lattice spacing in x,y,z = 3.98 3.98 3.98
region box block 0 5 0 5 0 5
create_box 1 box
Created orthogonal box = (0 0 0) to (19.9 19.9 19.9)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 500 atoms
pair_style meam/spline
pair_coeff * * Si_1.meam.spline Si
Reading potential file Si_1.meam.spline with DATE: 2012-02-01
mass * 28.085
velocity all create 500.0 44226611
fix 1 all nvt temp 500.0 500.0 1.0
thermo 50
run 500
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.5
ghost atom cutoff = 6.5
binsize = 3.25, bins = 7 7 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/spline, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/spline, 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) = 3.861 | 3.861 | 3.861 Mbytes
Step Temp E_pair E_mol TotEng Press
0 500 -1847.729 0 -1815.4786 1813162.7
50 1923.4262 -1940.0936 0 -1816.0311 -38700.835
100 2535.2542 -1982.6249 0 -1819.0989 10216.821
150 2592.8247 -1992.1569 0 -1824.9176 4839.3385
200 2484.7391 -1990.8452 0 -1830.5775 14040.141
250 2597.4401 -2003.1619 0 -1835.625 1261.5199
300 2513.0793 -2002.942 0 -1840.8463 6690.9815
350 2390.933 -2001.0761 0 -1846.859 -4880.1146
400 2269.0782 -1999.3441 0 -1852.9867 -4921.4391
450 2287.5096 -2006.8236 0 -1859.2774 -7313.6151
500 2303.0918 -2014.0693 0 -1865.518 -9995.1789
Loop time of 1.46588 on 4 procs for 500 steps with 500 atoms
Performance: 29.470 ns/day, 0.814 hours/ns, 341.093 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 | 1.4273 | 1.4292 | 1.432 | 0.1 | 97.50
Neigh | 0.0068567 | 0.0070301 | 0.0073655 | 0.2 | 0.48
Comm | 0.019111 | 0.022127 | 0.024148 | 1.2 | 1.51
Output | 0.00023174 | 0.00024784 | 0.00029206 | 0.0 | 0.02
Modify | 0.005043 | 0.0052016 | 0.0054417 | 0.2 | 0.35
Other | | 0.002066 | | | 0.14
Nlocal: 125 ave 131 max 118 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Nghost: 979.25 ave 986 max 975 min
Histogram: 1 1 0 1 0 0 0 0 0 1
Neighs: 4541.75 ave 4712 max 4362 min
Histogram: 1 1 0 0 0 0 0 0 0 2
FullNghs: 9083.5 ave 9485 max 8601 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Total # of neighbors = 36334
Ave neighs/atom = 72.668
Neighbor list builds = 14
Dangerous builds = 0
Total wall time: 0:00:01

248
examples/USER/misc/meam_spline/log.4May2017.meam-spline.TiO2.g++.1 Normal file
View File

@ -0,0 +1,248 @@
LAMMPS (13 Apr 2017)
using 1 OpenMP thread(s) per MPI task
#
variable T_depart equal 300
variable dt equal 0.0002
variable a equal 4.5937
variable c equal 2.9587
variable ca equal ${c}/${a}
variable ca equal 2.9587/${a}
variable ca equal 2.9587/4.5937
variable nx equal 6
variable ny equal 6
variable nz equal 11
variable bx equal ${a}*${nx}
variable bx equal 4.5937*${nx}
variable bx equal 4.5937*6
variable by equal ${a}*${ny}
variable by equal 4.5937*${ny}
variable by equal 4.5937*6
variable bz equal ${c}*${nz}
variable bz equal 2.9587*${nz}
variable bz equal 2.9587*11
# =======================================================================
units metal
atom_style atomic
dimension 3
boundary p p p
lattice sc 1.0
Lattice spacing in x,y,z = 1 1 1
region box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 32.5457 0.0 0.0 0.0
create_box 2 box_vide
Created triclinic box = (0 0 0) to (27.5622 27.5622 32.5457) with tilt (0 0 0)
1 by 1 by 1 MPI processor grid
#lattice sc 1.0
#region box_TiO2 block 0 ${bx} 0 ${by} 0 ${bz}
# titanium atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 2 region box_vide
Created 792 atoms
# Oxygen atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 1 region box_vide
Created 1584 atoms
mass 1 16.00
group Oxy type 1
1584 atoms in group Oxy
mass 2 47.867
group Ti type 2
792 atoms in group Ti
velocity all create ${T_depart} 277387
velocity all create 300 277387
pair_style meam/spline
pair_coeff * * TiO.meam.spline O Ti
Reading potential file TiO.meam.spline with DATE: 2016-06-05
neighbor 0.5 bin
neigh_modify every 2 delay 0 check yes
timestep ${dt}
timestep 0.0002
thermo_style custom step temp press pe ke etotal lx ly lz vol
thermo 10
#dump 5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z
fix 3 all nve
run 100
Neighbor list info ...
update every 2 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6
ghost atom cutoff = 6
binsize = 3, bins = 10 10 11
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/spline, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/spline, 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) = 5.146 | 5.146 | 5.146 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
0 300 22403.656 -14374.073 92.097853 -14281.975 27.5622 27.5622 32.5457 24724.15
10 301.41345 23612.297 -14374.507 92.531772 -14281.975 27.5622 27.5622 32.5457 24724.15
20 305.11674 25127.832 -14375.643 93.668657 -14281.974 27.5622 27.5622 32.5457 24724.15
30 313.28903 26655.89 -14378.151 96.17749 -14281.974 27.5622 27.5622 32.5457 24724.15
40 328.94567 26999.049 -14382.957 100.98397 -14281.974 27.5622 27.5622 32.5457 24724.15
50 354.05827 23023.294 -14390.667 108.69336 -14281.974 27.5622 27.5622 32.5457 24724.15
60 390.48404 13594.655 -14401.849 119.87581 -14281.973 27.5622 27.5622 32.5457 24724.15
70 442.69928 151.15709 -14417.877 135.90551 -14281.972 27.5622 27.5622 32.5457 24724.15
80 516.89551 -14984.124 -14440.654 158.68322 -14281.971 27.5622 27.5622 32.5457 24724.15
90 618.22135 -29948.066 -14471.76 189.78953 -14281.971 27.5622 27.5622 32.5457 24724.15
100 747.6193 -41964.291 -14511.487 229.51378 -14281.973 27.5622 27.5622 32.5457 24724.15
Loop time of 38.7948 on 1 procs for 100 steps with 2376 atoms
Performance: 0.045 ns/day, 538.817 hours/ns, 2.578 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 | 38.774 | 38.774 | 38.774 | 0.0 | 99.95
Neigh | 0.010751 | 0.010751 | 0.010751 | 0.0 | 0.03
Comm | 0.0039313 | 0.0039313 | 0.0039313 | 0.0 | 0.01
Output | 0.00048804 | 0.00048804 | 0.00048804 | 0.0 | 0.00
Modify | 0.0039241 | 0.0039241 | 0.0039241 | 0.0 | 0.01
Other | | 0.001809 | | | 0.00
Nlocal: 2376 ave 2376 max 2376 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4479 ave 4479 max 4479 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 106396 ave 106396 max 106396 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 212792 ave 212792 max 212792 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 212792
Ave neighs/atom = 89.5589
Neighbor list builds = 1
Dangerous builds = 0
unfix 3
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-3 1.0e-5 1000 10000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
Per MPI rank memory allocation (min/avg/max) = 6.271 | 6.271 | 6.271 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
100 747.6193 -41964.291 -14511.487 229.51378 -14281.973 27.5622 27.5622 32.5457 24724.15
101 747.6193 -39284.65 -14517.424 229.51378 -14287.91 27.569615 27.569695 32.513154 24712.789
Loop time of 0.814693 on 1 procs for 1 steps with 2376 atoms
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-14511.4866189 -14511.4866189 -14517.4235162
Force two-norm initial, final = 5602.25 5486.97
Force max component initial, final = 5232.05 5109.43
Final line search alpha, max atom move = 1.9113e-07 0.000976563
Iterations, force evaluations = 1 1
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.81429 | 0.81429 | 0.81429 | 0.0 | 99.95
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 6.485e-05 | 6.485e-05 | 6.485e-05 | 0.0 | 0.01
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.0003347 | | | 0.04
Nlocal: 2376 ave 2376 max 2376 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4449 ave 4449 max 4449 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 105639 ave 105639 max 105639 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 211278 ave 211278 max 211278 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 211278
Ave neighs/atom = 88.9217
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
reset_timestep 0
thermo 50
fix 3 all npt temp 300 300 0.1 aniso 1.0 1.0 1.0
run 500
Per MPI rank memory allocation (min/avg/max) = 5.162 | 5.162 | 5.162 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
0 747.6193 -39284.65 -14517.424 229.51378 -14287.91 27.569615 27.569695 32.513154 24712.789
50 1155.2849 30650.319 -14678.807 354.6642 -14324.143 27.608688 27.60914 32.375311 24678.15
100 790.03926 99869.991 -14678.858 242.5364 -14436.322 27.777994 27.77799 32.017001 24704.857
150 938.86463 -21488.442 -14803.782 288.22472 -14515.557 27.996584 27.995139 31.67008 24822.003
200 420.11331 -790.80799 -14671.687 128.97178 -14542.715 28.126911 28.125909 31.431033 24864.93
250 352.18149 -3244.2491 -14665.007 108.1172 -14556.889 28.222686 28.223673 31.238649 24883.078
300 622.91245 3657.7097 -14758.201 191.22967 -14566.972 28.301771 28.30503 31.07216 24891.363
350 888.25374 26274.358 -14852.568 272.68754 -14579.881 28.370312 28.375107 30.937051 24904.656
400 735.44163 63109.066 -14823.872 225.77532 -14598.097 28.446905 28.45227 30.838015 24959.642
450 804.81905 6221.0364 -14861.113 247.07369 -14614.039 28.543942 28.548719 30.775793 25078.977
500 628.19106 -33912.026 -14814.726 192.85016 -14621.876 28.611997 28.615169 30.74081 25168.642
Loop time of 176.167 on 1 procs for 500 steps with 2376 atoms
Performance: 0.049 ns/day, 489.353 hours/ns, 2.838 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 175.9 | 175.9 | 175.9 | 0.0 | 99.85
Neigh | 0.17043 | 0.17043 | 0.17043 | 0.0 | 0.10
Comm | 0.018243 | 0.018243 | 0.018243 | 0.0 | 0.01
Output | 0.00040984 | 0.00040984 | 0.00040984 | 0.0 | 0.00
Modify | 0.067142 | 0.067142 | 0.067142 | 0.0 | 0.04
Other | | 0.00828 | | | 0.00
Nlocal: 2376 ave 2376 max 2376 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4358 ave 4358 max 4358 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 102634 ave 102634 max 102634 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 205268 ave 205268 max 205268 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 205268
Ave neighs/atom = 86.3923
Neighbor list builds = 16
Dangerous builds = 0
Total wall time: 0:03:37

248
examples/USER/misc/meam_spline/log.4May2017.meam-spline.TiO2.g++.4 Normal file
View File

@ -0,0 +1,248 @@
LAMMPS (13 Apr 2017)
using 1 OpenMP thread(s) per MPI task
#
variable T_depart equal 300
variable dt equal 0.0002
variable a equal 4.5937
variable c equal 2.9587
variable ca equal ${c}/${a}
variable ca equal 2.9587/${a}
variable ca equal 2.9587/4.5937
variable nx equal 6
variable ny equal 6
variable nz equal 11
variable bx equal ${a}*${nx}
variable bx equal 4.5937*${nx}
variable bx equal 4.5937*6
variable by equal ${a}*${ny}
variable by equal 4.5937*${ny}
variable by equal 4.5937*6
variable bz equal ${c}*${nz}
variable bz equal 2.9587*${nz}
variable bz equal 2.9587*11
# =======================================================================
units metal
atom_style atomic
dimension 3
boundary p p p
lattice sc 1.0
Lattice spacing in x,y,z = 1 1 1
region box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 32.5457 0.0 0.0 0.0
create_box 2 box_vide
Created triclinic box = (0 0 0) to (27.5622 27.5622 32.5457) with tilt (0 0 0)
1 by 2 by 2 MPI processor grid
#lattice sc 1.0
#region box_TiO2 block 0 ${bx} 0 ${by} 0 ${bz}
# titanium atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 2 region box_vide
Created 792 atoms
# Oxygen atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 1 region box_vide
Created 1584 atoms
mass 1 16.00
group Oxy type 1
1584 atoms in group Oxy
mass 2 47.867
group Ti type 2
792 atoms in group Ti
velocity all create ${T_depart} 277387
velocity all create 300 277387
pair_style meam/spline
pair_coeff * * TiO.meam.spline O Ti
Reading potential file TiO.meam.spline with DATE: 2016-06-05
neighbor 0.5 bin
neigh_modify every 2 delay 0 check yes
timestep ${dt}
timestep 0.0002
thermo_style custom step temp press pe ke etotal lx ly lz vol
thermo 10
#dump 5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z
fix 3 all nve
run 100
Neighbor list info ...
update every 2 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6
ghost atom cutoff = 6
binsize = 3, bins = 10 10 11
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam/spline, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam/spline, 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) = 3.922 | 3.922 | 3.922 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
0 300 22403.656 -14374.073 92.097853 -14281.975 27.5622 27.5622 32.5457 24724.15
10 301.16725 23582.084 -14374.431 92.456192 -14281.975 27.5622 27.5622 32.5457 24724.15
20 304.58237 25059.749 -14375.479 93.504609 -14281.974 27.5622 27.5622 32.5457 24724.15
30 312.41477 26504.358 -14377.883 95.9091 -14281.974 27.5622 27.5622 32.5457 24724.15
40 327.67099 26687.057 -14382.566 100.59265 -14281.974 27.5622 27.5622 32.5457 24724.15
50 352.32125 22677.292 -14390.134 108.1601 -14281.974 27.5622 27.5622 32.5457 24724.15
60 388.40592 12472.705 -14401.211 119.23784 -14281.973 27.5622 27.5622 32.5457 24724.15
70 439.97199 -1520.4694 -14417.04 135.06825 -14281.972 27.5622 27.5622 32.5457 24724.15
80 513.34361 -16733.316 -14439.564 157.59282 -14281.971 27.5622 27.5622 32.5457 24724.15
90 613.3542 -31099.591 -14470.267 188.29535 -14281.971 27.5622 27.5622 32.5457 24724.15
100 741.02836 -42358.226 -14509.464 227.4904 -14281.973 27.5622 27.5622 32.5457 24724.15
Loop time of 8.92317 on 4 procs for 100 steps with 2376 atoms
Performance: 0.194 ns/day, 123.933 hours/ns, 11.207 timesteps/s
99.5% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 8.8912 | 8.9 | 8.9064 | 0.2 | 99.74
Neigh | 0.0027034 | 0.0028808 | 0.0032032 | 0.4 | 0.03
Comm | 0.010964 | 0.017648 | 0.026568 | 5.0 | 0.20
Output | 0.00037575 | 0.00047809 | 0.00053835 | 0.0 | 0.01
Modify | 0.00099134 | 0.001001 | 0.0010085 | 0.0 | 0.01
Other | | 0.001162 | | | 0.01
Nlocal: 594 ave 599 max 589 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 2290.25 ave 2296 max 2282 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Neighs: 26671.5 ave 26934 max 26495 min
Histogram: 1 0 0 2 0 0 0 0 0 1
FullNghs: 53343 ave 53828 max 52922 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Total # of neighbors = 213372
Ave neighs/atom = 89.803
Neighbor list builds = 1
Dangerous builds = 0
unfix 3
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-3 1.0e-5 1000 10000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
Per MPI rank memory allocation (min/avg/max) = 5.047 | 5.047 | 5.047 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
100 741.02836 -42358.226 -14509.464 227.4904 -14281.973 27.5622 27.5622 32.5457 24724.15
101 741.02836 -39686.588 -14515.398 227.4904 -14287.907 27.569587 27.569656 32.513154 24712.729
Loop time of 0.193516 on 4 procs for 1 steps with 2376 atoms
99.5% CPU use with 4 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-14509.46351 -14509.46351 -14515.3978891
Force two-norm initial, final = 5602.69 5487.77
Force max component initial, final = 5235.27 5113.06
Final line search alpha, max atom move = 1.91012e-07 0.000976657
Iterations, force evaluations = 1 1
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.19287 | 0.19299 | 0.19318 | 0.0 | 99.73
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00014043 | 0.00033247 | 0.00045896 | 0.0 | 0.17
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.0001886 | | | 0.10
Nlocal: 594 ave 601 max 586 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Nghost: 2263.25 ave 2271 max 2251 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Neighs: 26425.8 ave 26807 max 26121 min
Histogram: 1 0 0 1 1 0 0 0 0 1
FullNghs: 52851.5 ave 53580 max 52175 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Total # of neighbors = 211406
Ave neighs/atom = 88.9756
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
reset_timestep 0
thermo 50
fix 3 all npt temp 300 300 0.1 aniso 1.0 1.0 1.0
run 500
Per MPI rank memory allocation (min/avg/max) = 3.937 | 3.937 | 3.937 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
0 741.02836 -39686.588 -14515.398 227.4904 -14287.907 27.569587 27.569656 32.513154 24712.729
50 1157.347 29332.549 -14679.321 355.29725 -14324.024 27.60903 27.609325 32.375509 24678.772
100 777.55858 101883.12 -14674.854 238.70492 -14436.149 27.778518 27.777373 32.017262 24704.976
150 945.49014 -18305.383 -14806.687 290.25871 -14516.428 27.998313 27.99535 31.670225 24823.838
200 427.46608 -4045.0095 -14674.887 131.22903 -14543.658 28.130283 28.127147 31.431578 24869.438
250 362.82166 -7283.1332 -14669.07 111.38365 -14557.687 28.225232 28.222707 31.238451 24884.314
300 626.2858 7228.0309 -14760.128 192.26526 -14567.862 28.302384 28.299949 31.070038 24885.734
350 859.84293 30084.735 -14845.064 263.96563 -14581.099 28.372349 28.369334 30.934424 24899.261
400 755.26136 54745.408 -14830.701 231.85983 -14598.842 28.450301 28.448361 30.836159 24957.691
450 802.52344 5690.2863 -14860.193 246.36895 -14613.824 28.542311 28.541672 30.773339 25069.354
500 631.84734 -31473.795 -14816.101 193.97261 -14622.128 28.605857 28.605891 30.737955 25152.746
Loop time of 39.7881 on 4 procs for 500 steps with 2376 atoms
Performance: 0.217 ns/day, 110.522 hours/ns, 12.567 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 | 39.617 | 39.633 | 39.653 | 0.2 | 99.61
Neigh | 0.043624 | 0.046792 | 0.051708 | 1.4 | 0.12
Comm | 0.05215 | 0.072616 | 0.092142 | 5.6 | 0.18
Output | 0.00042915 | 0.00045079 | 0.00051546 | 0.0 | 0.00
Modify | 0.029836 | 0.030341 | 0.03094 | 0.2 | 0.08
Other | | 0.004489 | | | 0.01
Nlocal: 594 ave 606 max 582 min
Histogram: 1 0 0 0 1 1 0 0 0 1
Nghost: 2226 ave 2238 max 2214 min
Histogram: 1 0 0 0 1 1 0 0 0 1
Neighs: 25652.8 ave 26129 max 25153 min
Histogram: 1 0 0 0 1 1 0 0 0 1
FullNghs: 51305.5 ave 52398 max 50251 min
Histogram: 1 0 0 0 1 1 0 0 0 1
Total # of neighbors = 205222
Ave neighs/atom = 86.3729
Neighbor list builds = 16
Dangerous builds = 0
Total wall time: 0:00:49

130
potentials/TiO.meam.spline Normal file
View File

@ -0,0 +1,130 @@
# Ti-O cubic spline potential where O is in the dilute limit. DATE: 2016-06-05 CONTRIBUTOR: Pinchao Zhang, Dallas R. Trinkle
meam/spline 2 Ti O
spline3eq
13
-20 0
1.742692837 3.744277175966 99.4865081627958
2.05580176725 0.910839730906 10.8702523265355
2.3689106975 0.388045896634 -1.55322418749562
2.68201962775 -0.018840906533 2.43630041329215
2.995128558 -0.248098929639 2.67912713976835
3.30823748825 -0.264489550297 -0.125056384603077
3.6213464185 -0.227196189283 1.10662555360438
3.93445534875 -0.129293090176 -0.592053676745914
4.247564279 -0.059685366933 -0.470123414607672
4.56067320925 -0.031100025561 -0.0380739973059663
4.8737821395 -0.013847363202 -0.0711547960695406
5.18689106975 -0.003203412728 -0.081768292420175
5.5 0 -0.0571422964883619
spline3eq
5
0.155001355787331 0
1.9 0.533321679606674 0
2.8 0.456402081843862 -1.60311717015859
3.7 -0.324281383502201 1.19940299483249
4.6 -0.474029826906675 1.47909794595154
5.5 0 -2.49521499855605
spline3eq
13
0 0
1.742692837 0 0
2.05580176725 0 0
2.3689106975 0 0
2.68201962775 0 0
2.995128558 0 0
3.30823748825 0 0
3.6213464185 0 0
3.93445534875 0 0
4.247564279 0 0
4.56067320925 0 0
4.8737821395 0 0
5.18689106975 0 0
5.5 0 0
spline3eq
11
-1 0
2.055801767 1.7475279661 -525.869786904802
2.2912215903 -5.8677963945 252.796316927755
2.5266414136 -8.3376288737 71.7318388721015
2.7620612369 -5.8398712842 -1.93587742753693
2.9974810602 -3.1140648231 -39.2999192667503
3.2329008835 -1.7257245065 14.3424136002004
3.4683207068 -0.4428977017 -29.4925534559498
3.7037405301 -0.1466643003 -3.18010534572236
3.9391603534 -0.2095507945 3.33490838803603
4.1745801767 -0.1442384563 3.71918691359508
4.41 0 -9.66717019857564
spline3eq
5
-61.9827585211652 0
1.9 11.2293641315584 0
2.8 -27.9976343076148 122.648031332411
3.7 -8.32979773113248 -54.3340881766381
4.6 -1.00863195297399 3.23150064581724
5.5 0 -5.3514242228123
spline3eq
4
0.00776934946045395 0.105197706160344
-55.14233165 -0.29745568008 0.00152870603877451
-44.7409899033333 -0.15449458722 0.00038933722543571
-34.3396481566667 0.05098657168 0.00038124926922248
-23.93830641 0.57342694704 0.0156639264890892
spline3eq
5
-0.00676745157022662 -0.0159520381982146
-23.9928 0.297607384684645 0
-15.9241175 0.216691597077105 -0.0024248755353942
-7.855435 0.0637598673719069 0.00306245895013358
0.213247499999998 -0.00183450621970427 -0.00177588407633909
8.28193 -0.111277018874367 0
spline3eq
10
2.77327511656661 0
2.055801767 -0.1485215264 72.2010867146919
2.31737934844444 1.6845304918 -47.2744689053404
2.57895692988889 2.0113365977 -15.1859578405326
2.84053451133333 1.1444092747 3.33978204841873
3.10211209277778 0.2861606803 2.587867603808
3.36368967422222 -0.3459281126 6.14070694084556
3.62526725566667 -0.6257480601 3.7397696717154
3.88684483711111 -0.6119510826 4.64749084871402
4.14842241855556 -0.3112059651 2.83275746415936
4.41 0 -15.0612086827734
spline3eq
5
12.3315547862781 0
1.9 2.62105440156724 0
2.8 10.2850803058354 -25.439802988016
3.7 3.23933763743897 -7.20203673434025
4.6 -5.79049355858613 39.5509978688682
5.5 0 -41.221771373642
spline3eq
8
8.33642274810572 -60.4024574736564
-1 0.07651409193 -110.652321293778
-0.724509054371429 0.14155824541 44.8853405500508
-0.449018108742857 0.75788697341 -25.3065115342002
-0.173527163114286 0.63011570378 -2.48510144915082
0.101963782514286 0.09049597305 2.68769386908235
0.377454728142857 -0.35741586657 -1.01558570129633
0.652945673771428 -0.65293217647 13.4224786001212
0.9284366194 -6.00912190653 -452.752542694929
spline3eq
5
0.137191606537625 -1.55094230968985
-1 0.0513843442016519 0
-0.5 0.0179024412245673 -2.44986494990154
0 -0.260650876879273 3.91774583656401
0.5 -0.190163791764901 -4.84414871911743
1 -0.763795416646599 0
spline3eq
8
0 0
-1 0 0
-0.724509054371429 0 0
-0.449018108742857 0 0
-0.173527163114286 0 0
0.101963782514286 0 0
0.377454728142857 0 0
0.652945673771428 0 0
0.9284366194 0 0

627
src/USER-MISC/pair_meam_spline.cpp
View File

@ -13,6 +13,9 @@
/* ----------------------------------------------------------------------
Contributing author: Alexander Stukowski (LLNL), alex@stukowski.com
Will Tipton (Cornell), wwt26@cornell.edu
Dallas R. Trinkle (UIUC), dtrinkle@illinois.edu
Pinchao Zhang (UIUC)
see LLNL copyright notice at bottom of file
------------------------------------------------------------------------- */
@ -23,6 +26,9 @@
* 25-Mar-11 - AS: Fixed calculation of per-atom virial stress.
* 11-Apr-11 - AS: Adapted code to new memory management of LAMMPS.
* 24-Sep-11 - AS: Adapted code to new interface of Error::one() function.
* 20-Jun-13 - WT: Added support for multiple species types
* 25-Apr-17 - DRT/PZ: Modified format of multiple species type to
conform with pairing, updated to LAMMPS style
------------------------------------------------------------------------- */
#include <math.h>
@ -49,7 +55,6 @@ PairMEAMSpline::PairMEAMSpline(LAMMPS *lmp) : Pair(lmp)
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nelements = 0;
elements = NULL;
@ -77,6 +82,15 @@ PairMEAMSpline::~PairMEAMSpline()
if(allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete[] phis;
delete[] Us;
delete[] rhos;
delete[] fs;
delete[] gs;
delete[] zero_atom_energies;
delete [] map;
}
}
@ -85,11 +99,16 @@ PairMEAMSpline::~PairMEAMSpline()
void PairMEAMSpline::compute(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag, vflag);
else evflag = vflag_fdotr =
eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
const double* const * const x = atom->x;
double* const * const forces = atom->f;
const int ntypes = atom->ntypes;
double cutforcesq = cutoff*cutoff;
if (eflag || vflag) {
ev_setup(eflag, vflag);
} else {
evflag = vflag_fdotr = eflag_global = 0;
vflag_global = eflag_atom = vflag_atom = 0;
}
// Grow per-atom array if necessary
@ -99,22 +118,13 @@ void PairMEAMSpline::compute(int eflag, int vflag)
memory->create(Uprime_values,nmax,"pair:Uprime");
}
double** const x = atom->x;
double** forces = atom->f;
int nlocal = atom->nlocal;
bool newton_pair = force->newton_pair;
int inum_full = listfull->inum;
int* ilist_full = listfull->ilist;
int* numneigh_full = listfull->numneigh;
int** firstneigh_full = listfull->firstneigh;
// Determine the maximum number of neighbors a single atom has
int newMaxNeighbors = 0;
for(int ii = 0; ii < inum_full; ii++) {
int jnum = numneigh_full[ilist_full[ii]];
if(jnum > newMaxNeighbors) newMaxNeighbors = jnum;
for(int ii = 0; ii < listfull->inum; ii++) {
int jnum = listfull->numneigh[listfull->ilist[ii]];
if(jnum > newMaxNeighbors)
newMaxNeighbors = jnum;
}
// Allocate array for temporary bond info
@ -126,35 +136,35 @@ void PairMEAMSpline::compute(int eflag, int vflag)
}
// Sum three-body contributions to charge density and
// compute embedding energies
// the embedding energy
for(int ii = 0; ii < inum_full; ii++) {
int i = ilist_full[ii];
double xtmp = x[i][0];
double ytmp = x[i][1];
double ztmp = x[i][2];
int* jlist = firstneigh_full[i];
int jnum = numneigh_full[i];
double rho_value = 0;
for(int ii = 0; ii < listfull->inum; ii++) {
int i = listfull->ilist[ii];
int numBonds = 0;
MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
// compute charge density and numBonds
MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
double rho_value = 0;
const int ntypes = atom->ntypes;
const int itype = atom->type[i];
for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
int j = listfull->firstneigh[i][jj];
j &= NEIGHMASK;
double jdelx = x[j][0] - xtmp;
double jdely = x[j][1] - ytmp;
double jdelz = x[j][2] - ztmp;
double jdelx = x[j][0] - x[i][0];
double jdely = x[j][1] - x[i][1];
double jdelz = x[j][2] - x[i][2];
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
if(rij_sq < cutforcesq) {
if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq);
double partial_sum = 0;
const int jtype = atom->type[j];
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->f = fs[i_to_potl(jtype)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
@ -164,11 +174,11 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * g.eval(cos_theta);
partial_sum += bondk.f * gs[ij_to_potl(jtype,atom->type[bondk.tag],ntypes)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rho.eval(rij);
rho_value += rhos[i_to_potl(jtype)].eval(rij);
numBonds++;
nextTwoBodyInfo++;
@ -176,19 +186,20 @@ void PairMEAMSpline::compute(int eflag, int vflag)
}
// Compute embedding energy and its derivative
double Uprime_i;
double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy;
double embeddingEnergy = Us[i_to_potl(itype)].eval(rho_value, Uprime_i)
- zero_atom_energies[i_to_potl(itype)];
Uprime_values[i] = Uprime_i;
if(eflag) {
if(eflag_global) eng_vdwl += embeddingEnergy;
if(eflag_atom) eatom[i] += embeddingEnergy;
if(eflag_global)
eng_vdwl += embeddingEnergy;
if(eflag_atom)
eatom[i] += embeddingEnergy;
}
double forces_i[3] = {0, 0, 0};
// Compute three-body contributions to force
double forces_i[3] = {0, 0, 0};
for(int jj = 0; jj < numBonds; jj++) {
const MEAM2Body bondj = twoBodyInfo[jj];
double rij = bondj.r;
@ -198,6 +209,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double f_rij = bondj.f;
double forces_j[3] = {0, 0, 0};
const int jtype = atom->type[j];
MEAM2Body const* bondk = twoBodyInfo;
for(int kk = 0; kk < jj; kk++, ++bondk) {
@ -207,7 +219,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
bondj.del[1]*bondk->del[1] +
bondj.del[2]*bondk->del[2]);
double g_prime;
double g_value = g.eval(cos_theta, g_prime);
double g_value = gs[ij_to_potl(jtype,atom->type[bondk->tag],ntypes)].eval(cos_theta, g_prime);
double f_rik_prime = bondk->fprime;
double f_rik = bondk->f;
@ -271,40 +283,32 @@ void PairMEAMSpline::compute(int eflag, int vflag)
comm->forward_comm_pair(this);
int inum_half = listhalf->inum;
int* ilist_half = listhalf->ilist;
int* numneigh_half = listhalf->numneigh;
int** firstneigh_half = listhalf->firstneigh;
// Compute two-body pair interactions
for(int ii = 0; ii < listhalf->inum; ii++) {
int i = listhalf->ilist[ii];
const int itype = atom->type[i];
for(int ii = 0; ii < inum_half; ii++) {
int i = ilist_half[ii];
double xtmp = x[i][0];
double ytmp = x[i][1];
double ztmp = x[i][2];
int* jlist = firstneigh_half[i];
int jnum = numneigh_half[i];
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
for(int jj = 0; jj < listhalf->numneigh[i]; jj++) {
int j = listhalf->firstneigh[i][jj];
j &= NEIGHMASK;
double jdel[3];
jdel[0] = x[j][0] - xtmp;
jdel[1] = x[j][1] - ytmp;
jdel[2] = x[j][2] - ztmp;
jdel[0] = x[j][0] - x[i][0];
jdel[1] = x[j][1] - x[i][1];
jdel[2] = x[j][2] - x[i][2];
double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2];
if(rij_sq < cutforcesq) {
if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq);
const int jtype = atom->type[j];
double rho_prime;
rho.eval(rij, rho_prime);
double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]);
double rho_prime_i,rho_prime_j;
rhos[i_to_potl(itype)].eval(rij,rho_prime_i);
rhos[i_to_potl(jtype)].eval(rij,rho_prime_j);
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
double pair_pot_deriv;
double pair_pot = phi.eval(rij, pair_pot_deriv);
double pair_pot = phis[ij_to_potl(itype,jtype,ntypes)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient
@ -317,13 +321,14 @@ void PairMEAMSpline::compute(int eflag, int vflag)
forces[j][0] -= jdel[0]*fpair;
forces[j][1] -= jdel[1]*fpair;
forces[j][2] -= jdel[2]*fpair;
if (evflag) ev_tally(i, j, nlocal, newton_pair,
if (evflag) ev_tally(i, j, atom->nlocal, force->newton_pair,
pair_pot, 0.0, -fpair, jdel[0], jdel[1], jdel[2]);
}
}
}
if(vflag_fdotr) virial_fdotr_compute();
if(vflag_fdotr)
virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
@ -331,11 +336,23 @@ void PairMEAMSpline::compute(int eflag, int vflag)
void PairMEAMSpline::allocate()
{
allocated = 1;
int n = atom->ntypes;
int n = nelements;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
int nmultichoose2 = n*(n+1)/2;
//Change the functional form
//f_ij->f_i
//g_i(cos\theta_ijk)->g_jk(cos\theta_ijk)
phis = new SplineFunction[nmultichoose2];
Us = new SplineFunction[n];
rhos = new SplineFunction[n];
fs = new SplineFunction[n];
gs = new SplineFunction[nmultichoose2];
zero_atom_energies = new double[n];
map = new int[n+1];
}
@ -356,8 +373,6 @@ void PairMEAMSpline::coeff(int narg, char **arg)
{
int i,j,n;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
@ -366,45 +381,34 @@ void PairMEAMSpline::coeff(int narg, char **arg)
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read potential file: also sets the number of elements.
read_file(arg[2]);
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements
// elements = list of element names
if (elements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
elements = new char*[atom->ntypes];
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
nelements = 0;
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
map[i-2] = j;
if (j == nelements) {
n = strlen(arg[i]) + 1;
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
if ((nelements == 1) && (strlen(elements[0]) == 0)) {
// old style: we only have one species, so we're either "NULL" or we match.
for (i = 3; i < narg; i++)
if (strcmp(arg[i],"NULL") == 0)
map[i-2] = -1;
else
map[i-2] = 0;
} else {
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0)
break;
if (j < nelements) map[i-2] = j;
else error->all(FLERR,"No matching element in EAM potential file");
}
}
// for now, only allow single element
if (nelements > 1)
error->all(FLERR,
"Pair meam/spline only supports single element potentials");
// read potential file
read_file(arg[2]);
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
@ -425,65 +429,134 @@ void PairMEAMSpline::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
#define MAXLINE 1024
void PairMEAMSpline::read_file(const char* filename)
{
if(comm->me == 0) {
FILE *fp = force->open_potential(filename);
if(fp == NULL) {
char str[1024];
sprintf(str,"Cannot open spline MEAM potential file %s", filename);
error->one(FLERR,str);
}
int nmultichoose2; // = (n+1)*n/2;
// Skip first line of file.
char line[MAXLINE];
fgets(line, MAXLINE, fp);
if(comm->me == 0) {
FILE *fp = force->open_potential(filename);
if(fp == NULL) {
char str[1024];
sprintf(str,"Cannot open spline MEAM potential file %s", filename);
error->one(FLERR,str);
}
// Parse spline functions.
phi.parse(fp, error);
rho.parse(fp, error);
U.parse(fp, error);
f.parse(fp, error);
g.parse(fp, error);
// Skip first line of file. It's a comment.
char line[MAXLINE];
char *ptr;
fgets(line, MAXLINE, fp);
fclose(fp);
}
// Second line holds potential type ("meam/spline")
// in new potential format.
// Transfer spline functions from master processor to all other processors.
phi.communicate(world, comm->me);
rho.communicate(world, comm->me);
f.communicate(world, comm->me);
U.communicate(world, comm->me);
g.communicate(world, comm->me);
bool isNewFormat = false;
fgets(line, MAXLINE, fp);
ptr = strtok(line, " \t\n\r\f");
// Calculate 'zero-point energy' of single atom in vacuum.
zero_atom_energy = U.eval(0.0);
if (strcmp(ptr, "meam/spline") == 0) {
isNewFormat = true;
// parse the rest of the line!
ptr = strtok(NULL," \t\n\r\f");
if (ptr == NULL)
error->one(FLERR,"Need to include number of atomic species on"
" meam/spline line in multi-element potential file");
nelements = atoi(ptr);
if (nelements < 1)
error->one(FLERR, "Invalid number of atomic species on"
" meam/spline line in potential file");
elements = new char*[nelements];
for (int i=0; i<nelements; ++i) {
ptr = strtok(NULL," \t\n\r\f");
if (ptr == NULL)
error->one(FLERR, "Not enough atomic species in meam/spline"
" line of multi-element potential file");
elements[i] = new char[strlen(ptr)+1];
strcpy(elements[i], ptr);
}
} else {
isNewFormat = false;
nelements = 1; // old format only handles one species; (backwards compatibility)
elements = new char*[1];
elements[0] = new char[1];
strcpy(elements[0], "");
rewind(fp);
fgets(line, MAXLINE, fp);
}
// Determine maximum cutoff radius of all relevant spline functions.
cutoff = 0.0;
if(phi.cutoff() > cutoff) cutoff = phi.cutoff();
if(rho.cutoff() > cutoff) cutoff = rho.cutoff();
if(f.cutoff() > cutoff) cutoff = f.cutoff();
nmultichoose2 = ((nelements+1)*nelements)/2;
// allocate!!
allocate();
// Set LAMMPS pair interaction flags.
for(int i = 1; i <= atom->ntypes; i++) {
for(int j = 1; j <= atom->ntypes; j++) {
setflag[i][j] = 1;
cutsq[i][j] = cutoff;
}
}
// Parse spline functions.
for (int i = 0; i < nmultichoose2; i++)
phis[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nelements; i++)
rhos[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nelements; i++)
Us[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nelements; i++)
fs[i].parse(fp, error, isNewFormat);
for (int i = 0; i < nmultichoose2; i++)
gs[i].parse(fp, error, isNewFormat);
fclose(fp);
}
// Transfer spline functions from master processor to all other processors.
MPI_Bcast(&nelements, 1, MPI_INT, 0, world);
MPI_Bcast(&nmultichoose2, 1, MPI_INT, 0, world);
// allocate!!
if (comm->me != 0) {
allocate();
elements = new char*[nelements];
}
for (int i = 0; i < nelements; ++i) {
int n;
if (comm->me == 0)
n = strlen(elements[i]);
MPI_Bcast(&n, 1, MPI_INT, 0, world);
if (comm->me != 0)
elements[i] = new char[n+1];
MPI_Bcast(elements[i], n+1, MPI_CHAR, 0, world);
}
for (int i = 0; i < nmultichoose2; i++)
phis[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
rhos[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
fs[i].communicate(world, comm->me);
for (int i = 0; i < nelements; i++)
Us[i].communicate(world, comm->me);
for (int i = 0; i < nmultichoose2; i++)
gs[i].communicate(world, comm->me);
// Calculate 'zero-point energy' of single atom in vacuum.
for (int i = 0; i < nelements; i++)
zero_atom_energies[i] = Us[i].eval(0.0);
// Determine maximum cutoff radius of all relevant spline functions.
cutoff = 0.0;
for (int i = 0; i < nmultichoose2; i++)
if(phis[i].cutoff() > cutoff)
cutoff = phis[i].cutoff();
for (int i = 0; i < nelements; i++)
if(rhos[i].cutoff() > cutoff)
cutoff = rhos[i].cutoff();
for (int i = 0; i < nelements; i++)
if(fs[i].cutoff() > cutoff)
cutoff = fs[i].cutoff();
// Set LAMMPS pair interaction flags.
for(int i = 1; i <= atom->ntypes; i++) {
for(int j = 1; j <= atom->ntypes; j++) {
// setflag[i][j] = 1;
cutsq[i][j] = cutoff;
}
}
//phi.writeGnuplot("phi.gp", "Phi(r)");
//rho.writeGnuplot("rho.gp", "Rho(r)");
//f.writeGnuplot("f.gp", "f(r)");
//U.writeGnuplot("U.gp", "U(rho)");
//g.writeGnuplot("g.gp", "g(x)");
}
/* ----------------------------------------------------------------------
@ -491,16 +564,19 @@ void PairMEAMSpline::read_file(const char* filename)
------------------------------------------------------------------------- */
void PairMEAMSpline::init_style()
{
if(force->newton_pair == 0)
error->all(FLERR,"Pair style meam/spline requires newton pair on");
if(force->newton_pair == 0)
error->all(FLERR,"Pair style meam/spline requires newton pair on");
// Need both 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;
// Need both 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->requests[irequest_half]->half = 1;
// neighbor->requests[irequest_half]->halffull = 1;
// neighbor->requests[irequest_half]->halffulllist = irequest_full;
}
/* ----------------------------------------------------------------------
@ -509,8 +585,8 @@ void PairMEAMSpline::init_style()
------------------------------------------------------------------------- */
void PairMEAMSpline::init_list(int id, NeighList *ptr)
{
if(id == 1) listfull = ptr;
else if(id == 2) listhalf = ptr;
if(id == 1) listfull = ptr;
else if(id == 2) listhalf = ptr;
}
/* ----------------------------------------------------------------------
@ -518,33 +594,33 @@ void PairMEAMSpline::init_list(int id, NeighList *ptr)
------------------------------------------------------------------------- */
double PairMEAMSpline::init_one(int i, int j)
{
return cutoff;
return cutoff;
}
/* ---------------------------------------------------------------------- */
int PairMEAMSpline::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int* list_iter = list;
int* list_iter_end = list + n;
while(list_iter != list_iter_end)
*buf++ = Uprime_values[*list_iter++];
return n;
int* list_iter = list;
int* list_iter_end = list + n;
while(list_iter != list_iter_end)
*buf++ = Uprime_values[*list_iter++];
return n;
}
/* ---------------------------------------------------------------------- */
void PairMEAMSpline::unpack_forward_comm(int n, int first, double *buf)
{
memcpy(&Uprime_values[first], buf, n * sizeof(buf[0]));
memcpy(&Uprime_values[first], buf, n * sizeof(buf[0]));
}
/* ---------------------------------------------------------------------- */
int PairMEAMSpline::pack_reverse_comm(int n, int first, double *buf)
{
return 0;
return 0;
}
/* ---------------------------------------------------------------------- */
@ -558,141 +634,148 @@ void PairMEAMSpline::unpack_reverse_comm(int n, int *list, double *buf)
------------------------------------------------------------------------- */
double PairMEAMSpline::memory_usage()
{
return nmax * sizeof(double); // The Uprime_values array.
return nmax * sizeof(double); // The Uprime_values array.
}
/// Parses the spline knots from a text file.
void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error)
void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error,
bool isNewFormat)
{
char line[MAXLINE];
char line[MAXLINE];
// Parse number of spline knots.
fgets(line, MAXLINE, fp);
int n = atoi(line);
if(n < 2)
error->one(FLERR,"Invalid number of spline knots in MEAM potential file");
// If new format, read the spline format. Should always be "spline3eq" for now.
if (isNewFormat)
fgets(line, MAXLINE, fp);
// Parse first derivatives at beginning and end of spline.
fgets(line, MAXLINE, fp);
double d0 = atof(strtok(line, " \t\n\r\f"));
double dN = atof(strtok(NULL, " \t\n\r\f"));
init(n, d0, dN);
// Parse number of spline knots.
fgets(line, MAXLINE, fp);
int n = atoi(line);
if(n < 2)
error->one(FLERR,"Invalid number of spline knots in MEAM potential file");
// Skip line.
fgets(line, MAXLINE, fp);
// Parse first derivatives at beginning and end of spline.
fgets(line, MAXLINE, fp);
double d0 = atof(strtok(line, " \t\n\r\f"));
double dN = atof(strtok(NULL, " \t\n\r\f"));
init(n, d0, dN);
// Parse knot coordinates.
for(int i=0; i<n; i++) {
fgets(line, MAXLINE, fp);
double x, y, y2;
if(sscanf(line, "%lg %lg %lg", &x, &y, &y2) != 3) {
error->one(FLERR,"Invalid knot line in MEAM potential file");
}
setKnot(i, x, y);
}
// Skip line in old format
if (!isNewFormat)
fgets(line, MAXLINE, fp);
prepareSpline(error);
// Parse knot coordinates.
for(int i=0; i<n; i++) {
fgets(line, MAXLINE, fp);
double x, y, y2;
if(sscanf(line, "%lg %lg %lg", &x, &y, &y2) != 3) {
error->one(FLERR,"Invalid knot line in MEAM potential file");
}
setKnot(i, x, y);
}
prepareSpline(error);
}
/// Calculates the second derivatives at the knots of the cubic spline.
void PairMEAMSpline::SplineFunction::prepareSpline(Error* error)
{
xmin = X[0];
xmax = X[N-1];
xmin = X[0];
xmax = X[N-1];
isGridSpline = true;
h = (xmax-xmin)/(N-1);
hsq = h*h;
isGridSpline = true;
h = (xmax-xmin)/(N-1);
hsq = h*h;
double* u = new double[N];
Y2[0] = -0.5;
u[0] = (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - deriv0);
for(int i = 1; i <= N-2; i++) {
double sig = (X[i]-X[i-1]) / (X[i+1]-X[i-1]);
double p = sig * Y2[i-1] + 2.0;
Y2[i] = (sig - 1.0) / p;
u[i] = (Y[i+1]-Y[i]) / (X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
u[i] = (6.0 * u[i]/(X[i+1]-X[i-1]) - sig*u[i-1])/p;
double* u = new double[N];
Y2[0] = -0.5;
u[0] = (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - deriv0);
for(int i = 1; i <= N-2; i++) {
double sig = (X[i]-X[i-1]) / (X[i+1]-X[i-1]);
double p = sig * Y2[i-1] + 2.0;
Y2[i] = (sig - 1.0) / p;
u[i] = (Y[i+1]-Y[i]) / (X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
u[i] = (6.0 * u[i]/(X[i+1]-X[i-1]) - sig*u[i-1])/p;
if(fabs(h*i+xmin - X[i]) > 1e-8)
isGridSpline = false;
}
if(fabs(h*i+xmin - X[i]) > 1e-8)
isGridSpline = false;
}
double qn = 0.5;
double un = (3.0/(X[N-1]-X[N-2])) * (derivN - (Y[N-1]-Y[N-2])/(X[N-1]-X[N-2]));
Y2[N-1] = (un - qn*u[N-2]) / (qn * Y2[N-2] + 1.0);
for(int k = N-2; k >= 0; k--) {
Y2[k] = Y2[k] * Y2[k+1] + u[k];
}
double qn = 0.5;
double un = (3.0/(X[N-1]-X[N-2])) * (derivN - (Y[N-1]-Y[N-2])/(X[N-1]-X[N-2]));
Y2[N-1] = (un - qn*u[N-2]) / (qn * Y2[N-2] + 1.0);
for(int k = N-2; k >= 0; k--) {
Y2[k] = Y2[k] * Y2[k+1] + u[k];
}
delete[] u;
delete[] u;
#if !SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
if(!isGridSpline)
error->one(FLERR,"Support for MEAM potentials with non-uniform cubic splines has not been enabled in the MEAM potential code. Set SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES in pair_spline_meam.h to 1 to enable it");
if(!isGridSpline)
error->one(FLERR,"Support for MEAM potentials with non-uniform cubic splines has not been enabled in the MEAM potential code. Set SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES in pair_spline_meam.h to 1 to enable it");
#endif
// Shift the spline to X=0 to speed up interpolation.
for(int i = 0; i < N; i++) {
Xs[i] = X[i] - xmin;
// Shift the spline to X=0 to speed up interpolation.
for(int i = 0; i < N; i++) {
Xs[i] = X[i] - xmin;
#if !SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
if(i < N-1) Ydelta[i] = (Y[i+1]-Y[i])/h;
Y2[i] /= h*6.0;
if(i < N-1) Ydelta[i] = (Y[i+1]-Y[i])/h;
Y2[i] /= h*6.0;
#endif
}
xmax_shifted = xmax - xmin;
}
xmax_shifted = xmax - xmin;
}
/// Broadcasts the spline function parameters to all processors.
void PairMEAMSpline::SplineFunction::communicate(MPI_Comm& world, int me)
{
MPI_Bcast(&N, 1, MPI_INT, 0, world);
MPI_Bcast(&deriv0, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&derivN, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmin, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmax, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmax_shifted, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
if(me != 0) {
X = new double[N];
Xs = new double[N];
Y = new double[N];
Y2 = new double[N];
Ydelta = new double[N];
}
MPI_Bcast(X, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Xs, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Y, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Y2, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Ydelta, N, MPI_DOUBLE, 0, world);
MPI_Bcast(&N, 1, MPI_INT, 0, world);
MPI_Bcast(&deriv0, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&derivN, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmin, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmax, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&xmax_shifted, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
if(me != 0) {
X = new double[N];
Xs = new double[N];
Y = new double[N];
Y2 = new double[N];
Ydelta = new double[N];
}
MPI_Bcast(X, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Xs, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Y, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Y2, N, MPI_DOUBLE, 0, world);
MPI_Bcast(Ydelta, N, MPI_DOUBLE, 0, world);
}
/// Writes a Gnuplot script that plots the spline function.
///
/// This function is for debugging only!
void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename, const char* title) const
void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename,
const char* title) const
{
FILE* fp = fopen(filename, "w");
fprintf(fp, "#!/usr/bin/env gnuplot\n");
if(title) fprintf(fp, "set title \"%s\"\n", title);
double tmin = X[0] - (X[N-1] - X[0]) * 0.05;
double tmax = X[N-1] + (X[N-1] - X[0]) * 0.05;
double delta = (tmax - tmin) / (N*200);
fprintf(fp, "set xrange [%f:%f]\n", tmin, tmax);
fprintf(fp, "plot '-' with lines notitle, '-' with points notitle pt 3 lc 3\n");
for(double x = tmin; x <= tmax+1e-8; x += delta) {
double y = eval(x);
fprintf(fp, "%f %f\n", x, y);
}
fprintf(fp, "e\n");
for(int i = 0; i < N; i++) {
fprintf(fp, "%f %f\n", X[i], Y[i]);
}
fprintf(fp, "e\n");
fclose(fp);
FILE* fp = fopen(filename, "w");
fprintf(fp, "#!/usr/bin/env gnuplot\n");
if(title) fprintf(fp, "set title \"%s\"\n", title);
double tmin = X[0] - (X[N-1] - X[0]) * 0.05;
double tmax = X[N-1] + (X[N-1] - X[0]) * 0.05;
double delta = (tmax - tmin) / (N*200);
fprintf(fp, "set xrange [%f:%f]\n", tmin, tmax);
fprintf(fp, "plot '-' with lines notitle, '-' with points notitle pt 3 lc 3\n");
for(double x = tmin; x <= tmax+1e-8; x += delta) {
double y = eval(x);
fprintf(fp, "%f %f\n", x, y);
}
fprintf(fp, "e\n");
for(int i = 0; i < N; i++) {
fprintf(fp, "%f %f\n", X[i], Y[i]);
}
fprintf(fp, "e\n");
fclose(fp);
}
/* ----------------------------------------------------------------------
@ -734,3 +817,5 @@ void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename, const ch
* Lawrence Livermore National Security, LLC, and shall not be used for
* advertising or product endorsement purposes.
------------------------------------------------------------------------- */

357
src/USER-MISC/pair_meam_spline.h
View File

@ -1,4 +1,4 @@
/* -*- c++ -*- ----------------------------------------------------------
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -28,209 +28,230 @@ PairStyle(meam/spline,PairMEAMSpline)
namespace LAMMPS_NS {
/// Set this to 1 if you intend to use MEAM potentials with non-uniform spline knots.
/// Set this to 0 if you intend to use only MEAM potentials with spline knots on a uniform grid.
///
/// With SUPPORT_NON_GRID_SPLINES == 0, the code runs about 50% faster.
// Set this to 1 if you intend to use MEAM potentials with
// non-uniform spline knots.
// Set this to 0 if you intend to use only MEAM potentials with
// spline knots on a uniform grid.
//
// With SUPPORT_NON_GRID_SPLINES == 0, the code runs about 50% faster.
#define SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES 0
class PairMEAMSpline : public Pair
{
public:
PairMEAMSpline(class LAMMPS *);
virtual ~PairMEAMSpline();
virtual 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);
PairMEAMSpline(class LAMMPS *);
virtual ~PairMEAMSpline();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void get_coeff(double *, double *);
double pair_density(int );
double three_body_density(int );
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();
// helper functions for compute()
int ij_to_potl(const int itype, const int jtype, const int ntypes) const {
return jtype - 1 + (itype-1)*ntypes - (itype-1)*itype/2;
}
int i_to_potl(const int itype) const { return itype-1; }
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();
protected:
char **elements; // names of unique elements
int *map; // mapping from atom types to elements
int nelements; // # of unique elements
class SplineFunction {
public:
class SplineFunction {
public:
/// Default constructor.
SplineFunction() : X(NULL), Xs(NULL), Y(NULL), Y2(NULL), Ydelta(NULL), N(0) {}
/// Default constructor.
SplineFunction() : X(NULL), Xs(NULL), Y(NULL), Y2(NULL), Ydelta(NULL), N(0) {}
/// Destructor.
~SplineFunction() {
delete[] X;
delete[] Xs;
delete[] Y;
delete[] Y2;
delete[] Ydelta;
}
/// Destructor.
~SplineFunction() {
delete[] X;
delete[] Xs;
delete[] Y;
delete[] Y2;
delete[] Ydelta;
}
/// Initialization of spline function.
void init(int _N, double _deriv0, double _derivN) {
N = _N;
deriv0 = _deriv0;
derivN = _derivN;
// if (X) delete[] X;
// if (Xs) delete[] Xs;
// if (Y) delete[] Y;
// if (Y2) delete[] Y2;
// if (Ydelta) delete[] Ydelta;
X = new double[N];
Xs = new double[N];
Y = new double[N];
Y2 = new double[N];
Ydelta = new double[N];
}
/// Initialization of spline function.
void init(int _n, double _deriv0, double _derivN) {
N = _n;
deriv0 = _deriv0;
derivN = _derivN;
delete[] X;
delete[] Xs;
delete[] Y;
delete[] Y2;
delete[] Ydelta;
X = new double[N];
Xs = new double[N];
Y = new double[N];
Y2 = new double[N];
Ydelta = new double[N];
}
/// Adds a knot to the spline.
void setKnot(int n, double x, double y) { X[n] = x; Y[n] = y; }
/// Adds a knot to the spline.
void setKnot(int n, double x, double y) { X[n] = x; Y[n] = y; }
/// Returns the number of knots.
int numKnots() const { return N; }
/// Returns the number of knots.
int numKnots() const { return N; }
/// Parses the spline knots from a text file.
void parse(FILE* fp, Error* error, bool isNewFormat);
/// Parses the spline knots from a text file.
void parse(FILE* fp, Error* error);
/// Calculates the second derivatives of the cubic spline.
void prepareSpline(Error* error);
/// Calculates the second derivatives of the cubic spline.
void prepareSpline(Error* error);
/// Evaluates the spline function at position x.
inline double eval(double x) const
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
/// Evaluates the spline function at position x.
inline double eval(double x) const
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
#if SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
return a * Y[klo] + b * Y[khi] +
((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
#else
// For a spline with grid points, we can directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
return Y[khi] - a * Ydelta[klo] + ((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
// For a spline with regular grid, we directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
return Y[khi] - a * Ydelta[klo] +
((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
#endif
}
}
}
}
/// Evaluates the spline function and its first derivative at position x.
inline double eval(double x, double& deriv) const
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
deriv = deriv0;
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
deriv = derivN;
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
/// Evaluates the spline function and its first derivative at position x.
inline double eval(double x, double& deriv) const
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
deriv = deriv0;
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
deriv = derivN;
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
#if SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
deriv = (Y[khi] - Y[klo]) / h + ((3.0*b*b - 1.0) * Y2[khi] - (3.0*a*a - 1.0) * Y2[klo]) * h / 6.0;
return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
deriv = (Y[khi] - Y[klo]) / h +
((3.0*b*b - 1.0) * Y2[khi] -
(3.0*a*a - 1.0) * Y2[klo]) * h / 6.0;
return a * Y[klo] + b * Y[khi] +
((a*a*a - a) * Y2[klo] +
(b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
#else
// For a spline with grid points, we can directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
deriv = Ydelta[klo] + ((3.0*b*b - hsq) * Y2[khi] - (3.0*a*a - hsq) * Y2[klo]);
return Y[khi] - a * Ydelta[klo] + ((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
// For a spline with regular grid, we directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
deriv = Ydelta[klo] + ((3.0*b*b - hsq) * Y2[khi]
- (3.0*a*a - hsq) * Y2[klo]);
return Y[khi] - a * Ydelta[klo] +
((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
#endif
}
}
}
}
/// Returns the number of bytes used by this function object.
double memory_usage() const { return sizeof(*this) + sizeof(X[0]) * N * 3; }
/// Returns the number of bytes used by this function object.
double memory_usage() const { return sizeof(*this) + sizeof(X[0]) * N * 3; }
/// Returns the cutoff radius of this function.
double cutoff() const { return X[N-1]; }
/// Returns the cutoff radius of this function.
double cutoff() const { return X[N-1]; }
/// Writes a Gnuplot script that plots the spline function.
void writeGnuplot(const char* filename, const char* title = NULL) const;
/// Writes a Gnuplot script that plots the spline function.
void writeGnuplot(const char* filename, const char* title = NULL) const;
/// Broadcasts the spline function parameters to all processors.
void communicate(MPI_Comm& world, int me);
/// Broadcasts the spline function parameters to all processors.
void communicate(MPI_Comm& world, int me);
private:
double* X; // Positions of spline knots
double* Xs; // Shifted positions of spline knots
double* Y; // Function values at spline knots
double* Y2; // Second derivatives at spline knots
double* Ydelta; // If this is a grid spline, Ydelta[i] = (Y[i+1]-Y[i])/h
int N; // Number of spline knots
double deriv0; // First derivative at knot 0
double derivN; // First derivative at knot (N-1)
double xmin; // The beginning of the interval on which the spline function is defined.
double xmax; // The end of the interval on which the spline function is defined.
int isGridSpline; // Indicates that all spline knots are on a regular grid.
double h; // The distance between knots if this is a grid spline with equidistant knots.
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
};
private:
double* X; // Positions of spline knots
double* Xs; // Shifted positions of spline knots
double* Y; // Function values at spline knots
double* Y2; // Second derivatives at spline knots
double* Ydelta; // If this is a grid spline, Ydelta[i] = (Y[i+1]-Y[i])/h
int N; // Number of spline knots
double deriv0; // First derivative at knot 0
double derivN; // First derivative at knot (N-1)
double xmin; // The beginning of the interval on which the spline function is defined.
double xmax; // The end of the interval on which the spline function is defined.
int isGridSpline;// Indicates that all spline knots are on a regular grid.
double h; // The distance between knots if this is a grid spline with equidistant knots.
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
};
/// Helper data structure for potential routine.
struct MEAM2Body {
int tag;
double r;
double f, fprime;
double del[3];
};
/// Helper data structure for potential routine.
struct MEAM2Body {
int tag; // holds the index of the second atom (j)
double r;
double f, fprime;
double del[3];
};
SplineFunction phi; // Phi(r_ij)
SplineFunction rho; // Rho(r_ij)
SplineFunction f; // f(r_ij)
SplineFunction U; // U(rho)
SplineFunction g; // g(cos_theta)
double zero_atom_energy; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
SplineFunction* phis; // Phi_i(r_ij)
SplineFunction* rhos; // Rho_ij(r_ij)
SplineFunction* fs; // f_i(r_ij)
SplineFunction* Us; // U_i(rho)
SplineFunction* gs; // g_ij(cos_theta)
double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
double cutoff; // The cutoff radius
double cutoff; // The cutoff radius
double* Uprime_values; // Used for temporary storage of U'(rho) values
int nmax; // Size of temporary array.
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
MEAM2Body* twoBodyInfo; // Temporary array.
double* Uprime_values; // Used for temporary storage of U'(rho) values
int nmax; // Size of temporary array.
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
MEAM2Body* twoBodyInfo; // Temporary array.
void read_file(const char* filename);
void allocate();
void read_file(const char* filename);
void allocate();
};
}
@ -279,3 +300,5 @@ protected:
*
* See file 'pair_spline_meam.cpp' for history of changes.
------------------------------------------------------------------------- */

31
src/USER-OMP/pair_meam_spline_omp.cpp
View File

@ -110,6 +110,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
const int nthreads = comm->nthreads;
const int nlocal = atom->nlocal;
const int nall = nlocal + atom->nghost;
const int ntypes = atom->ntypes;
const double cutforcesq = cutoff*cutoff;
@ -135,33 +136,38 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
const double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
if (rij_sq < cutforcesq) {
const int jtype = atom->type[j];
const double rij = sqrt(rij_sq);
double partial_sum = 0;
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->f = fs[i_to_potl(jtype)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = myTwoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] + nextTwoBodyInfo->del[1]*bondk.del[1] + nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * g.eval(cos_theta);
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(jtype,atom->type[bondk.tag],ntypes)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rho.eval(rij);
rho_value += rhos[i_to_potl(jtype)].eval(rij);
numBonds++;
nextTwoBodyInfo++;
}
}
const int itype = atom->type[i];
// Compute embedding energy and its derivative.
double Uprime_i;
double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy;
double embeddingEnergy = Us[i_to_potl(itype)].eval(rho_value, Uprime_i)
- zero_atom_energies[i_to_potl(itype)];
Uprime_thr[i] = Uprime_i;
if (EFLAG)
e_tally_thr(this,i,i,nlocal,1/*newton_pair*/,embeddingEnergy,0.0,thr);
@ -173,6 +179,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
const MEAM2Body bondj = myTwoBodyInfo[jj];
const double rij = bondj.r;
const int j = bondj.tag;
const int jtype = atom->type[j];
const double f_rij_prime = bondj.fprime;
const double f_rij = bondj.f;
@ -187,7 +194,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
+ bondj.del[1]*bondk->del[1]
+ bondj.del[2]*bondk->del[2]);
double g_prime;
double g_value = g.eval(cos_theta, g_prime);
double g_value = gs[ij_to_potl(jtype,atom->type[bondk->tag],ntypes)].eval(cos_theta, g_prime);
const double f_rik_prime = bondk->fprime;
const double f_rik = bondk->f;
@ -279,6 +286,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
const double ztmp = x[i][2];
const int* const jlist = firstneigh_half[i];
const int jnum = numneigh_half[i];
const int itype = atom->type[i];
for(int jj = 0; jj < jnum; jj++) {
const int j = jlist[jj] & NEIGHMASK;
@ -291,13 +299,16 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
if(rij_sq < cutforcesq) {
double rij = sqrt(rij_sq);
const int jtype = atom->type[j];
double rho_prime;
rho.eval(rij, rho_prime);
double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]);
double rho_prime_i,rho_prime_j;
rhos[i_to_potl(itype)].eval(rij,rho_prime_i);
rhos[i_to_potl(jtype)].eval(rij,rho_prime_j);
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
double pair_pot_deriv;
double pair_pot = phi.eval(rij, pair_pot_deriv);
double pair_pot = phis[ij_to_potl(itype,jtype,ntypes)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient.